CUT&Tag 数据处理与分析教程 二：质控（不需要修剪 reads！不需要修剪 reads！不需要修剪reads！）和数据比对

第一部分我们介绍了数据的下载相关内容：

此部分需要重点注意的内容：

reads 开始时序列内容不一致是 CUT&Tag reads 中常见的现象。没有通过每个基本的 sequnence 内容不意味着你的数据失败。

• 可能由于是 Tn5 的偏好性
• 可能检测到的是 10 bp 的周期性，该周期性在长度分布中显示为锯齿状。如果是这样，这是正常现象，不会影响校准或 Peak calling。无论如何，我们都不建议您修剪 reads，因为我们列出的 bowtie2 参数将提供准确的 mapping 信息而无需修剪。

II. 数据预处理

使用 FASTQC 进行质量控制（可选）

这一步不是必需的。如果用户正在生成自己的数据，并且 FastQC 是用户组中的例行检查过程之一，那么我们提供这个步骤作为故障排除说明。

FASTQC 安装

##== linux 命令 ==##
mkdir -p $projPath/tools
wget -P $projPath/tools https://www.bioinformatics.babraham.ac.uk/projects/fastqc/fastqc_v0.11.9.zip
cd $projPath/tools
unzip fastqc_v0.11.9.zip

运行 FASTQC 进行质量检查

##== linux 命令 ==##
mkdir -p ${projPath}/fastqFileQC/${histName}

$projPath/tools/FastQC/fastqc -o ${projPath}/fastqFileQC/${histName} -f fastq ${projPath}/fastq/${histName}_R1.fastq.gz
$projPath/tools/FastQC/fastqc -o ${projPath}/fastqFileQC/${histName} -f fastq ${projPath}/fastq/${histName}_R2.fastq.gz

对质量检查结果进行记录

FASTQC 结果参考链接: https://www.bioinformatics.babraham.ac.uk/projects/fastqc/bad_sequence_fastqc.html^[1]

“> 图 3：每个碱基序列内容未通过 FastQC 质量检查。

reads 开始时序列内容不一致是 CUT&Tag reads 中常见的现象。 没有通过每个基本的 sequnence 内容不意味着你的数据失败。

• 可能由于是 Tn5 的偏好性
• 您可能检测到的是 10 bp 的周期性，该周期性在长度分布中显示为锯齿状。如果是这样，这是正常现象，不会影响校准或 Peak calling。无论如何，我们都不建议您修剪 reads，因为我们列出的 bowtie2 参数将提供准确的 mapping 信息而无需修剪。

如果有需要，合并技术重复或者一个 lanes

有时，为了提高效率，样本常常跨多个通道进行排序，可以在回帖到基因组之前进行汇总。如果您希望检查相同样本的不同行序列之间的重现性，可以跳过此步骤，并分别对每个排序文件 (fastq 文件) 进行排列。

##== linux 命令 ==##
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

---

III. 比对

Bowtie2 比对

带有 Tn5 接头和条形码 barcoded PCR 引物的 CUT&Tag 插入文库的结构如下所示：

“> 图 4：CUT&Tag insert libraries with the sequence of adapters.

我们的标准流程是对一个 HiSeq 2500 flowcell 上的 90 个汇总样本进行单指数 25x25 PE Illumina 测序，每个样本都有一个独特的 PCR 引物条形码。每个文库的数量调整为提供约 500 万对双端 reads，这为丰富的染色质特征提供了高质量的分析，具有特异性和高产的抗体。较少丰富的特征通常需要较少的 reads 量，而较低质量的抗体可能会增加产生稳健染色质谱所需的 reads 量。一个深入讨论的 feature recall 和测序深度的 CUT&Tag 已经发表（Kaya-Okur et al 2020）。

比对到 hg38 基因组上

##== linux 命令 ==##
cores=8
ref="/path/to/bowtie2Index/hg38"

mkdir -p ${projPath}/alignment/sam/bowtie2_summary
mkdir -p ${projPath}/alignment/bam
mkdir -p ${projPath}/alignment/bed
mkdir -p ${projPath}/alignment/bedgraph

## Build the bowtie2 reference genome index if needed:
## bowtie2-build path/to/hg38/fasta/hg38.fa /path/to/bowtie2Index/hg38
bowtie2 --end-to-end --very-sensitive --no-mixed --no-discordant --phred33 -I 10 -X 700 -p ${cores} -x ${ref} -1 ${projPath}/fastq/${histName}_R1.fastq.gz -2 ${projPath}/fastq/${histName}_R2.fastq.gz -S ${projPath}/alignment/sam/${histName}_bowtie2.sam &> ${projPath}/alignment/sam/bowtie2_summary/${histName}_bowtie2.txt

