本周學(xué)習(xí)一下CHIP-seq。 并根據(jù)網(wǎng)上的教程，自己實踐一下， 一方面是要為了弄清楚什么是chip-seq, 這個技術(shù)有什么用，另一個方面是想學(xué)習(xí)一下這個技術(shù)如何來實踐， 本文參考的文章主要來自生信技能樹，以及簡書上的其他作者寫的教程，由于每個人在做分析時，使用的操作系統(tǒng)不一樣，所以網(wǎng)上的代碼在自己的電腦上進(jìn)行運行的時候經(jīng)常出現(xiàn)問題，這需要每個人針對自己的情況進(jìn)行分析和總結(jié)。 本次分析采用的系統(tǒng)是MACos 10.12.6。

先不談CHIP-seq 的原理， 因為我正在學(xué)習(xí)中，后期將這部分內(nèi)容補上。先根據(jù)文獻(xiàn)提供的數(shù)據(jù)和網(wǎng)上的代碼進(jìn)行實際操作。

CHIP-seq能干什么？

明確每一類組蛋白或者轉(zhuǎn)錄因子在整個基因組上結(jié)合基因的位置
如果比較多個組蛋白在亞基，可以看這些亞基之間在基因組上結(jié)合的基因的包含關(guān)系，即用韋恩圖展示這些組蛋白結(jié)合基因相互之間是否包含。
檢查每一類組蛋白結(jié)合基因在TSS上的位置。
檢查每一組（不同組蛋白之間結(jié)合相同的基因）在TSS上的位置。（這樣可以看出缺少某一類組蛋白之后，基因是否表達(dá)，驗證這個組蛋白具有的功能和意義）
不同組蛋白結(jié)合基因的功能（GO），及參與的代謝通路（KEGG）
可以研究每一個組蛋白targets 的基因的表達(dá)
第一步：從NCBI下載數(shù)據(jù)，并解壓到本地電腦。 從NCBI的GEO dataset 中輸入作者上傳的GEO號碼GSE42466,如下圖所示：

在本文文章的哪里打?qū)?#xff0c;并進(jìn)入，之后看一看到文章的具體信息，包括作者的信息，以及實驗方法，實驗設(shè)計，實驗上傳的具體數(shù)據(jù)編號。如下所示：

在最后的一欄中找到本數(shù)據(jù)的SRA號碼，點擊進(jìn)入，如下：

在上面作者提供的6個數(shù)據(jù)上打上對勾，并在右上側(cè)的send to 這個框中選擇file, 在format 中選擇 runinfo. 點擊生成文件，即可生成下載的文件，這個是個excel文件，包含了具體run的信息，我們需要的是run 的ID號碼，打開EXCEL文件，并在download 那一列獲取SRA號SRR620204。 這個號碼用于下載。代碼如下：

prefetch SRR620204
如果批量下載，則將上面的文件編號存放到一個文本文章中，如 sradata.txt ,下載代碼如下：

prefetch --option-file sradata.txt #install prefetch first before run the code
之后用fastq-dump 軟件進(jìn)行解壓，如下：

ls *sra |while read id; do fastq-dump –split-3 $id;done # --split-3的目的是如果文件包含兩個以上文件，則分別命名，如果是一個文件則直接是文件原名，而不是根據(jù)reads 來進(jìn)行拆分
以上可以獲得本實驗的的所有數(shù)據(jù)，但是在實際操作中，SRR620207一致下載不下來，原因不知道。 我暫且不比對此文件，在后期使用作者提供的文件。

第二步： 對原始數(shù)據(jù)進(jìn)行質(zhì)控，采用fastqc, multiqc 兩個軟件， fastqc對每個文件進(jìn)行質(zhì)控分析， multiqc對fastqc的結(jié)果進(jìn)行整合，方便讀者從總體上對數(shù)據(jù)質(zhì)量進(jìn)行把控。

ls *.fastq |while read id; do fastqc -t 3 $id; done
multiqc *fastqc.zip --pdf # multiqc 在我的電腦山無法生成PDF，雖然我按照說明書安裝了pandoc, latex. 但依然不行.
根據(jù)multiqc的質(zhì)控顯示，本測序數(shù)據(jù)在3’端的數(shù)據(jù)質(zhì)量不高，fastqc 也顯示警告，原因是平均的堿基質(zhì)量有點低，如下圖所示：

