CUT&Tag数据分析笔记（1）
2020-12-26 本文已影响0人
生信start_site
前几天学习了文献CUT&Tag，了解了原理以后，可以跟着官网来学习如何分析CUT&Tag数据了。官网地址：这里。需要说明的是：（1）笔记里的代码不完全和官网的一样，例如我的代码加了loop，用于批量处理（官网的代码只写了loop里的内容，而不是完整代码）。（2）再例如在比对的可视化部分，由于我的比对txt文件的内容与官网展示的不同，所涉及到的行数和列数也就有所不同，请阅读的同学注意这一点。根据自己的电脑运行出来的文件为准，适当调整。（3）就是官网的代码也出现了错误，有的地方缺少诸如“>”这样的符号导致代码报错。
NOTE：再次强调，即便是官网的代码，也要自己亲自运行一遍，你不会知道这些代码在你的电脑上是否报错。
第一节 前言
（一）CUT&Tag概述
发生在真核细胞核DNA上的所有动态过程都发生在染色质 landscape的背景下，染色质 landscape由核小体及其修饰、转录因子和染色质相关复合物组成。不同的染色质特征标志具有激活和沉默转录调控元件的位点和染色质结构域，它们在不同的细胞类型和发育过程中发生变化。
染色质全基因组特征的绘制传统上是使用染色质免疫沉淀(ChIP)进行的，其中染色质交联和溶解，并使用蛋白抗体或感兴趣的组蛋白修饰蛋白抗体来免疫沉淀结合的DNA(图1a)。自从35年前ChIP首次被描述以来，它的一般操作方法几乎没有改变，仍然有信噪比的问题和artifacts。另一种染色质分析策略是酶原位栓系，通过抗体或融合蛋白靶向染色质蛋白或修饰。然后，DNA被标记或切割，在过去的二十年里，一系列的酶系方法被引入。CUT&Tag是一种使用蛋白A - Tn5 (pA-Tn5)转座体融合蛋白的栓系方法(图1b)。在CUT&Tag中，渗透化的细胞或细胞核用针对特定染色质蛋白的抗体孵育，然后负载有嵌合末端接头的pA-Tn5被成功地连接到抗体结合位点。通过加入镁离子激活转座体，导致接头与附近的DNA结合。然后将它们扩增，产生测序文库。基于抗体栓系Tn5的方法由于pA-Tn5栓系后对样品进行严格的清洗和adaptor整合效率高，因此具有较高的灵敏度。与ChIP-seq相比，改进的信噪比转化为一个数量级的减少了绘制染色质特征所需的测序量，允许样本池(通常多达90个样本)通过文库的 barcode PCR在Illumina NGS测序仪上进行paired-end测序。
图1
（二）目标
本教程是为了按照Bench top CUT&Tag V.3 protocol处理和分析CUT&Tag数据而设计的。本教程中使用的数据是人类淋巴瘤K562细胞系中组蛋白修饰的图谱，但本教程一般适用于任何染色质蛋白，包括转录因子、RNA聚合酶II和表位标记蛋白。
（三）CUT&Tag数据处理和分析流程
（四）数据下载
我们使用的数据是Kaya-Okur et al. (2020)文章里的数据。
我们先从EBI上把下载地址都存到一个txt文件里：
$ cat download_list
ftp.sra.ebi.ac.uk/vol1/fastq/SRR122/017/SRR12246717/SRR12246717_1.fastq.gz
ftp.sra.ebi.ac.uk/vol1/fastq/SRR122/017/SRR12246717/SRR12246717_2.fastq.gz
ftp.sra.ebi.ac.uk/vol1/fastq/SRR110/040/SRR11074240/SRR11074240_1.fastq.gz
ftp.sra.ebi.ac.uk/vol1/fastq/SRR110/040/SRR11074240/SRR11074240_2.fastq.gz
ftp.sra.ebi.ac.uk/vol1/fastq/SRR110/054/SRR11074254/SRR11074254_1.fastq.gz
ftp.sra.ebi.ac.uk/vol1/fastq/SRR110/054/SRR11074254/SRR11074254_2.fastq.gz
ftp.sra.ebi.ac.uk/vol1/fastq/SRR110/058/SRR11074258/SRR11074258_1.fastq.gz
ftp.sra.ebi.ac.uk/vol1/fastq/SRR110/058/SRR11074258/SRR11074258_2.fastq.gz
ftp.sra.ebi.ac.uk/vol1/fastq/SRR119/024/SRR11923224/SRR11923224_1.fastq.gz
ftp.sra.ebi.ac.uk/vol1/fastq/SRR119/024/SRR11923224/SRR11923224_2.fastq.gz
ftp.sra.ebi.ac.uk/vol1/fastq/SRR875/001/SRR8754611/SRR8754611_1.fastq.gz
ftp.sra.ebi.ac.uk/vol1/fastq/SRR875/001/SRR8754611/SRR8754611_2.fastq.gz
ftp.sra.ebi.ac.uk/vol1/fastq/SRR875/002/SRR8754612/SRR8754612_1.fastq.gz
ftp.sra.ebi.ac.uk/vol1/fastq/SRR875/002/SRR8754612/SRR8754612_2.fastq.gz
调整到ascp下载的格式：
$ sed -e 's/ftp.sra.ebi.ac.uk//g' download_list >> download_list_2
$ cat download_list_2
/vol1/fastq/SRR122/017/SRR12246717/SRR12246717_1.fastq.gz
/vol1/fastq/SRR122/017/SRR12246717/SRR12246717_2.fastq.gz
/vol1/fastq/SRR110/040/SRR11074240/SRR11074240_1.fastq.gz
/vol1/fastq/SRR110/040/SRR11074240/SRR11074240_2.fastq.gz
/vol1/fastq/SRR110/054/SRR11074254/SRR11074254_1.fastq.gz
/vol1/fastq/SRR110/054/SRR11074254/SRR11074254_2.fastq.gz
/vol1/fastq/SRR110/058/SRR11074258/SRR11074258_1.fastq.gz
/vol1/fastq/SRR110/058/SRR11074258/SRR11074258_2.fastq.gz
/vol1/fastq/SRR119/024/SRR11923224/SRR11923224_1.fastq.gz
/vol1/fastq/SRR119/024/SRR11923224/SRR11923224_2.fastq.gz
/vol1/fastq/SRR875/001/SRR8754611/SRR8754611_1.fastq.gz
/vol1/fastq/SRR875/001/SRR8754611/SRR8754611_2.fastq.gz
/vol1/fastq/SRR875/002/SRR8754612/SRR8754612_1.fastq.gz
/vol1/fastq/SRR875/002/SRR8754612/SRR8754612_2.fastq.gz
使用aspera下载（参数解释请参考文章：ChIP-seq实践(H3K27Ac,enhancer的筛选和enhancer相关基因的GO分析)）：
$ ascp -QT -l 300m -P33001 -k1 -i ~/.aspera/connect/etc/asperaweb_id_dsa.openssh --mode recv --host fasp.sra.ebi.ac.uk --user era-fasp --file-list download_list_2 ./
下载后，把文件名改成样品名（我这里是手动的，有兴趣的同学可以自己搜索一下批量改文件名）：
$ ll
total 583040
-rw------- 1 fangy04 fangy04 41887306 Dec 22 19:53 SH_Hs_IgG_rep1_1.fastq.gz
-rw------- 1 fangy04 fangy04 43398133 Dec 22 19:53 SH_Hs_IgG_rep1_2.fastq.gz
-rw------- 1 fangy04 fangy04 44682456 Dec 22 19:55 SH_Hs_IgG_rep2_R1_1.fastq.gz
-rw------- 1 fangy04 fangy04 45488548 Dec 22 19:56 SH_Hs_IgG_rep2_R1_2.fastq.gz
-rw------- 1 fangy04 fangy04 23285383 Dec 22 19:56 SH_Hs_IgG_rep2_R2_1.fastq.gz
-rw------- 1 fangy04 fangy04 24206889 Dec 22 19:57 SH_Hs_IgG_rep2_R2_2.fastq.gz
-rw------- 1 fangy04 fangy04 57747126 Dec 22 19:54 SH_Hs_K27m3_rep1_1.fastq.gz
-rw------- 1 fangy04 fangy04 60076447 Dec 22 19:55 SH_Hs_K27m3_rep1_2.fastq.gz
-rw------- 1 fangy04 fangy04 55390834 Dec 22 19:50 SH_Hs_K27m3_rep2_1.fastq.gz
-rw------- 1 fangy04 fangy04 56801290 Dec 22 19:50 SH_Hs_K27m3_rep2_2.fastq.gz
-rw------- 1 fangy04 fangy04 30945039 Dec 22 19:51 SH_Hs_K4m3_rep1_1.fastq.gz
-rw------- 1 fangy04 fangy04 32532245 Dec 22 19:51 SH_Hs_K4m3_rep1_2.fastq.gz
-rw------- 1 fangy04 fangy04 39599278 Dec 22 19:52 SH_Hs_K4m3_rep2_1.fastq.gz
-rw------- 1 fangy04 fangy04 40380295 Dec 22 19:52 SH_Hs_K4m3_rep2_2.fastq.gz
第二节 数据预处理
（一）FastQC【可选步骤】
这一步不是必需的。如果用户生成自己的数据，而FastQC是用户的常规检查程序之一，我们提供此步骤作为例子。（这里我就用脚本来批量进行fastqc分析了）
#!/bin/bash
cat /gpfs/home/fangy04/cut_tag/filenames | while read i
do
xargs fastqc -t 20 -o /gpfs/home/fangy04/cut_tag/fastqc/
done
我们可以随便打开一个结果查看：
需要注意的是：read开始时序列内容不一致是CUT&Tag reads的常见现象。没有通过fastqc测试并不意味着数据失败。
（1）这可能是由于Tn5偏好。
（2）你可能检测到的是10 bp的周期性，在长度分布中显示为锯齿模式（如果不理解这句话，可以阅读文献阅读笔记：CUT&Tag）。如果是，这是正常的，不会影响比对或peak calling。在任何情况下，我们都不建议进行trimming!!!因为我们列出的bowtie2参数将在不进行trimming的情况下给出准确的映射信息。
（二）合并技术性重复/lines【可选步骤】
有时，为了提高效率，样本通常需要跨多个通道进行测序，并且可以在比对前将其放在一起。如果想检查同一样品的不同lane序列之间的重现性，可以跳过这一步，分别对各测序文件(fastq文件)进行比对。
#代码来自官网
##== linux command ==##
histName="K27me3_rep1"
mkdir -p ${projPath}/fastq
cat ${projPath}/data/${histName}/*_R1_*.fastq.gz >${projPath}/fastq/${histName}_R1.fastq.gz
cat ${projPath}/data/${histName}/*_R2_*.fastq.gz >${projPath}/fastq/${histName}_R2.fastq.gz
第三节 比对
（一）Bowtie2比对【必需步骤】
结合Tn5 adaptors和barcode PCR引物的CUT&Tag插入文库结构如下图所示:
我们的标准流程是在HiSeq 2500 flowcell上对多达90个混合样本进行单index 25x25 PE Illumina测序，每个样本都有一个独特的PCR引物barcode。通过调整每个文库的数量，可提供约500万个paired-end reads，从而利用一种特异性和高产抗体为丰富的染色质特征提供高质量的分析。较少富集的特征通常需要较少的reads，而较低质量的抗体可能会增加reads的数量，以产生稳定的染色质谱。关于CUT&Tag的特征recall和测序深度的深入讨论已经发表(Kaya-Okur et al 2020)。
（1）比对到HG38
首先写个文件名的txt：
$ cat filenames
SH_Hs_IgG_rep1
SH_Hs_IgG_rep2_R1
SH_Hs_IgG_rep2_R2
SH_Hs_K27m3_rep1
SH_Hs_K27m3_rep2
SH_Hs_K4m3_rep1
SH_Hs_K4m3_rep2
写个脚本进行比对：
#!/bin/bash
bowtie2_ref=/gpfs/share/apps/iGenomes/Homo_sapiens/UCSC/hg38/Sequence/Bowtie2Index/genome
cat filenames | while read i
do
bowtie2 --end-to-end --very-sensitive --no-mixed --phred33 -I 10 -X 700 -p 8 -x $bowtie2_ref -1 ${i}_1.fastq.gz -2 ${i}_2.fastq.gz -S ./alignment/${i}.sam 2> ./alignment/${i}_bowtie2.txt
done
（2）比对到spike-in基因组进行校准【可选步骤/推荐】
这一步是可选的，但是推荐你进行这一步，但要根据你的实验方法来进行选择。
大肠杆菌的DNA与细菌产生的pA-Tn5蛋白一起携带，并在反应过程中受到非特异性标记。定位到大肠杆菌基因组的总reads的比例取决于表位靶向的CUT&Tag的产量，因此也取决于所使用的细胞数量和染色质中表位的丰度。由于在CUT&Tag反应中加入一定量的pA-Tn5并带来一定量的大肠杆菌DNA，因此大肠杆菌reads可用于一系列实验中表位丰度的标准化。更多讨论请见第五节。
首先构建Ecoli的基因组索引（参考：Create bowtie2 index）：
$ bowtie2-build /gpfs/share/apps/iGenomes/Escherichia_coli_K_12_MG1655/NCBI/2001-10-15/Sequence/WholeGenomeFasta/genome.fa /gpfs/home/fangy04/Ecoli_bowtie2index/ecoli
写个脚本批量处理：
#!/bin/bash
spikeInRef="/gpfs/home/fangy04/Ecoli_bowtie2index/ecoli" #这里要注意的是这是Ecoli的索引
chromSize="/gpfs/share/apps/iGenomes/Homo_sapiens/UCSC/hg38/Sequence/WholeGenomeFasta/GenomeSize.xml" #这里是Hg38的基因组大小
## bowtie2 to Ecoli
cat filenames | while read i
do
bowtie2 --end-to-end --very-sensitive --no-mixed --no-discordant --phred33 -I 10 -X 700 -p 8 -x ${spikeInRef} -1 ${i}_1.fastq.gz -2 ${i}_2.fastq.gz -S ./alignment/${i}_spikein.sam 2> ./alignment/${i}_spikein.txt
seqDepthDouble=`samtools view -F 0x04 ./alignment/${i}_spikein.sam | wc -l`
seqDepth=$((seqDepthDouble/2))
echo $seqDepth >./alignment/${i}_bowtie2_spikeIn.seqDepth
done
生成文件：
-rw------- 1 fangy04 fangy04 5 Dec 23 10:38 SH_Hs_K4m3_rep2_bowtie2_spikeIn.seqDepth
-rw------- 1 fangy04 fangy04 366130966 Dec 23 10:38 SH_Hs_K4m3_rep2_spikein.sam
-rw------- 1 fangy04 fangy04 503 Dec 23 10:38 SH_Hs_K4m3_rep2_spikein.txt
-rw------- 1 fangy04 fangy04 4 Dec 23 10:37 SH_Hs_K4m3_rep1_bowtie2_spikeIn.seqDepth
-rw------- 1 fangy04 fangy04 306316978 Dec 23 10:37 SH_Hs_K4m3_rep1_spikein.sam
-rw------- 1 fangy04 fangy04 501 Dec 23 10:37 SH_Hs_K4m3_rep1_spikein.txt
-rw------- 1 fangy04 fangy04 4 Dec 23 10:37 SH_Hs_K27m3_rep2_bowtie2_spikeIn.seqDepth
-rw------- 1 fangy04 fangy04 524919105 Dec 23 10:37 SH_Hs_K27m3_rep2_spikein.sam
-rw------- 1 fangy04 fangy04 501 Dec 23 10:37 SH_Hs_K27m3_rep2_spikein.txt
-rw------- 1 fangy04 fangy04 4 Dec 23 10:36 SH_Hs_K27m3_rep1_bowtie2_spikeIn.seqDepth
-rw------- 1 fangy04 fangy04 579907784 Dec 23 10:36 SH_Hs_K27m3_rep1_spikein.sam
-rw------- 1 fangy04 fangy04 501 Dec 23 10:36 SH_Hs_K27m3_rep1_spikein.txt
-rw------- 1 fangy04 fangy04 6 Dec 23 10:35 SH_Hs_IgG_rep2_R2_bowtie2_spikeIn.seqDepth
-rw------- 1 fangy04 fangy04 229807348 Dec 23 10:35 SH_Hs_IgG_rep2_R2_spikein.sam
-rw------- 1 fangy04 fangy04 505 Dec 23 10:35 SH_Hs_IgG_rep2_R2_spikein.txt
-rw------- 1 fangy04 fangy04 6 Dec 23 10:35 SH_Hs_IgG_rep2_R1_bowtie2_spikeIn.seqDepth
-rw------- 1 fangy04 fangy04 433174854 Dec 23 10:35 SH_Hs_IgG_rep2_R1_spikein.sam
-rw------- 1 fangy04 fangy04 505 Dec 23 10:35 SH_Hs_IgG_rep2_R1_spikein.txt
-rw------- 1 fangy04 fangy04 6 Dec 23 10:34 SH_Hs_IgG_rep1_bowtie2_spikeIn.seqDepth
-rw------- 1 fangy04 fangy04 424349914 Dec 23 10:34 SH_Hs_IgG_rep1_spikein.sam
-rw------- 1 fangy04 fangy04 505 Dec 23 10:34 SH_Hs_IgG_rep1_spikein.txt
为了进行spike-in标准化，Reads比对到E.coli基因组要多用2个参数--no-overlap 和--no-dovetail防止交叉比对到human基因组。
（3）比对总结
挑一个比对结果查看：
$ cat SH_Hs_K4m3_rep1_bowtie2.txt #比对到Hg38的情况
1581710 reads; of these:
1581710 (100.00%) were paired; of these:
87391 (5.53%) aligned concordantly 0 times
1360015 (85.98%) aligned concordantly exactly 1 time
134304 (8.49%) aligned concordantly >1 times
----
87391 pairs aligned concordantly 0 times; of these:
6102 (6.98%) aligned discordantly 1 time
94.86% overall alignment rate
$ cat SH_Hs_K4m3_rep1_spikein.txt #比对到ecoli的情况
1581710 reads; of these:
1581710 (100.00%) were paired; of these:
1581335 (99.98%) aligned concordantly 0 times
365 (0.02%) aligned concordantly exactly 1 time
10 (0.00%) aligned concordantly >1 times
0.02% overall alignment rate
$ cat SH_Hs_K4m3_rep1_bowtie2_spikeIn.seqDepth
375
上面的比对结果可以解读为（Hg38 ）：
1581710 是测序深度，比如说：总的paired reads数。
87391是没有比对上的read pairs。
1360015 + 134304是成功比对上的read pairs数。
94.86%是总比对率。
（二）报告测序比对总结【必需步骤】
总结原始reads和唯一比对reads，以报告比对的效率。对于高质量数据，比对频率预计为> 80%。一般来说，CUT&Tag数据的背景非常低，因此，只需100万个比对上的片段就可以为人类基因组中的组蛋白修饰提供可靠的profiles。低丰度的转录因子和染色质蛋白的谱分析可能需要10倍于下游分析的图谱片段。
我们可以评估以下指标:
测序深度
比对率
比对上片段的数量
重复率
unique文库大小
片段大小分布
（1）测序深度
这一部分是在R里进行的，因为我们要画列表。
> library(magrittr)
#Sequencing depth
> sampleList = c("SH_Hs_K27m3_rep1", "SH_Hs_K27m3_rep2", "SH_Hs_K4m3_rep1", "SH_Hs_K4m3_rep2", "SH_Hs_IgG_rep1","SH_Hs_IgG_rep2_R1","SH_Hs_IgG_rep2_R2")
> histList = c("K27m3", "K4m3", "IgG")
## Collect the alignment results from the bowtie2 alignment summary files
> alignResult = c()
> for(hist in sampleList){
alignRes = read.table(paste0(hist, "_bowtie2.txt"), header = FALSE, fill = TRUE,skip = 7)
alignRate = substr(alignRes$V1[10], 1, nchar(as.character(alignRes$V1[10]))-1)
histInfo = strsplit(hist, "_")[[1]]
alignResult = data.frame(Histone = histInfo[3], Replicate = histInfo[4],
SequencingDepth = alignRes$V1[1] %>% as.character %>% as.numeric,
MappedFragNum_hg38 = alignRes$V1[4] %>% as.character %>% as.numeric + alignRes$V1[5] %>% as.character %>% as.numeric,
AlignmentRate_hg38 = alignRate %>% as.numeric) %>% rbind(alignResult, .)
}
> alignResult$Histone = factor(alignResult$Histone, levels = histList)
> alignResult %>% mutate(AlignmentRate_hg38 = paste0(AlignmentRate_hg38, "%"))
（2）Spike-in比对
> spikeAlign = c()
> for(hist in sampleList){
spikeRes = read.table(paste0(hist, "_spikein.txt"), header = FALSE, fill = TRUE,skip = 7)
alignRate = substr(spikeRes$V1[6], 1, nchar(as.character(spikeRes$V1[6]))-1)
histInfo = strsplit(hist, "_")[[1]]
spikeAlign = data.frame(Histone = histInfo[3], Replicate = histInfo[4],
SequencingDepth = spikeRes$V1[1] %>% as.character %>% as.numeric,
MappedFragNum_spikeIn = spikeRes$V1[4] %>% as.character %>% as.numeric + spikeRes$V1[5] %>% as.character %>% as.numeric,
AlignmentRate_spikeIn = alignRate %>% as.numeric) %>% rbind(spikeAlign, .)
}
> spikeAlign$Histone = factor(spikeAlign$Histone, levels = histList)
> spikeAlign %>% mutate(AlignmentRate_spikeIn = paste0(AlignmentRate_spikeIn, "%"))
（3）比对到hg38和E.coli的总结
> alignSummary = left_join(alignResult, spikeAlign, by = c("Histone", "Replicate", "SequencingDepth")) %>%
mutate(AlignmentRate_hg38 = paste0(AlignmentRate_hg38, "%"),
AlignmentRate_spikeIn = paste0(AlignmentRate_spikeIn, "%"))
> alignSummary
（4）可视化测序深度和比对结果
> library(ggridges)
> library(ggplot2)
> library(viridis)
> library(ggpubr)
> fig3A = alignResult %>% ggplot(aes(x = Histone, y = SequencingDepth/1000000, fill = Histone)) +
geom_boxplot() +
geom_jitter(aes(color = Replicate), position = position_jitter(0.15)) +
scale_fill_viridis(discrete = TRUE, begin = 0.1, end = 0.9, option = "magma", alpha = 0.8) +
scale_color_viridis(discrete = TRUE, begin = 0.1, end = 0.9) +
theme_bw(base_size = 18) +
ylab("Sequencing Depth per Million") +
xlab("") +
ggtitle("A. Sequencing Depth")
> fig3B = alignResult %>% ggplot(aes(x = Histone, y = MappedFragNum_hg38/1000000, fill = Histone)) +
geom_boxplot() +
geom_jitter(aes(color = Replicate), position = position_jitter(0.15)) +
scale_fill_viridis(discrete = TRUE, begin = 0.1, end = 0.9, option = "magma", alpha = 0.8) +
scale_color_viridis(discrete = TRUE, begin = 0.1, end = 0.9) +
theme_bw(base_size = 18) +
ylab("Mapped Fragments per Million") +
xlab("") +
ggtitle("B. Alignable Fragment (hg38)")
> fig3C = alignResult %>% ggplot(aes(x = Histone, y = AlignmentRate_hg38, fill = Histone)) +
geom_boxplot() +
geom_jitter(aes(color = Replicate), position = position_jitter(0.15)) +
scale_fill_viridis(discrete = TRUE, begin = 0.1, end = 0.9, option = "magma", alpha = 0.8) +
scale_color_viridis(discrete = TRUE, begin = 0.1, end = 0.9) +
theme_bw(base_size = 18) +
ylab("% of Mapped Fragments") +
xlab("") +
ggtitle("C. Alignment Rate (hg38)")
> fig3D = spikeAlign %>% ggplot(aes(x = Histone, y = AlignmentRate_spikeIn, fill = Histone)) +
geom_boxplot() +
geom_jitter(aes(color = Replicate), position = position_jitter(0.15)) +
scale_fill_viridis(discrete = TRUE, begin = 0.1, end = 0.9, option = "magma", alpha = 0.8) +
scale_color_viridis(discrete = TRUE, begin = 0.1, end = 0.9) +
theme_bw(base_size = 18) +
ylab("Spike-in Alignment Rate") +
xlab("") +
ggtitle("D. Alignment Rate (E.coli)")
#把图组合起来
> ggarrange(fig3A, fig3B, fig3C, fig3D, ncol = 2, nrow=2, common.legend = TRUE, legend="bottom")
在一项针对65,000个 K562细胞中富集H3K27me3组蛋白修饰的CUT&Tag实验中，E.coli的reads比例在~0.01%至10%之间。如果细胞数量少或表位数量少，E.coli的reads可占总reads的70%。对于IgG对照，E.coli的reads比例通常比组蛋白修饰的要高得多。
（三）去重【可选/推荐】
CUT&Tag将adaptors整合到抗体栓系pA-Tn5附近的DNA中，并且整合的准确位置受周围DNA可及性的影响。由于这个原因，共享精确起始和结束位置的片段是常见的，而这种“duplicates”可能不是由于PCR过程中的重复。在实践中，我们发现对于高质量的CUT&Tag数据集，重复率是很低的，甚至“重复”片段也很可能是真实的片段。因此，我们不建议删除duplicates。在用非常少量的材料或怀疑是PCR重复的时候，duplicates可以被去除。下面的命令展示了如何使用Picard检查重复率。
这里要说明：运行picard最好要在配置高一点的电脑运行，或者在服务器上，请参看我之前的文章：ATAC-seq分析练习，我之前尝试过用自己的电脑运行，一个样品跑了过夜也没跑完，电脑带不动。如果没有服务器，但实在想过一遍CUT&tag分析流程，可以用别的去重软件，在ATAC这篇文章里也有提到。
NOTE：这里我手动把比对到hg38的sam文件名都加了*_bowtie2，用来区别spikein.sam。另外这里的gatk安装途径需要改，不要直接copy。
#!/bin/bash
file_path="/gpfs/home/fangy04/cut_tag"
cat filenames | while read i
do
## Sort by coordinate
java -jar /gpfs/share/apps/gatk/4.1.7.0/gatk-package-4.1.7.0-local.jar SortSam -I ${file_path}/alignment/${i}_bowtie2.sam -O ${file_path}/deduplicates/${i}_bowtie2.sorted.sam --SORT_ORDER coordinate
## mark duplicates
java -jar /gpfs/share/apps/gatk/4.1.7.0/gatk-package-4.1.7.0-local.jar MarkDuplicates -I ${file_path}/deduplicates/${i}_bowtie2.sorted.sam -O ${file_path}/deduplicates/${i}_bowtie2.sorted.dupMarked.sam --METRICS_FILE ${file_path}/deduplicates/${i}_picard.dupMark.txt
## remove duplicates
java -jar /gpfs/share/apps/gatk/4.1.7.0/gatk-package-4.1.7.0-local.jar MarkDuplicates -I ${file_path}/deduplicates/${i}_bowtie2.sorted.sam -O ${file_path}/deduplicates/${i}_bowtie2.sorted.rmDup.sam --REMOVE_DUPLICATES true --METRICS_FILE ${file_path}/deduplicates/${i}_picard.rmDup.txt
done
运行完后，再回到R里进行可视化：
> library(magrittr)
> library(dplyr)
> library(tidyverse)
> sampleList = c("SH_Hs_K27m3_rep1", "SH_Hs_K27m3_rep2", "SH_Hs_K4m3_rep1", "SH_Hs_K4m3_rep2", "SH_Hs_IgG_rep1","SH_Hs_IgG_rep2_R1","SH_Hs_IgG_rep2_R2")
> histList = c("K27m3", "K4m3", "IgG")
> dupResult = c()
> for(hist in sampleList){
dupRes = read.table(paste0(hist, "_picard.rmDup.txt"), header = TRUE, fill = TRUE,skip = 6)
histInfo = strsplit(hist, "_")[[1]]
dupResult = data.frame(Histone = histInfo[3], Replicate = histInfo[4], MappedFragNum_hg38 = dupRes$READ_PAIRS_EXAMINED[1] %>% as.character %>% as.numeric, DuplicationRate = dupRes$PERCENT_DUPLICATION[1] %>% as.character %>% as.numeric * 100, EstimatedLibrarySize = dupRes$ESTIMATED_LIBRARY_SIZE[1] %>% as.character %>% as.numeric) %>% mutate(UniqueFragNum = MappedFragNum_hg38 * (1-DuplicationRate/100)) %>% rbind(dupResult, .)
}
> dupResult$Histone = factor(dupResult$Histone, levels = histList)
> dupResult$percent = c("%")
> dupResult$DuplicationRate <- str_c(dupResult$DuplicationRate,dupResult$percent)
> dupResult = dupResult[,-7]
> alignDupSummary = cbind(alignSummary,dupResult[,c(4,5,6)])
> alignDupSummary
（1）在这些样本数据集中，IgG对照样本有相对较高的重复率，因为该样本中的reads来自于CUT&Tag反应中的非特异性标记。因此，在进行下游分析之前，应该将重复的IgG数据集删除。
（2）文库的大小是Picard根据PE重复计算出的文库中唯一分子的数量。
（3）估计的文库大小与靶标表位的丰度和所用抗体的质量成正比，而IgG样本的文库估计大小非常低。
（4）唯一的片段数由MappedFragNum_hg38 * (1-DuplicationRate/100)计算。
## generate boxplot figure for the duplication rate
> library(ggridges)
> library(ggplot2)
> library(viridis)
> library(ggpubr)
> m = substr(dupResult$DuplicationRate,1,6)
> m = as.numeric(m)
> fig4A = dupResult %>% ggplot(aes(x = Histone, y = m, fill = Histone)) +
geom_boxplot() +
geom_jitter(aes(color = Replicate), position = position_jitter(0.15)) +
scale_fill_viridis(discrete = TRUE, begin = 0.1, end = 0.9, option = "magma", alpha = 0.8) +
scale_color_viridis(discrete = TRUE, begin = 0.1, end = 0.9) +
theme_bw(base_size = 18) +
ylab("Duplication Rate (*100%)") +
xlab("")
> fig4B = dupResult %>% ggplot(aes(x = Histone, y = EstimatedLibrarySize, fill = Histone)) +
geom_boxplot() +
geom_jitter(aes(color = Replicate), position = position_jitter(0.15)) +
scale_fill_viridis(discrete = TRUE, begin = 0.1, end = 0.9, option = "magma", alpha = 0.8) +
scale_color_viridis(discrete = TRUE, begin = 0.1, end = 0.9) +
theme_bw(base_size = 16) +
ylab("Estimated Library Size") +
xlab("")
> fig4C = dupResult %>% ggplot(aes(x = Histone, y = UniqueFragNum, fill = Histone)) +
geom_boxplot() +
geom_jitter(aes(color = Replicate), position = position_jitter(0.15)) +
scale_fill_viridis(discrete = TRUE, begin = 0.1, end = 0.9, option = "magma", alpha = 0.8) +
scale_color_viridis(discrete = TRUE, begin = 0.1, end = 0.9) +
theme_bw(base_size = 16) +
ylab("# of Unique Fragments") +
xlab("")
> p <- ggarrange(fig4A, fig4B, fig4C, ncol = 3, common.legend = TRUE, legend="bottom")
> p
（四）评估比对上的片段大小分布【必需步骤】
CUT&Tag在被酶栓附近的染色质颗粒的两侧插入adaptors，虽然染色质颗粒内的标记也可能发生。因此，针对组蛋白修饰的CUT&Tag反应主要导致核小体长度(~180 bp)或数倍于该长度的片段。CUT&Tag靶向转录因子主要产生核小体大小的片段和不同数量的短片段，分别来自邻近的核小体和因子结合位点。DNA的标记在核小体表面的也发生，用单碱基对分辨率绘制片段长度显示了10 bp的锯齿状周期性，这是成功的CUT&Tag实验的典型特征。
#!/bin/bash
projPath="/gpfs/home/fangy04/cut_tag/alignment/"
cat filenames | while read i
do
## Extract the 9th column from the alignment sam file which is the fragment length
samtools view -F 0x04 $projPath/${i}_bowtie2.sam | awk -F'\t' 'function abs(x){return ((x < 0.0) ? -x : x)} {print abs($9)}' | sort | uniq -c | awk -v OFS="\t" '{print $2, $1/2}' >$projPath/${i}_fragmentLen.txt
done
生成文件：
-rw------- 1 fangy04 fangy04 6033 Dec 23 19:55 SH_Hs_K4m3_rep2_fragmentLen.txt
-rw------- 1 fangy04 fangy04 5959 Dec 23 19:55 SH_Hs_K4m3_rep1_fragmentLen.txt
-rw------- 1 fangy04 fangy04 6025 Dec 23 19:55 SH_Hs_K27m3_rep2_fragmentLen.txt
-rw------- 1 fangy04 fangy04 5988 Dec 23 19:55 SH_Hs_K27m3_rep1_fragmentLen.txt
-rw------- 1 fangy04 fangy04 5763 Dec 23 19:54 SH_Hs_IgG_rep2_R2_fragmentLen.txt
-rw------- 1 fangy04 fangy04 6012 Dec 23 19:54 SH_Hs_IgG_rep2_R1_fragmentLen.txt
-rw------- 1 fangy04 fangy04 6076 Dec 23 19:54 SH_Hs_IgG_rep1_fragmentLen.txt
在R里进行可视化：
## Collect the fragment size information
> library(magrittr)
> library(dplyr)
> library(tidyverse)
> sampleList = c("SH_Hs_K27m3_rep1", "SH_Hs_K27m3_rep2", "SH_Hs_K4m3_rep1", "SH_Hs_K4m3_rep2", "SH_Hs_IgG_rep1","SH_Hs_IgG_rep2_R1","SH_Hs_IgG_rep2_R2")
> histList = c("K27m3", "K4m3", "IgG")
> fragLen = c()
> for(hist in sampleList){
histInfo = strsplit(hist, "_")[[1]]
fragLen = read.table(paste0(hist, "_fragmentLen.txt"), header = FALSE) %>% mutate(fragLen = V1 %>% as.numeric, fragCount = V2 %>% as.numeric, Weight = as.numeric(V2)/sum(as.numeric(V2)), Histone = histInfo[3], Replicate = histInfo[4], sampleInfo = hist) %>% rbind(fragLen, .)
}
> fragLen$sampleInfo = factor(fragLen$sampleInfo, levels = sampleList)
> fragLen$Histone = factor(fragLen$Histone, levels = histList)
## Generate the fragment size density plot (violin plot)
> fig5A = fragLen %>% ggplot(aes(x = sampleInfo, y = fragLen, weight = Weight, fill = Histone)) +
geom_violin(bw = 5) +
scale_y_continuous(breaks = seq(0, 800, 50)) +
scale_fill_viridis(discrete = TRUE, begin = 0.1, end = 0.9, option = "magma", alpha = 0.8) +
scale_color_viridis(discrete = TRUE, begin = 0.1, end = 0.9) +
theme_bw(base_size = 20) +
ggpubr::rotate_x_text(angle = 30) +
ylab("Fragment Length") +
xlab("")
> fig5B = fragLen %>% ggplot(aes(x = fragLen, y = fragCount, color = Histone, group = sampleInfo, linetype = Replicate)) +
geom_line(size = 1) +
scale_color_viridis(discrete = TRUE, begin = 0.1, end = 0.9, option = "magma") +
theme_bw(base_size = 20) +
xlab("Fragment Length") +
ylab("Count") +
coord_cartesian(xlim = c(0, 500))
> ggarrange(fig5A, fig5B, ncol = 2)
较小的片段(50-100 bp)可能是由于Tn5被栓系在核小体表面以及linker区域，因此小片段可能不是背景。
（五）评估重复样品的重现率
重复之间的数据再现性是通过跨基因组的比对上的read counts的相关性分析来评估的。为了简单化，当文件格式被转换成bed文件时，我们将在第四节之后进行这个分析。
