相変わらずいろいろ画策中。
やっぱり何が面倒ってab1ファイルを開いて2つ重なったピークを分離するところなわけで。
ピークコールの自動化ができるととても楽ちんになる。
BiopythonモジュールでもSeqデータを見られるらしい。
from Bio import SeqIO
from Bio.SeqIO import AbiIO
test = SeqIO.read('/Volumes/Documents/works/testseq.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)こんな感じにプロットを出せるのね。
ということはこれを4塩基で重ねていけば、ヘテロで重なったsecondary peakも取れるんだろうか。
とりあえずtestオブジェクトのアトリビュートを確認すると
>>> dir(test)
['__add__', '__bool__', '__class__', '__contains__', '__delattr__', '__dict__', '__dir__', '__doc__', '__eq__', '__format__', '__ge__', '__getattribute__', '__getitem__', '__gt__', '__hash__', '__init__', '__init_subclass__', '__iter__', '__le__', '__le___', '__len__', '__lt__', '__module__', '__ne__', '__new__', '__nonzero__', '__radd__', '__reduce__', '__reduce_ex__', '__repr__', '__setattr__', '__sizeof__', '__str__', '__subclasshook__', '__weakref__', '_per_letter_annotations', '_seq', '_set_per_letter_annotations', '_set_seq', 'annotations', 'dbxrefs', 'description', 'features', 'format', 'id', 'letter_annotations', 'lower', 'name', 'reverse_complement', 'seq', 'upper']
こんな感じ。
>>> 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=[])( ちなみにここで見ているシークエンスファイルはネットに落ちていたもので何かは知らない）
さらにannotationsは
>>> test.annotations.keys()
dict_keys(['sample_well', 'dye', 'polymer', 'machine_model', 'run_start', 'run_finish', 'abif_raw'])
>>> for k,v in test.annotations.items():
... print('key:'+k+' , values:'+v)
...
key:sample_well , values:A5
key:dye , values:Z-BigDyeV3
key:polymer , values:POP7
key:machine_model , values:3730
key:run_start , values:2018-05-25 18:23:59
key:run_finish , values:2018-05-25 20:38:32
Traceback (most recent call last):
File "<stdin>", line 2, in <module>
TypeError: must be str, not dictabif_rawキーに対するvalueには
'DATA1': 'Channel 1 raw data',
'DATA2': 'Channel 2 raw data',
'DATA3': 'Channel 3 raw data',
'DATA4': 'Channel 4 raw data',
'DATA5': 'Short Array holding measured volts/10 (EP voltage) during run',
'DATA6': 'Short Array holding measured milliAmps trace (EP current) during run',
'DATA7': 'Short Array holding measured milliWatts trace (Laser EP Power) during run',
'DATA8': 'Short Array holding measured oven Temperature (polymer temperature) trace during run',
'DATA9': 'Channel 9 processed data',
'DATA10': 'Channel 10 processed data',
'DATA11': 'Channel 11 processed data',
'DATA12': 'Channel 12 processed data',となっているようで、
>>> print(test.annotations['abif_raw']['DATA9'][100])
196
というふうに数値を取り出せる。
ピークを抽出するのはどうする？
9~12のデータの並び順は
>>> print(test.annotations['abif_raw']['FWO_1'])
GATC
からGATCの順のようだ。
と、いろいろ調べてきたが、結局2番めに高いピークをひろうにはすでにあるab1 fileではどうにもならないようだ。
ベースコールを行うphredで-dオプションを付けて再解析する必要があるらしい。
しかしPhredのダウンロードサイトが閉じているようで、入手ができそうにない。
Mac,Windowsのデスクトップアプリはあるようなのだけど、それだと自前のwebアプリにビルトインできないからなあ。
やっぱりコマンドラインで動くUnix版がほしいところ。
というわけで代わりにTraceTunerを試してみる。
sourceforge.net
これまたとっつきにくいアプリケーションの雰囲気を漂わせている。CVSだけど、気にせずsrc/に入ってmake
めちゃめちゃワーニングが出るな。
とりあえず、
$ cd ../rel/Linux
$ ./ttuner -h
-h (Help) This message
-Q (Quiet) Turn off status messages
-V (Verbose) Output more status messages
-nocall Disable base recalling and just use the original
called bases read from the input sample file
-recalln Disable adding bases to or deleting from the
original called sequence. Only recall Ns
-het Call call hetezygotes
-mix Call mixed bases
-min_ratio <ratio> Override the default threshold ratio of heights of
-trim_window <size> Set the trimming window size for averaging quality
to the specified value. The default is 10.
-trim_threshold <qv> Set the average quality value used in trimming to
-C <consensusfile> Specify the name of the FASTA file which contains
the consensus sequence
-edited_bases Start base recalling from the ABI's edited bases
-t <table> Use specified lookup table. This option overrides
the default (automatic choice of the lookup table)
as well as the options -3700pop5, -3700pop6, -3100,
and -mbace. To get a message showing
which table was used, specify -V option
-3730 Use the built-in ABI 3730-pop7 lookup table
-3700pop5 Use the built-in ABI 3700-pop5 lookup table
-3700pop6 Use the built-in ABI 3700-pop6 lookup table
-3100 Use the built-in ABI 3100-pop6 lookup table
-mbace Use the built-in MegaBACE lookup table
-c Output SCF file(s), in the current directory
-cd <dir> Output SCF file(s), in the specified directory
-cv3 Use version 3 for output SCF file. Default is
version 2.
-o <dir> Output multi-fasta files of bases (tt.seq),
their locations (tt.pos), quality values (tt.qual)
and status reports (tt.status) to directory <dir>
-p Output .phd.1 file(s), in the current directory
-pd <dir> Output .phd.1 file(s), in the specified directory
-q Output .qual file(s), in the current directory
-qa <file> Append .qual file(s) to <file>
-qd <dir> Output .qual file(s), in the specified directory
-s Output .seq file(s) in FASTA format, in the current
directory
-sa <file> Append .seq file(s) in FASTA format to <file>
-sd <dir> Output .seq file(s) in FASTA format, in the specified
directory
-qr <file> Output a quality report that gives data for a
histogram on the number of reads with quality
values >= 20, to the specified file
-if <file> Read the input sample filenames from the specified
file
-id <dir> Read the input sample files from specified directory
-tab Call heterozygotes or mixed bases and output .tab
file(s) in the current directory
-tabd <dir> Call mixed bases and output .tab file(s), in the
specified directory
-tal Output .tal file(s),in the current directory
-tald <dir> Output .tal file(s),in the specified directory
-hpr Output a homopolymer runs file in current directory
-hprd <dir> Output a homopolymer runs file(s),in the specified directory
-d Output .poly file(s),in the current directory
-dd <dir> Output .poly file(s),in the specified directory
-ipd <dir> Input the original bases and peak locations from a
.phd file in the specified directory.
さてオプションがいっぱいあるが
$ /Applications/tracetuner_3.0.6beta/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
t 26 414
上に上げたKB basecallerで読まれたシークエンス（チャート）と比べると
CNGNNNCGAGTCTTTGATGCAGTTGCGC
CGGTCTCGAGTCTTTGATGCAGTTGCGC
とまあほぼ一致している。
$ /Applications/tracetuner_3.0.6beta/rel/Linux/ttuner -mix -tab ./1I1_F_P1815443_047.ab1
こんな感じに投げてみると
# CHROMAT_FILE: 1I1_F_P1815443_047.ab1
# SOFTWARE_VERSION: TT_3.0.4beta
# NUM_ABC: 236
# NUM_SUBSTITUTIONS: 236
# NUM_DELETIONS: 0
# base2 qv2 pos2 ind2
Y 11 115 0
T 7 116 0
S 13 133 2
C 7 131 2
A 7 154 4
C 7 154 4
M 10 176 6
A 7 176 6
K 8 187 7
T 7 192 7
M 12 198 8
C 10 198 8
W 11 254 14
T 7 255 14
R 14 262 15
S 9 260 15
A 9 265 15
C 11 261 15
M 9 277 16
C 7 279 16
W 11 290 17
A 7 292 17
W 16 319 20
T 9 319 20
T 9 332 21
G 25 330 21
W 18 356 23
A 11 355 23
Y 11 374 25
T 9 373 25
K 24 385 26
T 9 384 26
W 23 412 28
A 9 411 28
・・・・・・・・・・
・・・・・・・・・・
こういうファイルが出来上がるわけだが、これでいいのかな？
CNGNNNCGAGTCTTTGATGCAGTTGCGC
CGGTCTCGAGTCTTTGATGCAGTTGCGC
YTSCACMAKTMCWTRSACMCWAWTTGWAY
なんかおかしいね。
$ /Applications/tracetuner_3.0.6beta/rel/Linux/ttuner -min_ratio 0.1 -d ./1I1_F_P1815443_047.ab1
こうすると、最も高いピークに加えて10%以上の高さのピークを次に選んでくれそう
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
これかな。
$ /Applications/tracetuner_3.0.6beta/rel/Linux/ttuner -mix -s ./1I1_F_P1815443_047.ab1
こっちでもいいかも
>1I1_F_P1815443_047.ab1 1195 16 579
CGGTMTCGAGTCTTAGATGCAKTTGCGCTCGASGCCATTASSTTGAGAGC
ACGTCTGTTTGGGCGTCATGCCTTGCGTCATTCTAGCCATCCATCTACTC
TCTGTGGGTGATGGGGGATGTGGAGATTGACCTTCCGTGCTTTAATTGTA
なお
R A or G
Y C or T
S G or C
W A or T
K G or T
M A or C
B C or G or T
D A or G or T
H A or C or T
V A or C or G
N any base
