69. 《Bioinformatics Data Skills》

使用pysam包读取BAM/SAM文件返回的AlignmentFile对象包含丰富的属性。

header属性

header以字典形式存储，所有的键为：

>>> bamfile.header.keys()
['HD', 'SQ', 'RG', 'PG']

例如RG（read group）的内容为：

>>> bamfile.header['RG'][0]
{'LB': 'g1k-sc-NA12891-CEU-1', 'CN': 'SC', 'DS': 'SRP000032', 'SM': 'NA12891', 'PI': '200', 'ID': 'ERR001776', 'PL': 'ILLUMINA'}

染色体序列为：

>>> bamfile.header['SQ'][0]
{'LN': 249250621, 'M5': '1b22b98cdeb4a9304cb5d48026a85128', 'AS': 'NCBI37', 'SN': '1', 'UR': 'file:/lustre/scratch102/projects/g1k/ref/main_project/human_g1k_v37.fasta'}

方便起见，不通过header也可直接得到参考染色体及其长度：

>>> bamfile.references
('1', '2', '3', '4', '5', '6', '7', '8', '9', '10', '11', '12', '13', '14', '15', '16', '17', '18', '19', '20', '21', '22', 'X', 'Y', 'MT', ...)

>>> bamfile.lengths
(249250621, 243199373, 198022430, 191154276, 180915260, 171115067, 159138663, 146364022, 141213431, 135534747, 135006516, 133851895, 115169878, ...)

read属性

遍历AlignmentFile返回的read对象也有自己的属性，比如上节提到的qname。

>>> read = bamfile.next()

tid或者reference_id属性代表此read比对的参考序列ID，这个值为非负整数。以第一个read为例，首先确定read是能够比对到参考序列的（否则没有意义）：

>>> read.is_unmapped
False

靶向染色体ID为：

>>> read.reference_id
0
>>> read.tid
0

通过AlignmentFile对象方法getrname可以得到靶向id对应的参考基因组名称：

>>> bamfile.getrname(read.reference_id)
'1'

我们也可以将参考基因组名称转为ID：

>>> bamfile.gettid("MT")
24

如果名称出错，将返回不可能存在的-1，基于此可以检验自己的输入是否正确。

>>> bamfile.gettid("Mt")
-1