因此在接下來比對的時候，我們要去掉3‘端質(zhì)量不好的堿基，用bowtie2自動去掉。

第三步： 利用bowtie2（本次分析采用的是bowtie2-2.3.5）對數(shù)據(jù)進(jìn)行mapping， 數(shù)據(jù)索引下載如下：

wget -4 -q ftp://ftp.ccb.jhu.edu/pub/data/bowtie2_indexes/mm10.zip
unzip mm10.zip
比對程序如下：

ls *.fastq|while read id;
do bowtie2 -p 3 -3 5 --local -x ~/src/GSE42466/data/chip-se/reference/mm10 -U $id | /opt/biosoft/samtools-1.3.1/samtools sort -O bam -o $id.bam;
done
比對后的總體結(jié)果如下：

SRR620204

19687027 reads; of these:

19687027 (100.00%) were unpaired; of these:

7960096 (40.43%) aligned 0 times8614711 (43.76%) aligned exactly 1 time3112220 (15.81%) aligned >1 times

59.57% overall alignment rate

SRR620205

20710168 reads; of these:

20710168 (100.00%) were unpaired; of these:

2716388 (13.12%) aligned 0 times12169372 (58.76%) aligned exactly 1 time5824408 (28.12%) aligned >1 times

86.88% overall alignment rate

SRR620206

21551927 reads; of these:

21551927 (100.00%) were unpaired; of these:

7316904 (33.95%) aligned 0 times10534163 (48.88%) aligned exactly 1 time3700860 (17.17%) aligned >1 times

66.05% overall alignment rate

SRR620208

13291429 reads; of these:

13291429 (100.00%) were unpaired; of these:

5731176 (43.12%) aligned 0 times5235529 (39.39%) aligned exactly 1 time2324724 (17.49%) aligned >1 times

56.88% overall alignment rate

SRR620209

31218866 reads; of these:

31218866 (100.00%) were unpaired; of these:

5699289 (18.26%) aligned 0 times17025381 (54.54%) aligned exactly 1 time8494196 (27.21%) aligned >1 times

81.74% overall alignment rate

#############

比對完之后，利用IGV進(jìn)行可視化，先用samtools進(jìn)行index, 之后在IGV中導(dǎo)入，導(dǎo)入的時候選擇BAM文件即可，INDEX文件也必須在bam文件相同的文件夾。

samtools index FGENESH_C15+B10.sorted.bam FGENESH_C15+B10.sorted.bai
從IGV上確實可以發(fā)現(xiàn)， chip-seq 測序完在有target 的地方有峰值，沒有target的地方對照和實驗都沒有，如下所示：

之后，我們使用macs2進(jìn)行peal calling， 在使用前請先安裝macs2.

macs2 callpeak -c igG_old.sorted.bam -t Suz12.sorted.bam -q 0.05 -f BAM -g mm -n suz12 2> suz12.macs2.log
macs2 callpeak -c igG_old.sorted.bam -t Ring1B.sorted.bam -q 0.05 -f BAM -g mm -n Ring1B 2> Ring1B.macs2.log
macs2 callpeak -c igG.sorted.bam -t Cbx7.sorted.bam -q 0.05 -f BAM -g mm -n cbx7 2> cbx7.macs2.log
比對之后進(jìn)行數(shù)據(jù)的作圖和可視化，代碼

ls bam |while read id;
do file=$(basename $id);sample=${file%%.};
echo $sample;
bamCoverage -b $id -o $sample.bw; ##–normalizeUsing RPKM -p 5
computeMatrix reference-point --referencePoint TSS -b 10000 -a 10000 -R mm10_ensenble_ref.bed -S KaTeX parse error: Expected group after '_' at position 33: …eros -o matrix1_?{sample}TSS.gz --outFileSortedRegions regions1KaTeX parse error: Expected 'EOF', got '#' at position 28: …es.bed -p 20; #?計算matrix 的時候特別慢…{sample}_TSS.gz -out ${sample}.png;
done
在這個時候，在我的macos 系統(tǒng)上出現(xiàn)了一個小插曲， 使用bamCoverage 的時候報錯，原因如下：

