生命情報科学

pysam

pysam は、TopHat などのマッピングプログラムが出力 SAM/BAM フォーマットファイルを取り扱うための Python ライブラリーである。詳細なドキュメントは pysam のウェブサイトに記載されている。

ファイルの読み込み

pysam を利用して SAM/BAM を読み込む際は pysam.AlignmentFile を利用してファイルを開き、for 文で一行ずつ回す方法で行う。pysam.AlignmentFile を利用してファイルを開くとき、SAM ならば読み込みモード "r" で開けばよく、BAM ならばバイナリーの読み込みモード "rb" で開けばよい。

import pysam
bamfile = pysam.AlignmentFile("SRR616268.bam", "rb")

for read in bamfile:
    print(read)

## SRR616268.11783	163	21	15722	40	93M	21	15758	93	CGTAGCCGAACTGAGAGGTTGATCGGCCACATTGGGACTGAGACACGGCCCAAACTCCTACGGGAGGCAGCAGTAGGGAATCTTCCACAATGG	array('B', [68, 70, 70, 70, 70, 70, 72, 72, 72, 72, 72, 72, 72, 69, 71, 69, 71, 72, 67, 70, 71, 72, 72, 72, 72, 71, 72, 72, 72, 72, 72, 71, 72, 72, 72, 72, 71, 72, 72, 72, 72, 72, 72, 70, 71, 72, 72, 72, 70, 70, 68, 68, 66, 67, 66, 67, 66, 66, 66, 66, 66, 66, 66, 66, 66, 66, 66, 66, 64, 66, 64, 66, 66, 66, 66, 66, 66, 66, 65, 64, 64, 65, 66, 66, 67, 66, 66, 66, 66, 66, 66, 66, 65])	[('AS', 0), ('XN', 0), ('XM', 0), ('XO', 0), ('XG', 0), ('NM', 0), ('MD', '93'), ('YS', -29), ('YT', 'CP')]
## SRR616268.11784	83	21	18352	42	4M3I100M	21	18275	107	ATTTCGGCAGATAAACTCCGAATACCATCGATCTTACTCCGGGAGTCAGACAGTGAGCGATAAGGTCCATTGTCGAAAGGGGAACAGCCCAGATCACCAGTTAAGGT	array('B', [68, 68, 68, 68, 65, 65, 65, 56, 66, 66, 66, 66, 65, 65, 61, 64, 61, 66, 66, 66, 67, 66, 65, 63, 66, 66, 64, 62, 66, 66, 66, 65, 63, 60, 67, 63, 66, 66, 70, 70, 70, 72, 71, 71, 70, 65, 72, 72, 72, 71, 72, 72, 72, 72, 72, 72, 72, 72, 72, 71, 69, 71, 72, 71, 70, 64, 71, 71, 70, 68, 69, 72, 72, 72, 72, 72, 72, 72, 72, 72, 72, 72, 72, 72, 71, 69, 72, 71, 71, 70, 69, 69, 70, 67, 68, 63, 70, 70, 70, 68, 68, 68, 68, 68, 65, 65, 65])	[('AS', -20), ('XN', 0), ('XM', 1), ('XO', 1), ('XG', 3), ('NM', 4), ('MD', '1A102'), ('YS', 0), ('YT', 'CP')]
## SRR616268.11784	163	21	18275	42	93M	21	18352	93	AGCTTTAGGGTTAGCCTCGGAGGATGGATCATGGAGGTAGAGCACTGTTTGAACTAGGGGCCCATCAAGGGTTACTGAATTCAGATAAACTCC	array('B', [68, 70, 70, 70, 70, 70, 71, 72, 72, 72, 72, 72, 72, 72, 72, 72, 72, 72, 72, 72, 72, 72, 72, 70, 71, 72, 72, 71, 72, 72, 72, 72, 72, 72, 69, 71, 72, 66, 70, 69, 70, 71, 72, 70, 71, 72, 72, 72, 72, 72, 72, 69, 71, 72, 70, 70, 70, 70, 70, 68, 68, 66, 66, 66, 66, 66, 66, 66, 66, 66, 66, 58, 63, 64, 66, 66, 66, 66, 66, 66, 67, 66, 66, 66, 66, 66, 66, 66, 66, 65, 66, 65, 65])	[('AS', 0), ('XN', 0), ('XM', 0), ('XO', 0), ('XG', 0), ('NM', 0), ('MD', '93'), ('YS', -20), ('YT', 'CP')]
## SRR616268.11785	99	8	1416	40	107M	8	1486	107	CTTCATTGTTTCAAAAAGATTATCTAACGTGACATGGTTGGTGTTCTTACGTTCTTCTTCCTGTGCCCGCTTGGCCCGTTCTTCACCAGCAGCACGCGCGTTGAAAT	array('B', [65, 65, 65, 68, 68, 68, 68, 68, 70, 70, 70, 70, 70, 72, 72, 72, 72, 71, 69, 72, 72, 72, 72, 72, 72, 72, 72, 71, 71, 70, 72, 72, 72, 72, 72, 72, 72, 70, 71, 72, 72, 72, 72, 71, 72, 72, 72, 72, 72, 72, 72, 72, 72, 72, 72, 72, 72, 72, 71, 71, 72, 72, 72, 72, 72, 72, 72, 72, 70, 70, 70, 68, 64, 66, 67, 62, 65, 66, 66, 66, 66, 66, 66, 66, 67, 66, 66, 65, 64, 64, 65, 64, 66, 66, 66, 66, 66, 66, 66, 66, 64, 62, 64, 66, 68, 68, 68])	[('AS', -24), ('XN', 0), ('XM', 4), ('XO', 0), ('XG', 0), ('NM', 4), ('MD', '100G2C0G0T1'), ('YS', 0), ('YT', 'CP')]

