CUT&Tag數(shù)據(jù)分析筆記（3）

這是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)的圖很有可能是這樣的：

參考：https://www.biostars.org/p/322414/

運(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é)果。

?著作權(quán)歸作者所有,轉(zhuǎn)載或內(nèi)容合作請(qǐng)聯(lián)系作者
【社區(qū)內(nèi)容提示】社區(qū)部分內(nèi)容疑似由AI輔助生成，瀏覽時(shí)請(qǐng)結(jié)合常識(shí)與多方信息審慎甄別。
平臺(tái)聲明：文章內(nèi)容（如有圖片或視頻亦包括在內(nèi)）由作者上傳并發(fā)布，文章內(nèi)容僅代表作者本人觀點(diǎn)，簡(jiǎn)書系信息發(fā)布平臺(tái)，僅提供信息存儲(chǔ)服務(wù)。
禁止轉(zhuǎn)載，如需轉(zhuǎn)載請(qǐng)通過(guò)簡(jiǎn)信或評(píng)論聯(lián)系作者。

相關(guān)閱讀更多精彩內(nèi)容

友情鏈接更多精彩內(nèi)容