使用 Bowtie2 进行双端比对参数应设置为：--end-to-end --very-sensitive --no-mixed --no-discordant --phred33 -I 10 -X 700，用于回帖长度为 10-700 bp 的插入。
关键步骤：没有必要修剪从标准 25x25 PE 测序输出的测序结果，因为接头序列将不包括在 ** 插入 >25 bp** 的 reads。然而，对于执行 ** 较长测序 的用户， reads** 将需要通过 Cutadapt 进行修剪，并且在比对步骤中，比对参数应该为 --local --very-sensitive --no-mixed --no-discordant --phred33 -I 10 -X 700 以便忽略 reads 的 3' 末端 的任何剩余接头序列。

对 spike-in 基因组进行校准，以进行 spike-in 校准（可选 / 推荐）

“> 这部分是可选的，但建议取决于您的实验方案。

大肠杆菌的 DNA 与细菌产生的 pA-Tn5 蛋白一起携带，在反应过程中非特异性地被标记。比对到大肠杆菌的 reads 量占总 reads 量的比例。大肠杆菌的基因组依赖于 epitope-targeted CUT&Tag 的产量，因此也依赖于使用的细胞数量和染色质中表位的丰度。由于在 CUT&Tag 反应中加入了一定量的 pATn5，并带来了一定量的大肠杆菌 DNA，大肠杆菌 reads 可用于一系列实验中表位丰度的标准化。更多讨论请见第五部分。

##== linux command ==##
spikeInRef="/shared/ngs/illumina/henikoff/Bowtie2/Ecoli"
chromSize="/fh/fast/gottardo_r/yezheng_working/SupplementaryData/hg38/chromSize/hg38.chrom.size"

## bowtie2-build path/to/Ecoli/fasta/Ecoli.fa /path/to/bowtie2Index/Ecoli
bowtie2 --end-to-end --very-sensitive --no-mixed --no-discordant --phred33 -I 10 -X 700 -p ${cores} -x ${spikeInRef} -1 ${projPath}/fastq/${histName}_R1.fastq.gz -2 ${projPath}/fastq/${histName}_R2.fastq.gz -S $projPath/alignment/sam/${histName}_bowtie2_spikeIn.sam &> $projPath/alignment/sam/bowtie2_summary/${histName}_bowtie2_spikeIn.txt

seqDepthDouble=`samtools view -F 0x04`
seqDepth=$((seqDepthDouble/2))

echo $seqDepth >$projPath/alignment/sam/bowtie2_summary/${histName}_bowtie2_spikeIn.seqDepth

为了进行 **spike-in** 标准化，reads 与大肠杆菌基因组 U00096.3 再加上两个参数对齐 `--no-overlap` 和 `--no-dovetail` (`--end-to-end --very-sensitive --no-overlap --no-dovetail --no-mixed --no-discordant --phred33 -I 10 -X 700`) 避免可能的实验基因组与用于校准的携带的大肠杆菌 DNA 的交叉比对。

本内容疑问解答：

1、上面代码 samtools  中有一行 -F 0x04 不禁会想，这是什么意思？
首先自然是从 samtools 的官方说明书开始（请记住，绝多大数情况下，用一个软件，从软件本身所发表的文章和说明入手是最好的，不要为了偷懒，在网上到处看别人教程，比如此文，哈哈。。）， samtools 官网 samtools 官方说明书 ^[2] 。

你可以看到，非常详细的说明，包括 samtools 一些常用的命令。那么我们回到正题：命令中中参数 -F 表示 什么？我们可以看到此行命令前面是使用的 samtools view 功能，那么我们就需要查看 samtools view 的说明书：samtools view 说明书 ^[3] 。我们可以看到其中 -F INT 表示 **Do not output alignments with any bits set in INT present in the FLAG field. INT can be specified in hex by beginning with 0x' (i.e. /^0x[0-9A-F]+/) or in octal by beginning with0' (i.e./^0 [0-7]+/) [0]。其中 ****INT** 为 BAM/SAM 文件的 flag 值，至于什么是 flag 值，稍后再说。大致意思就是不出书符合 flag 条件对应的信息， flag 可以用 0x 开头的一些字符，也可以是 0 开始的八进制数值。与之相反的参数 -f 即输出符合 flag 条件的信息。

好了接着上面话题，那么 -F 后面跟着的 0x04 表示什么？那么我们先大致的从 flag 值说起。首先这个 flag 值表示我们测序 reads 回帖基因组后的状态信息。其次我们可以从 samtools 官网 samtools 官方说明书 ^[4] 页面可以看到其中有这样的一个内容：

flags
samtools flags INT|STR[,...]

Convert between textual and numeric flag representation.