bamfile.close()

BAM または SAM ファイルに保存されているその他の情報を取得するメソッドは以下のようなものがある。

import pysam

bamfile = pysam.AlignmentFile("SRR616268.bam", "rb")

bamfile.filename
## b'SRR616268.bam'

bamfile.header
## {'SQ': [{'LN': 899, 'SN': 'ctg09_00003'}, {'LN': 2417, 'SN': 'ctg03_00032'}, {'LN': 4695, 'SN': 'ctg10_00023'}, {'LN': 4886, 'SN': 'ctg20_00036'}, {'LN': 11173, 'SN': 'ctg16_00033'}, {'LN': 14785, 'SN': 'ctg19_00007'}, {'LN': 15226, 'SN': 'ctg08_00002'}, {'LN': 16183, 'SN': 'ctg25_00029'}, {'LN': 20788, 'SN': 'ctg14_00035'}, {'LN': 21848, 'SN': 'ctg26_00057'}, {'LN': 25877, 'SN': 'ctg21_00025'}, {'LN': 34382, 'SN': 'ctg23_00054'}, {'LN': 46078, 'SN': 'ctg07_00012'}, {'LN': 51669, 'SN': 'ctg05_00050'}, {'LN': 64425, 'SN': 'ctg22_00009'}, {'LN': 64885, 'SN': 'ctg18_00021'}, {'LN': 91565, 'SN': 'ctg13_00018'}, {'LN': 109464, 'SN': 'ctg06_00026'}, {'LN': 122299, 'SN': 'ctg17_00030'}, {'LN': 125041, 'SN': 'ctg15_00005'}, {'LN': 131336, 'SN': 'ctg02_00047'}, {'LN': 134838, 'SN': 'ctg24_00011'}, {'LN': 141486, 'SN': 'ctg04_00046'}, {'LN': 222389, 'SN': 'ctg27_00027'}, {'LN': 231011, 'SN': 'ctg11_00044'}, {'LN': 237614, 'SN': 'ctg01_00022'}, {'LN': 465659, 'SN': 'ctg12_00024'}, {'LN': 472701, 'SN': 'ctg28_00010'}], 'PG': [{'PN': 'bowtie2', 'ID': 'bowtie2', 'VN': '2.2.4', 'CL': '"/usr/local/bin/../Cellar/bowtie2/2.2.4/bin/bowtie2-align-s --wrapper basic-0 -x L12Av1 -S SRR616268sub.sam -1 SRR616268sub.qc2_1.fastq -2 SRR616268sub.qc2_2.fastq"'}], 'HD': {'SO': 'coordinate', 'VN': '1.0'}}

bamfile.nreferences
## 28

bamfile.references
## ('ctg09_00003', 'ctg03_00032', 'ctg10_00023', 'ctg20_00036', 'ctg16_00033', 'ctg19_00007', 'ctg08_00002', 'ctg25_00029', 'ctg14_00035', 'ctg26_00057', 'ctg21_00025', 'ctg23_00054', 'ctg07_00012', 'ctg05_00050', 'ctg22_00009', 'ctg18_00021', 'ctg13_00018', 'ctg06_00026', 'ctg17_00030', 'ctg15_00005', 'ctg02_00047', 'ctg24_00011', 'ctg04_00046', 'ctg27_00027', 'ctg11_00044', 'ctg01_00022', 'ctg12_00024', 'ctg28_00010')

