如何高效地从BAM文件中提取fastq

Posted

tags:

篇首语:本文由小常识网(cha138.com)小编为大家整理,主要介绍了如何高效地从BAM文件中提取fastq相关的知识,希望对你有一定的参考价值。

参考技术A

在一年前,我写过一篇文章,叫做 如何从BAM文件中提取fastq ,之前也发现了从BAM里面提取Fastq是有些麻烦,只不过最后通过 samtools 的子命令实现了数据提取,实现功能之后也没有再去思考如何提高效率。

最近读到 每周文献-190419-植物单细胞BAM重比对以及假基因研究 时,发现里面提到了一个工具叫做 bazam , 功能就是提取Fastq文件,文章发表在 Genome Biology 上。

软件地址为 https://github.com/ssadedin/bazam ,因为他是个Groovy工具(据说Groovy比Java好写)所以安装很方便,

可以通过如下命令检查是否安装成功

假如你原来的BAM文件是双端测序,想提取成两个文件,那么可以用如下命令。 :这里的your.bam 指的是你实际BAM文件

此外,bazam还支持使用 -L 指定区间来提取某个区间的数据。

这里补充下为啥从BAM文件中提取双端序列那么麻烦,这是因为一般而言BAM文件都是按照位置信息排序,想要找到配对的reads,要么是根据read的编号进行排序(这个方法要求额外的内存和存储空间),或者就是在提取的时候记录当前的read的ID,再找到另一端ID后释放内存空间。

Pysam 处理bam文件

Pysam可用来处理bam文件

安装:

用 pip 或者 conda即可

 

使用:

Pysam的函数有很多,主要的读取函数有:

  • AlignmentFile:读取BAM/CRAM/SAM文件

  • VariantFile:读取变异数据(VCF或者BCF)

  • TabixFile:读取由tabix索引的文件;

  • FastaFile:读取fasta序列文件;

  • FastqFile:读取fastq测序序列文件

一般常用的是第一个和第二个。

 

例子:

import pysam

bf = pysam.AlignmentFile("in.bam","rb");  其中r = read, b:binary.  二进制文件。   bam文件index

bf是一个迭代器,可以next()或者for读取

for i in bf:

    print i.reference_name,i.pos,i.mapq,i.isize

结果:

ctg000331_np121 144935 27 -284
ctg000331_np121 144940 48 291
ctg000331_np121 144941 48 309
ctg000331_np121 144944 48 255
ctg000331_np121 144946 27 -370
ctg000331_np121 144947 27 -346

  • i.reference_name代表read比对到的参考序列染色体id;

  • i.pos代表read比对的位置;

  • i.mapq代表read的比对质量值;

  • i.isize代表PE read直接的插入片段长度,有时也称Fragment长度;

 

一些功能:

  • check_index() 

          检测index文件是否存在存在即为true

  • close()

          用完记得关闭

  • count(self,contig=None, start=None, stop=None, region=None, until_eof=False, read_callback=‘nofilter‘, reference=None,end=None)
             计算目标区域内比对上的reads数目

           bf.count(contig="ctg000331_np121", start=1, stop=6000)
           24

  • count_coverage(self, contig=None, start=None, stop=None, region=None, quality_threshold=15, read_callback=‘all‘, reference=None, end=None)
         计算目标区域内的覆盖度。返回1个4维的array,代表ACGT的覆盖度,而每个维度的array长度为100,里面的数字代表该碱基在各个位置上的覆盖度。

       bf.count_coverage(contig="ctg000331_np121",start=1,stop=100)
  • fetch(self, contig=None, start=None, stop=None, region=None, tid=None, until_eof=False, multiple_iterators=False, reference=None, end=None)
          提取出比对到目标区域内的全部reads。返回的是一个迭代器,可以通过for循环或者next函数从中取出reads,我们使用next()函数取出第一条reads,reads是用         AlignedSegment对象表示,可以通过该对象的内置方法再对这条reads进行一些查询操作。
          allreads=bf.fetch(contig="ctg000331_np121",start=1,stop=10000)
   是一个迭代器,可以用for循环获得
  • get_index_statistics(self)
    通过index统计该BAM文件中在各个染色体上mapped/unmapped的reads个数。
bf.get_index_statistics()
 

以上是关于如何高效地从BAM文件中提取fastq的主要内容,如果未能解决你的问题,请参考以下文章

perl从2个数组中提取常见元素(fastq文件中的常见序列)

2021-04-22 bam文件中提取比对(mapped)或未比对上(unmapped)reads

python 此代码将从执行它的文件夹中的所有fastq的前1000行中提取标题信息。然后它需要

如何正确地从主仓库中提取,更新本地文件,然后将更改推送到自己的仓库?

sh 提取所有FASTQ标头

bam比对flag说明以及提取未必对上的reads