FLAGS:

0x1 PAIRED paired-end (or multiple-segment) sequencing technology
0x2 PROPER_PAIR each segment properly aligned according to the aligner
0x4 UNMAP segment unmapped
0x8 MUNMAP next segment in the template unmapped
0x10 REVERSE SEQ is reverse complemented
0x20 MREVERSE SEQ of the next segment in the template is reverse complemented
0x40 READ1 the first segment in the template
0x80 READ2 the last segment in the template
0x100 SECONDARY secondary alignment
0x200 QCFAIL not passing quality controls
0x400 DUP PCR or optical duplicate
0x800 SUPPLEMENTARY supplementary alignment

额，好吧，遗憾这里没有介绍（其实就是 0x4 ）。哈哈，那么我们通过查看 flag 值对应解释的 在线网页 explain SAM flags^[5]  。我们输入 0x04 可以看到，表示没有回帖到基因组上的 reads 。结合参数 -F 来看 -F 0X04 表示的是不输出没有比对上的 reads ，换句话说，只输出比对上的 reads 。好了，这里应该都懂了。


总之 samtools 作为生信分析相关人员不可或缺的一个软件，我们有必要花点时间去看看它的原始文档。

最后贴一下官方文档中介绍的 samtools 常用的一些命令：

# 将压缩的 sam 文件转换为 bam 文件
samtools view -bt ref_list.txt -o aln.bam aln.sam.gz

# 将 bam 文件进行排序
samtools sort -T /tmp/aln.sorted -o aln.sorted.bam aln.bam

# 对 bam 文件建立索引，比如 IGV 查看 bam 文件等操作所必须的步骤；注意建立索引前先进行对 bam 文件排序
samtools index aln.sorted.bam

# 统计 bam 文件信息
samtools idxstats aln.sorted.bam

# 统计 bam 文件中 reads 的状态信息，比如唯一比对有多少，多重比对有多少，匹配率有多少
samtools flagstat aln.sorted.bam

# 
samtools stats aln.sorted.bam

samtools bedcov aln.sorted.bam

# 每个碱基位点的测序深度
samtools depth aln.sorted.bam

# 查看指定区域的 reads 信息，一般可以通过此命令来获得我们感兴趣的区域或者基因上的 reads 信息
samtools view aln.sorted.bam chr2:20,100,000-20,200,000

# 合并多个 bam 文件为一个
samtools merge out.bam in1.bam in2.bam in3.bam


samtools faidx ref.fasta

samtools fqidx ref.fastq

samtools tview aln.sorted.bam ref.fasta

samtools split merged.bam

# 检查 bam 文件是否有问题；如果任何输入文件没有有效的头文件或缺少 EOF 块，此命令将以非零的退出代码退出。否则它将成功退出(使用零退出码)。
samtools quickcheck in1.bam in2.cram

samtools dict -a GRCh38 -s "Homo sapiens" ref.fasta

samtools fixmate in.namesorted.sam out.bam

# 此命令对于重测序等分析的很重要，用来生成 bcf 文件，然后对接 bcftools 进行 SNP 和 INDEL 的分析
samtools mpileup -C50 -f ref.fasta -r chr3:1,000-2,000 in1.bam in2.bam

# 输出每条染色体覆盖的直方图或表格
samtools coverage aln.sorted.bam

# 在文本和数字标志表示之间切换，比如 0x1 → PAIRED
samtools flags PAIRED,UNMAP,MUNMAP

# 将 bam 文件转换为 fastq 文件
samtools fastq input.bam > output.fastq

# 同上，只不过是将 bam 文件转换为 fasta 文件
samtools fasta input.bam > output.fasta

samtools addreplacerg -r 'ID:fish' -r 'LB:1334' -r 'SM:alpha' -o output.bam input.bam

samtools collate -o aln.name_collated.bam aln.sorted.bam

samtools depad input.bam

# 标记重复序列用
samtools markdup in.algnsorted.bam out.bam

# 去除重复序列
samtools rmdup in.algnsorted.bam out.bam

好了，就这些吧，更多的内容，详细查看官方文档。

参考资料

[1]

https://www.bioinformatics.babraham.ac.uk/projects/fastqc/bad_sequence_fastqc.html: https://www.bioinformatics.babraham.ac.uk/projects/fastqc/bad_sequence_fastqc.html

[2]

samtools 官方说明书: http://www.htslib.org/doc/samtools.html

[3]

samtools view 说明书: http://www.htslib.org/doc/samtools-view.html

[4]

samtools 官方说明书: http://www.htslib.org/doc/samtools.html

[5]

在线网页 explain SAM flags: https://broadinstitute.github.io/picard/explain-flags.htm