kuroの覚え書き

96の個人的覚え書き

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

なかなか先に進まないが、とりあえずちょとでもいじってみるか。
f:id:k-kuro:20190615223410p:plain
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

こんな感じでターミナル上にデータが送られてくることを確認。
f:id:k-kuro:20190522202019p:plain
これを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

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

f:id:k-kuro:20190522200859p:plain

#!/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のテストで1EPOCHあたりにかかった時間は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)を機械的に作成する(その2)

ちょっと調べたらいけそうな気がしてきた。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)
    if match:
        print(">" + record.id)
        print(match)

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

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)
    if match:
        seq = Seq(match, IUPAC.ambiguous_dna)
        print(">" + record.id)
        print(seq.translate())

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

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)

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

multi FASTA (DNA)からmulti FASTA (Amino Acid)を機械的に作成する

やりたいことは
複数の遺伝子のcDNA情報をまとめて記載したFASTA形式のファイルがあったとして、それをアミノ酸に翻訳し、clustalw等でアライメントを作成する。
cDNA情報はUTRを含んでいたりいなかったりまちまちである。
フレームを3フレームともチェックし、最も長いORFが作れるものを選択。
MetからStopまでを取り出してFASTA形式で保存する。
というもの。
これまで、そういうことをやってくれるウェブサービス
http://shigen.nig.ac.jp/tools/translatorV2/
を利用していたのだが、あろうことかサービスが閉じられてしまった。
研究ツールは公開するならやっぱりソースも公開しておいてほしいなあ。そうすればローカルで使い続けられるのに。

ということで自前でそういう処理をできないか模索中
DNA の翻訳 | Python を利用して DNA をアミノ酸配列に翻訳する方法

スマニュー砲キター

f:id:k-kuro:20190512154452p:plain

これが噂のスマートニュース?
全然アクセスが上がっている感じはないですがね。

で、スマホにアプリを入れてみたけど、どこからリンクされているのか見つけられんかった。

f:id:k-kuro:20190512172810p:plain
肝心のアクセスは・・・ぜんぜん変わりません。あれえ?

ま、アフィリエイトとかやってるわけじゃないしどうでもいいっちゃあどうでもいいんだけど。