数据分析-cuttag分析流程分享1-linux代码流程分析
原创数据分析-cuttag分析流程分享1-linux代码流程分析
原创
老板最近比较痴迷于各种seq，由于俩师姐外加一师妹的chip-seq建库老不成功，于是改成了CUTTAG建库，其实读了文献，发现都是相似的原理，只不过是CUTTAG用的细胞起始量远低于chip-seq，同时用的酶不同，建库的时间相对少很多。具体可以参考一下这篇文献CUT&Tag for efficient epigenomic profiling of small samples and single cell（doi:10.1038/s41467-019-09982-5）。
文章里面对CUTTAG与其他几种seq建库的结果进行了比较，并且在文末进行了单细胞cuttag建库尝试，获得了相对比较好的结果。
鉴于近期分析了大批量的3个不同物种的cuttag数据，准备分享一下相关的代码。
数据完整性检测
首先是需要对测序公司给的测序数据的完整性进行test，如果md5值不吻合，后面还需要让公司发一份。
//##multiqc查看数据质量
multiqc *.fq.gz -o multiqc ##multiqc会自动寻找指定文件下面fq的结果文件，-o 指定输出目录
##检查双端测序的行数是否相等
ls *.fq.gz
zcat *.fq.gz | wc -l
md5sum -c MD5.txt或者可以选用将一个大的文件目录下面的fq文件都进行md5值校正
// ##将文件下面的fq文件全部计算MD5值
find ./file/* -type f -print0 |xargs -0 md5sum>md5.txt参考基因组下载与构建索引
一般是去自己比较喜欢的数据库上下载自己需要的参考基因组，因为我是做植物的，所以喜欢去phytozome网站(https://phytozome-next.jgi.doe.gov/)上下载基因组和注释文件，主要以拟南芥基因组构建索引为例：(由于个人比较懒，很多代码都是在后台运行的，而且一般都是test完了，全部都是在后台跑循环，鉴于是为了对这个分析的坑更好的说清楚，就对每个shell脚本拆开运行。如果有友友们不喜欢后台运行，可以把nohup和&去掉，就可以在虚拟机界面看到啦！)
// ##参考基因组的构建主要选用的是bowtie2软件进行构建
nohup bowtie2-build --threads 8 TAIR10.fa TAIR10_index > TAIR10_bowtie2_index.log 2>&1 &
##同时官网推荐要用大肠杆菌的基因组进行校准，所以我在分析中也加入了相关的校准，但是发现比对率很低，通常是低于0.1%。
nohup bowtie2-build --threads 8 Ecoli.fna Ecoli &数据指控与标准化
大部分我所选用的代码都是cuttag文章分析流程推荐的代码（https://yezhengstat.github.io/CUTTag_tutorial/#VII_Visualization），但是官网上有部分内容不太适合我们分析的内容，所以做了一些微小的改动。
//## -I 是最短序列 -X 是最长序列 -p 线程 -x 指定参考基因组的index -1 序列1文件 -2 序列2文件 -S 生成的sam文件
##由于文件比较多，可以写一个简单的for循环的shell脚本放到文件目录下面（ for i in `cat file_list` ; do done）
##比对到研究物种的参考基因组上，其中*-rep*是加了研究抗体的，IgG_rep1是IgG，用来call峰的时候去除背景的，因为目前只做了一个样本，4个重复，因此后面是没有办法做差异peak的，所以在后面用基因注释的时候选用的是全部的peak做的基因注释。
for i in D_rep1 D_rep2 D_rep3 D_rep4 IgG_rep1;do
{
cores=8
ref="/nnj/TAIR10"
bowtie2 --end-to-end --very-sensitive --no-mixed \
--no-discordant --phred33 -I 10 -X 700 -p ${cores} \
-x ${ref} -1 ${i}_1.clean.fq.gz -2 ${i}_2.clean.fq.gz \
-S ${i}_bowtie2.sam &> ${i}_bowtie2.txt
}&
done
##比对到大肠杆菌的基因组上
for i in D_rep1 D_rep2 D_rep3 D_rep4 IgG_rep1;do
{
cores=8
spikeInRef="/Ecoli/Ecoli"
bowtie2 --end-to-end --very-sensitive --no-mixed \
--no-discordant --phred33 -I 10 -X 700 -p 8 \
-x ${spikeInRef} -1 ${i}_1.clean.fq.gz \
-2 ${i}_2.clean.fq.gz -S ${i}_bowtie2_spikeIn.sam &> ${i}_bowtie2_spikeIn.txt ##进行比对
seqDepthDouble=`samtools view -F 0x04 ${i}_bowtie2_spikeIn.sam | wc -l`
seqDepth=$((seqDepthDouble/2))
echo $seqDepth >${i}_bowtie2_spikeIn.seqDepth
}&
done随后给予权限，一般都是chmod 775 bidui.sh
// ##后台运行shell脚本
nohup ./bidui.sh &>bidui.out 2>&1 &官网这个时候就准备对测序的一些数据质量进行可视化了，但是由于我比较懒，所以准备在后面贴R语言的代码，前面的话可以通过shell脚本运行出来的几个txt文件，大致的知道自己测序文件的mapping率及数据质量，一般是mapping率大于80%，可以算是比较不错的测序数据，一般建议可以接着往下分析。
文件格式的转换，获得bed文件
// ##文件格式转换
## Filter and keep the mapped read pairs
## 筛选和保留比对上的双端 reads
for i in D_rep1 D_rep2 D_rep3 D_rep4 IgG_rep1; do
{
samtools view -bS -F 0x04 ${i}_bowtie2.sam > ${i}_bowtie2.mapped.bam
## Convert into bed file format
## 将 BAM 文件转换为 bed 文件格式
bedtools bamtobed -i ${i}_bowtie2.mapped.bam -bedpe > ${i}_bowtie2.bed
## Keep the read pairs that are on the same chromosome and fragment length less than 1000bp.
## 保留那些在同一条染色体且片段长度小于 1000bp 的双端 reads
awk '$1==$4 && $6-$2 < 1000 {print $0}' ${i}_bowtie2.bed > ${i}_bowtie2.clean.bed
## Only extract the fragment related columns
## 仅提取片段相关的列
cut -f 1,2,6 ${i}_bowtie2.clean.bed | sort -k1,1 -k2,2n -k3,3n > ${i}_bowtie2.fragments.bed
}&
donechmod 775 zhuanhuan.sh
// ##后台运行shell脚本
nohup ./zhuanhuan.sh &>zhuanhuan.out 2>&1 &对重复之间的相关性进行比较
// ##评估重复性
##为了研究重复之间和不同条件下的可重复性，基因组被分成 500 bp bin，每个 bin reads 计数的 log2 转换值的皮尔逊相关性在重复数据集之间计算。多个重复和 IgG 对照数据集显示在层次化聚类关联矩阵中。
for i in D_rep1 D_rep2 D_rep3 D_rep4 IgG_rep1; do
{
binLen=500
awk -v w=$binLen '{print $1, int(($2 + $3)/(2*w))*w + w/2}' ${i}_bowtie2.fragments.bed |\
sort -k1,1V -k2,2n |\
uniq -c |\
awk -v OFS="\t" '{print $2, $3, $1}' |\
sort -k1,1V -k2,2n >${i}_bowtie2.fragmentsCount.bin$binLen.bed
}&
donechmod 775 chongfuxing.sh
// ##后台运行脚本
nohup ./chongfuxing.sh &>chongfuxing.out 2>&1 &加入测序深度对bed文件进行标准化处理
// #标准化
samtools faidx ref.genome ##后面的bedtools -g的内容通过samtools的faidx来进行构建获得索引文件
##这里的seqDepth是前面的数据比对处理当中获得测序深度值，可以调取txt文件，来进行后续的循环处理
for i in D_rep1 D_rep2 D_rep3 D_rep4 IgG_rep1; do
{
if [[ "$seqDepth" -gt "1" ]]; then
scale_factor=`echo "10000 / $seqDepth" | bc -l`
echo "Scaling factor for $histName is: $scale_factor!"
bedtools genomecov -scale $scale_factor -i {i}_rep1_bowtie2.fragments.bed -g /nnj/ref.genome > {i}.fragments.normalized.bedgraph
fi
}&
donechmod 775 biaozhunhua.sh
// 后台运行shell脚本
nohup ./biaozhunhua.sh &>biaozhunhua.out 2>&1 &现在已经通过{i}.fragments.normalized.bedgraph文件在IGV上对call峰的结果进行可视化了，不推荐bam文件，主要原因是bam文件相比于bedgraph文件较大，而且bedgraph文件还可以用notepad打开，查看相对应峰的具体位置，减少对windows的运行空间。如果是用的linux系统的IGV，那就不需要考虑这些问题啦，最近发现了一个远程控制服务器的软件，MobaXterm，充分满足了我这个懒人，不想配置X11的想法，下载完了，就可以可视化了，后续用R出可视化的图，也不用在加载ggplot2来保存pdf文件。
call peak
call peak是选用的官网上推荐的SEACR_1.3.sh，需要提前对这个文件进行下载，并加入环境变量中，来调取。
// ##Peak calling
##== linux 命令==##
for i in D_rep1 D_rep2 D_rep3 D_rep4 ; do
{
seacr="/biotools/SEACR-master/SEACR_1.3.sh"
histControl=$4
bash $seacr ${i}_bowtie2.fragments.normalized.bedgraph \
IgG_rep1_bowtie2.fragments.normalized.bedgraph \
non stringent ${i}_seacr_control.peaks
bash $seacr ${i}_bowtie2.fragments.normalized.bedgraph \
0.01 non stringent ${i}_seacr_top0.01.peaks
}&
done其中获得两个结果，一个是与IgG进行对比后，去除了IgG的背景（${i}_seacr_control.peaks），一个文件峰值的数量在top0.01的文件（${i}_seacr_top0.01.peaks）
chmod 775 peakcalling.sh
// 后台运行shell脚本
nohup ./peakcalling.sh &>peakcalling.out 2>&1 &总结
目前是通过这些shell脚本获得了可以可视化的bedgraph文件，可以初步确定自己做的蛋白与DNA的结合位置是否准确，同时还获得了一些峰值的结果，可以大致看到cuttag建库的背景情况，那下一篇对R语言可视化（数据分析-cuttag分析流程分享2-R代码可视化流程处理）的代码进行相关的整理，可以更清楚的看到相关的结果。
其实可以发现以上的流程都是可以放到一个大的for循环当中进行后台流程分析，主要是需要在前期把需要的软件下载好，并给予环境变量，在进行shell命令的时候，可以找到软件的路径，同时还要提前将各个基因组的索引文件构建好，不至于在跑到一半的时候命令卡住，有利于后续的下游分析。目前按照我的经验来看，如果测序获得fq文件不大，上游的这些流程分析大约可以在2d内拿到相应的结果，来去做后面的个性化分析。
也有相关的流程建议是跑完bam文件的时候进行picard的去除，但是在我第一次test文件后，发现有很多的峰都去掉了，怀疑是cuttag数据不需要去除，因此查看了几个博主的文章，有推荐的，也有不推荐的，主要还是考虑我们的测序数据的结果，同时也问了几个做分析很厉害的师兄，目前也是不建议去重的，所以目前这个去重的流程分析还是需要针对项目的数据质量来进行后续分析的。
原创声明：本文系作者授权腾讯云开发者社区发表，未经许可，不得转载。
如有侵权，请联系 cloudcommunity@tencent.com 删除。
原创声明：本文系作者授权腾讯云开发者社区发表，未经许可，不得转载。
如有侵权，请联系 cloudcommunity@tencent.com 删除。