测序数据饱和度分析

TIME_@

于 2022-07-13 19:44:22 发布

阅读量2.5k

点赞数 2

分类专栏：生物信息文章标签：大数据

本文链接：https://blog.csdn.net/geekfocus/article/details/125724446

版权

生物信息专栏收录该内容

50 篇文章 27 订阅

订阅专栏

ref1
ref2

二代测序获得数据后，正式分析前通常需要：reads的饱和度分析，评估RNA-seq测序的数据与mRNA真实表达水平之间的一致性。
测序深度和基因组覆盖率为正相关关系。
影响：

DNA序列：测序不饱和，影响拼接
RNA-seq：定量不准
宏基因组：不能反映物种组成
chip，atac…
评估方式一：gene number 与reads number关系

抽样bam

测序数据量对于NGS分析非常重要，测序数据量过低，不能有效覆盖基因组完整信息，测序数据量过高，会造成冗余，不够经济。为了验证当前测序量能否满足需求，或者说加大测序量是否能够进一步挖掘信息，通常需要进行饱和度分析。
ATACseqQC 作为一个专门针对ATAC文库进行QC的R包
如上图所示，ATACseqQC的饱和度分析思路如下，对effectice fragment进行随机抽样，比如10%，20%到100%，对于每个梯度分别进行peak calling, 统计其peak个数，然后以effectice fragment数目为横坐标，peak个数为纵坐标绘制散点图，并进行拟合，如果拟合曲线的末端斜率上升平缓，表明继续加大测序数据量，新识别到的peak个数也不会新增太多，就可以认为当前的测序数据量是比较饱和的，不比再加测数据量了。如果曲线末端仍处于一个斜率急剧上升的阶段，说明测序数据量不够，还需加大测序数据量。
在这里插入图片描述

理解了分析思路，再来看一下具体的操作过程。第一步是对effectice fragment进行抽样。对原始bam文件过滤，去除MAPQ较低, 去除PCR重复，去除比对到线粒体的序列，filter后的bam文件。

对effectice fragment进行抽样，本质就是操作这个bam文件。通过linux下的shuf命令对sam文件随机抽取对应的行，从而实现对effectice fragment抽样的目的，基本思路如下

提取sam文件的header部分作为一个单独的文件；
使用linux内置的shuf命令对sam文件的正文部分随机抽取一定比例的行数，追加到header文件中，生成一个新的sam文件
在文章的附件中提供了一个bash脚本，用于对bam文件进行抽样，用法如下

利用排序好的bam文件可以直接进行peak calling, 然后统计每个抽样百分比对应的effectice fragment数目和peak个数，就可以绘制测序饱和度图了。在ATACseqQC中，提供了saturationPlot函数，可以直接读取一系列peak文件，然后绘制测序饱和度图，具体的用法可以查看该函数的帮助文档。
结果

GL1.bam.name.sorted.0.1.sam
GL1.bam.name.sorted.0.1.coord.sorted.bam
GL1.bam.name.sorted.0.1.coord.sorted.bam.bai

ATACseqQC给我们提供了一个很好的思路，来进行ATAC文库的测序饱和度分析，类似的，这个做法也可以推广到chip_seq, m6A_seq等IP类型的实验中去。

shuf命令 – 产生随机的排列

The shuf command in Linux writes a random permutation of the input lines to standard output. It pseudo randomizes an input in the same way as the cards are shuffled. It is a part of GNU Coreutils and is not a part of POSIX. This command reads either from a file or standard input in bash and randomizes those input lines and displays the output. Like other Linux commands, shuf command comes with –help option.

Syntax:

shuf [OPTION] [FILE] //file shuf
shuf -i LO-HI [OPTION] // range shuf
shuf -e [OPTION]… [ARG] //list shuf

https://www.linuxcool.com/shuf
模拟体彩超级大乐透：

shuf -i 1-35 -n 5|sort -n && shuf -i 1-12 -n 2|sort -n
4
17
20
29
31
6
11

shuf -n  [number of sub reads] sam  >>  sub.sam

((i=$j+$k))    等价于 i=`expr $j + $k`
((i=$j-$k))     等价于   i=`expr $j -$k`
((i=$j*$k))     等价于   i=`expr $j \*$k`
((i=$j/$k))     等价于   i=`expr $j /$k`