CHIP-seq 分析笔记
本周学习一下CHIP-seq。 并根据网上的教程，自己实践一下， 一方面是要为了弄清楚什么是chip-seq, 这个技术有什么用，另一个方面是想学习一下这个技术如何来实践， 本文参考的文章主要来自生信技能树，以及简书上的其他作者写的教程，由于每个人在做分析时，使用的操作系统不一样，所以网上的代码在自己的电脑上进行运行的时候经常出现问题，这需要每个人针对自己的情况进行分析和总结。 本次分析采用的系统是MACos 10.12.6。
先不谈CHIP-seq 的原理， 因为我正在学习中，后期将这部分内容补上。先根据文献提供的数据和网上的代码进行实际操作。
CHIP-seq能干什么？
明确每一类组蛋白或者转录因子在整个基因组上结合基因的位置
如果比较多个组蛋白在亚基，可以看这些亚基之间在基因组上结合的基因的包含关系，即用韦恩图展示这些组蛋白结合基因相互之间是否包含。
检查每一类组蛋白结合基因在TSS上的位置。
检查每一组（不同组蛋白之间结合相同的基因）在TSS上的位置。（这样可以看出缺少某一类组蛋白之后，基因是否表达，验证这个组蛋白具有的功能和意义）
不同组蛋白结合基因的功能（GO），及参与的代谢通路（KEGG）
可以研究每一个组蛋白targets 的基因的表达
第一步：从NCBI下载数据，并解压到本地电脑。 从NCBI的GEO dataset 中输入作者上传的GEO号码GSE42466,如下图所示：
在本文文章的哪里打对勾，并进入，之后看一看到文章的具体信息，包括作者的信息，以及实验方法，实验设计，实验上传的具体数据编号。如下所示：
在最后的一栏中找到本数据的SRA号码，点击进入，如下：
在上面作者提供的6个数据上打上对勾，并在右上侧的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 软件进行解压，如下：
ls *sra |while read id; do fastq-dump –split-3 $id;done # --split-3的目的是如果文件包含两个以上文件，则分别命名，如果是一个文件则直接是文件原名，而不是根据reads 来进行拆分
以上可以获得本实验的的所有数据，但是在实际操作中，SRR620207一致下载不下来，原因不知道。 我暂且不比对此文件，在后期使用作者提供的文件。
第二步： 对原始数据进行质控，采用fastqc, multiqc 两个软件， fastqc对每个文件进行质控分析， multiqc对fastqc的结果进行整合，方便读者从总体上对数据质量进行把控。
ls *.fastq |while read id; do fastqc -t 3 $id; done
multiqc *fastqc.zip --pdf # multiqc 在我的电脑山无法生成PDF，虽然我按照说明书安装了pandoc, latex. 但依然不行.
根据multiqc的质控显示，本测序数据在3’端的数据质量不高，fastqc 也显示警告，原因是平均的碱基质量有点低，如下图所示：
因此在接下来比对的时候，我们要去掉3‘端质量不好的碱基，用bowtie2自动去掉。
第三步： 利用bowtie2（本次分析采用的是bowtie2-2.3.5）对数据进行mapping， 数据索引下载如下：
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
比对后的总体结果如下：
SRR620204
19687027 reads; of these:
19687027 (100.00%) were unpaired; of these:
7960096 (40.43%) aligned 0 times
8614711 (43.76%) aligned exactly 1 time
3112220 (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 times
12169372 (58.76%) aligned exactly 1 time
5824408 (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 times
10534163 (48.88%) aligned exactly 1 time
3700860 (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 times
5235529 (39.39%) aligned exactly 1 time
2324724 (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 times
17025381 (54.54%) aligned exactly 1 time
8494196 (27.21%) aligned >1 times
81.74% overall alignment rate
#############
比对完之后，利用IGV进行可视化，先用samtools进行index, 之后在IGV中导入，导入的时候选择BAM文件即可，INDEX文件也必须在bam文件相同的文件夹。
samtools index FGENESH_C15+B10.sorted.bam FGENESH_C15+B10.sorted.bai
从IGV上确实可以发现， chip-seq 测序完在有target 的地方有峰值，没有target的地方对照和实验都没有，如下所示：
之后，我们使用macs2进行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
比对之后进行数据的作图和可视化，代码
ls bam |while read id;
do file=$(basename
i
d
)
;
s
a
m
p
l
e
=
id ); sample=
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 系统上出现了一个小插曲， 使用bamCoverage 的时候报错，原因如下：
可能的原因是没有找到相应的库，于是我重新安装了libssh2. 结果显示可以运行了。
此外在运行computeMatrix 的时候也报错，错误的原因是提供的bed文件有问题，于是，我有重新用gff文件制作了一个bed文件（注意：这个地方的gff文件必须和参考基因组是一直的，必须是针对这个参考基因组的gff3文件，在R中如果要用CHIPSeeker包进行注释的时候，当上面没有注释的数据库时，需要用GFF3文件来自己制作一个注释数据库，这个时候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
至此，所有要画图的文件都准备好了，用plotHeatmap进行画图，代码如下： 在这里我只画了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上下游画的，那可不可以对某一部分感兴趣的基因在所有的样本中进行展示，这是由于在注释的时候，发现不同类型的样本直接有相同的target基因，这些基因在TSS的上下游有什么变化，可以用deeptools展示吗？ 在CHIPseeker 包中的是可以的。
用deeptools 展示相同的traget 的一部分基因可以采用以下策略结合CHIPseeker 和deeptools两种工具，首先找共有的基因，这个地方先的用CHIPseeker进行注释，对每个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)
g
e
n
e
I
d
,
r
i
n
g
1
B
=
a
s
.
d
a
t
a
.
f
r
a
m
e
(
r
i
n
g
1
B
p
e
a
k
A
n
n
o
)
geneId, ring1B=as.data.frame(ring1B_peakAnno)
geneId,ring1B=as.data.frame(ring1BpeakAnno)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产生的peak文件，但最好是改名加bed后缀，不然在R中转化为数据框的时候，第一列的抬头有错位，目前还不知道是什么原因。
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
————————————————
版权声明：本文为CSDN博主「samhuairen」的原创文章，遵循CC 4.0 BY-SA版权协议，转载请附上原文出处链接及本声明。
原文链接：https://blog.csdn.net/samhuairen/article/details/88718431