bamfile.lengths
## (899, 2417, 4695, 4886, 11173, 14785, 15226, 16183, 20788, 21848, 25877, 34382, 46078, 51669, 64425, 64885, 91565, 109464, 122299, 125041, 131336, 134838, 141486, 222389, 231011, 237614, 465659, 472701)

bamfile.mapped
## 1960816

bamfile.unmapped
## 14626

bamfile.nocoordinate
## 24516

bamfile.close()

リードのソート

ソートされた BAM ファイルからインデックスを作成する際に pysam.index を利用する。この機能は samtools index コマンドに等しい。

import pysam

pysam.sort("SRR616268.bam", "SRR616268.sorted")
## samtools sort SRR616268.bam SRR616268.sorted

インデックスの作成

BAM ファイルのソートは pysam.sort メソッドで行う。この機能は samtools sort コマンドに等しい。

import pysam

pysam.index("SRR616268.sorted.bam")
## samtools index SRR616268.sorted.bam

条件を満たすリードを出力

pysam.AlignmentFile を利用して開いたファイルからある範囲にあるリードを取得する際に fetch メソッドを利用する。ただし、paired-end リードの場合は予めソートする必要がある。

import pysam
bamfile = pysam.AlignmentFile("SRR616268.sorted.bam", "rb")

bamfile.references
## ('ctg09_00003', 'ctg03_00032', 'ctg10_00023', 'ctg20_00036', 'ctg16_00033', 'ctg19_00007', 'ctg08_00002', 'ctg25_00029', 'ctg14_00035', 'ctg26_00057', 'ctg21_00025', 'ctg23_00054', 'ctg07_00012', 'ctg05_00050', 'ctg22_00009', 'ctg18_00021', 'ctg13_00018', 'ctg06_00026', 'ctg17_00030', 'ctg15_00005', 'ctg02_00047', 'ctg24_00011', 'ctg04_00046', 'ctg27_00027', 'ctg11_00044', 'ctg01_00022', 'ctg12_00024', 'ctg28_00010')

bamfile.count("ctg09_00003", 1, 100)
## 1647

for read in bamfile.fetch("ctg09_00003", 1, 100):
    print(read)

