學(xué)習(xí)linux常用命令麦锯，并且實戰(zhàn)ChIP-seq已經(jīng)一個月了√成疲現(xiàn)在需要用R包ChIPseeker進(jìn)行注釋，發(fā)現(xiàn)R已經(jīng)忘了大半搞乏。補(bǔ)課R語言波桩。。對peak的注釋分為兩個部分——結(jié)構(gòu)注釋和功能注釋
結(jié)構(gòu)注釋：會將peak所落在基因組上的區(qū)域結(jié)構(gòu)注釋出來请敦，比如說啟動子區(qū)域镐躲，UTR區(qū)域，內(nèi)含子區(qū)域等等侍筛。同時萤皂，也會將與peak最臨近的基因注釋出來，非常好用匣椰。
功能注釋：其實并不是對peak進(jìn)行注釋的裆熙，而是對peak最臨近的基因注釋的，因此這部分就和傳統(tǒng)的GO/KEGG等分析如出一轍啦禽笑。
一入录、安裝R包ChIPseeker，各種報錯及解決
按照教程在console輸入：BiocManager::install("ChIPseeker")
報錯：Error in loadNamespace(x) : there is no package called ‘BiocManager’
于是輸入：install.packages("BiocManager") 先安裝BiocManager佳镜。
再輸入BiocManager::install("ChIPseeker")
運(yùn)行了一會兒出現(xiàn)如下報錯：
Warning messages:
1: In install.packages(...) :
installation of package ‘DO.db’ had non-zero exit status
2: In install.packages(...) :
installation of package ‘GO.db’ had non-zero exit status
3: In install.packages(...) :
installation of package ‘igraph’ had non-zero exit status
4: In install.packages(...) :
installation of package ‘gtools’ had non-zero exit status
還需要安裝一個GO和DO包,和其他包僚稿。
BiocManager::install("GO.db")
BiocManager::install("DO.db")
BiocManager::install("topGO")
BiocManager::install("GSEABase")
BiocManager::install("Rgraphviz")
install.packages("igraph")
install.packages("gtools")
在安裝igraph和gtools時候，出現(xiàn)提示：Do you want to install from sources the packages which need compilation?
【解決辦法】
選擇“NO”蟀伸，而不是“YES”蚀同，便可正確安裝igraph和gtools缅刽。
二、注釋庫
BiocManager::install("org.Hs.eg.db")
BiocManager::install("org.Mm.eg.db")
BiocManager::install("TxDb.Mmusculus.UCSC.mm10.knownGene")
BiocManager::install("TxDb.Hsapiens.UCSC.hg38.knownGene")
這些是基因組注釋庫唤崭，最常用的人和小鼠的拷恨。
當(dāng)然也可以用gff/gtf構(gòu)建一個TxDb（注釋庫），略谢肾。
三、bed數(shù)據(jù)
把數(shù)據(jù)從服務(wù)器下載到Mac電腦上小泉。一定要先連上VPN芦疏，在本電腦的Terminal中進(jìn)行操作。
scp -r zds209@222.28.163.113:~/ChIP-seqtest/bed/ /Users/meraner/Desktop
-r是整個文件夾里的文件都下載到bed里面微姊。速度能到3.7M/s酸茴。
四、使用ChIPseeker
參考教程
http://events.jianshu.io/p/9aa719faa4b5
http://www.reibang.com/p/c76e83e6fa57
https://www.modb.pro/db/401713
http://www.bioconductor.org/packages/devel/bioc/vignettes/ChIPseeker/inst/doc/ChIPseeker.html #官網(wǎng)
https://blog.csdn.net/qq_39306047/article/details/106601095
1. 先看看整體分布吧兢交。
library(ChIPseeker) #加載ChIPseeker包
library(ggplot2) #加載ggplot2包
library(org.Mm.eg.db) #加載參考基因組數(shù)據(jù)
library(TxDb.Mmusculus.UCSC.mm10.knownGene) #加載參考基因組注釋數(shù)據(jù)包
library(clusterProfiler) #加載clusterProfiler包
library(GenomicFeatures) #加載GenomicFeatures包
eithtWG16_1<-readPeakFile("8WG16-1.sorted.bam_peaks.narrowPeak") #讀取narrowPeak數(shù)據(jù)
covplot(eithtWG16_1,weightCol=5) #可視化一下數(shù)據(jù)全景
這樣的圖片薪捍，指示的是在每條染色體上的富集程度
還能進(jìn)行多樣本，及制定區(qū)域的展示配喳。
peak=GenomicRanges::GRangesList(WG1=readPeakFile("8WG16-1.sorted.bam_peaks.narrowPeak"),WG2=readPeakFile("8WG16-2.sorted.bam_peaks.narrowPeak"))#（先指定的在右側(cè)酪穿，這個比較奇怪，所以以后還是要改成2號在前) 使用ChIPseeker包中的readPeakFile函數(shù)將bed文件讀入到R晴裹，并存儲為GRanges對象
covplot(peak, weightCol="V5") + facet_grid(chr ~ .id)#畫圖多樣本展示每條染色體的富集程度被济，也可以理解為coverage plot 展示峰在染色體上分布的圖。
covplot(peak, weightCol="V5", chrs=c("chr17", "chr18"), xlim=c(4e7, 5e7)) + facet_grid(chr ~ .id)#截取某些片段進(jìn)行展示涧团。
2.在啟動子區(qū)(TSS)附近的 ChIP peaks 的可視化
（注：這個也可以用deeptools在linux服務(wù)器里面去跑只磷。day26都在學(xué)deeptools。）
peak <- readPeakFile("8WG16-1.sorted.bam_summits.bed")#讀取bed文件
txdb <- TxDb.Mmusculus.UCSC.mm10.knownGene #指定txdb
promoter <- getPromoters(TxDb=txdb, upstream=3000, downstream=3000) #getPromoters函數(shù)定義啟動子區(qū)域的范圍
tagMatrix <- getTagMatrix(peak, windows=promoter) # getTagMatrix函數(shù)將Peaks匹配到啟動子區(qū)域泌绣，然后計算出tagMatrix用于畫圖
tagHeatmap(tagMatrix, xlim=c(-3000, 3000), color="red") #熱圖
plotAvgProf(tagMatrix, xlim=c(-3000, 3000),
xlab="Genomic Region (5'->3')", ylab = "Read Count Frequency") #折線圖
plotAvgProf(tagMatrix, xlim=c(-3000, 3000),conf = 0.95, resample = 1000,
xlab="Genomic Region (5'->3')", ylab = "Read Count Frequency") #折線圖上加置信區(qū)間
多個樣本钮追，在啟動子區(qū)(TSS)附近的 ChIP peaks 的可視化
peak1 <- "/Users/meraner/Desktop/bed/8WG16-1.sorted.bam_summits.bed"
peak2 <- "/Users/meraner/Desktop/bed/8WG16-1.sorted.bam_summits.bed"
files <-list(peak1=peak1,peak2=peak2)#把兩個bed文件指定給files, 就是一次讀入多個summits文件，使用list存儲阿迈，然后使用lapply注釋元媚。
tagMatrixList <- lapply(files, getTagMatrix, windows=promoter)
plotAvgProf(tagMatrixList, xlim=c(-3000, 3000),
conf=0.95,resample=500, facet="row") # 添加置信區(qū)間并分面
3. 結(jié)構(gòu)注釋：
library(ChIPseeker)
library(ggplot2)
library(TxDb.Mmusculus.UCSC.mm10.knownGene)
library(org.Mm.eg.db)
peak <- readPeakFile("8WG16-1.sorted.bam_summits.bed")#讀取bed文件
peakAnno <- annotatePeak(peak,TxDb = TxDb.Mmusculus.UCSC.mm10.knownGene, annoDb ="org.Mm.eg.db" ,tssRegion=c(-3000, 3000)) # 進(jìn)行結(jié)構(gòu)注釋
注釋完，進(jìn)行可視化仿滔，多種圖可供選擇
plotAnnoBar(peakAnno) #柱狀圖
plotDistToTSS(peakAnno, title="Distribution of transcription factor-binding loci \n relative to TSS") #展示轉(zhuǎn)錄因子(TF)與基因啟動子區(qū)的相對位置vennpie(peakAnno)
plotAnnoPie(peakAnno) #餅圖
上面這幾個圖都是能看出peak的分布結(jié)構(gòu)的
還有另外軟件可以對圖進(jìn)行優(yōu)化惠毁。安裝下面兩個包。
BiocManager::install("ggupset")
install.packages("ggimage")
library(ggupset)
library(ggimage)
然后可以運(yùn)行
upsetplot(peakAnno)
upsetplot(peakAnno, vennpie=TRUE)
4.多個樣本結(jié)構(gòu)注釋
library(ChIPseeker) #加載ChIPseeker包
library(ggplot2) #加載ggplot2包
library(org.Mm.eg.db) #加載參考基因組數(shù)據(jù)
library(TxDb.Mmusculus.UCSC.mm10.knownGene) #加載參考基因組注釋數(shù)據(jù)包
library(clusterProfiler) #加載clusterProfiler包
library(GenomicFeatures) #加載GenomicFeatures包
txdb <- TxDb.Mmusculus.UCSC.mm10.knownGene
setwd("/Users/meraner/Desktop/bed") #制定當(dāng)前工作目錄崎页，為數(shù)據(jù)存儲目錄
peak1 <- "/Users/meraner/Desktop/bed/8WG16-1.sorted.bam_summits.bed"
peak2 <- "/Users/meraner/Desktop/bed/8WG16-1.sorted.bam_summits.bed"
files <-list(peak1=peak1,peak2=peak2)#把兩個bed文件指定給files, 就是一次讀入多個summits文件鞠绰，使用list存儲，然后使用lapply注釋飒焦。
peakAnnoList <- lapply(files, annotatePeak, TxDb=TxDb.Mmusculus.UCSC.mm10.knownGene,
tssRegion=c(-3000, 3000),
verbose=FALSE) #注釋
plotAnnoBar(peakAnnoList) #畫圖bar
plotDistToTSS(peakAnnoList) #畫圖距離
#不能做多樣本的餅圖和韋恩圖蜈膨。
在看完數(shù)據(jù)全景之后屿笼，就會想知道這些peaks和什么類型的基因有關(guān)。
5.功能注釋：
基因注釋+富集分析
利用ChIPseeker的seq2gene 將peak的位置與所有的基因關(guān)聯(lián)起來【包括 host gene (exon/intron), promoter region and flanking gene from intergenomic region】翁巍，然后用clusterProfiler拿這些基因跑ORA驴一，做富集
下面是根據(jù)寫好的R腳本，執(zhí)行順利灶壶。
6. 單樣本ChIPseeker 腳本
#工作目錄肝断，樣本種屬來源需要根據(jù)情況調(diào)整
setwd("/Users/meraner/Desktop/bed") #制定當(dāng)前工作目錄，為數(shù)據(jù)存儲目錄
library(ChIPseeker) #加載ChIPseeker包
library(ggplot2) #加載ggplot2包
library(org.Mm.eg.db) #加載參考基因組數(shù)據(jù)
library(TxDb.Mmusculus.UCSC.mm10.knownGene) #加載參考基因組注釋數(shù)據(jù)包
library(clusterProfiler) #加載clusterProfiler包
library(GenomicFeatures) #加載GenomicFeatures包
txdb <- TxDb.Mmusculus.UCSC.mm10.knownGene
#peak的整體分布
peak<-readPeakFile("8WG16-1.sorted.bam_peaks.narrowPeak") #讀取narrowPeak數(shù)據(jù)
covplot(peak, weightCol=5) #可視化一下數(shù)據(jù)全景
covplot(peak, weightCol=5, chrs=c("chr17", "chr18"), xlim=c(4e7, 5e7)) #截取某些染色體驰凛，某些片段進(jìn)行展示
#單一樣本胸懈，在啟動子區(qū)(TSS)附近的 ChIP peaks 的可視化
peak <- readPeakFile("8WG16-1.sorted.bam_summits.bed") #讀取bed文件
promoter <- getPromoters(TxDb=txdb, upstream=3000, downstream=3000) #getPromoters函數(shù)定義啟動子區(qū)域的范圍
tagMatrix <- getTagMatrix(peak, windows=promoter) # getTagMatrix函數(shù)將Peaks匹配到啟動子區(qū)域，然后計算出tagMatrix用于畫圖
tagHeatmap(tagMatrix, xlim=c(-3000, 3000), color="red") #一個樣本在TSS附近結(jié)合區(qū)的熱圖
plotAvgProf(tagMatrix, xlim=c(-3000, 3000),
xlab="Genomic Region (5'->3')", ylab = "Read Count Frequency") #折線圖恰响，
plotAvgProf(tagMatrix, xlim=c(-3000, 3000),conf = 0.95, resample = 1000,
xlab="Genomic Region (5'->3')", ylab = "Read Count Frequency") #前面的折線圖上加置信區(qū)間
#注釋趣钱，結(jié)構(gòu)可視化
peakAnno <- annotatePeak(peak,TxDb = txdb, annoDb ="org.Mm.eg.db" ,tssRegion=c(-3000, 3000)) # 進(jìn)行結(jié)構(gòu)注釋
plotAnnoBar(peakAnno) #柱狀圖
plotDistToTSS(peakAnno, title="Distribution of transcription factor-binding loci \n relative to TSS") #展示轉(zhuǎn)錄因子(TF)與基因啟動子區(qū)的相對位置vennpie(peakAnno)
plotAnnoPie(peakAnno) #餅圖
library(ggupset)
library(ggimage)
upsetplot(peakAnno)
upsetplot(peakAnno, vennpie=TRUE)
#功能可視化-基因注釋 + 富集分析
library(DOSE)
library(GO.db)
library(topGO)
library(GSEABase)
library(clusterProfiler)
require(clusterProfiler)
gene <- seq2gene(peak, tssRegion = c(-1000, 1000), flankDistance = 3000, TxDb=txdb)
enk <- enrichKEGG(gene=gene,
organism = 'mmu',
pvalueCutoff = 0.05,
pAdjustMethod = "BH",
qvalueCutoff = 0.05) #KEGG通路富集
dotplot(enk)#可視化KEGG富集
browseKEGG(enk, 'mmu04110') # 會打開瀏覽器調(diào)到KEGG數(shù)據(jù)庫的這條通路上，并且富集的基因會以不同的顏色標(biāo)出
engo_MF <- enrichGO(gene = gene,
OrgDb = org.Mm.eg.db,
ont = "MF",
pAdjustMethod = "BH",
pvalueCutoff = 0.05,
qvalueCutoff = 0.05,
readable = TRUE) #GO富集分析MF
dotplot(engo_MF) #可視化GO富集
barplot(engo_MF, showCategory=15) #可視化GO富集-bar圖
plotGOgraph(engo_MF) #對Go通路樹狀圖
engo_BP <- enrichGO(gene = gene,
OrgDb = org.Mm.eg.db,
ont = "BP",
pAdjustMethod = "BH",
pvalueCutoff = 0.05,
qvalueCutoff = 0.05,
readable = TRUE) #GO富集分析BP
dotplot(engo_BP) #可視化GO富集-點(diǎn)圖
barplot(engo_BP, showCategory=15) #可視化GO富集-bar圖
enrichMap(engo_BP) #這個命令有問題,報錯could not find function "enrichMap"胚宦？沒搞懂首有，也許是R版本問題。
plotGOgraph(engo_BP) #對Go通路樹狀圖
#輸出注釋信息
peakAnnotable<-as.data.frame(peakAnno)
write.table(peakAnnotable,file="result_genes",sep="\t",quote=F) #把peak相關(guān)基因?qū)С霰砀?
enktable<-as.data.frame(enk)
write.table(enktable,file="result_kegg.txt",sep="\t",quote=F)#把KEGG富集結(jié)果導(dǎo)出表格
MFtable<-as.data.frame(engo_MF)
write.table(MFtable,file="result_GO_MF.txt",sep="\t",quote=F)#把GO_MF富集結(jié)果導(dǎo)出表格
BPtable<-as.data.frame(engo_BP)
write.table(BPtable,file="result_GO_MF.txt",sep="\t",quote=F)#把GO_BP富集結(jié)果導(dǎo)出表格
7.多樣本ChIPseeker分析腳本
#工作目錄枢劝，樣本種屬來源需要根據(jù)情況調(diào)整
setwd("/Users/meraner/Desktop/bed") #制定當(dāng)前工作目錄井联，為數(shù)據(jù)存儲目錄
library(ChIPseeker) #加載ChIPseeker包
library(ggplot2) #加載ggplot2包
library(org.Mm.eg.db) #加載參考基因組數(shù)據(jù)
library(TxDb.Mmusculus.UCSC.mm10.knownGene) #加載參考基因組注釋數(shù)據(jù)包
library(clusterProfiler) #加載clusterProfiler包
library(GenomicFeatures) #加載GenomicFeatures包
txdb <- TxDb.Mmusculus.UCSC.mm10.knownGene
#peak的整體分布
peak=GenomicRanges::GRangesList("8WG16_1"=readPeakFile("8WG16-1.sorted.bam_peaks.narrowPeak"),"8WG16_2"=readPeakFile("8WG16-2.sorted.bam_peaks.narrowPeak"))
covplot(peak, weightCol="V5") + facet_grid(chr ~ .id)#coverageplot
covplot(peak, weightCol="V5", chrs=c("chr17", "chr18"), xlim=c(4e7, 5e7)) + facet_grid(chr ~ .id)#截取某些染色體，某些片段進(jìn)行展示
#多個樣本呈野，在啟動子區(qū)(TSS)附近的 ChIP peaks 的可視化
peak1 <- "/Users/meraner/Desktop/bed/8WG16-1.sorted.bam_summits.bed"
peak2 <- "/Users/meraner/Desktop/bed/8WG16-1.sorted.bam_summits.bed"
files <-list(peak1=peak1,peak2=peak2)#把兩個bed文件指定給files, 就是一次讀入多個summits文件低矮，使用list存儲，然后使用lapply注釋被冒。
promoter <- getPromoters (TxDb=txdb, upstream=3000, downstream=3000) #getPromoters函數(shù)定義啟動子區(qū)域的范圍
files <-list(peak1=peak1,peak2=peak2)#把兩個bed文件指定給files, 就是一次讀入多個summits文件军掂，使用list存儲，然后使用lapply注釋昨悼。
tagMatrixList <- lapply(files, getTagMatrix, windows=promoter)
peakHeatmap(files, TxDb=txdb,
upstream=3000, downstream=3000,
color=rainbow(length(files))) #多個樣本在TSS附近結(jié)合區(qū)的熱圖
plotAvgProf(tagMatrixList, xlim=c(-3000, 3000),
conf=0.95,resample=500, facet="row") # 添加置信區(qū)間并分面
#注釋蝗锥，結(jié)構(gòu)可視化
peakAnnoList <- lapply(files, annotatePeak, TxDb=txdb,
tssRegion=c(-3000, 3000),
verbose=FALSE) #注釋
plotAnnoBar(peakAnnoList) #畫圖bar
plotDistToTSS(peakAnnoList) #畫圖距離
#單樣本可以做餅圖和韋恩圖，多樣本的不行率触。如果需要终议，只能用單樣本流程。
#功能可視化-基因注釋 + 富集分析
require(clusterProfiler)
seq=lapply(files, readPeakFile)
genes <-lapply(seq, function(i)
seq2gene(i, c(-1000, 3000), 3000, TxDb=txdb))
cc <- compareCluster(geneClusters = genes, organism="mmu",
fun="enrichKEGG")
dotplot(cc, showCategory=10)
#輸出注釋信息
peakAnnotable<-as.data.frame(peakAnnoList)
write.table(peakAnnotable,file="result_multisample_gene",sep="\t",quote=F) #導(dǎo)出基因信息
cctable<-as.data.frame(cc)
write.table(cctable,file="result_multisample_kegg.txt",sep="\t",quote=F)#把KEGG富集結(jié)果導(dǎo)出表格