可能的原因是沒有找到相應(yīng)的庫，于是我重新安裝了libssh2. 結(jié)果顯示可以運行了。

此外在運行computeMatrix 的時候也報錯，錯誤的原因是提供的bed文件有問題，于是，我有重新用gff文件制作了一個bed文件（注意：這個地方的gff文件必須和參考基因組是一直的，必須是針對這個參考基因組的gff3文件，在R中如果要用CHIPSeeker包進(jìn)行注釋的時候，當(dāng)上面沒有注釋的數(shù)據(jù)庫時，需要用GFF3文件來自己制作一個注釋數(shù)據(jù)庫，這個時候gff3也必須是同mapping時用的的參考基因組的注釋文件），首先下載bedops軟件，之后采用如下代碼運行：

convert2bed --input=gff [–output=bed] <ensembl.mm10.gff3> mm10_ensembl_GFF3_genes.bed
#取上述bed文件的前三行即可
cut -f 1-3 mm10_ensembl_GFF3_genes.bed > mm10_ensenble_ref.bed
至此，所有要畫圖的文件都準(zhǔn)備好了，用plotHeatmap進(jìn)行畫圖，代碼如下： 在這里我只畫了cbx7的，由于計算comuteMatrix 太慢了。

plotHeatmap -m matrix.gz -out plotHeatmap.png --plotTitle “Heatmap”

也可以對多個文件整合在一起畫，代碼如下：

computeMatrix reference-point -p 28 --referencePoint TSS -b 2500 -a 2500 -S *bw -R mm10_ensembl_GFF3_genes.bed --skipZeros -o tmp4.mat.gz --outFileNameMatrix alBamFile
plotHeatmap -m tmp4.mat.gz -out tmp4.merge.png
plotProfile --dpi 720 -m tmp4.mat.gz -out tmp4.profile.pdf --plotFileFormat pdf --perGroup
plotHeatmap --dpi 720 -m tmp4.mat.gz -out tmp4.merge.pdf --plotFileFormat pdf
效果如下：

上面的圖都是針對整個基因組上，所有基因在選定的TSS上下游畫的，那可不可以對某一部分感興趣的基因在所有的樣本中進(jìn)行展示，這是由于在注釋的時候，發(fā)現(xiàn)不同類型的樣本直接有相同的target基因，這些基因在TSS的上下游有什么變化，可以用deeptools展示嗎？ 在CHIPseeker 包中的是可以的。

用deeptools 展示相同的traget 的一部分基因可以采用以下策略結(jié)合CHIPseeker 和deeptools兩種工具，首先找共有的基因，這個地方先的用CHIPseeker進(jìn)行注釋，對每個peak文件讀取，代碼如下：

library(ChIPseeker)
library(org.Mm.eg.db)
library(TxDb.Mmusculus.UCSC.mm10.knownGene)
library(VennDiagram)

cbx7 <- readPeakFile(“cbx7_peaks.narrowPeak.bed”)
ring1B <- readPeakFile(“ring1B_peaks.narrowPeak.bed”)
suz12 <- readPeakFile(“suz12_peaks.narrowPeak.bed”)

cbx7_peakAnno <- annotatePeak(cbx7, tssRegion = c(-2500, 2500), TxDb = TxDb.Mmusculus.UCSC.mm10.knownGene, annoDb=“org.Mm.eg.db”)
ring1B_peakAnno <- annotatePeak(ring1B, tssRegion = c(-2500, 2500), TxDb = TxDb.Mmusculus.UCSC.mm10.knownGene, annoDb=“org.Mm.eg.db”)
suz12_peakAnno <- annotatePeak(suz12, tssRegion = c(-2500, 2500), TxDb = TxDb.Mmusculus.UCSC.mm10.knownGene, annoDb=“org.Mm.eg.db”)

