pythonでseq data

相変わらずいろいろ画策中。
やっぱり何が面倒って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 dict

abif_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

CRISPRの編集を調べるプログラム

以前、CRISPRによって編集された遺伝子配列を解析するプログラムを書いたが、100％マッチするときしか検出できないのは不便だな、ということで、アライメントを取るツールを利用できないだろうかと考え中。

pairwise2 | BioPython の pairwise2 ライブラリーを利用したペアワイズアラインメント
こちらの情報を参考に、biopythonのpairwise2を使ってみる。

Python 3.6.5 (default, Apr 25 2018, 14:26:36) 
[GCC 4.2.1 Compatible Apple LLVM 9.0.0 (clang-900.0.39.2)] on darwin
Type "help", "copyright", "credits" or "license" for more information.
>>> from Bio import pairwise2
>>> from Bio.pairwise2 import format_alignment
>>> for a in pairwise2.align.localms("ACCGTCC", "ACG", 2, -1, -0.5, -0.1):
...     print(format_alignment(*a))
... 
ACCGTCC
| ||
A-CG---
  Score=5.5

ACCGTCC
|| |
AC-G---
  Score=5.5

なるほど。これで比べる2つ目の配列をregexで混合した配列から使えるとよいのだけれど

例えば

AACTCTGTTGTACTCTGGTCATAATTCCCTCGGCTAATTTAACCCGCT

こんな配列があったとして

シークエンス結果が

AACTCTGTTGTACTCTaatTCtTAccTCCgTCGatTtATT
AACTCTGTTGTACTCTGGGaaAtcATtcgCctaGCtAaac

こんなだったとする。
正解は

AACTCTGTTGTACTCTAGGTCATAATTCCCTCGGCTAATT
                ^
AACTCTGTTGTACTCTG    ATAATTCCCTCGGCTAATTTAAC

こういうふうになっているのだけれど

import regex
>>> w = 'AACTCTGTTGTACTCTGGTCATAATTCCCTCGGCTAATTTAACCCGCT'
>>> s1 = 'AACTCTGTTGTACTCTAATTCTTACCTCCGTCGATTTATT'
>>> s2 = 'AACTCTGTTGTACTCTGGGAAATCATTCGCCTAGCTAAAC'
>>> s = "".join('['+i+j+']' for i, j in zip(s1, s2))
>>> m = regex.finditer(s,w,overlapped=True)
>>> for match in m:
...     print(match.start(), match.group(), match.end())
... 
>>> 

こんなふうに以前作ったアプリではうまくヒットしない。一方はミスマッチのinsertionがあり、もう一方はそれとずれてdeletionが入っているからだ
とはいえpairwise2になんとなく組み込んでみても

>>> for a in pairwise2.align.localms(w, s, 2, -1, -0.5, -0.1):
...     print(format_alignment(*a))
... 
-----AA--C---T---C---T---G---TT------G---T---A---C---T---C---T----G---G--T-------CA--TAA--TT---C--C---CT----------CG--GC--T--------A--A-------TT--TA--A--------CCCGCT
     ||  |   |   |   |   |   ||      |   |   |   |   |   |   |    |   |  |       ||  | |  ||   |  |   ||          ||  ||  |        |  |       ||  ||  |        |
[AA][AA][CC][TT][CC][TT][GG][TT][TT][GG][TT][AA][CC][TT][CC][TT][AG][AG][TG][TA][CA][T-A][TT][AC][CA][CT][TT][CC][CG][GC][TC][CT][GA][AG][TC][TT][TA][AA][TA][TC]----
  Score=59.9
................

これじゃだめだね。
もうちょっと考えよう。

ちなみに

>>> s1 = 'ATTCTTACCTCCGTCGATTTATT'
>>> s2 = 'GGAAATCATTCGCCTAGCTAAAC'
>>> s = "".join('['+i+j+']' for i, j in zip(s1, s2))
>>> m = regex.finditer(s,w,overlapped=True)
>>> for match in m:
...     print(match.start(), match.group(), match.end())
... 
16 GGTCATAATTCCCTCGGCTAATT 39
20 ATAATTCCCTCGGCTAATTTAAC 43
>>> 

