ABIシークエンサーのデータを解析するには当然ABI謹製のアプリケーションで行えばいいのだが、シークエンサーに出向くのは面倒だし、最近はそもそも外注に出してデータだけメールで送られてくるというのが普通になってきたので、データ解析をローカルのPCで行う必要がある。普通にシークエンスチェックするだけならまあ、外注業者が返してくれるfastaファイルを見るだけで事足りるのだろうが、生のデータを見る必要があることもまあまあある。
とりあえずシークエンスデータ.ab1ファイルを開くpython moduleを探すと
https://pypi.org/project/abifpy/
このあたりが見つかる。
早速pipでインストールしてみる
$ pip3 install abifpy
$ python3
Python 3.6.1 (v3.6.1:69c0db5050, Mar 21 2017, 01:21:04)
[GCC 4.2.1 (Apple Inc. build 5666) (dot 3)] on darwin
Type "help", "copyright", "credits" or "license" for more information.
>>> from abifpy import Trace
>>> test = Trace('good_example_control.ab1')
>>> test.seq
'NNNNNNGGGGCATCCTGTGTTCTACCTGGCACCTGTCCCCATAGAAATGAGCGTGAGTGCCCGGGATCTGCTGCGGGGCTGTGCTGGGCTCTTTCTCAGCCTGGCCCGAAGTTTCCAGATCTGATTGAGCGAGAGAGCAGCAGGACCTGCCCCTCTGCTGGGCTCTTACCTTCGCGGCACTCGCCACTGCCCAGCAGCAGGTGAGGCCCAACACAACCAGTTGCAGGCGCCCCATGGTGAGCATCAGCCTCTGGGTGGCCCTCCCTCTGGGCCTCGGGTATTTATGGAGCTGGATCCAAGGTCACATGCTTGTTCATGAGCTCTCAGGCA'
こんな感じにシークエンスデータを読み出すことができる。
>>> test.qual_val
[5, 4, 5, 6, 6, 6, 9, 10, 9, 16, 16, 13, 19, 41, 42, 40, 47, 29, 62, 52, 46, 47, 32, 33, 42, 56, 47, 49, 35, 32, 30, 47, 62, 59, 62, 62, 59, 62, 62, 43, 25, 54, 59, 43, 59, 62, 39, 32, 54, 59, 49, 59, 59, 46, 54, 62,・・・・
>>> test.data['well']
'E10'
>>> test.data['model']
'3730'
>>> test.data['run start date']
datetime.date(2017, 9, 6)
シークエンスクオリティや諸々の情報を引き出すことも可能。
また、Biopythonモジュールを利用すると
from Bio import SeqIO
from Bio.SeqIO import AbiIO
test = SeqIO.read('good_example_control.ab1', 'abi')
print(test.annotations['abif_raw']['DATA1'])
(-7, 6, -7, -4, -8, 1, -4, -1, -3, 0, -2, -12, -3, -4, 6, 3, -4, -4, -1, 2, -3, -15, -6, -5, -8, -4, 1, -14, -8, -4, -3, -7, -6, -8, -3, 0, -3, -6, -12, -10, 0, 0, -15, -7, -14, -8, 0, -10, -2, -13, -3, -10, -8, -8, -14, -3, -9, -7, -10, -7, -8, -5, -5, -11, 0, -7, -1, -8, -6, -6, -4, -2, -13, -10, -5, -6, -9, -8, -4, -6, -19, -11, -11, -8, -6, 0, -2, -6, -8, -13, -5, -8, -4, -7, -6, -4, -10, -4,
・・・・・・
4, 488, 501, 529, 523, 530, 531, 550, 566, 577, 583, 579, 592, 577, 584, 593, 610, 597, 631, 633, 612, 612, 623, 645, 660, 649, 665, 688, 666, 640, 652, 658, 652, 660, 650, 635, 654, 667, 616, 609, 600, 588, 568, 550, 517, 524, 517, 510, 494, 470, 457, 439, 417, 408, 401, 396, 360, 334, 334, 317, 306, 285, 270, 241, 255, 230, 224, 203, 215, 188, 191, 192, 172, 151, 148, 146, 150, 137, 151, 144, 134, 136, 150, 145, 124, 132, 137, 145, 138, 143, 132, 126, 140, 134)
こんなふうにシークエンス波形プロットのデータも取れて、matplotlibを使うと波形データを描かせることも可能なようだ。
>>> test.id
'1I1_F'
>>> test.name
'1I1_F_P1815443_047'
>>> test.seq
Seq('CNGNNNCGAGTCTTTGATGCAGTTGCGCTCGAGGCCATTANGTTGAGAGCACGT...NNN', IUPACUnambiguousDNA())
>>> test.reverse_complement()
SeqRecord(seq=Seq('NNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNN...CNG', IUPACUnambiguousDNA()), id='<unknown id>', name='<unknown name>', description='<unknown description>', dbxrefs=[])
なお、ここに上げているシークエンスデータはネット上から拾ってきたデータなので、何者かは知らない。（上のabifpyでみたデータとは別）
さて、基本的なシークエンスを読むだけならこれでいいのだが、hetero genusなサンプルのシークエンスデータを読んで、どういう配列でヘテロになっているかを調べたい場合はどうすればよいだろうか。ピークコールは一番高いピークデータだけを拾うため、特に塩基欠損mutationとのヘテロ接合のシークエンスを読むと、ずれたシークエンスが混ざり合って全然データベースにマッチしない結果となる。故に、２番めに高いピークを抽出する方法がないかと模索する。
シークエンスデータからピークコールするソフトウェアといえばPhredが王道なのだけど、どうやら配布していたサイトが閉じてしまい、ダウンロードできないようになってしまっている。
そこで他にピークコールできるソフトウェアがないかと探したところ
https://sourceforge.net/projects/tracetuner/
こちらのソフトウェアでできそうなことがわかった。
安定版らしい3.0.1をダウンロードしsrcに入ってmakeを実施した。
上で解析した.ab1データを使って解析してみると
$ cd ./src
$ make
$ cd ../rel/Linux
$ ttuner -p ./1I1_F_P1815443_047.ab1
BEGIN_SEQUENCE 1I1_F_P1815443_047.ab1
BEGIN_COMMENT
CHROMAT_FILE: 1I1_F_P1815443_047.ab1
ABI_THUMBPRINT: 0
PHRED_VERSION: TT_3.0.4beta
CALL_METHOD: ttuner
QUALITY_LEVELS: 45
TIME: Sat Jun 29 21:04:35 2019
TRACE_ARRAY_MIN_INDEX: 0
TRACE_ARRAY_MAX_INDEX: 14037
TRIM: 16 622 0.010000
CHEM: unknown
DYE: big
END_COMMENT
BEGIN_DNA
c 14 114
g 16 118
g 13 135
t 17 138
c 7 154
t 7 168
c 13 176
g 11 183
a 13 199
g 18 203
t 17 219
c 17 223
t 20 235
t 14 245
t 14 255
g 14 259
a 15 276
t 15 289
g 16 295
c 16 306
a 17 319
g 24 330
t 25 346
t 26 357
g 26 364
c 25 376
g 26 387
c 25 399
........
このように解析が実施される。
上で読まれたシークエンスと比べると
CNGNNNCGAGTCTTTGATGCAGTTGCGC
CGGTCTCGAGTCTTTGATGCAGTTGCGC
こちらではNがない。シークエンスクオリティーのカットオフの基準がデフォルトでは緩いのかもしれない。
さて、このソフトウェアではピークコールだけでなく2nd ピークを抽出することも可能となっている模様。これが欲しかった。
以下のようにすると、最も高いピークに加えて10%以上の高さのピークを次に選んでコールしてくれるようだ。
$ ttuner -min_ratio 0.1 -d ./1I1_F_P1815443_047.ab1
1I1_F_P1815443_047.ab1 1.0 1.0 1.0 1.0 1.0
C 114 2990.000000 1.207398 N -1 -1.000000 -1.000000 51.000000 196.000000 155.000000 79.000000
G 118 3210.000000 1.296236 T 115 912.500000 0.368478 62.000000 140.000000 170.000000 72.000000
G 135 2444.500000 0.987118 C 132 935.500000 0.377766 59.000000 61.000000 208.000000 133.000000
T 138 2012.500000 0.812672 N -1 -1.000000 -1.000000 37.000000 35.000000 173.000000 154.000000
C 154 3376.500000 1.363471 A 154 3609.500000 1.457559 215.000000 265.000000 0.000000 0.000000
T 168 1846.500000 0.745639 N -1 -1.000000 -1.000000 25.000000 76.000000 0.000000 133.000000
C 176 2879.000000 1.162575 A 176 268.500000 0.108424 29.000000 224.000000 17.000000 102.000000
G 183 1223.000000 0.493862 N -1 -1.000000 -1.000000 5.000000 85.000000 122.000000 5.000000
A 199 1633.000000 0.659425 C 198 852.500000 0.344250 136.000000 89.000000 199.000000 0.000000
G 203 3149.000000 1.271604 N -1 -1.000000 -1.000000 110.000000 44.000000 286.000000 0.000000
T 219 4653.500000 1.879139 N -1 -1.000000 -1.000000 0.000000 178.000000 4.000000 425.000000
C 223 3368.000000 1.274430 N -1 -1.000000 -1.000000 0.000000 329.000000 0.000000 299.000000
T 235 2353.000000 0.885069 N -1 -1.000000 -1.000000 0.000000 0.000000 0.000000 351.000000
T 245 7217.500000 2.724202 N -1 -1.000000 -1.000000 22.000000 0.000000 0.000000 574.000000
T 255 5279.000000 1.665352 A 252 912.000000 0.287706 94.000000 2.000000 132.000000 536.000000
G 259 2757.000000 0.820499 A 262 726.500000 0.216211 77.000000 36.000000 273.000000 317.000000
A 276 6011.000000 1.741713 C 280 659.000000 0.190948 411.000000 39.000000 0.000000 0.000000
T 289 3652.000000 0.970141 N -1 -1.000000 -1.000000 101.000000 3.000000 233.000000 404.000000
G 295 8538.000000 2.130612 A 293 841.500000 0.209992 84.000000 0.000000 831.000000 144.000000
C 306 3019.000000 0.642641 N -1 -1.000000 -1.000000 14.000000 315.000000 11.000000 8.000000
A 319 4293.000000 0.916368 T 318 794.500000 0.169591 312.000000 0.000000 9.000000 73.000000
G 330 4400.000000 0.946491 T 332 857.000000 0.184351 23.000000 0.000000 391.000000 101.000000
T 346 5086.500000 1.070403 N -1 -1.000000 -1.000000 14.000000 14.000000 1.000000 550.000000
T 357 4582.000000 0.911786 C 352 601.000000 0.119595 46.000000 17.000000 83.000000 510.000000
G 364 6482.000000 1.361264 N -1 -1.000000 -1.000000 2.000000 0.000000 679.000000 62.000000
C 376 5179.000000 1.060825 T 373 475.000000 0.097295 0.000000 455.000000 0.000000 58.000000
G 387 4755.000000 0.927941 T 383 813.000000 0.158657 0.000000 90.000000 488.000000 46.000000
C 399 6513.500000 1.303052 N -1 -1.000000 -1.000000 0.000000 592.000000 0.000000 0.000000
T 414 4645.000000 0.878936 A 412 989.000000 0.187140 94.000000 73.000000 0.000000 544.000000
C 421 4935.500000 1.008171 A 421 501.000000 0.102339 57.000000 513.000000 0.000000 66.000000
........
2nd peakを並べて書くと
CGGTCTCGAGTCTTTGATGCAGTTGCGC
NTCNANANCNNNNNAACNANTTNCNTTN
となる。
もしこれがhetero genusなシークエンス配列を読んだものだとするとNはおそらく1st peakと同一な配列である可能性が高く、
CGGTCTCGAGTCTTTGATGCAGTTGCGC
CTCTATAGCGTCTTAACTACTTTCGTTC
ということになるのだが、ではこの２本のシークエンスがhetero接合しているのかというとさにあらず、上と下の塩基はランダムに入れ替わっている可能性がある。それを正しく読みほぐすにはまた別の操作が必要となる。その件はまた今度書く。