|
轉(zhuǎn)自 簡書 生信start_site 在實(shí)戰(zhàn)之前，首先對CHIP-seq分析做一些了解，以下是從各個(gè)地方copy過來綜合起來的，有些散亂，是我認(rèn)為重要以及應(yīng)該了解的知識點(diǎn)。
chip-seq 實(shí)戰(zhàn)就跟在轉(zhuǎn)錄組和外顯子上的實(shí)戰(zhàn)是沒有本質(zhì)區(qū)別的，只是它多了一些分析，轉(zhuǎn)錄組只找表達(dá)量做差異分析，外顯子只找變異做變異注釋，那chip-seq 找的是peaks和motif 。
這里biobabble公眾號有比較完整的chip-seq的介紹：
CS1: ChIPseq簡介
CS2: BED文件
CS3: peak注釋
CS4：關(guān)于ChIPseq注釋的幾個(gè)問題
CS5: 吃著火鍋，唱著歌，還把分析給做了
CS6: ChIPseeker的可視化方法
CS7：Genomic coordination的富集性分析（1）
CS8：Genomic coordination的富集性分析（2）
CS9: GEO數(shù)據(jù)挖掘
CS番外1：如何定義下游的區(qū)間？ ChIP-seq 數(shù)據(jù)分析流程
1.得到原始序列文件后，第一步是執(zhí)行標(biāo)準(zhǔn)的高通量數(shù)據(jù)質(zhì)量控制。
2.原始的 reads 回比到研究物種的參考基因組上
3.特異性 ChIP-seq 的質(zhì)量控制，檢查 ChIP-seq 樣品中的富集情況，并排除過度碎片
4.【核心】鑒定全基因上感興趣的因子富集的區(qū)域，也叫 Peak calling
5.一旦確定了峰值區(qū)域，還可以評估需要峰值位置的特定 ChIP-seq 的質(zhì)量控制措施
6.在 Peak calling 之后，要確保預(yù)測的峰值準(zhǔn)確地捕捉到binding pattern。
7.Peak calling后的分析取決于實(shí)驗(yàn)設(shè)計(jì)與生物學(xué)問題。比如，如果在實(shí)驗(yàn)中有重復(fù)和/或多個(gè)條件，那么下一步可能是樣本之間的比較（見第 8 章）。生物學(xué)重復(fù)將提供有關(guān)數(shù)據(jù)(噪聲)中的重現(xiàn)性和內(nèi)在生物和技術(shù)變異性的信息。差異結(jié)合分析解決了哪一Peak區(qū)域在兩種情況下顯示出明顯不同的豐度（例如：與觀察到的相同條件的重復(fù)之間的變化相比，差異比預(yù)期的要大得多）。
8.Peak區(qū)間或差異Peak做下游分析，比如基因組注釋、GO分析、Pathway 分析、motif 查找、與其他基因組數(shù)據(jù)聯(lián)合分析。（參考文章：https://www.jianshu.com/p/19b178ea406e） 
image.png 怎樣識別富集區(qū)域？(https://www.jianshu.com/p/df33bcb68548)
從DNA測序的角度來說，因?yàn)闇y序都是5'端的reads，對于一個(gè)DNA序列來說（有正負(fù)鏈的），它mapping的位置正負(fù)鏈都有的（也就是紅色和藍(lán)色的reads都有），對這些reads位置進(jìn)行統(tǒng)計(jì)畫圖可以看到一個(gè)紅色的peak，一個(gè)藍(lán)色的peak。這兩個(gè)peak說明的是一個(gè)事情，就是這個(gè)地方有富集。最后對這兩個(gè)peak進(jìn)行merge，最后變成了一個(gè)富集區(qū)域?；疑膒eak! 
image ChIP-seq的常見峰型
不同類型的蛋白或者組蛋白修飾會(huì)得到不同的峰形：
sharp binding sites, CTCF (red);
a mixture of shapes, RNA Polymerase II (orange);
medium size, H3K36me3 (green);
large domains, H3K27me3 (blue) 
image 富集倍數(shù)
一般來說富集倍數(shù)要在5以上才算是顯著，如何計(jì)算富集比率呢？看圖~ 
image Chip-Seq實(shí)踐需要用的軟件：
fastqc, sra-tools, bowtie2, samtools, MACS2, deeptools
上述軟件之前在RNA-seq實(shí)踐中已經(jīng)安裝了大部分，需要安裝MACS2和deeptools。
MACS2比較特殊，因?yàn)樗枰狿ython2.7的環(huán)境才能安裝，而我的電腦安裝的是python3環(huán)境。所以首先創(chuàng)建一個(gè)虛擬環(huán)境，在python2.7的虛擬環(huán)境中使用（https://mp.weixin.qq.com/s/mEwl7-f2WTNfmK-FkNkP0A）： conda create -n py27 python=2
source activate py27
pip install MACS2
本次實(shí)踐的文章是：
RYBP and Cbx7 define specific biological functions of polycomb complexes in mouse embryonic stem cells
（探索PRC1，PCR2蛋白復(fù)合物，不是轉(zhuǎn)錄因子或者組蛋白的CHIP-seq）
本文主要講了PRC1(Polycomb repressive complex 1)在小鼠的胚胎干細(xì)胞中有兩類亞型Cbx7-PRC1和RYBP-PRC1，但是Cbx7和RYBP這兩種是不能共存的，盡管兩者的全基因組定位在某些gene上交叉，也就是說在基因組上結(jié)合在了同一個(gè)位置。在分子水平上，Cbx7用于招募Ring1B結(jié)合到染色質(zhì)上，然而RYBP可以增強(qiáng)PRC1的酶活性。RYBP結(jié)合的基因，其有著了較低水平的Ring1B和H2AK119ub，但其相比結(jié)合Cbx7有較高的表達(dá)量。在功能上，RYBP結(jié)合的基因主要涉及代謝調(diào)節(jié)和細(xì)胞周期，Cbx7結(jié)合的基因主要涉及early-lineage commitment of ESCs。（這段來源：http://www./archives/358） 1.下載原始數(shù)據(jù)：
可以寫一個(gè)腳本，批量下載： #!/bin/bash
for ((i=204;i<=209;i++))
do wget https://sra-download.ncbi.nlm./traces/sra6/SRR/000605/SRR620$i
done
2.下載之后提取fastq文件： #!/bin/bash
for i in SRR62020*
do
echo $i
fastq-dump --gzip --split-files $i
done
3.fastqc檢查測序質(zhì)量
fastqc發(fā)現(xiàn)Ring1B、Suz12,IgG_old，IgG樣本的3’端有大約5個(gè)bp長度的堿基質(zhì)量不是太好，可以在之后的bowtie2比對的時(shí)候設(shè)置參數(shù)，把質(zhì)量差的幾個(gè)堿基去掉。 4.bowtie2比對
首先下載bowtie2的mm10基因組索引： wget ftp://ftp.ccb.jhu.edu/pub/data/bowtie2_indexes/mm10.zip
還需要mm10 refseq注釋bed文件，進(jìn)行網(wǎng)址下載http://genome./cgi-bin/hgTables，其中track和table選擇RefSeq，然后輸出格式為bed即可，最后下載（可以先在自己電腦上下載再傳到服務(wù)器上，也可以像下面直接在服務(wù)器上下載）（參考http://www./archives/358） curl 'http://genome./cgi-bin/hgTables?hgsid=646311755_P0RcOBvAQnWZSzQz2fQfBiPPSBen&boolshad.hgta_printCustomTrackHeaders=0&hgta_ctName=tb_ncbiRefSeq&hgta_ctDesc=table+browser+query+on+ncbiRefSeq&hgta_ctVis=pack&hgta_ctUrl=&fbQual=whole&fbUpBases=200&fbExonBases=0&fbIntronBases=0&fbDownBases=200&hgta_doGetBed=get+BED' >mm10.refseq.bed
用bowtie2進(jìn)行比對（用腳本進(jìn)行）： #!/bin/bash
bowtie2 -p 6 -3 5 --local -x /media/yanfang/FYWD/CHIPSEQ/bowtie2mm10/mm10 -U /media/yanfang/FYWD/CHIPSEQ/fqfile/SRR620204_1.fastq.gz | samtools sort -O bam -o ring1B.bam
bowtie2 -p 6 -3 5 --local -x /media/yanfang/FYWD/CHIPSEQ/bowtie2mm10/mm10 -U /media/yanfang/FYWD/CHIPSEQ/fqfile/SRR620205_1.fastq.gz | samtools sort -O bam -o cbx7.bam
bowtie2 -p 6 -3 5 --local -x /media/yanfang/FYWD/CHIPSEQ/bowtie2mm10/mm10 -U /media/yanfang/FYWD/CHIPSEQ/fqfile/SRR620206_1.fastq.gz | samtools sort -O bam -o Suz12.bam
bowtie2 -p 6 -3 5 --local -x /media/yanfang/FYWD/CHIPSEQ/bowtie2mm10/mm10 -U /media/yanfang/FYWD/CHIPSEQ/fqfile/SRR620207_1.fastq.gz | samtools sort -O bam -o RYBP.bam
bowtie2 -p 6 -3 5 --local -x /media/yanfang/FYWD/CHIPSEQ/bowtie2mm10/mm10 -U /media/yanfang/FYWD/CHIPSEQ/fqfile/SRR620208_1.fastq.gz | samtools sort -O bam -o IgGold.bam
bowtie2 -p 6 -3 5 --local -x /media/yanfang/FYWD/CHIPSEQ/bowtie2mm10/mm10 -U /media/yanfang/FYWD/CHIPSEQ/fqfile/SRR620209_1.fastq.gz | samtools sort -O bam -o IgG.bam
~
但是貌似對比率不高： 19687027 reads; of these:
19687027 (100.00%) were unpaired; of these:
7778437 (39.51%) aligned 0 times
7766495 (39.45%) aligned exactly 1 time
4142095 (21.04%) aligned >1 times
60.49% overall alignment rate
[bam_sort_core] merging from 4 files and 1 in-memory blocks...
20710168 reads; of these:
20710168 (100.00%) were unpaired; of these:
2584870 (12.48%) aligned 0 times
10914615 (52.70%) aligned exactly 1 time
7210683 (34.82%) aligned >1 times
87.52% overall alignment rate
[bam_sort_core] merging from 4 files and 1 in-memory blocks...
21551927 reads; of these:
21551927 (100.00%) were unpaired; of these:
7109852 (32.99%) aligned 0 times
9455003 (43.87%) aligned exactly 1 time
4987072 (23.14%) aligned >1 times
67.01% overall alignment rate
[bam_sort_core] merging from 4 files and 1 in-memory blocks...
39870646 reads; of these:
39870646 (100.00%) were unpaired; of these:
6312903 (15.83%) aligned 0 times
20702502 (51.92%) aligned exactly 1 time
12855241 (32.24%) aligned >1 times
84.17% overall alignment rate
[bam_sort_core] merging from 8 files and 1 in-memory blocks...
13291429 reads; of these:
13291429 (100.00%) were unpaired; of these:
5608796 (42.20%) aligned 0 times
4663584 (35.09%) aligned exactly 1 time
3019049 (22.71%) aligned >1 times
57.80% overall alignment rate
[bam_sort_core] merging from 2 files and 1 in-memory blocks...
31218866 reads; of these:
31218866 (100.00%) were unpaired; of these:
5370101 (17.20%) aligned 0 times
15180536 (48.63%) aligned exactly 1 time
10668229 (34.17%) aligned >1 times
82.80% overall alignment rate
[bam_sort_core] merging from 6 files and 1 in-memory blocks...
5.用MACS2進(jìn)行call peak:
首先要進(jìn)入python2.7的環(huán)境： source activate py27
然后寫一個(gè)腳本： #!/bin/bash
macs2 callpeak -c /media/yanfang/FYWD/CHIPSEQ/sambam/IgGold.bam -t /media/yanfang/FYWD/CHIPSEQ/sambam/Suz12.bam -m 10 30 -p 1e-5 -f BAM -g mm -n Suz12 2>Suz12.masc2.log
macs2 callpeak -c /media/yanfang/FYWD/CHIPSEQ/sambam/IgGold.bam -t /media/yanfang/FYWD/CHIPSEQ/sambam/ring1B.bam -m 10 30 -p 1e-5 -f BAM -g mm -n ring1B 2>ring1B.masc2.log
macs2 callpeak -c /media/yanfang/FYWD/CHIPSEQ/sambam/IgGold.bam -t /media/yanfang/FYWD/CHIPSEQ/sambam/cbx7.bam -m 10 30 -p 1e-5 -f BAM -g mm -n cxb7 2>cxb7.masc2.log
macs2 callpeak -c /media/yanfang/FYWD/CHIPSEQ/sambam/IgGold.bam -t /media/yanfang/FYWD/CHIPSEQ/sambam/RYBP.bam -m 10 30 -p 1e-5 -f BAM -g mm -n RYBP 2>RYBP.masc2.log
-t/–treatment 輸入文件，支持很多格式，搭配-f/–format使用
-f/–format 設(shè)定輸入文件的格式，這里選擇bam格式
-c/–control 對照組（或者模擬的）數(shù)據(jù)，需要跟-t的輸出文件在相同目錄下
-n/–name 輸出文件（有很多個(gè)文件）的前綴
–outdir 輸出文件所在目錄（不設(shè)定的話就默認(rèn)當(dāng)前目錄了）
-g/–gsize 提供基因組的大小，程序有默認(rèn)的幾個(gè)物種可以選hs,mm,ce,dm
-q/–qvalue 設(shè)定FDR的閾值，默認(rèn)是0.05
-p/–pvalue 設(shè)定pvalue的閾值，如果參數(shù)設(shè)定了-p，則其會(huì)覆蓋參數(shù)-q的作用
-B/–bdg If this flag is on, MACS will store the fragment pileup, control lambda, -log10pvalue and -log10qvalue scores in bedGraph files. 參考文章（http://www./archives/363） 運(yùn)行后得到的文件不只是上述提到的一類，還有如下格式：
1.NAME_peaks.xls: 以表格形式存放peak信息，雖然后綴是xls，但其實(shí)能用文本編輯器打開，和bed格式類似，但是以1為基，而bed文件是以0為基.也就是說xls的坐標(biāo)都要減一才是bed文件的坐標(biāo)。
2.NAME_peaks.narrowPeak： NAME_peaks.broadPeak 類似。后面4列表示為， integer score for display， fold-change，-log10 pvalue，-log10 qvalue，relative summit position to peak start。內(nèi)容和NAME_peaks.xls基本一致，適合用于導(dǎo)入R進(jìn)行分析。
3.NAME_summits.bed：記錄每個(gè)peak的peak summits，話句話說就是記錄極值點(diǎn)的位置。MACS建議用該文件尋找結(jié)合位點(diǎn)的motif。
4.NAME_model.r： 能通過$ Rscript NAME_model.r作圖，得到是基于你提供數(shù)據(jù)的peak模型。 運(yùn)行之后，看一下每個(gè)樣本找到多少個(gè)peak（bed的格式，儲(chǔ)存了每個(gè)peaks的位置信息）: $ wc -l *.bed 2754 cxb7_summits.bed 7193 ring1B_summits.bed 473 RYBP_summits.bed 6696 Suz12_summits.bed
本文作者在GEO給的PEAKS個(gè)數(shù)如下：
2754 GSE42466_Cbx7_peaks_10.txt
6982 GSE42466_Ring1b_peaks_10.txt
6872 GSE42466_RYBP_peaks_5.txt
8054 GSE42466_Suz12_peaks_10.txt
然而我查了幾篇對該文獻(xiàn)重復(fù)的文章，call peak的數(shù)目和文章也有較大差距。（https://cloud.tencent.com/developer/article/1055109
http://www./archives/363
https://mp.weixin.qq.com/s/_A0rHldzEgVk7bgwt457qQ?
https://mp.weixin.qq.com/s/mEwl7-f2WTNfmK-FkNkP0A） 6.結(jié)果注釋與可視化
可視化這里我查閱了很多文章，分成兩種，有些人用的是chpseeker包進(jìn)行peak的可視化和注釋；而還有一些人用了deeptool軟件進(jìn)行分析。我兩種都試了試。 首先我用的是Chipseeker包。
ChIPseeker的功能分為三類(https://www.jianshu.com/p/c83a38915afc):
(1) 注釋：提取peak附近最近的基因， 注釋peak所在區(qū)域
(2)比較：估計(jì)ChIP peak數(shù)據(jù)集中重疊部分的顯著性；整合GEO數(shù)據(jù)集，以便于將當(dāng)前結(jié)果和已知結(jié)果比較
(3)可視化： peak的覆蓋情況；TSS區(qū)域結(jié)合的peak的平均表達(dá)譜和熱圖；基因組注釋；TSS距離；peak和基因的重疊。 Chip peaks coverage plot #調(diào)用要用的包
library("ChIPseeker")
library("org.Mm.eg.db")
library("TxDb.Mmusculus.UCSC.mm10.knownGene")
library("clusterProfiler")
#定位文件位置
setwd("/media/yanfang/FYWD/CHIPSEQ/bed")
#讀入bed文件
cbx7 <- readPeakFile("cxb7_peaks.narrowPeak.bed")
ring1B <- readPeakFile("ring1B_peaks.narrowPeak.bed")
Suz12 <- readPeakFile("Suz12_peaks.narrowPeak.bed")
#查看不同的蛋白的peak在染色體上的分布，因?yàn)镽YBP的peak幾乎沒有，所以沒有用這個(gè)蛋白做后續(xù)的畫圖處理
covplot(cbx7, weightCol="V5")
covplot(ring1B, weightCol="V5")
covplot(Suz12, weightCol="V5")
例：cbx7的peak在染色體上的分布 #查看在指定的染色體上的分布
covplot(cbx7,chrs=c("chr17", "chr18"))
cbx7在第17和18號染色體上的分布 熱圖（Heatmap of ChIP binding to TSS regions） # 定義TSS上下游的距離
> promoter <- promoters(TxDb.Mmusculus.UCSC.mm10.knownGene, upstream = 3000, downstream = 3000)> list_peak <- list(cbx7, ring1B, Suz12)> tagMatrix <- lapply(list_peak, getTagMatrix, window=promoter)> tagHeatmap(tagMatrix, xlim=c(-3000, 3000), color="red")
image.png Average Profile of ChIP peaks binding to TSS region
用譜圖的形式來展示不同蛋白結(jié)合的強(qiáng)度： > library(org.Mm.eg.db)
> library(TxDb.Mmusculus.UCSC.mm10.knownGene)
> library(VennDiagram)
> setwd("/media/yanfang/FYWD/CHIPSEQ/bed")
> cbx7 <- readPeakFile("cxb7_peaks.narrowPeak.bed")> ring1B <- readPeakFile("ring1B_peaks.narrowPeak.bed")> Suz12 <- readPeakFile("Suz12_peaks.narrowPeak.bed")> peaks <- list(cbx7=cbx7,ring1B=ring1B,Suz12=Suz12)
> txdb <- TxDb.Mmusculus.UCSC.mm10.knownGene> promoter <- getPromoters(TxDb=txdb, upstream=3000, downstream=3000)> tagMatrixList <- lapply(peaks, getTagMatrix, windows=promoter)> plotAvgProf(tagMatrixList, xlim=c(-3000, 3000))
image.png 還支持使用bootstrap估計(jì)置信區(qū)間： plotAvgProf(tagMatrixList, xlim=c(-3000, 3000), conf = 0.95, resample = 1000)
image.png 注釋peak
ChIPseeker負(fù)責(zé)注釋的就是annotatePeak #對于非向量的循環(huán)，推薦構(gòu)建list，然后用lapply> peaks <- list(ring1B, cbx7, Suz12)#注釋。啟動(dòng)子區(qū)域是沒有明確的定義的，所以你可能需要指定tssRegion，把基因起始轉(zhuǎn)錄位點(diǎn)的上下游區(qū)域來做為啟動(dòng)子區(qū)域。> peakAnnoList <- lapply(peaks, annotatePeak, tssRegion=c(-2500,2500), TxDb=TxDb.Mmusculus.UCSC.mm10.knownGene,
addFlankGeneInfo=TRUE, flankDistance=5000)#結(jié)果可以用as.GRanges或者as.data.frame查看#Granges是R語言中的數(shù)據(jù)結(jié)構(gòu)，它主要用來存儲(chǔ)基因組中的坐標(biāo)區(qū)間。> as.GRanges(peakAnnoList[[1]])> as.data.frame(peakAnnoList[[1]])#畫圖，不同蛋白的靶位點(diǎn)的注釋情況> plotDistToTSS(peakAnnoList)
可視化TSS區(qū)域的TF binding loci 韋恩圖 genes= lapply(peakAnnoList, function(i) as.data.frame(i)$geneId)vennplot(genes[1:2:])
韋恩圖 注釋的可視化： #先注釋
peakAnno <- annotatePeak(cbx7, tssRegion=c(-3000, 3000),
TxDb=txdb, annoDb="org.Mm.eg.db")
#柱形圖
plotAnnoBar(peakAnno)
#vennpie
vennpie(peakAnno)
#餅圖
plotAnnoPie(peakAnno)
柱形圖
vennpie
餅圖 # 展示注釋信息的交集$ upsetplot(peakAnno)
image.png 7.deeptools的使用
BAM文件是SAM的二進(jìn)制轉(zhuǎn)換版，應(yīng)該都知道。那么bigWig格式是什么？bigWig是wig或bedGraph的二進(jìn)制版，存放區(qū)間的坐標(biāo)軸信息和相關(guān)計(jì)分(score)，主要用于在基因組瀏覽器上查看數(shù)據(jù)的連續(xù)密度圖，可用wigToBigWig從wiggle進(jìn)行轉(zhuǎn)換。為什么要用bigWig? 主要是因?yàn)锽AM文件比較大，直接用于展示時(shí)對服務(wù)器要求較大。因此在GEO上僅會(huì)提供bw,即bigWig下載，便于下載和查看。
（一）將bam文件轉(zhuǎn)換成bw(bigwig)文件
這一步需要用到的是deeptool里的bamCoverage工具： #!/bin/bash
bamCoverage -b cbx7.bam -o cbx7.bw -p 10 --normalizeUsing RPKM
bamCoverage -b ring1B.bam -o ring1B.bw -p 10 --normalizeUsing RPKM
bamCoverage -b Suz12.bam -o Suz12.bw -p 10 --normalizeUsing RPKM
bamCoverage -b IgGold.bam -o IgGold.bw -p 10 --normalizeUsing RPKM
-b:指要進(jìn)行操作的bam文件，后面跟的是bam文件名。
-o:指輸出的文件，后面跟的是要輸出的文件名。
-p: --numberOfProcessors，進(jìn)程的數(shù)量。
--normalizeUsing: 標(biāo)準(zhǔn)化的方法。這個(gè)工具有RPKM, CPM, BPM, RPGC, None幾個(gè)不同的選擇，這里我選了RPKM。 （二）peak分布可視化
也是deeptool里的工具。為了統(tǒng)計(jì)全基因組范圍的peak在基因特征的分布情況，需要用到computeMatrix計(jì)算，用plotHeatmap以熱圖的方式對覆蓋進(jìn)行可視化，用plotProfile以折線圖的方式展示覆蓋情況。
computeMatrix有兩種模式，scale-regions mode和reference-point mode（參考文章：https://www.jianshu.com/p/2b386dd437d3）
reference-point(relative to a point): 計(jì)算某個(gè)點(diǎn)的信號豐度。給定一個(gè)bed file，以某個(gè)點(diǎn)為中心開始統(tǒng)計(jì)信號（TSS/TES/center）。
scale-regions(over a set of regions): 把所有基因組區(qū)段縮放至同樣大小，然后計(jì)算其信號豐度。簡單來說會(huì)將給定bed file范圍內(nèi)的結(jié)合信號做一個(gè)統(tǒng)計(jì)（指的是一段長度），并將基因長度統(tǒng)一scale到設(shè)定regionBdoyLength的長度，加上統(tǒng)計(jì)基因上游和下游Xbp的信號（beforeRegionStartLength參數(shù)和afterRegionStartLength參數(shù)）
不管是哪一種，有兩個(gè)參數(shù)是必須的，-S提供bigwig文件，-R提供基因的注釋信息。
image 如果想得到單獨(dú)的文件，就是說每個(gè)樣品都有單獨(dú)的matrix文件，可以用下面的代碼: #!/bin/bash
computeMatrix reference-point --referencePoint TSS -a 2500 -b 2500 -S /media/yanfang/FYWD/CHIPSEQ/bw/cbx7.bw -R /media/yanfang/FYWD/CHIPSEQ/bw/mm10.bed --skipZeros -o cbx7_matrix_TSS.gz
computeMatrix reference-point --referencePoint TSS -a 2500 -b 2500 -S /media/yanfang/FYWD/CHIPSEQ/bw/ring1B.bw -R /media/yanfang/FYWD/CHIPSEQ/bw/mm10.bed --skipZeros -o ring1B_matrix_TSS.gz
computeMatrix reference-point --referencePoint TSS -a 2500 -b 2500 -S /media/yanfang/FYWD/CHIPSEQ/bw/Suz12.bw -R /media/yanfang/FYWD/CHIPSEQ/bw/mm10.bed --skipZeros -o Suz12_matrix_TSS.gz
computeMatrix reference-point --referencePoint TSS -a 2500 -b 2500 -S /media/yanfang/FYWD/CHIPSEQ/bw/IgGold.bw -R /media/yanfang/FYWD/CHIPSEQ/bw/mm10.bed --skipZeros -o IgGold_matrix_TSS.gz
reference-point: 選擇模式，這里我選擇了第二種。
-p 10:指10個(gè)進(jìn)程。我的電腦比較渣，所以沒有加這個(gè)參數(shù)。
--referencePoint TSS：選擇參考點(diǎn)，這里選擇轉(zhuǎn)錄起始位點(diǎn)。
-a 3000:下游3000bp
-b 3000:上游3000bp
-S:指bw文件，后面跟的是文件的位置。
-R:參考基因組的注釋文件，這里要求是bed格式的。后面跟的是文件的位置。
--skipZeros：Whether regions with only scores of zero should be included or not. Default is to include them.
-o:輸出文件，后面跟的是輸出文件的名字 如果想把四個(gè)樣品合在一起，后續(xù)畫圖也是merge在一起的，可以用下面的代碼： $ computeMatrix reference-point --referencePoint TSS -a 2500 -b 2500 -S /media/yanfang/FYWD/CHIPSEQ/bw/*.bw -R /media/yanfang/FYWD/CHIPSEQ/bw/mm10.bed --skipZeros -o merge_matrix_TSS.gz
畫基因的TSS附近的profile和heatmap圖
image.png #heatmap圖$ plotHeatmap -m merge_matrix_TSS.gz -out merge_heatmap.png#將heatmap圖輸出為pdf格式$ plotHeatmap --dpi 720 -m merge_matrix_TSS.gz -out merge_heatmap.pdf --plotFileFormat pdf#輪廓圖plotProfile -m merge_matrix_TSS.gz -out merge_profile.png#將輪廓圖輸出為pdf格式$ plotProfile --dpi 720 -m merge_matrix_TSS.gz -out merge_profile.pdf --plotFileFormat pdf#如果想把每個(gè)樣品的輪廓圖merge上在一個(gè)圖片里，在代碼后面加上--perGroup$ plotProfile -m merge_matrix_TSS.gz -out merge_profile.png --perGroup
$ plotProfile --dpi 720 -m merge_matrix_TSS.gz -out merge_profile.pdf --plotFileFormat pdf --perGroup
熱圖
把幾個(gè)樣品的輪廓圖merge到同一張圖里 也可以整合所有的chipseq的bw文件，畫基因的genebody附近的profile和heatmap圖（圖就不放上來了）: $ computeMatrix scale-regions -a 3000 -b 3000 -m 5000 -S /media/yanfang/FYWD/CHIPSEQ/bw/*.bw -R /media/yanfang/FYWD/CHIPSEQ/bw/mm10.bed --skipZeros -o merge_scaleregion_matrix_TSS.gz
$ plotHeatmap -m merge_scaleregion_matrix_TSS.gz -out merge_scaleregion_matrix_TSS.png
$ plotProfile -m merge_scaleregion_profile_TSS.gz -out merge_scaleregion_profile_TSS.png
|