こういうふうに投げる配列をうまく調整すると答えにだいたいたどり着けるのだけど、どれくらいの位置にどれくらいdeletionが入っているのかというのが運次第というのではちょっとね。

>>> s1 = 'AACTCTGTTGTACTCTAATTCTTACCTCCGTCGATTTATT'
>>> s2 = 'AACTCTGTTGTACTCTGGGAAATCATTCGCCTAGCTAAAC'
>>> s = "".join('['+i+j+']' for i, j in zip(s1, s2))
>>> n = 1
>>> while s:
...     m = regex.finditer(s,w,overlapped=True)
...     print(str(n)+': '+s)
...     for match in m:
...             print(match.start(), match.group(), match.end())
...     s = s[4:]
...     n += 1
... 
1: [AA][AA][CC][TT][CC][TT][GG][TT][TT][GG][TT][AA][CC][TT][CC][TT][AG][AG][TG][TA][CA][TA][TT][AC][CA][CT][TT][CC][CG][GC][TC][CT][GA][AG][TC][TT][TA][AA][TA][TC]
2: [AA][CC][TT][CC][TT][GG][TT][TT][GG][TT][AA][CC][TT][CC][TT][AG][AG][TG][TA][CA][TA][TT][AC][CA][CT][TT][CC][CG][GC][TC][CT][GA][AG][TC][TT][TA][AA][TA][TC]
3: [CC][TT][CC][TT][GG][TT][TT][GG][TT][AA][CC][TT][CC][TT][AG][AG][TG][TA][CA][TA][TT][AC][CA][CT][TT][CC][CG][GC][TC][CT][GA][AG][TC][TT][TA][AA][TA][TC]
4: [TT][CC][TT][GG][TT][TT][GG][TT][AA][CC][TT][CC][TT][AG][AG][TG][TA][CA][TA][TT][AC][CA][CT][TT][CC][CG][GC][TC][CT][GA][AG][TC][TT][TA][AA][TA][TC]
5: [CC][TT][GG][TT][TT][GG][TT][AA][CC][TT][CC][TT][AG][AG][TG][TA][CA][TA][TT][AC][CA][CT][TT][CC][CG][GC][TC][CT][GA][AG][TC][TT][TA][AA][TA][TC]
6: [TT][GG][TT][TT][GG][TT][AA][CC][TT][CC][TT][AG][AG][TG][TA][CA][TA][TT][AC][CA][CT][TT][CC][CG][GC][TC][CT][GA][AG][TC][TT][TA][AA][TA][TC]
7: [GG][TT][TT][GG][TT][AA][CC][TT][CC][TT][AG][AG][TG][TA][CA][TA][TT][AC][CA][CT][TT][CC][CG][GC][TC][CT][GA][AG][TC][TT][TA][AA][TA][TC]
8: [TT][TT][GG][TT][AA][CC][TT][CC][TT][AG][AG][TG][TA][CA][TA][TT][AC][CA][CT][TT][CC][CG][GC][TC][CT][GA][AG][TC][TT][TA][AA][TA][TC]
9: [TT][GG][TT][AA][CC][TT][CC][TT][AG][AG][TG][TA][CA][TA][TT][AC][CA][CT][TT][CC][CG][GC][TC][CT][GA][AG][TC][TT][TA][AA][TA][TC]
10: [GG][TT][AA][CC][TT][CC][TT][AG][AG][TG][TA][CA][TA][TT][AC][CA][CT][TT][CC][CG][GC][TC][CT][GA][AG][TC][TT][TA][AA][TA][TC]
11: [TT][AA][CC][TT][CC][TT][AG][AG][TG][TA][CA][TA][TT][AC][CA][CT][TT][CC][CG][GC][TC][CT][GA][AG][TC][TT][TA][AA][TA][TC]
12: [AA][CC][TT][CC][TT][AG][AG][TG][TA][CA][TA][TT][AC][CA][CT][TT][CC][CG][GC][TC][CT][GA][AG][TC][TT][TA][AA][TA][TC]
13: [CC][TT][CC][TT][AG][AG][TG][TA][CA][TA][TT][AC][CA][CT][TT][CC][CG][GC][TC][CT][GA][AG][TC][TT][TA][AA][TA][TC]
14: [TT][CC][TT][AG][AG][TG][TA][CA][TA][TT][AC][CA][CT][TT][CC][CG][GC][TC][CT][GA][AG][TC][TT][TA][AA][TA][TC]
15: [CC][TT][AG][AG][TG][TA][CA][TA][TT][AC][CA][CT][TT][CC][CG][GC][TC][CT][GA][AG][TC][TT][TA][AA][TA][TC]
16: [TT][AG][AG][TG][TA][CA][TA][TT][AC][CA][CT][TT][CC][CG][GC][TC][CT][GA][AG][TC][TT][TA][AA][TA][TC]
17: [AG][AG][TG][TA][CA][TA][TT][AC][CA][CT][TT][CC][CG][GC][TC][CT][GA][AG][TC][TT][TA][AA][TA][TC]
18: [AG][TG][TA][CA][TA][TT][AC][CA][CT][TT][CC][CG][GC][TC][CT][GA][AG][TC][TT][TA][AA][TA][TC]
16 GGTCATAATTCCCTCGGCTAATT 39
20 ATAATTCCCTCGGCTAATTTAAC 43
19: [TG][TA][CA][TA][TT][AC][CA][CT][TT][CC][CG][GC][TC][CT][GA][AG][TC][TT][TA][AA][TA][TC]
17 GTCATAATTCCCTCGGCTAATT 39
21 TAATTCCCTCGGCTAATTTAAC 43
20: [TA][CA][TA][TT][AC][CA][CT][TT][CC][CG][GC][TC][CT][GA][AG][TC][TT][TA][AA][TA][TC]
18 TCATAATTCCCTCGGCTAATT 39
22 AATTCCCTCGGCTAATTTAAC 43
21: [CA][TA][TT][AC][CA][CT][TT][CC][CG][GC][TC][CT][GA][AG][TC][TT][TA][AA][TA][TC]
19 CATAATTCCCTCGGCTAATT 39
23 ATTCCCTCGGCTAATTTAAC 43
22: [TA][TT][AC][CA][CT][TT][CC][CG][GC][TC][CT][GA][AG][TC][TT][TA][AA][TA][TC]
20 ATAATTCCCTCGGCTAATT 39
24 TTCCCTCGGCTAATTTAAC 43
23: [TT][AC][CA][CT][TT][CC][CG][GC][TC][CT][GA][AG][TC][TT][TA][AA][TA][TC]
21 TAATTCCCTCGGCTAATT 39
25 TCCCTCGGCTAATTTAAC 43
24: [AC][CA][CT][TT][CC][CG][GC][TC][CT][GA][AG][TC][TT][TA][AA][TA][TC]
22 AATTCCCTCGGCTAATT 39
26 CCCTCGGCTAATTTAAC 43
25: [CA][CT][TT][CC][CG][GC][TC][CT][GA][AG][TC][TT][TA][AA][TA][TC]
23 ATTCCCTCGGCTAATT 39
27 CCTCGGCTAATTTAAC 43
26: [CT][TT][CC][CG][GC][TC][CT][GA][AG][TC][TT][TA][AA][TA][TC]
24 TTCCCTCGGCTAATT 39
28 CTCGGCTAATTTAAC 43
27: [TT][CC][CG][GC][TC][CT][GA][AG][TC][TT][TA][AA][TA][TC]
25 TCCCTCGGCTAATT 39
29 TCGGCTAATTTAAC 43
28: [CC][CG][GC][TC][CT][GA][AG][TC][TT][TA][AA][TA][TC]
26 CCCTCGGCTAATT 39
30 CGGCTAATTTAAC 43
29: [CG][GC][TC][CT][GA][AG][TC][TT][TA][AA][TA][TC]
27 CCTCGGCTAATT 39
31 GGCTAATTTAAC 43
30: [GC][TC][CT][GA][AG][TC][TT][TA][AA][TA][TC]
28 CTCGGCTAATT 39
32 GCTAATTTAAC 43
31: [TC][CT][GA][AG][TC][TT][TA][AA][TA][TC]
29 TCGGCTAATT 39
33 CTAATTTAAC 43
32: [CT][GA][AG][TC][TT][TA][AA][TA][TC]
30 CGGCTAATT 39
34 TAATTTAAC 43
33: [GA][AG][TC][TT][TA][AA][TA][TC]
31 GGCTAATT 39
35 AATTTAAC 43
34: [AG][TC][TT][TA][AA][TA][TC]
32 GCTAATT 39
36 ATTTAAC 43
35: [TC][TT][TA][AA][TA][TC]
33 CTAATT 39
37 TTTAAC 43
36: [TT][TA][AA][TA][TC]
21 TAATT 26
34 TAATT 39
38 TTAAC 43
37: [TA][AA][TA][TC]
21 TAAT 25
22 AATT 26
34 TAAT 38
35 AATT 39
39 TAAC 43
38: [AA][TA][TC]
0 AAC 3
22 AAT 25
23 ATT 26
35 AAT 38
36 ATT 39
40 AAC 43
39: [TA][TC]
1 AC 3
3 TC 5
7 TT 9
11 AC 13
13 TC 15
18 TC 20
20 AT 22
23 AT 25
24 TT 26
25 TC 27
29 TC 31
36 AT 38
37 TT 39
38 TT 40
41 AC 43
40: [TC]
2 C 3
3 T 4
4 C 5
5 T 6
7 T 8
8 T 9
10 T 11
12 C 13
13 T 14
14 C 15
15 T 16
18 T 19
19 C 20
21 T 22
24 T 25
25 T 26
26 C 27
27 C 28
28 C 29
29 T 30
30 C 31
33 C 34
34 T 35
37 T 38
38 T 39
39 T 40
42 C 43
43 C 44
44 C 45
46 C 47
47 T 48

