新年快乐~
1、chipseeker包的介绍：
①、南科大Y叔写的包，感谢
②、主要是对peaks的注释和可视化
2、这次用的数据，是chipseeker包内置的数据。
"GSM1174480_ARmo_0M_peaks.bed.gz" ---雄激素受体的过表达增强了前列腺癌中chip的数据
"GSM1174481_ARmo_1nM_peaks.bed.gz"
"GSM1174482_ARmo_100nM_peaks.bed.gz" "GSM1295076_CBX6_BF_ChipSeq_mergedReps_peaks.bed.gz"
"GSM1295077_CBX7_BF_ChipSeq_mergedReps_peaks.bed.gz"---是乳腺癌的chip的分析的数据
3、bed文件介绍：
BED 文件格式
chrom - 染色体号; 例如，chr1，chrX。。。。。。。
chromStart - feature在染色体上起始位置. 从0开始算，染色体上第一个碱基位置标记为0。
chromEnd - feature在染色体上终止位置。
name
strand - 正负链标记 这里没有
thickStart - feature起始位置
thickEnd - feature编码终止位置
blockSizes - blocks (exons)大小列表，逗号分隔，对应于blockCount.
blockStarts -blocks (exons)起始位置列表，逗号分隔，对应于blockCount.；这个起始位置是与chromStart的一个相对位置。
4、开始流程
①、chip在整个染色体上的分布
BiocManager::install("ChIPseeker")
library(ChIPseeker)
library(ggplot2)
files <- getSampleFiles()
basename(unlist(files))
tiff(filename = "chip_peaks_over_chromosomes",res = 600, width = 4800, height = 4800, compression = "lzw") lzw是压缩的模式
covplot(files[[4]])
dev.off()
input：
可以看出，这个数据“GSM1295077_CBX7_BF_ChipSeq_mergedReps_peaks.bed.gz”的chip在染色体上的分布。其中是1和2号染色体的chip最多。
②、peaks注释：
介绍：peaks注释，主要是利用TxDb文件注释、注释信息一般要包含基因的起始终止，基因的外显子、内含子及它们的起始终止、非编码区域位置、功能元件的位置等。
Bioconductor包提供了30个TxDb包
如果没有的话，好像是可以自己做一个TxDb文件的。
注释流程
library(TxDb.Hsapiens.UCSC.hg19.knownGene)
txdb = TxDb.Hsapiens.UCSC.hg19.knownGene
f = getSampleFiles()[[4]]
x = annotatePeak(f, tssRegion=c(-1000, 1000), TxDb=txdb)
peakAnno <- annotatePeak(files[[4]],
tssRegion=c(-3000, 3000),
TxDb=txdb, annoDb="org.Hs.eg.db")
③、对注释的结果可视化
#pie图
plotAnnoPie(peakAnno)
vennpie(peakAnno)
#TSS的强度分布图
peakHeatmap(f, weightCol="V5", TxDb=txdb,
upstream=3000, downstream=3000,
color=rainbow(length(f)))
#几个数据的chip的交集的韦恩图
peakAnnoList <- lapply(files, annotatePeak,
TxDb=txdb,tssRegion=c(-3000, 3000))
genes <- lapply(peakAnnoList, function(i)
as.data.frame(i)$geneId)
library(Vennerable)
vennplot(genes[2:4], by='Vennerable')
input：
可以看出，chip主要是集中在promoter启动子，distal intergenic远端基因间。
转录起始点（TSS）的强度分布。
韦恩图看，几个数据集的交集。？？作用是？
④、KEGG分析
require(clusterProfiler)
bedfile=getSampleFiles() # 将bed文件读入（readPeakFile是利用read.delim读取，然后转为GRanges对象）
seq=lapply(bedfile, readPeakFile)
genes=lapply(seq, function(i)
seq2gene(i, c(-1000, 3000), 3000, TxDb=txdb))
cc = compareCluster(geneClusters = genes,
fun="enrichKEGG", organism="hsa")
dotplot(cc, showCategory=10)
input：
可以看到chip主要富集的几个通路：
Arrhythmogenic right ventricular cardiomyopathy --致心律失常性右心室心肌病
Glutamatergic synapse--Glutamatergic--突触
Axon guidance---轴突的指导
Transcriptional misregulation in cancer--癌症中的转录调控失调
