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からはアミノ酸を読み出せないことかな。

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 をアミノ酸配列に翻訳する方法

スマニュー砲キター

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

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

肝心のアクセスは・・・ぜんぜん変わりません。あれえ？

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

文字列とQRコードを並べて合体した画像を作成する。

テプラには文字列からQRコードを生成して、テープに印字する機能があるのだが、このQRコード生成と文字情報の印字の入力は別系統な操作が必要で結構面倒臭い。
なので、QRコードと文字情報をまとめて画像にしてしまって画像挿入機能一操作だけで完結できるようにすると良いのではないかな。
ということでQR生成と合わせて文字も画像化してこれらを合体してやりたい。
QR生成と文字の画像化まではできているのであとは合体だな。

これにはpillowモジュールが使えるようだ。
要するに新規画像を作成して、その画像の適当な座標にパーツとなる画像をprintするという手順のようだ。
つなげ合わせる画像のサイズにぴったりの余白のない合成画像を作りたければ、それぞれの画像のサイズを調べてそれらを足し算した新規画像を用意すると良いというわけ。

    from PIL import ImageFont
    from PIL import Image
    from PIL import ImageDraw

    fontFile = 'ipaexg.ttf'
    font = ImageFont.truetype(fontFile, 64, encoding='utf-8')

    p_id2 = 'ID: ' + p_id
    w, h = font.getsize(p_id2)
    print('%s w:%d h:%d' % (p_id2, w, h))
    im = Image.new('RGB', (w, h), (255,255,255))
    draw = ImageDraw.Draw(im)
    draw.text((0, 0), p_id2, fill=(0,0,0), font=font)
    data = dir + '/id_' + p_id + '_id.png'
    im.save(data)

    plant_id = 'Plant: ' + plant_id
    w, h = font.getsize(plant_id)
    print('%s w:%d h:%d' % (plant_id, w, h))
    im = Image.new('RGB', (w, h), (255,255,255))
    draw = ImageDraw.Draw(im)
    draw.text((0, 0), plant_id, fill=(0,0,0), font=font)
    data = dir + '/id_' + p_id + '_pid.png'
    im.save(data)

    p_date ='Date: ' + p_date
    w, h = font.getsize(p_date)
    print('%s w:%d h:%d' % (p_date, w, h))
    im = Image.new('RGB', (w, h), (255,255,255))
    draw = ImageDraw.Draw(im)
    draw.text((0, 0), p_date, fill=(0,0,0), font=font)
    data = dir + '/id_' + p_id + '_pdt.png'
    im.save(data)

    owner = 'Owner: ' + owner
    w, h = font.getsize(owner)
    print('%s w:%d h:%d' % (owner, w, h))
    im = Image.new('RGB', (w, h), (255,255,255))
    draw = ImageDraw.Draw(im)
    draw.text((0, 0), owner, fill=(0,0,0), font=font)
    data = dir + '/id_' + p_id + '_own.png'
    im.save(data)

ここまでで４つの文字データを画像データに変換している。

    im0 = Image.open(dir + '/id_' + p_id + '_id.png')
    im1 = Image.open(dir + '/id_' + p_id + '_pid.png')
    im2 = Image.open(dir + '/id_' + p_id + '_pdt.png')
    im3 = Image.open(dir + '/id_' + p_id + '_own.png')
    w0, h0 = im0.size
    w1, h1 = im1.size
    w2, h2 = im2.size
    w3, h3 = im3.size

ここで書く画像のサイズを調べている

    w = max(w0, w1, w2, w3)

こうやって一番大きい画像に合わせて画像を作成すれば無駄が出ないはずなんだけど、なぜかエラーが出てうまくいかない・・・
maxに必要なモジュールってなんかいるのかな？
仕方がないので、とりあえず横幅は十分なサイズで一律固定として処理する。

    w = 800
    h = h0+h1+h2+h3
    im = Image.new('RGB', (w, h), (255,255,255))
    im.paste(im0, (0, 0))
    im.paste(im1, (0, h0))
    im.paste(im2, (0, h0 + h1))
    im.paste(im3, (0, h0 + h1 + h2))
    data = dir + '/id_' + p_id + '_merge.png'
    im.save(data)

    im1 = Image.open(dir + '/id_' + p_id + '.png')
    im2 = Image.open(dir + '/id_' + p_id + '_merge.png')
    w1, h1 = im1.size
    w2, h2 = im2.size

    im = Image.new('RGB', (w1 + w2, h1), (255,255,255))
    im.paste(im1, (0, 0))
    im.paste(im2, (w1, 30))
    data = dir + '/id_tag.png'
    im.save(data)