というわけでちょっと改良して、前から1文字分ずつ削っていって、合致する配列が現れるポイントを探る感じにしてみた。
下の方の数文字しかないところは全く意味がないので削る文字数の制限とかつけるともっと見やすい感じになるかも。

というわけでかなり作り込みました。

これはかなりいい感じじゃなかろうか

またまたサーバの構成をいじる

サーバ　　　CPU　　　　　　　　　　　　　 　メモリ　　　　　　　　　　　　理論性能
RX1330M3　E3-1230v6(4C/8T, 3.50GHz)　　　2400 UDIMM　64GB　　224.0GFLOPS
RX300S7　　E5-2667(6C/12T, 2.90GHz)x2　　1600 LV-RDIMM　32GB　278.4GFLOPS
RX300S7　　E5-2643(4C/8T, 3.30GHz)x2　　1600 LV-UDIMM　24GB　　　　211.2GFLOPS
RX200S7　　E5-2630(6C/12T, 2.30GHz)x2　　1333 RDIMM　40GB　　220.8GFLOPS
RX300S7　　E5-2620(6C/12T, 2.00GHz)x2　　1600 LV-RDIMM　24GB　　　192.0GFLOPS　　　　

RX200S6　　E5630(4C/8T, 2.50GHz)x2　　　1333 RDIMM　12GB　　　　　80.0GFLOPS         ＃GPU test機
RX200S6　　E5630(4C/8T, 2.50GHz)x2　　　1333 UDIMM　4GB　　　　　80.0GFLOPS　　　#バックアップ機

