如何高效地从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)
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)
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)
- get_index_statistics(self)
通过index统计该BAM文件中在各个染色体上mapped/unmapped的reads个数。
以上是关于如何高效地从BAM文件中提取fastq的主要内容,如果未能解决你的问题,请参考以下文章
perl从2个数组中提取常见元素(fastq文件中的常见序列)
2021-04-22 bam文件中提取比对(mapped)或未比对上(unmapped)reads
python 此代码将从执行它的文件夹中的所有fastq的前1000行中提取标题信息。然后它需要