common_gene <- Reduce(intersect, list(cbx7=as.data.frame(cbx7_peakAnno)$geneId,ring1B=as.data.frame(ring1B_{p}eakAnno)$geneId,
suz12=as.data.frame(suz12_peakAnno)$geneId))

cbx7_common <- cbx7_df[cbx7_df$geneId %in% common_gene, 1:12][,1:3] ## 這些common gene 在cbx7中的 peak位置，基因只有2955個，但是這些peak位置很多。
head(cbx7_common)
seqnames start end
1 chr1 3062894 3063062
2 chr1 3671191 3671347
3 chr1 3671842 3671960
4 chr1 4491542 4493590
5 chr1 4493643 4493888
6 chr1 4494395 4494544

write.table(cbx7_common,file=“cbx7_commen.bed”,row.names = F,quote = F,col.names = F,sep="\t")
將上面保存的bed文件用于計算矩陣，用deeptools，如下：
computeMatrix reference-point -p 3 --referencePoint TSS -b 2500 -a 2500 -S Cbx7.bw -R cbx7_commen.bed --skipZeros -o
commongene_cbx7.mat.gz --outFileNameMatrix commengenecbx7
plotHeatmap -m commongene_cbx7.mat.gz -out commongene.png

這樣就能將cbx7上的commeon gene 做出來heatmap

上邊左邊是2500多個共有的基因的peak, 右邊的是所有基因的peak做的

chipseek 注釋peak

library(ChIPseeker)
library(GenomicFeatures)
library(GenomicRanges)
library(rtracklayer)
gffRangdata<-import.gff("~/src/GSE42466/data/chip-se/reference/ensembl.mm10.gff3") #獲取本地下載的參考基因組的gff3文件，自己生成txdb
myRanges<-as(gffRangdata,“GRanges”) ##獲得GRanges對象
txdb<-makeTxDbFromGRanges(myRanges)
cbx7<-readPeakFile("~/ncbi/public/sra/igv_check/cbx7_peaks.narrowPeak.bed")#這個地方是macs2產(chǎn)生的peak文件，但最好是改名加bed后綴，不然在R中轉(zhuǎn)化為數(shù)據(jù)框的時候，第一列的抬頭有錯位，目前還不知道是什么原因。
peak_cbx7<-annotatePeak(cbx7,tssRegion = c(-2500,2500),TxDb = txdb,addFlankGeneInfo = T,flankDistance = 5000)
Long coding RNA 注釋

##載入lincRNA注釋文件
library(ChIPseeker)
#source(“https://bioconductor.org/biocLite.R”)
#biocLite(“TxDb.Hsapiens.UCSC.hg19.lincRNAsTranscripts”)
library(“TxDb.Hsapiens.UCSC.hg19.lincRNAsTranscripts”)
lincRNA_txdb=TxDb.Hsapiens.UCSC.hg19.lincRNAsTranscripts

##讀入bed文件
HSC <- readPeakFile(“hsc478_summits.bed”)
HSC_lincRNA <- annotatePeak(HSC, TxDb=lincRNA_txdb)

visualazition

plotAnnoBar(HSC_lincRNA)

Visualize Genomic Annotation

upsetplot(HSC_lincRNA, vennpie=TRUE)
參考文章

https://www.jianshu.com/p/af3e492a84a9

https://mp.weixin.qq.com/s/_A0rHldzEgVk7bgwt457qQ

https://www.jianshu.com/p/a7b6ce208f98

http://www.bioinfo-scrounger.com/archives/365

https://mp.weixin.qq.com/s/vWTf59KDs1lp_sQXjEhI_g
————————————————
版權(quán)聲明：本文為CSDN博主「samhuairen」的原創(chuàng)文章，遵循CC 4.0 BY-SA版權(quán)協(xié)議，轉(zhuǎn)載請附上原文出處鏈接及本聲明。
原文鏈接：https://blog.csdn.net/samhuairen/article/details/88718431

總結(jié)

以上是生活随笔為你收集整理的CHIP-seq 分析笔记的全部內(nèi)容，希望文章能夠幫你解決所遇到的問題。

如果覺得生活随笔網(wǎng)站內(nèi)容還不錯，歡迎將生活随笔推薦給好友。