ちょいとCPUとメモリを追加して、構成をいじってみた。

ヒートシンクを調達すればE5-2630L→E5-2620 x2　192.0GFLOPSにアップできるので、これもそのうち。
ただし、いまあんまり重いデータ解析をやってないので、ほぼ趣味の範疇に入るかも。
それよりはGPUのバージョンアップのほうが優先かも。
（それよりは深層学習をちゃんと理解できるレベルまで勉強するほうが先決）
（ていうか論文書けよ）

追記
ヒートシンクを調達し、RX300S7をE5-2620 x 2CPUにした。クラスタはこれで十分かな。RX200S6は基本テスト用途かな。
バックアップ機はファンの音が結構うるさいので連続通電運転は辛い。
現状余っている機材はE5-2630L、E5603

ディープラーニングちょっとずつ

なかなか先に進まないが、とりあえずちょとでもいじってみるか。

jupyter notebookの使い方を確かめながらmnistのデータを使った練習をやってみる。

いろいろわからないまま言われるままに入力し、その出力をまずは眺めてみる。

zeissのlsmファイルから画像を取り出して重ね合わせる

zeissの共焦点レーザー顕微鏡で撮影したマルチチャンネルな画像ファイルの各チャンネルをバラバラにしたファイルを出力し、それらをstackではなく1枚の画像にmergeしたものを作成したい。使うのはImageJ。
とりあえずImageメニューの中のツールでできることは確認したので、それをマクロ登録してちゃっちゃと連続で処理したい。

