跟着官网学习CUT&Tag数据分析，参考：protocol。用的自己的数据：
A_rep1, A_rep2, WT_rep1, WT_rep2。
1 前言
1.3 CUT&Tag 数据处理和分析大纲
2 数据预处理
2.1 质控【选做】
公司返回数据已经质检，这里就不做，也不做 trimming。
在任何情况下，我们不建议修剪，因为后面的 bowtie2 参数将提供准确的 mapping 信息而无需修剪。
3 比对
3.1 Bowtie2比对【必做】
3.1.1 比对到 mm10
本实验用的是小鼠细胞，所以参考基因组选 mm10。UCSC下载 latest/
使用 Bowtie2 进行双端比对，参数设置：
--end-to-end --very-sensitive --no-mixed --no-discordant --phred33 -I 10 -X 700， mapping 长度为 10-700 bp 的插入片段。关键步骤：无需修剪标准 25x25 PE 测序的 reads，因为接头序列不会包含在插入片段 >25 bp 的reads。但是，对于执行较长测序的用户，需要通过 Cutadapt 修剪，比对时，参数设置为--local --very-sensitive --no-mixed --no-discordant --phred33 -I 10 -X 700，以忽略 3' 端的任何剩余接头序列。
3.1.2 比对到 spike-in【根据实验方案选做】个人感觉不用做，后面用CPM校准好像差不多
因为试剂盒提供的是三段 spike-in 序列，所以这里另外分别计算了 spike-in1&2&3 的 reads，用以验证理论占比。
projPath="~/0908_test"
config1=~/0908_test/config1 #用于后续批量处理。第一列是样本名，第二、三列分别是双端测序原始数据的名称——R1.fq.gz和R2.fq.gz
ref="~/INDEX/bowtie2Index/mouse/mm10"
spikeInRef="~/INDEX/bowtie2Index/ecoli3/ecoli"
spikeInRef1="~/INDEX/bowtie2Index/ecoli3/ecoli1"
spikeInRef2="~/INDEX/bowtie2Index/ecoli3/ecoli2"
spikeInRef3="~/INDEX/bowtie2Index/ecoli3/ecoli3"
mkdir -p ${projPath}/alignment/sam/bowtie2_summary
mkdir -p ${projPath}/alignment/bam/picard_summary
mkdir -p ${projPath}/alignment/bed
module load tools/bowtie2-2.3.5.1
## Build the bowtie2 reference genome index if needed:
bowtie2-build -c GTTCCACTCCTGAAGTGTCAAGTACATCGCAAAGTCTCCGCAATTACACGCAAGAAAAAACCGCCATCAGGCGGCTTGGTGTTCTTTCAGTTCTTCAATTCGAATATTGGTTACGTCTGCATGTGCTATCTGCGCCCATATCATCCAGTGGTCGTAGCAGTCGTTGATGTTCTCCGCTTCGATAACTCTGTTGAATGGCTCTCCATTCCATTCTCCTGTGACTCGGAAGT ~/INDEX/bowtie2Index/ecoli3/ecoli1
bowtie2-build -c AGTCGGTGTGAATCCCATCAGCGTTACCGTTTCGCGGTGCTTCTTCAGTACGCTACGGCAAATGTCATCGACGTTTTTATCCGGAAACTGCTGTCTGGCTTTTTTTGATTTCAGAATTAGCCTGACGGGCAATGCTGCGAAGGGCGTTTTCCTGCTGAGGTGTCATTGAACAAGTCCCATGTCGGCAAGCATAAGCACACAGAATATGAAGCCCGCTGCCAGAAAAATGCATTCCGTGGTTGTCATACCT ~/INDEX/bowtie2Index/ecoli3/ecoli2
bowtie2-build -c CCTTACTGGAATCGATGGTGTCTCCGGTGTGAAAGAACACCAACAGGGGTGTTACCACTACCGCAGGAAAAGGAGGACGTGTGGCGAGACAGCGACGAAGTATCACCGACATAATCTGCGAAAACTGCAAATACCTTCCAACGAAACGCACCAGAAATAAACCCAAGCCAATCCCAAAAGAATCTGACGTAAAAACCTTCAACTACACGGCTCACCTGTGGGATATCCGGTGGCTAAGACGTCGTGCGAGGAAAACAAGGTGATTGACCAAAATCGAAGTTACGAACAAGAAAGCGTCGA ~/INDEX/bowtie2Index/ecoli3/ecoli3
bowtie2-build ~/GENOME/ecoli3/ecoli3.fasta ~/INDEX/bowtie2Index/ecoli3/ecoli #ecoli3.fasta是三段序列的汇总
bowtie2-build ~/GENOME/mouse/mm10.fa ~/INDEX/bowtie2Index/mouse/mm10
cat $config1 | while read line
do
arr=($line)
fq1=${arr[1]}
fq2=${arr[2]}
sample=${arr[0]}
bowtie2 --end-to-end --very-sensitive --no-mixed --no-discordant --phred33 -I 10 -X 700 \
-p 32 -x ${ref} -1 $fq1 -2 $fq2 \
-S ${projPath}/alignment/sam/$sample.sam 2> ${projPath}/alignment/sam/bowtie2_summary/$sample.txt
bowtie2 --end-to-end --very-sensitive --no-mixed --no-discordant --phred33 -I 10 -X 700 \
-p 32 -x ${spikeInRef} -1 $fq1 -2 $fq2 \
-S $projPath/alignment/sam/$sample.spikein.sam 2> $projPath/alignment/sam/bowtie2_summary/$sample.spikein.txt
bowtie2 --end-to-end --very-sensitive --no-mixed --no-discordant --phred33 -I 10 -X 700 \
-p 32 -x ${spikeInRef1} -1 $fq1 -2 $fq2 \
-S $projPath/alignment/sam/$sample.spikein1.sam 2> $projPath/alignment/sam/bowtie2_summary/$sample.spikein1.txt
bowtie2 --end-to-end --very-sensitive --no-mixed --no-discordant --phred33 -I 10 -X 700 \
-p 32 -x ${spikeInRef2} -1 $fq1 -2 $fq2 \
-S $projPath/alignment/sam/$sample.spikein2.sam 2> $projPath/alignment/sam/bowtie2_summary/$sample.spikein2.txt
bowtie2 --end-to-end --very-sensitive --no-mixed --no-discordant --phred33 -I 10 -X 700 \
-p 32 -x ${spikeInRef3} -1 $fq1 -2 $fq2 \
-S $projPath/alignment/sam/$sample.spikein3.sam 2> $projPath/alignment/sam/bowtie2_summary/$sample.spikein3.txt
done
为了进行 spike-in 标准化，可以增加两个参数：
--no-overlap --no-dovetail，以避免实验基因组与用于校准的大肠杆菌DNA的交叉比对。
用这两个参数，虽然三段 spike-in 之和有 reads 数，但 spike-in1号 直接被筛归零，这里就不用了。
3.1.3 比对结果
有关更详细的参数说明，用户可以参考bowite2手册。
${projPath}/alignment/sam/bowtie2_summary/$sample.txt
- 30319290 reads; of these:
- 30319290 (100.00%) were paired; of these:
- 4850287 (16.00%) aligned concordantly 0 times
- 20208420 (66.65%) aligned concordantly exactly 1 time
- 5260583 (17.35%) aligned concordantly >1 times
- 84.00% overall alignment rate
———
30319290 是测序深度，即 paired reads 的总数。
4850287 是没比对上的 read pairs 数。
20208420 + 5260583 是成功比对上的 read paris 数。
84.00% 是总体比对率。
也可以用 samtools flagstat ${projPath}/alignment/bam/$sample.mapped.bam 计算 reads 数：
- 50938006 + 0 in total (QC-passed reads + QC-failed reads)
- 0 + 0 secondary
- 0 + 0 supplementary
- 0 + 0 duplicates
- 50938006 + 0 mapped (100.00% : N/A)
- 50938006 + 0 paired in sequencing
- 25469003 + 0 read1
- 25469003 + 0 read2
- 50938006 + 0 properly paired (100.00% : N/A)
- 50938006 + 0 with itself and mate mapped
- 0 + 0 singletons (0.00% : N/A)
- 0 + 0 with mate mapped to a different chr
-
0 + 0 with mate mapped to a different chr (mapQ>=5)
——
25469003 才是跟上面一样，成功比对上的 read paris 数：20208420 + 5260583。
3.4 评估比对上的片段大小分布【必做】
mkdir -p $projPath/alignment/sam/fragmentLen
cat $config1 | while read line
do
arr=($line)
sample=${arr[0]}
## Extract the 9th column from the alignment sam file which is the fragment length
samtools view -F 0x04 $projPath/alignment/sam/$sample.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/alignment/sam/fragmentLen/$sample.fragmentLen.txt
done
4 比对结果筛选和文件格式转换
4.1 通过 mapping 质量筛选 mapped reads [可选]
某些项目可能需要对比对质量分数进行更严格的筛选。这篇博客详细讨论了 bowtie 如何计算质量分数的。
MAPQ(x) = -10 * log10(P(x is mapped wrongly)) = -10 * log10(p)，MAPQ范围从 0 到 37、40 或 42。
samtools view -q minQualityScore将过滤低于用户定义的minQualityScore的所有比对结果。
4.2 文件格式转换【必做】
3.3 去除重复【选做/必做】
CUT&Tag 将接头整合到 antibody-tethered pA-Tn5 附近的 DNA 中，整合的确切位点受到周围 DNA 可及性的影响。出于这个原因，具有相同的起始和结束位置的片段应该是常见的，并且这种“重复”可能不是 PCR 过程中产生的重复。在实践中，我们发现高质量的 CUT&Tag 数据集的 apparent 重复率很低，即使是 apparent “重复”片段也可能是真正的片段。因此，我们不建议删除重复 reads。在使用非常少材料或者怀疑是PCR重复的时候，可以去除重复。以下命令显示如何使用 Picard 检查重复率。
官网不建议 rmdup，但这里因为样本原因是跑了15个循环，而且我们很明确是pcr重复，直接用 mapped 数据跑出来的图后面有很多“杂质”。这里转化成 bam 文件后再做的去重，方便。一定要记得两种排序的转化！不然就会一直报错，得不到 PE 形式的结果（那还做什么 fragment！我猜的）。还有之前傻子似的把spike-in做了去重，千万不要！本来就是pcr产物，去了就没了！！
picardCMD="java -jar ～/.conda/envs/cuttagg/share/picard-2.18.29-0/picard.jar"
cat $config1 | while read line
do
arr=($line)
sample=${arr[0]}
## Filter and keep the mapped read pairs
samtools view -bS -F 0x04 $projPath/alignment/sam/$sample.sam >$projPath/alignment/bam/$sample.mapped.bam #这里加-q参数，完成步骤4.1
## 坐标排序为 rmDup 做准备
samtools sort -o $projPath/alignment/bam/$sample.sorted.bam $projPath/alignment/bam/$sample.mapped.bam
## remove duplicates
$picardCMD MarkDuplicates \
REMOVE_DUPLICATES=true \
I=$projPath/alignment/bam/$sample.sorted.bam \
O=$projPath/alignment/bam/$sample.sorted.rmDup.bam \
M=$projPath/alignment/bam/picard_summary/$sample.sorted.rmDup.txt
## 这里一定一定要再按名称排序一次，不然会报错几个G！因为后面 -bedpe 参数只能识别这种排序！！！
samtools sort -n -o $projPath/alignment/bam/$sample.name.rmDup.bam $projPath/alignment/bam/$sample.sorted.rmDup.bam
## Convert into bed file format
bedtools bamtobed -i $projPath/alignment/bam/$sample.name.rmDup.bam -bedpe >$projPath/alignment/bed/$sample.bed
## Keep the read pairs that are on the same chromosome and fragment length less than 1000bp.
awk '$1==$4 && $6-$2 < 1000 {print $0}' $projPath/alignment/bed/$sample.bed >$projPath/alignment/bed/$sample.clean.bed
## Only extract the fragment related columns
cut -f 1,2,6 $projPath/alignment/bed/$sample.clean.bed | sort -k1,1 -k2,2n -k3,3n >$projPath/alignment/bed/$sample.fragments.bed
done
4.3 评估重复性
为了研究重复之间和不同条件下的可重复性，将基因组分成 500 个 bp 的 bin【取中点，再算上前后250bp】，计算每个 bin reads计数的 log2 转换值的 （在重复数据集之间）Pearson 相关性。
binLen=500
mkdir -p $projPath/alignment/bed/binLen
cat $config1 | while read line
do
arr=($line)
sample=${arr[0]}
## We use the mid point of each fragment to infer which 500bp bins does this fragment belong to.
awk -v w=$binLen '{print $1, int(($2 + $3)/(2*w))*w + w/2}' $projPath/alignment/bed/$sample.fragments.bed | sort -k1,1V -k2,2n | uniq -c | awk -v OFS="\t" '{print $2, $3, $1}' | sort -k1,1V -k2,2n >$projPath/alignment/bed/binLen/$sample.fragmentsCount.bin$binLen.bed
done
5 Spike-in 校正
使用常量 C 来避免归一化数据中出现小分数，我们将标准化系数 scaling factor: S 定义为
S = C / (fragments mapped to E. coli genome) 这个常数是一个任意乘数，通常为 10,000。
归一化覆盖率的计算公式为：Normalized coverage = (primary_genome_coverage) * S
chromSize="/lustre/home/yfxie/GENOME/mouse/mm10.chrom.sizes.txt" # UCSC latest/下载
mkdir -p ${projPath}/alignment/bedgraph
mkdir -p $projPath/alignment/bigwig
cat $config1 | while read line
do
arr=($line)
sample=${arr[0]}
samtools sort -o $projPath/alignment/bam/$sample.sortedbw.rmDup.bam $projPath/alignment/bam/$sample.sorted.rmDup.bam
samtools index $projPath/alignment/bam/$sample.sortedbw.rmDup.bam
seqDepthDouble=`samtools view -F 0x04 $projPath/alignment/sam/$sample.spikein.sam | wc -l`
seqDepth=$((seqDepthDouble/2))
echo $seqDepth >$projPath/alignment/sam/bowtie2_summary/$sample.spikein.seqDepth
if [[ "$seqDepth" -gt "1" ]]; then
scale_factor=`echo "10000 / $seqDepth" | bc -l`
echo "Scaling factor for $sample is: $scale_factor!"
bedtools genomecov -bg -scale $scale_factor -i $projPath/alignment/bed/$sample.fragments.bed -g $chromSize > $projPath/alignment/bedgraph/$sample.fragments.normalized.bedgraph
bamCoverage --scaleFactor $scale_factor -b $projPath/alignment/bam/$sample.sortedbw.rmDup.bam -o $projPath/alignment/bigwig/$sample.normalized.bw
fi
done
这里直接通过 spike-in 校正因子转换成 bw 文件（--bamCoverage）也 ok!
得到 bedgraph 文件（--genomecov）了，下回见！Peak calling。