こんな感じに座標を指定して貼り付けてやればいいのだね。

ちなみに12ｍｍ幅のテープに印刷したらQRのドットが抜けてしまうのか読み取れず18ｍｍのテープが必要であった。
テプラのソフトで生成した QRコードだと読めたんだけどな。まあ４段に組んだ文字列は12mmじゃ虫眼鏡で見ないと読めないから18mmでいいんだけど。

QRコードをpythonで作る（その２）

QRコード生成をFlaskのwebアプリから行ってみる。
2つのパターンを想定。
前もってデータベースのIDごとのQRコードを生成。
データベースの1エントリーからQRコードを生成。

from functools import wraps
from flask import request, redirect, url_for, render_template, session
from flask_wtf import FlaskForm
from werkzeug import secure_filename
from wtforms import StringField, IntegerField, SubmitField, BooleanField
from ..models import Test
from testDB import Session
import os, shutil,, qrcode

class qrForm(FlaskForm):
    id = StringField('ID')
    qr1 = IntegerField('min')
    qr2 = IntegerField('max')
    submit = SubmitField('Submit')

@app.route('/qr/', methods=['GET', 'POST'])
def qr_index():
    form = qrForm()
    p_id = request.args.get('id', type = str)
    session = Session()
    q = session.query(Test)

# qr1~qr2の連番でQRコードをまとめて生成する。
    if form.qr1.data:
        dir = "./qr"
        if os.path.exists(dir):
            shutil.rmtree(dir)
        if not os.path.exists(dir):
            os.makedirs(dir)
        min = form.qr1.data
        max = form.qr2.data + 1
        for qr in range(min, max):
            qr = str(qr)
            img = qrcode.make('http://k-kuro.hatenadiary.jp/test_edit/?id=' + qr)
            data = dir + "/id_" + qr + '.png'
            img.save(data)
        flash('QR codes generated')


#データベースから１エントリーを選んでQRコードを生成。
    if p_id:
        # q = q.filter(Test.id==p_id).one()
        # plant_id = q[2]
        # date = q[4]
        # owner = q[8]   #この辺はまだ実装していない
        dir = "./qr"
        if os.path.exists(dir):
            shutil.rmtree(dir)
        if not os.path.exists(dir):
            os.makedirs(dir)
        img = qrcode.make('http://k-kuro.hatenadiary.jp/test_edit/?id=' + p_id)
        data = dir + "/id_" + p_id + '.png'
        img.save(data)
        flash('QR code generated')
    return render_template('qr.html', form=form, contents=q)

だいたいこんな感じで./qrディレクトリに画像が作られることを確認。
これをダウンロードする仕組みとテプラなり、タックシールなりに印刷する手段を考えないと。

＃今回のアドレスは架空。

for loopの中でintとstrを連結しようとしてハマったのはお約束。
やっぱりちゃんとローカルでテストしながらやらないと遠回りになるな

次は文字を画像化してQRコードと連結してテプラで印刷しやすいように加工だ。

こんな感じ

Pillow(PIL)で文字列を画像として保存(位置・サイズ自動調整) - cloverrose's blog
指定のフォントファイルを使って描画した文字列を画像ファイルに保存 - 強火で進め
Python, Pillowで画像を縦・横に連結（結合） | note.nkmk.me
この辺を参考に

fontとしてIPAフォントを使ってみる
IPAexフォント/IPAフォント | IPAフォントのダウンロードサイトです

import sys
from PIL import ImageFont
from PIL import Image
from PIL import ImageDraw

fontFile = 'ipaexg.ttf'
font = ImageFont.truetype(fontFile, 64, encoding='utf-8')

text = u'謹賀新年'
w, h = font.getsize(text)

print('%s w:%d h:%d' % (text, w, h))

im = Image.new('RGB', (w, h), (255,255,255))
draw = ImageDraw.Draw(im)
draw.text((0, 0), text, fill=(0,0,0), font=font)

im.show() 
im.save('image.png')

載っていたコードをちょいといじってpython3で実行してみた。

なぜ謹賀新年ｗ
あ、今日から令和。おめでとう。