dir = getDirectory("image");
name = replace(getInfo("image.filename"), ".lsm", "");
run("Stack to Images");
run("Merge Channels...", "c2=Ch2 c4=ChD keep");
for (i=0;i<nImages;i++) {
        selectImage(i+1);
        title = getTitle;
        print(title);
        saveAs("png", dir+"/"+name+"_"+title+".png");
}
while (nImages>0) {
          selectImage(nImages);
          close();
} 

これでよし。

VNCとpyenv

deep learning machineの構成をあれこれいじっているうちに動作がおかしくなってきたので、一旦リセットしてOSインストールからやり直すことにした。

ここまでの手順では

CentOS7をデフォルトの最小構成でインストール
ネットワーク設定
一般ユーザー追加
gnome desktopインストール
vncサーバインストール
NVIDIAドライバインストール
CUDAインストール
pyenvインストール
anacondaインストール
python3.6にダウングレード
tensorflow-gpuインストール
kerasインストール

という感じであったのだが、今回最小構成ではなく最初からgnome desktop、開発環境付きでインストールしたところ、これがマズったらしく、かなりツボにはまり未だ脱出できていない。

問題はvncサーバのインストールを忘れていて最後にこれを入れようとして起こった。

どうやらpyenvとvncサーバの組み合わせには難があるらしく、pyenvの環境変数を.bashrcに加えるとvncサーバを起動しても画面が表示できないらしい。

現状解決策としては.bashrcに入るはずの

export PYENV_ROOT="$HOME/.pyenv"
export PATH="$PYENV_ROOT/bin:$PATH"
eval "$(pyenv init -)"
export PATH="$PYENV_ROOT/versions/anaconda3-5.3.1/bin/:$PATH"

をコメントアウトし、systemctlでvncserverを起動してからログイン後にこれら環境変数をexportする、というものだが、メンドウクサス。

どうにかならないものか。

もう一回最初から手順通りにインストールしたほうが早そうだが、そのためにはサーバに直接でむかなくてはならず、これもまた面倒。

ここで気がついた。以前の環境じゃpyenvとvncが共存していたと思っていたが、単にインストールしたあと一回もrebootしてなかっただけなんじゃね？
どうインストールしてもpyenvとは両立しないのかも。

苦肉の解決策
まず、systemdによる自動起動は諦める。
.bashrcには上記環境変数を記載。
.bash_profileには

 User specific environment and startup programs

PATH=$PATH:$HOME/.local/bin:$HOME/bin

export PATH

vncserver :1 -geometry 1280x1024 -depth 24

# .bash_profile

# Get the aliases and functions
if [ -f ~/.bashrc ]; then
        . ~/.bashrc
fi

のように.bashrcを読み込む前にvncserver起動を行うように順番を入れ替える。

ただ、この方法だとvncでアクセスする前にsshでログインをしておく必要があり、sshの方でログアウトしちゃうとvncも終了してしまうのでちょっとなあ。

Raspberry piで温度ロガー

ラズパイを温度監視用ロガーとして使えないかなと。

