是用pip下载的pysam,为了做index文件用brew install下载了samtools。
import pysam
读取bam文件如下。
n1 = pysam.AlignmentFile('./atac_bam/scN1.pe.q10.sort.rmdup.bam', "rb")
为了能够读取特定位置的read,用samtools做index文件。
samtools index scN1.pe.q10.sort.rmdup.bam
然后用fetch的方法读取行。fetch返回的是一个迭代器(iterator),可以迭代读取内容。
fetch的API如下:
fetch(self, reference=None, start=None, end=None, region=None, tid=None, until_eof=False, multiple_iterators=False)
for read in n1.fetch('chr6'):
print read
n1.close()