## SRR616268.843433	163	0	99	6	93M	0	142	93	ACGAGTTGCGAGACCGCGAGGTTAAGCTAATCTCTTAAAGCCATTCTCAGTTCGGACTGTAGGCTGCAACTCGCCTACACGAAGTCGGAATC	array('B', [68, 70, 70, 70, 70, 70, 71, 71, 72, 72, 72, 71, 72, 72, 72, 72, 72, 72, 71, 72, 72, 70, 70, 71, 72, 71, 72, 72, 72, 72, 71, 71, 72, 72, 72, 72, 72, 71, 70, 70, 70, 70, 70, 66, 67, 68, 68, 68, 67, 68, 67, 67, 67, 66, 66, 66, 66, 66, 66, 65, 67, 66, 66, 66, 66, 66, 66, 66, 66, 66, 66, 66, 66, 66, 66, 64, 66, 66, 66, 66, 66, 62, 66, 66, 66, 66, 66, 66, 66, 64, 66, 66, 47])	[('AS', -1), ('XS', -6), ('XN', 1), ('XM', 1), ('XO', 0), ('XG', 0), ('NM', 1), ('MD', '22N70'), ('YS', -30), ('YT', 'CP')]
## SRR616268.864982	163	0	99	6	93M	0	147	93	ACGAGTTGCGAGACCGCGAGGTTAAGCTAATCTCTTAAAGCCATTCTCAGTTCGGACTGTAGGCTGCAACTCGCCTACACGAAGTCGGAATC	array('B', [68, 70, 70, 70, 70, 70, 72, 72, 72, 72, 72, 72, 71, 72, 72, 72, 72, 72, 71, 72, 72, 72, 72, 72, 72, 71, 72, 72, 72, 72, 72, 71, 71, 72, 72, 72, 72, 71, 70, 70, 70, 70, 70, 68, 68, 68, 68, 68, 68, 68, 67, 68, 67, 66, 66, 66, 66, 66, 66, 65, 66, 66, 66, 66, 66, 66, 66, 66, 65, 65, 66, 66, 66, 64, 66, 66, 66, 66, 65, 66, 66, 66, 66, 64, 64, 66, 66, 66, 64, 66, 66, 66, 63])	[('AS', -1), ('XS', -6), ('XN', 1), ('XM', 1), ('XO', 0), ('XG', 0), ('NM', 1), ('MD', '22N70'), ('YS', -26), ('YT', 'CP')]
## SRR616268.411741	83	0	99	2	4M2I101M	0	42	107	GTCATCCTTGCGAGACCGCGAGGTTAAGCTAATCTCTTAAAGCCATTCTCAGTTCGGACTGTAGGCTGCAACTCGCCTACACGAAGTCGGAATCGCTAGTAATCGCG	array('B', [68, 68, 68, 68, 65, 65, 65, 60, 66, 66, 66, 66, 66, 66, 64, 64, 66, 66, 64, 66, 66, 66, 65, 66, 66, 66, 66, 65, 66, 66, 66, 65, 66, 66, 66, 65, 66, 66, 66, 66, 66, 66, 67, 67, 67, 67, 67, 67, 68, 68, 68, 68, 70, 70, 69, 71, 72, 72, 72, 71, 72, 72, 72, 71, 71, 68, 72, 71, 71, 71, 70, 71, 70, 72, 72, 72, 72, 71, 70, 72, 72, 71, 72, 72, 72, 72, 72, 72, 72, 72, 72, 72, 71, 72, 70, 69, 70, 70, 70, 68, 68, 68, 68, 68, 65, 65, 65])	[('AS', -36), ('XS', -41), ('XN', 1), ('XM', 5), ('XO', 1), ('XG', 2), ('NM', 7), ('MD', '0A0C0G1G17N82'), ('YS', -9), ('YT', 'CP')]
## SRR616268.446068	83	0	99	6	5M1D102M	0	81	107	ACGTGTGCGAGACCGCGAGGTTAAGCTAATCTCTTAAAGCCATTCTCAGTTCGGACTGTAGGCTGCAACTCGCCTACACGAAGTCGGAATCGCTAGTAATCGCGGAC	array('B', [67, 68, 66, 68, 65, 65, 65, 49, 65, 66, 66, 66, 66, 66, 66, 66, 66, 66, 66, 66, 65, 65, 66, 66, 66, 66, 66, 66, 65, 65, 65, 63, 63, 66, 66, 66, 64, 65, 61, 66, 66, 66, 65, 65, 66, 66, 65, 61, 66, 63, 64, 63, 62, 57, 63, 68, 68, 70, 63, 69, 68, 70, 72, 71, 70, 67, 67, 63, 62, 57, 63, 55, 70, 71, 69, 71, 68, 66, 68, 61, 69, 69, 71, 69, 72, 71, 72, 72, 72, 72, 71, 70, 64, 71, 70, 70, 70, 70, 70, 68, 68, 66, 68, 68, 65, 65, 65])	[('AS', -21), ('XS', -26), ('XN', 1), ('XM', 3), ('XO', 1), ('XG', 1), ('NM', 4), ('MD', '3A1^T16N84T0'), ('YS', -1), ('YT', 'CP')]
## SRR616268.490935	83	0	99	6	4M3I100M	0	33	107	ACCTTCCGTTGCGAGACCGCGAGGTTAAGCTAATCTCTTAAAGCCATTCTCAGTTCGGACTGTAGGCTGCAACTCGCCTACACGAAGTCGGAATCGCTAGTAATCGC	array('B', [66, 68, 68, 68, 65, 65, 64, 64, 66, 66, 66, 65, 66, 66, 66, 66, 66, 66, 66, 66, 66, 66, 66, 66, 66, 66, 66, 66, 66, 66, 65, 65, 66, 66, 66, 66, 65, 66, 66, 66, 66, 66, 66, 66, 67, 67, 67, 68, 68, 68, 68, 68, 70, 70, 70, 72, 71, 71, 70, 72, 72, 72, 72, 72, 71, 71, 72, 72, 72, 72, 71, 71, 72, 72, 71, 71, 71, 71, 69, 70, 72, 71, 72, 71, 72, 72, 72, 72, 72, 72, 72, 72, 72, 71, 70, 70, 70, 70, 70, 68, 68, 68, 68, 68, 65, 65, 65])	[('AS', -27), ('XS', -32), ('XN', 1), ('XM', 3), ('XO', 1), ('XG', 3), ('NM', 6), ('MD', '2G0A18N81'), ('YS', -9), ('YT', 'CP')]
## ...

bamfile.close()

リードの情報の取得

ここで SAM / BAM ファイルに保存されている各リード（アラインメント）の情報を取得する方法を紹介する。

import pysam
bamfile = pysam.AlignmentFile("SRR616268.bam", "rb")

