kuroの覚え書き

96の個人的覚え書き

pythonでABIのシークエンスデータをゴニョゴニョする

シークエンスファイルとかfastaファイルとかMacのApEとかで開いてどうにかするのがだんだん億劫になってきた。
pythonでどうにかあんなことやこんなことができないかと調査中

abifpy · PyPI

まあこんなモジュールでも使えばどうにかなりそうな感じ。
引き続き使い方を調査だ。

ところで企業再編の嵐が吹き荒れた結果、Applied BiosystemsもThermoFisherグループに組み込まれちゃったんだな。
キャピラリーシークエンサーが登場したときはすげーもんが出てきたもんだと・・・年がバレる。
そういえば昔は"ABI"って呼んでたよな。ABI Prism 310とかって。Iってなんだったんだ?Industory?
Incか。今どきの若者はABIとか呼ばないのかな。だってIncとしては存在してないからな。


さて、モジュールを早速使ってみよう。

$ 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'

なるほどabifpyからTraceを読み込んでおいてtestオブジェクトににab1ファイルのデータを並べて記述、という感じなのかな。

>>> test.qual
"&%&'''*+*11.4JKIP>_UOPABKYPRDA?P_\\__\\__L:W\\L\\_HAW\\R\\\\OW_\\\\\\\\_\\WWY\\M_\\\\\\\\WT\\___\\_________\\__WRRW_\\\\ORRRRK___\\_YS____\\\\Y_\\________\\___\\_\\_\\__\\___RBW______W___________KWWRW_RR_\\___\\_\\__\\__\\\\\\\\\\______\\__\\_\\_WWKW___\\_Y_\\\\W_\\_____\\\\_YS__W_\\\\_\\OW_Y_\\__\\CW_YWW__WY____\\Y833)2-_______\\\\YDW_____W__\\____W___\\_____\\_\\__WW____RLY>G<(2R\\_\\3A=3"

これの意味するところはなんだろか。
asciiコードなのか。

           51  3      71  G      91  [     111  o
           52  4      72  H      92  \     112  p
33  !      53  5      73  I      93  ]     113  q
34  "      54  6      74  J      94  ^     114  r
35  #      55  7      75  K      95  _     115  s
36  $      56  8      76  L      96  `     116  t
37  %      57  9      77  M      97  a     117  u
38  &      58  :      78  N      98  b     118  v
39  '      59  ;      79  O      99  c     119  w
40  (      60  <      80  P     100  d     120  x
41  )      61  =      81  Q     101  e     121  y
42  *      62  >      82  R     102  f     122  z
43  +      63  ?      83  S     103  g     123  {
44  ,      64  @      84  T     104  h     124  |
45  -      65  A      85  U     105  i     125  }
46  .      66  B      86  V     106  j     126  ~
47  /      67  C      87  W     107  k
48  0      68  D      88  X     108  l 
49  1      69  E      89  Y     109  m 
50  2      70  F      90  Z     110  n

で、ASCII code=Q+33なので、

>>> 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)

この辺はシークエンサーの情報や、作業日時など