m6A图文复现05-比对结果以及可视化

2022-01-06 本文已影响0人信你个鬼

在上一期中我们得到hisat2的比对结果，现在一起来看下比对情况。

一，比对结果解读

比对结果文件如下，每一个样本一个排序后的bam文件，bam的索引bai以及存有比对信息的log日志。

image

选取一个样本的log查看比对详细情况：

cat SRR866991.log

58186286 reads; of these:
  58186286 (100.00%) were unpaired; of these:
    4153900 (7.14%) aligned 0 times
    35063517 (60.26%) aligned exactly 1 time
    18968869 (32.60%) aligned >1 times
92.86% overall alignment rate

我们可以看到这个样本的总比对率为92.86%，这是非常重要的一个指标，一般来说，样本的总比对率需要达到80%以上，70%多也能接受。如果你的总比对率太低如20%以下，有可能有以下几种情况：

1.参考基因组序列与fastq文件中的reads不是来源于同一个物种，要么是参考基因组用错了，要么是fastq文件搞错了。
2.序列中含有其他序列如细菌等污染。

这个时候可以对你的fastq序列抽取其中一部分reads做NCBI的nr/nt库的blast比对来确定fastq中的序列来源。

比对log中其他信息：

第一行：58186286 reads; of these，统计了这个样本的总reads数
第二行：58186286 (100.00%) were unpaired; of these:，read是单端
第三行：4153900 (7.14%) aligned 0 times，没有比对上参考基因组序列的read数和比例
第四行：35063517 (60.26%) aligned exactly 1 time，比对上参考基因组并且只比对上一个位置的序列。
第五行：18968869 (32.60%) aligned >1 times，比对上参考基因组且比对到多个位置。
第六行：92.86% overall alignment rate，总比对率

清楚了每一行的信息之后，我们来算一下总比对率是如何算的：

>((58186286-4153900)/58186286)*100 = 92.86%

也就是说总比对率其实包含了多比对read的部分。

将我们的log日志进行整理得到总比对的表格或者一个html文档

># multiqc可以将log日志中的数据进行可视化方便一起查看多个样本的比对结果。
multiqc *log

image

所有样本的总比对，有一个样本的比对率有些低，可以回去查查这个样的的质控等情况。

image

详细比对结果：

image

综合以上的情况，SRR866993号样本在这里的结果不是特别好，后面分析的有效数据会比较少。

二，bam转bw文件并绘热图

现在我们来将bam文件转成bw文件，可以使用IGV来查看比对结果，以及peak可视化（IGV教程参考隔壁公众号：保姆级 IGV 基因组浏览器使用指南（图文详解）https://mp.weixin.qq.com/s/Dtin-r85QRBexR2Rik9p2g）。

bam转bw有代码资料的那篇文献中代码使用的是最简单的参数：

image

这里我们加上--centerReads：可获得更陡峭的峰

># bam2bw
ls ../mapping/*bam | while read id
do
sam=${id##*/};sam1=${sam%%.*};
nohup bamCoverage -p 6 -b ${id} --centerReads --ignoreDuplicates -o ${sam1}.bw &
done

# 制作gene的bed文件
cat Mus_musculus.GRCm38.104.gtf |awk 'BEGIN{ FS=" ";OFS="\t" }{if($3 ~/gene/) print $1,$4-1,$5-1,$10,$6,$7}'|sed -e s'/"//'g -e s'/;//'g >Mus_musculus.GRCm38.104.bed

# computeMatrix
computeMatrix scale-regions -p 8 \
-b 3000 -a 3000  \
-R Mus_musculus.GRCm38.101.bed \
--regionBodyLength 6000 \
-S SRR*.bw \
-o All_Sample.TSS-body-TES_distribute.gz 

# plot heatmap
plotHeatmap --dpi 200 --heatmapHeight 15 \
--matrixFile All_Sample.TSS-body-TES_distribute.gz \
--outFileName All_Sample.TSS-body-TES_Heatmap.png