這是CUT&Tag數(shù)據(jù)分析的最后一部分，這一部分基本就是進(jìn)行各種的可視化了报账。
第七節(jié) 可視化
通常研底，我們感興趣的是使用基因組瀏覽器在各個(gè)區(qū)域可視化染色質(zhì)landscape。Integrative Genomic Viewer提供了一個(gè)易于使用的網(wǎng)頁(yè)應(yīng)用程序版本和本地桌面版本透罢。UCSC Genome Browser提供了最全面的補(bǔ)充基因組信息榜晦。
（一）查看標(biāo)準(zhǔn)化的bedgraph文件
根據(jù)官網(wǎng)展示的染色體位置區(qū)域，我們可以使用IGV也查看一下：hg38 羽圃，chr7:131,000,000-134,000,000
下面是我運(yùn)行后得到的IGV峰圖：
再看一下官網(wǎng)的圖：
基本上是一樣的乾胶。
（二）熱圖可視化特異區(qū)域
我們還感興趣的是在一系列的注釋位點(diǎn)上觀察染色質(zhì)特征，例如基因啟動(dòng)子上的組蛋白修飾信號(hào)朽寞。我們將使用deepTools的computeMatrix和plotHeatmap函數(shù)來(lái)生成熱圖识窿。
首先，我們要先把bam文件sort脑融，生成其索引喻频，在生成bw文件：
#!/bin/bash
projPath="/gpfs/home/fangy04/cut_tag"
cat /gpfs/home/fangy04/cut_tag/filenames | while read i
do
samtools sort -o $projPath/bam_files/${i}.sorted.bam $projPath/bam_files/${i}_bowtie2.mapped.bam
samtools index $projPath/bam_files/${i}.sorted.bam
bamCoverage -b $projPath/bam_files/${i}.sorted.bam -o $projPath/bigwig/${i}_raw.bw
done
（1）Heatmap over transcription units
這一步首先你要先下載hg38的基因注釋文件，可以參考我之前的文章：ChIP-seq實(shí)踐（非轉(zhuǎn)錄因子吨掌，非組蛋白）進(jìn)行下載。然后運(yùn)行（這一步非常非常耗時(shí)以及耗內(nèi)存脓恕，我用了服務(wù)器上120G膜宋，運(yùn)行了9個(gè)小時(shí)。炼幔。秋茫。）：
#!/bin/bash
projPath="/gpfs/home/fangy04/cut_tag"
cores=8
computeMatrix scale-regions -S $projPath/bigwig/SH_Hs_K27m3_rep1_raw.bw \
$projPath/bigwig/SH_Hs_K27m3_rep2_raw.bw \
$projPath/bigwig/SH_Hs_K4m3_rep1_raw.bw \
$projPath/bigwig/SH_Hs_K4m3_rep2_raw.bw \
-R $projPath/data/hg38_genes.bed \ #這是我們下載的基因組注釋文件，bed格式的
--beforeRegionStartLength 3000 \
--regionBodyLength 5000 \
--afterRegionStartLength 3000 \
--skipZeros -o $projPath/data/matrix_gene.mat.gz -p $cores
plotHeatmap -m $projPath/data/matrix_gene.mat.gz -out $projPath/data/Histone_gene.png --sortUsing sum --missingDataColor "#cc2200"
這里你可以看到上面的plotheatmap代碼里乃秀，我使用了--missingDataColor "#cc2200"這個(gè)參數(shù)肛著，是因?yàn)槟阋付ó媹D的時(shí)候圆兵，你的missing data的顏色，如果你不指定枢贿，并且在computeMatrix代碼里也沒(méi)有指定--missingDataAsZero殉农，那么你畫出來(lái)的圖很有可能是這樣的：
運(yùn)行完成后生成兩個(gè)文件：
-rw------- 1 fangy04 fangy04 2069525 Dec 25 03:08 Histone_gene.png
-rw------- 1 fangy04 fangy04 134316380 Dec 25 03:04 matrix_gene.mat.gz
（2）在CUT&Tag peaks上的熱圖
我們使用從SEACR返回的信號(hào)塊的中點(diǎn)來(lái)對(duì)齊熱圖中的信號(hào)。SEACR輸出的第六列是一個(gè)chr:start-end形式的條目局荚，表示該區(qū)域信號(hào)最大的第一個(gè)和最后一個(gè)堿基的位置超凳。我們首先在第6列中生成一個(gè)包含中點(diǎn)信息的bed文件，并使用deeptools進(jìn)行熱圖可視化耀态。
#!/bin/bash
projPath="/gpfs/home/fangy04/cut_tag"
cores=12
cat /gpfs/home/fangy04/cut_tag/filenames2 | while read i
do
awk '{split($6, summit, ":"); split(summit[2], region, "-"); print summit[1]"\t"region[1]"\t"region[2]}' $projPath/peak_calling/${i}_seacr_control.peaks.stringent.bed > $projPath/peak_calling/${i}_seacr_control.peaks.summitRegion.bed
computeMatrix reference-point -S $projPath/bigwig/${i}_raw.bw \
-R $projPath/peak_calling/${i}_seacr_control.peaks.summitRegion.bed \
--skipZeros -o $projPath/data/${i}_SEACR.mat.gz -p $cores \
-a 3000 -b 3000 --referencePoint center \
--missingDataAsZero #這次我加了這個(gè)參數(shù)轮傍，因?yàn)榈谝淮伟凑展倬W(wǎng)的代碼沒(méi)有加這個(gè)參數(shù)，畫出來(lái)的熱圖里很多黑色的線首装，很難看
plotHeatmap -m $projPath/data/${i}_SEACR.mat.gz -out $projPath/data/${i}_SEACR_heatmap.png --sortUsing sum --startLabel "Peak Start" --endLabel "Peak End" --xAxisLabel "" --regionsLabel "Peaks" --samplesLabel "${i}"
done
第八節(jié) 差異分析
評(píng)估來(lái)自高通量測(cè)序分析count數(shù)據(jù)的方差-均值依賴性和基于使用負(fù)二項(xiàng)分布的模型的差異表達(dá)測(cè)試创夜。
（一）創(chuàng)建peak x sample矩陣
通常，差異檢測(cè)比較相同組蛋白修飾的兩種或多種條件仙逻。在本教程中驰吓，受演示數(shù)據(jù)的限制，我們將通過(guò)比較兩個(gè)重復(fù)的H3K27me3和兩個(gè)重復(fù)的H3K4me3來(lái)說(shuō)明差異檢測(cè)桨醋。我們將使用DESeq2 (complete tutorial)作為說(shuō)明棚瘟。
（1）創(chuàng)建主要peak列表，合并每個(gè)樣品的peaks
> library(GenomicRanges)
> library(magrittr)
> mPeak = GRanges()
## overlap with bam file to get count
> histL = c("SH_Hs_K27m3", "SH_Hs_K4m3")
> repL = c("rep1","rep2")
> for(hist in histL){
for(rep in repL){
peakRes = read.table(paste0(hist, "_", rep, "_seacr_control.peaks.stringent.bed"), header = FALSE, fill = TRUE)
mPeak = GRanges(seqnames = peakRes$V1, IRanges(start = peakRes$V2, end = peakRes$V3), strand = "*") %>% append(mPeak, .)
}
}
> masterPeak = reduce(mPeak)
（2）在主要peak列表里獲取每個(gè)peak的片段counts
> library(DESeq2)
# 這里官網(wǎng)的代碼仍然是有問(wèn)題的喜最，不能直接用偎蘸，需要自己加for語(yǔ)句
> for (hist in histL){
for(rep in repL){
countMat = matrix(NA, length(masterPeak), length(histL)*length(repL))
}
}
## overlap with bam file to get count
> library(chromVAR)
> i = 1
> for(hist in histL){
for(rep in repL){
bamFile = paste0(hist, "_", rep, "_bowtie2.mapped.bam")
fragment_counts <- getCounts(bamFile, masterPeak, paired = TRUE, by_rg = FALSE, format = "bam")
countMat[, i] = counts(fragment_counts)[,1]
i = i + 1
}
}
> colnames(countMat) = paste(rep(histL, 2), rep(repL, each = 2), sep = "_")
（二）測(cè)序深度標(biāo)準(zhǔn)化和差異富集peaks檢測(cè)
> selectR = which(rowSums(countMat) > 5) ## remove low count genes
> dataS = countMat[selectR,]
> condition = factor(rep(histL, each = length(repL)))
> dds = DESeqDataSetFromMatrix(countData = dataS,
colData = DataFrame(condition),
design = ~ condition)
> DDS = DESeq(dds)
> normDDS = counts(DDS, normalized = TRUE) ## normalization with respect to the sequencing depth
> colnames(normDDS) = paste0(colnames(normDDS), "_norm")
> res = results(DDS, independentFiltering = FALSE, altHypothesis = "greaterAbs")
> countMatDiff = cbind(dataS, normDDS, res)
> head(countMatDiff)
DataFrame with 6 rows and 14 columns
SH_Hs_K27m3_rep1 SH_Hs_K4m3_rep1 SH_Hs_K27m3_rep2 SH_Hs_K4m3_rep2 SH_Hs_K27m3_rep1_norm SH_Hs_K4m3_rep1_norm
<numeric> <numeric> <numeric> <numeric> <numeric> <numeric>
1 3 4 5 4 0.748764 1.248831
2 0 0 298 198 0.000000 0.000000
3 0 1 173 91 0.000000 0.312208
4 0 0 266 206 0.000000 0.000000
5 4 3 0 2 0.998352 0.936623
6 4 2 0 1 0.998352 0.624415
SH_Hs_K27m3_rep2_norm SH_Hs_K4m3_rep2_norm baseMean log2FoldChange lfcSE stat pvalue padj
<numeric> <numeric> <numeric> <numeric> <numeric> <numeric> <numeric> <numeric>
1 19.9532 11.66190 8.40317 3.97867 1.06773 3.726287 1.94321e-04 3.08812e-03
2 1189.2085 577.26424 441.61819 14.06841 1.55908 9.023519 1.82142e-19 1.31189e-17
3 690.3795 265.30831 238.99999 11.73399 1.54086 7.615238 2.63205e-14 5.63249e-13
4 1061.5083 600.58805 415.52408 13.98064 1.54744 9.034672 1.64495e-19 1.20258e-17
5 0.0000 5.83095 1.94148 1.70735 1.82556 0.935247 3.49661e-01 9.53197e-01
6 0.0000 2.91548 1.13456 1.02468 2.28264 0.448903 6.53502e-01 9.53197e-01
（1）DESeq2要求輸入矩陣應(yīng)該是非標(biāo)準(zhǔn)化counts或測(cè)序reads的估計(jì)counts。
（2）DESeq2模型在內(nèi)部修正文庫(kù)大小瞬内。
（3）countMatDiff總結(jié)了差異分析結(jié)果:
前4列:在過(guò)濾低counts的峰區(qū)域后迷雪，原始reads計(jì)數(shù)
第二個(gè)4列:消除文庫(kù)大小差異的標(biāo)準(zhǔn)化reads計(jì)數(shù)。
剩余列:差示檢測(cè)結(jié)果虫蝶。