ここまでラズパイはマイクロLinuxボックスとしてしか使っておらず、IoT的な電子工作はやっていない。
で、どうなんよ、と今更ながらにIO関係を調べてみたところ、AD変換とかはない。
抵抗とコンデンサを使って、あとはソフトウェアでゴリ押しでなんとかするとかできなくもないだろうけど、今どきそんなことをするやつはおらず、みんな安いAD変換モジュールとかを買ってきてつなぐだけなんだ。つまんね。

とまあ、話はそれたが、結局ラズパイはただのミニPCなので、やはりセンサ関係をコントロールするならArduinoとかマイコンがまだまだ活躍するわけだ。
ラズパイにはユーザインターフェースに専念してもらうことになる。

で、そうなると結局10年ほど前にごにょごにょやっていたArduino温度計に戻っちゃう。あの頃はラズパイなんてなかったから改造FONルータにopen-WRTを入れて作ってた。
k-kuro.hatenadiary.jp

結局同じことをすることになる。

温度センサは
温室の温度警報装置 - kuroの覚え書き
これでも使ったLM-35DZが残っていたのでそれを利用。

Arduinoのスケッチはサンプルをちょいとだけいじって

#include <TimeLib.h>
// These constants won't change.  They're used to give names
// to the pins used:
const int analogInPin = A0;  // Analog input pin that the potentiometer is attached to
const int analogOutPin = 9; // Analog output pin that the LED is attached to

int sensorValue = 0;        // value read from the pot
int outputValue = 0;        // value output to the PWM (analog out)

void setup() {
  // initialize serial communications at 9600 bps:
  Serial.begin(9600);
}

void loop() {
  // read the analog in value:
  sensorValue = analogRead(analogInPin);
  // map it to the range of the analog out:
  outputValue = map(sensorValue, 0, 1023, 0, 500);
  // change the analog out value:
  analogWrite(analogOutPin, outputValue);

  // print the results to the serial monitor:
  Serial.print(now());
  Serial.print("\t sensor = " );
  Serial.print(sensorValue);
  Serial.print("\t temp = ");
  Serial.println(outputValue);

  // wait 2 milliseconds before the next loop
  // for the analog-to-digital converter to settle
  // after the last reading:
  delay(60000);
}

こんな感じ。
結局Arduinoは時計を持っていないのでリアルタイムは残せない。
なのでシリアルを受信した時間から換算して、ラズパイ側で日時を割り出す方式でお茶を濁す。
RTCもいるな、やっぱり。そもそもラズパイにも時計入れたいし。

さて、ラズパイの側でシリアルを受信するプログラムであるがcuコマンドをインストールして利用する。

$ cu -l /dev/ttyUSB0 -s 9600

こんな感じでターミナル上にデータが送られてくることを確認。

これをpythonでファイルに落とし込むところを作る。

import datetime
import subprocess

date1 = datetime.datetime.now()
date2 = "{0:%Y%m%d_%H%M%S}".format(date1)
file1 = "./" + date2 + ".txt"
cmd = "cu -l /dev/ttyUSB0 -s 9600 > %s" % file1
subprocess.call(cmd, shell=True)

こんだけ。

このスクリプトをどうにかラズパイ起動と同時に走らせたいのだが、案外うまくいかない。

そういえばシリアルを受信するだけだったら

$ cat /dev/ttyUSB0 >> log.txt

みたいにするだけでできたような気がする。送信しないし、こっちのほうが簡潔かも。

#!/bin/sh
date1=`date +%Y%m%d_%H%M%S`
file1='./temp/'$date1'.txt'
cat /dev/ttyUSB0>>$file1

わざわざpython持ち出すよりこのほうがいいな。
これをcronに

crontab -e

@reboot /home/pi/temp.sh&

と追加してやったところ、問題なく記録が開始された。

primer設計でblast

適当に選んだ配列がoff targetを増幅しないか調べるためにblastnで検索するとき、普通にデフォルトでやっても何も引っかからない。
そんなときは

$ blastn -db ath -query ~/Desktop/act1_f.txt -word_size 7

のように-word_sizeオプションを付けると良い。

Deep learning マシンを変更

手持ちのサーバの構成を考え、Deep learningの環境テストを行うノードをPRIMERGY RX200S6に移した。

