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()
原文地址:https://www.cnblogs.com/zhanmaomao/p/11990448.html