for read in bamfile:
    read.cigarstring
    ## '4M3I100M'
    read.cigartuples
    ## [(0, 4), (1, 3), (0, 100)]
    read.flag
    ## 89
    read.get_aligned_pairs()
    ## [(0, 19710), (1, 19711), (2, 19712), (3, 19713), (4, None), (5, None), (6, None), (7, 19714), (8, 19715), (9, 19716), (10, 19717), (11, 19718), (12, 19719), (13, 19720), (14, 19721), (15, 19722), (16, 19723), (17, 19724), (18, 19725), (19, 19726), (20, 19727), (21, 19728), (22, 19729), (23, 19730), (24, 19731), (25, 19732), (26, 19733), (27, 19734), (28, 19735), (29, 19736), (30, 19737), (31, 19738), (32, 19739), (33, 19740), (34, 19741), (35, 19742), (36, 19743), (37, 19744), (38, 19745), (39, 19746), (40, 19747), (41, 19748), (42, 19749), (43, 19750), (44, 19751), (45, 19752), (46, 19753), (47, 19754), (48, 19755), (49, 19756), (50, 19757), (51, 19758), (52, 19759), (53, 19760), (54, 19761), (55, 19762), (56, 19763), (57, 19764), (58, 19765), (59, 19766), (60, 19767), (61, 19768), (62, 19769), (63, 19770), (64, 19771), (65, 19772), (66, 19773), (67, 19774), (68, 19775), (69, 19776), (70, 19777), (71, 19778), (72, 19779), (73, 19780), (74, 19781), (75, 19782), (76, 19783), (77, 19784), (78, 19785), (79, 19786), (80, 19787), (81, 19788), (82, 19789), (83, 19790), (84, 19791), (85, 19792), (86, 19793), (87, 19794), (88, 19795), (89, 19796), (90, 19797), (91, 19798), (92, 19799), (93, 19800), (94, 19801), (95, 19802), (96, 19803), (97, 19804), (98, 19805), (99, 19806), (100, 19807), (101, 19808), (102, 19809), (103, 19810), (104, 19811), (105, 19812), (106, 19813)]
    read.get_reference_positions()
    ## [19710, 19711, 19712, 19713, 19714, 19715, 19716, 19717, 19718, 19719, 19720, 19721, 19722, 19723, 19724, 19725, 19726, 19727, 19728, 19729, 19730, 19731, 19732, 19733, 19734, 19735, 19736, 19737, 19738, 19739, 19740, 19741, 19742, 19743, 19744, 19745, 19746, 19747, 19748, 19749, 19750, 19751, 19752, 19753, 19754, 19755, 19756, 19757, 19758, 19759, 19760, 19761, 19762, 19763, 19764, 19765, 19766, 19767, 19768, 19769, 19770, 19771, 19772, 19773, 19774, 19775, 19776, 19777, 19778, 19779, 19780, 19781, 19782, 19783, 19784, 19785, 19786, 19787, 19788, 19789, 19790, 19791, 19792, 19793, 19794, 19795, 19796, 19797, 19798, 19799, 19800, 19801, 19802, 19803, 19804, 19805, 19806, 19807, 19808, 19809, 19810, 19811, 19812, 19813]
    read.get_tags()
    ## [('AS', -32), ('XN', 0), ('XM', 3), ('XO', 1), ('XG', 3), ('NM', 6), ('MD', '0C0C0C101'), ('YT', 'UP')]
    read.is_read1
    ## True
    read.is_read2
    ## False
    read.is_reverse
    ## True
    read.is_secondary
    ## False
    read.is_unmapped
    ## False
    read.mate_is_reverse
    ## False
    read.mate_is_unmapped
    ## True
    break

bamfile.close()

ファイルの保存

ある BAM / SAM ファイル中の特定のリードを取得して別のファイルに保存することも pysam を利用すれば簡単に行える。次の例では、SRR616268.sorted.bam ファイル中にあるリードのうち、ペアエンドリードの自分自身がマップされていて、かつ相手もマップされているリードだけを取得して別のファイルに保存する方法である。

import pysam

in_bamfile = pysam.AlignmentFile("SRR616268.sorted.bam", "rb")
out_bamfile = pysam.AlignmentFile("SRR616268.mapped.bam", "wb", template = in_bamfile)

for read in bamfile:
    if (not read.is_unmapped) and (not read.mate_is_unmapped):
        out_bamfile.write(read)


in_bamfile.close()
out_bamfile.close()

このように作成した SRR616268.mapped.bam ファイルのサイズは SRR616268.sorted.bam よりも大きくなる場合がある。一度 SRR616268.mapped.bam をソートすれば、サイズが小さくなる。