测序数据饱和度分析

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`
  • 2
    点赞
  • 0
    收藏
    觉得还不错? 一键收藏
  • 0
    评论

“相关推荐”对你有帮助么?

  • 非常没帮助
  • 没帮助
  • 一般
  • 有帮助
  • 非常有帮助
提交
评论
添加红包

请填写红包祝福语或标题

红包个数最小为10个

红包金额最低5元

当前余额3.43前往充值 >
需支付:10.00
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包
实付
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明:

1.余额是钱包充值的虚拟货币,按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载,可以购买VIP、付费专栏及课程。

余额充值