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接合しているのかというとさにあらず、上と下の塩基はランダムに入れ替わっている可能性がある。それを正しく読みほぐすにはまた別の操作が必要となる。その件はまた今度書く。