手順は

• クラスタから切り離す。
• 一般ユーザを作成。
• グラボを移し替え。
• NVIDIAのドライバインストール。
• CUDAのインストール。
• pyenvインストール。
• anacondaインストール。
• pythonのバージョンを3.6に落とす。
• tensorflow-gpuインストール。
• kerasインストール。

以上。
以前やったとおりで特に問題なく構築完了。
CPUが1世代前になって、理論速度では96GFLOPSから80GFLOPSに落ちたが、mnist_cnn.pyのテストで１EPOCHあたりにかかった時間は76秒と全く同じ。CPUの速度はほぼ関係ないということだろう。ということでRX300S7の方をクラスタに組み込むことにする。

EPOCH=1のトータルタイムは
GPUあり85.33秒
GPUなし138.65秒
sandy bridgeの96GFLOPSマシンで108秒だったのが80GFLOPSになったのでGPUなしは確実に遅くなったが、GPUありだとほとんど一緒だな。

multi FASTA (DNA)からmulti FASTA (Amino Acid)を機械的に作成する（その２）

ちょっと調べたらいけそうな気がしてきた。Biopythonを使うといろいろ簡単にできる模様。

まずはmultifastaを開いて配列を順番に読み込む

import sys
from Bio import SeqIO

fasta_file = sys.argv[1]

for record in SeqIO.parse(fasta_file, 'fasta'):
    ids = record.id
    desc = record.description
    seq = record.seq

    print('id:', ids)
    print('desc:', desc)
    print('seq:', seq)

とりあえずこんなスクリプトをfasta_in.pyという名前で作れば

$ python3 fasta_in.py <fasta_file.fasta>

という感じに投げるとidと説明書きとシークエンスを取り出すことができる。

bioinformatics - Python find longest ORF in DNA sequence - Stack Overflow
そんでもって最長ORFはここに出ている例を使って取り出せるので

import sys
from Bio import SeqIO
import re

fasta_file = sys.argv[1]

for record in SeqIO.parse(fasta_file, 'fasta'):
    for strand, seq in (1, record.seq), (-1, record.seq.reverse_complement()):
        for frame in range(3):
            index = frame
            while index < len(record) - 6:
                match = re.match('(ATG(?:\S{3})*?T(?:AG|AA|GA))', str(seq[index:]))
                if match:
                    orf = match.group()
                    index += len(orf)
                    if len(orf) > 1300:
                        pos = str(record.seq).find(orf) + 1
                        print(">{}, pos {}, length {}, strand {}, frame {}".format\
                           (record.id, pos, len(orf), strand, frame ))
                        print(orf)
                else: index += 3

これでいいかな？と思ったが、なんかおかしい。

そもそもcDNAのfastaを見ているのでstrandは1方向固定でいいから

for strand, seq in (1, record.seq), (-1, record.seq.reverse_complement()):

このforループはいらないな。
matchじゃなくてfindallのほうがいいかな？
もうちょっと調べる。

import sys
from Bio import SeqIO
import re

fasta_file = sys.argv[1]

for record in SeqIO.parse(fasta_file, 'fasta'):
    match = max(re.findall('(ATG(?:\S{3})*?T(?:AG|AA|GA))', str(record.seq)), key = len,  default=0)
    if match:
        print(">" + record.id)
        print(match)

これで最長ORFを取り出せた。
あとはこれをアミノ酸に置き換える。

#!/usr/bin/python3
import sys, re
from Bio import SeqIO
from Bio.Alphabet import IUPAC
from Bio.Seq import Seq

fasta_file = sys.argv[1]

for record in SeqIO.parse(fasta_file, 'fasta'):
    match = max(re.findall('(ATG(?:\S{3})*?T(?:AG|AA|GA))', str(record.seq)), key = len,  default=0)
    if match:
        seq = Seq(match, IUPAC.ambiguous_dna)
        print(">" + record.id)
        print(seq.translate())

できた。これで一気に処理できる。
問題があるとすると、stopコドンがなくて尻切れトンボのORFを持つcDNAからはアミノ酸を読み出せないことかな。