本文来源于panxiaoguang.github.io
samtools是我们比较常用的处理bam文件的软件，但很多关于bam的具体信息，samtools不能给出，比如想获得某个ref位置，每条reads的碱基情况，及对应的readsID；获得比对到参考具体区间的碱基分布矩阵等。
Pysam内置了较多的函数，可以非常灵活的满足，各种个性化的关于bam文件信息的提取。
- AlignmentFile对象内置方法
- 2.AlignedSegmennt对象内置方法
- PileupColumn对象内置方法
- PileupRead对象内置方法
1 获得pysam.Alignment对象
import pysam
samfile = pysam.AlignmentFile("test.bam","rb")
#核查索引
samfile.check_index() #返回true
2 AlignmentFile对象内置方法
-
计算目标区域内比对上的reads数目
count(self，contig=None, start=None, stop=None, region=None, until_eof=False, read_callback=’nofilter’, reference=None，end=None)
#计算目标区域内比对上的reads数目
samfile.count(contig='chr1', start=500, stop=600) #注意是stop，而是end
-
计算目标区域里的覆盖度
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，里面的数字代表该碱基在各个位置上的覆盖度。
#计算目标区域里的覆盖度
samfile.count_coverage(contig='chr1', start=500, end=600)
-
提取出比对到目标区域内的全部reads
fetch(self, contig=None, start=None, stop=None, region=None, tid=None, until_eof=False, multiple_iterators=False, reference=None, end=None)
提取比对到目标区域的全部reads，返回的是迭代器，可使用next、for循环访问。reads用的是AlignedSegment对象表示，以下主要说reads的AlignedSengment对象的内置方法
allreads=sam.fetch(contig='chr1', start=500, end=600)
reads1=next(allreads)
find_introns(self, read_iterator)返回一个字典{(start, stop): count}，列出了reads中的intronic sites（cigar字符串中的’N’）
get_index_statistics(self)<u>通过index统计该BAM文件中在各个染色体上mapped/unmapped的reads个数</u>。
sam.get_index_statistics()
[IndexStats(contig='chr1', mapped=310538, unmapped=940, total=311478), IndexStats(contig='chr2', mapped=0, unmapped=0, total=0),
IndexStats(contig='chrM', mapped=0, unmapped=0, total=0)]
get_reference_name(self, tid)
-
根据染色体顺序，得到染色体名称
得到tid（数字）相对应的reference name（字符串），例如在这个BAM中tid为0代表chr1，tid为1代表chr2，tid为24代表chrM。
samfile.get_reference_name(0)'chr1'
'chr1'
与之对应的是<u>将染色体名称转为数字</u>，即将reference转换成tid
sam.get_tid('chr2')
1
pileup(self, contig=None, start=None, stop=None, region=None, reference=None, end=None, **kwargs)
类似于samtools的pileup操作，返回一个迭代器，每次迭代返回一个位点的PileupColumn对象,该对象的操作在后面进行详细介绍。
pileupcolumn = sam.pileup("chr1", 999863, 999865,truncate=True)
pc=next(pileupcolumn)
#stepper='nofilter', min_base_quality=0, min_mapping_quality=0
sam.references
('chr1', 'chr2', 'chr3', 'chr4', 'chr5', 'chr6', 'chr7', 'chr8', 'chr9', 'chr10', 'chr11', 'chr12', 'chr13', 'chr14', 'chr15', 'chr16', 'chr17', 'chr18', 'chr19', 'chr20', 'chr21', 'chr22', 'chrX', 'chrY', 'chrM')
#获得reads的cigar值
reads1.cigarstring
'2S151M'
cigar记录的元组格式，下面的例子就代表了151M（0 代表 M）
reads1.cigartuples
[(0, 151)]
Image.png
- 谨记reads在pysam包内部的对象名称是AlignedSegment
get_aligned_pairs(self, matches_only=False, with_seq=False)
展示了query和reference之间的比对情况，如果遇到indel或者skipping的情况，相应的query或者reference会变为None
# 比较两条reads的在ref上坐标的差值
reads.compare(reads2)
reads1.get_aligned_pairs(matches_only=False, with_seq=True)
#[(0, 999850, 'C'), (1, 999851, 'G'), (2, 999852, 'G')]
read1.get_blocks() #返回该reads比对到ref的起始与终止位置
#返回该reads比对到的ref区域的reference碱基
reads1.get_reference_sequence()
来源：
https://panxiaoguang.github.io/2019/03/06/Pysam%E5%A4%84%E7%90%86bam%E6%96%87%E4%BB%B6/
网友评论