kuroの覚え書き

96の個人的覚え書き

SRAToolkit

公開されているSRAファイルをダウンロードして再解析がしたいとき、データのダウンロードが結構面倒なので(サイトの構造が複雑すぎてなかなかファイル本体にたどり着けない)専用のツールを利用したい。
とおもってインストールをbrewからやってみた。

$ brew install sratoolkit

ところがだ。

$ fastq-dump --gzip --sprit-files SRRxxxxxxxx

localeがおかしいとかのwarningはまあいいとして、どうもperlのライブラリがぶつかって動かないらしい。


というわけで、バイナリを手動でインストールしてみた。

$ wget https://ftp-trace.ncbi.nlm.nih.gov/sra/sdk/current/sratoolkit.current-centos_linux64.tar.gz
$ tar xvzf sratoolkit.current-centos_linux64.tar.gz
$ mv sratoolkit.2.10.7-centos_linux64 ~/src/
$ ln -s ~/src/sratoolkit.2.10.7-centos_linux64/bin/* ~/bin

homebrewでpython3.6(前のバージョン)をインストールする

$ cd /home/linuxbrew/.linuxbrew/Homebrew/Library/Taps/homebrew/homebrew-core/Formula
$ git log python.rb
commit b02eee23d9ac8e1027a212f269da6edd4c1f0f5f
Author: BrewTestBot <homebrew-test-bot@lists.sfconservancy.org>
Date:   Mon Apr 20 06:02:17 2020 +0000

    python: update 3.7.7_1 bottle.

commit 8bcd19011df9b5d42f5286bd5f3ec65e2aafc8c7
Author: Michka Popoff <michkapopoff@gmail.com>
Date:   Mon Apr 20 08:01:09 2020 +0200

    python: bump revision for libffi (#20152)

commit 4f455d5f8c2d45682448416d53d0c47355445177
Author: Michka Popoff <michkapopoff@gmail.com>
Date:   Sat Apr 18 18:17:50 2020 +0200

    python: fix style for Linuxbrew

・・・・・

こんな感じにgit logを開く
とにかく3.7を入れちゃうとモジュール関係も最新を入れてしまうのでGLIBC_2.18が必要となってしまうので、3.6に戻ってもらう。

commit ac64d11bca23f0d221967fef3b89e1d30053e067
Author: LinuxbrewTestBot <testbot@linuxbrew.sh>
Date:   Sun Jun 17 15:53:29 2018 +0000

    python: update 3.6.5_1 bottle for Linuxbrew.

commit 3ec730ba2ea8fcbd8fb853778b07236c3edb02c1
Merge: d07628b798 597da37dcc
Date:   Sun Jun 17 17:23:04 2018 +0200

    Merge branch homebrew/master into linuxbrew/master
    
     Conflicts:
            Formula/gdbm.rb
            Formula/pipenv.rb
            Formula/pre-commit.rb
            Formula/python.rb
            Formula/python@2.rb
            Formula/sip.rb
            Formula/vim@7.4.rb

commit f2a764ef944b1080be64bd88dca9a1d80130c558
Author: BrewTestBot <brew-test-bot@googlegroups.com>
Date:   Sun Jun 17 13:24:46 2018 +0000

    python: update 3.6.5_1 bottle.

commit b60c0325575ed4c3c90ba061880a79bc0c208cf0
Author: commitay <commitay@users.noreply.github.com>
Date:   Sun Jun 17 18:01:44 2018 +1000

    python: revision for gdbm
    
    Also, update `setuptools` and `wheel` resources

commit 5f5331172328b5d21592518fece5c332c42bc770
Merge: dbc1429eab c1f561a772
Author: Michka Popoff <michkapopoff@gmail.com>
Date:   Wed Jun 13 23:07:55 2018 +0200

    Merge branch homebrew/master into linuxbrew/master

commit c1f561a772a747d711657afa2644287aa0cc9a40
Author: commitay <commitay@users.noreply.github.com>
Date:   Wed Jun 13 19:26:30 2018 +1000

    python 3.7.0rc1 (devel) (#28975)

ここでやっと3.6登場

commit 0036460b7cf3f16e5e02880d893717ee86404ef3
Author: Chuancong Gao <chuanconggao@users.noreply.github.com>
Date:   Thu Mar 29 19:40:41 2018 -0700

    python 3.6.5

このへんで行くか。

$ git checkout 0036460b7cf3f16e5e02880d893717ee86404ef3 python.rb
$ HOMEBREW_NO_AUTO_UPDATE=1 brew install python
Warning: Calling 'devel' blocks in formulae is deprecated! Use 'head' blocks or @-versioned formulae instead.
Please report this issue to the homebrew/core tap (not Homebrew/brew or Homebrew/core), or even better, submit a PR to fix it:
  /home/linuxbrew/.linuxbrew/Homebrew/Library/Taps/homebrew/homebrew-core/Formula/python.rb:15

==> Downloading https://files.pythonhosted.org/packages/72/c2/c09362ab29338413ab
######################################################################## 100.0%
==> Downloading https://files.pythonhosted.org/packages/c4/44/e6b8056b6c8f2bfd14
######################################################################## 100.0%
==> Downloading https://files.pythonhosted.org/packages/fa/b4/f9886517624a4dcb81
######################################################################## 100.0%
==> Downloading https://www.python.org/ftp/python/3.6.5/Python-3.6.5.tar.xz
##############################                                   [  938.835230] serial8250: too much work for irq4
######################################################################## 100.0%
Warning: Calling 'devel' blocks in formulae is deprecated! Use 'head' blocks or @-versioned formulae instead.
Please report this issue to the homebrew/core tap (not Homebrew/brew or Homebrew/core), or even better, submit a PR to fix it:
  /home/linuxbrew/.linuxbrew/Homebrew/Library/Taps/homebrew/homebrew-core/Formula/python.rb:15

Error: An exception occurred within a child process:
  NoMethodError: undefined method `installed?' for OS::Mac::CLT:Module
Did you mean?  instance_of?

うーん、やっぱりだめか?

次はanacondaでトライ

Python3.7.7/CentOS7でエラー

Traceback (most recent call last):
  File "manage.py", line 2, in <module>
    from app import app, manager
  File "/home/kuro/database/app/__init__.py", line 23, in <module>
    from app.views import samples, genes, anno, fasta, blast, blast_nb, crispr, clustalw, seeds, plasmid, stat
  File "/home/kuro/database/app/views/crispr.py", line 8, in <module>
    import matplotlib
  File "/home/linuxbrew/.linuxbrew/Cellar/python/3.7.7_1/lib/python3.7/site-packages/matplotlib/__init__.py", line 205, in <module>
    _check_versions()
  File "/home/linuxbrew/.linuxbrew/Cellar/python/3.7.7_1/lib/python3.7/site-packages/matplotlib/__init__.py", line 190, in _check_versions
    from . import ft2font
ImportError: /lib64/libc.so.6: version `GLIBC_2.18' not found (required by /home/linuxbrew/.linuxbrew/lib/libstdc++.so.6)

matplotlibをインポートしようとするとGLIBC_2.18が無いと言って落ちる。
CentOS7ではGLIBC_2.18はサポートしないらしいので、そろそろCentOS7も8に上げるべき時期が来たのかもしれない。

KVM仮想環境でテストを開始するか。
それともbrewで入れずanacondaにしてみるか?

brewで新規にインストールすると3.7.7がインストールされてしまうので、brewの環境を一旦アンインストールし、yumで標準リポジトリから3.6系を入れてみることにした。そうすると問題はなくなることは確認できたが、yumで入れるといちいちsudoでパッケージを追加したりと若干使い勝手が悪い。

次はbrewで明示的に3.6系をインストールする方法を試してみる。
参照
macOS Mojave に Python 3.6 環境を構築する - Satoshi Oikawa - Medium

cDNA FASTAファイルから最長ORFを抽出し、5'UTR/CDS/3'UTRに分割してそれぞれのFASTAファイルを作成する(改訂版)

#DNA FASTAファイルから最長のORFを抽出し5UTR,CDS,3UTRに分割して保存する。

import sys, os, re
from Bio import SeqIO
from Bio.Alphabet import IUPAC
from Bio.Seq import Seq

fasta_file = sys.argv[1]
faname = os.path.basename(fasta_file)
fdir = os.path.dirname(fasta_file)
fname = os.path.splitext(faname)[0]
fext = os.path.splitext(fasta_file)[1]
utr_5_fasta = os.path.splitext(fasta_file)[0] + '_5utr' + fext
cds_fasta = os.path.splitext(fasta_file)[0] + '_cds' + fext
utr_3_fasta = os.path.splitext(fasta_file)[0] + '_3utr' + fext
for record in SeqIO.parse(fasta_file, 'fasta'):
    cds = re.findall('(ATG(?:\S{3})*?T(?:AG|AA|GA))', str(record.seq).upper())

    if cds:
        match = max(cds, key = len)
        seq = Seq(match, IUPAC.ambiguous_dna)
        utr_5 = re.sub(str(seq.strip()+'[ACGT]*'), '', str(record.seq).upper())
        utr_3 = re.sub(str(utr_5 + seq.strip()), '', str(record.seq).upper())
        with open(utr_5_fasta, 'a') as f:
            f.write(">" + record.id +'\n')
            f.write(str(utr_5) + '\n')
        with open(cds_fasta, 'a') as f:
            f.write(">" + record.id +'\n')
            f.write(str(seq) + '\n')
        with open(utr_3_fasta, 'a') as f:
            f.write(">" + record.id +'\n')
            f.write(str(utr_3) + '\n')
    else:
        t_cds = re.findall('(ATG[ACGT]*)', str(record.seq).upper())
        if t_cds:
            match = max(t_cds, key = len)
            seq = Seq(match, IUPAC.ambiguous_dna)
            utr_5 = re.sub(str(seq.strip()+'[ACGT]*'), '', str(record.seq).upper())
            utr_3 = re.sub(str(utr_5 + seq.strip()), '', str(record.seq).upper())
            with open(utr_5_fasta, 'a') as f:
                f.write(">" + record.id +'\n')
                f.write(str(utr_5) + '\n')
            with open(cds_fasta, 'a') as f:
                f.write(">" + record.id +'\n')
                f.write(str(seq) + '\n')
            with open(utr_3_fasta, 'a') as f:
                f.write(">" + record.id +'\n')
                f.write(str(utr_3) + '\n')

CDSがstopまで入っていないとき5UTRも切り出せていなかったところを改善してみた。

Pythonのライブラリアップデート

$ python3 -m pip install --upgrade <ライブラリ>

が基本書式。
インストールしているライブラリを一括でアップデートしたいときは

$ python3 -m pip install pip-review
$ pip-review --auto

とする。
アップデートする前に

$ python3 -m pip freeze > lib.txt

とやってリストを作っておくほうが良い。アップデートして不具合が生じるようなら

$ python3 -m pip install <パッケージ>==<バージョン>

という書式でバージョン指定してインストールし直す。

OSXで改行をまとめて消す。

普通にテキストエディタの置換でやろうとすると1万行とかあると時間がかかって仕方がない。
そこでsedで置換してやりたいのだがOSXsedGNU sedではなくBSD sedなのでお作法が違う。

$ sed -e :loop -e 'N; $!b loop' -e 's/\n/* /g' gene.txt >gene1.txt

こう書くと改行を消してアスタリスクとスペースを付け足してくれる。
1万行でも一瞬で完了。

seqkitでfastaファイルから一部を取り出す

multi fastaファイルから一部の遺伝子だけ取り出したサブfastaファイルを作るには

samtools faidx TAIR10_cDNA.fasta AT1G01010.1 AT1G01020.1 AT1G01030.1 AT1G01040.1 AT1G01050.1 > subset.fasta

のようにsamtoolsを使えばいいのだけれど、fastaファイルがcDNAで遺伝子のリストがtranscript IDではなくてgene IDだったとしたら、これでは抽出できない。
そんなときは

seqkit faidx TAIR10_cDNA.fasta -r AT1G01010* AT1G01020* AT1G01030* AT1G01040* AT1G01050* > subset.fasta

という感じでseqkitを使って正規表現で抽出すると良い。

更に遺伝子のリストがlist.txtなどテキストファイルに列記されているなら。

list=`cat list.txt`
seqkit faidx TAIR10_cDNA.fasta -r $list > subset.fasta

みたいにテキストを読み込んで渡してやる。

さらにさらに、遺伝子のリストが10000遺伝子とかめちゃめちゃ膨大だと引数の文字制限にぶち当たって

Argument list too long

とか文句言われるので、

list=`cat list.txt`
echo "$list" | xargs seqkit faidx TAIR10_cDNA.fa -r  >> subset.fasta

とxargsを使うと良いみたい。

>>ではなくて>でも内容は同じだが、何故か順番は変わってくる模様。

`はバックスラッシュなので注意。

Raspberry pi にarduinoをのっける

Raspberry piにはピンが出ていて、デジタル信号を操作した工作などができるようになっている。ただし、そのピンは基本デジタルであってアナログな信号は直接やり取りできない。アナログ(電圧値など)の入出力をコントロールするならワンチップマイコンの出番なわけだ。ArduinoならばUSB経由でPCと通信できるので、工作するまでもなく、単純にUSBケーブルでつなぐだけでよいのだけれど、それじゃあんまりおもしろくない。てことでワンチップArduinoRaspberry piのピンヘッダに直接乗っけてやろうと言うことになった。

ワンチップと言いながら、シリアル通信をすることもあり、ある程度精度はほしいのでクリスタルとかコンデンサは最低限つなぐためにユニバーサル基板の切れ端でちょっとしたIOボードを作ってやる。猛者は亀の子のような工作をしたりするがAtMegaを使っている時点である程度の大きさはあるため、実用上は意味はない。

f:id:k-kuro:20200602162355j:plainf:id:k-kuro:20200602162425j:plain

Piの1pinが電源なのでこれをAtmagaの7pinに繋ぎ、6pin→8pin (GND)、9pin→10pF→1pin と1Pin→10kΩ→1pinでリセット回路、後は8pin→1kΩ→2pin、10pin→1kΩ→3pinでシリアル通信。クリスタルは9,10ピンに繋いでやる。10pFコンデンサをそれぞれの足からGNDに落として完成。

Raspberry Piには設定>Add/Remove softwareでArduinoを検索して出てくるIDEを入れてやる。

プログラムの書き込みはシリアルポートを設定してやらないとならない。
シリアルポートはRaspbianの標準ではシリアルコンソールに使う設定となっているので、シリアル通信ポートに設定を変えてやっておく必要がある。

$ sudo raspi-confit

で5 Interfacing Option→P6 Serial→No→Yes→Ok→Finish
/etc/crazy.localの1. Exit 0の前に

/user/bin/gpio mode 0 alt5

と入れて、11pin をRTS信号に割り付けておく。

今回、AtMega168pに前もってArduinoファームウェアをインストール済みだったので、そのまま繋ぐとプログラムは動き出す。
Arduino IDEでシリアルポートを指定して、ボードをArduino Diecimila or Duemilanove w/ ATmega168pを選択する。
これでプログラムコンパイル→書き込みしてやればプログラムを即実行できる。


LM35DZを使った温度ロガープログラムの例

#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 ledPin =  13;      // the number of the LED pin
int sensorValue = 0;        // value read from the pot
int outputValue = 0;        // value output to the PWM (analog out)
int ledState = LOW;             // ledState used to set the LED
long previousMillis = 0;        // will store last time LED was updated
long interval = 10000;
void setup() {
  // initialize serial communications at 9600 bps:
  Serial.begin(9600);
   pinMode(ledPin, OUTPUT); 
}

void loop() {
  unsigned long currentMillis = millis();
  if(currentMillis - previousMillis > interval) {
    previousMillis = currentMillis;
    if (ledState == LOW)
      ledState = HIGH;
    else
      ledState = LOW;
     digitalWrite(ledPin, ledState);
  // read the analog in value:
      sensorValue = analogRead(analogInPin);
      // map it to the range of the analog out:
      outputValue = map(sensorValue, 0, 1023, 0, 330); //3.3V駆動としているので330
      // 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);
    }
}

シリアルコンソールを開くと10秒ごとにアナログポートの値とそれを温度変換した数値が表示されるはずだ。
温度計測と同時にLEDがついたり消えたりするので動作も確認できる。
なお、LM35DZを写真のような基板上にはんだ付けしてしまうとATmegaやRaspberry Piの熱の影響がもろに伝わってしまい、常に30度以上になってしまうので、センサーはツイストケーブルとかを使って本体から離したほうがよさそうだ。

matplotlibで自動的に振られる色が10色しか無い点を変更する

marplotlibは複数の要素をグラフ表示すると、自動的に色を分けてくれるわけだが、これが10色しか用意されていない。まあ実際10本以上折れ線が重なっていたらかなり見にくいことになるのは確かなのだが、それでもたくさんの色で表現したいこともあるわけで。そこでこの色の自動振り当てをどうにかする方法を考えた。

color_codes={3:'#E60012',7:'#F39800',11:'#FFF100',4:'#8FC31F',8:'#009944',12:'#009E96',1:'#00A0E9',5:'#0068B7',9:'#1D2088',2:'#920783',6:'#E4007F',10:'#E5004F',15:'#EB6100',19:'#FCC800',23:'#CFDB00',16:'#22AC38',20:'#009B6B',24:'#00A0C1',13:'#0086D1',17:'#00479D',21:'#601986',14:'#BE0081',18:'#E5006A',22:'#E60033',25:'#848484',0:'#000000'}

まずはこういうカラーチャートに番号を振っておく。番号が順番どおりでないのはグラフに表示させる順番で近い色が隣り合わせにならないように調整したため。

c_code = 0
#と初期化しておいて複数の要素を描かせるためにforループを回す。
for y in q1:
    y0 = [y.data1, y.data2, y.data3]
    y2 = [y.data10,y.data11,y.data12]
    y3 = [y.data13,y.data14,y.data15]
    y5 = [y.data19,y.data20,y.data21]
    y0m = statistics.mean(y0)
    y2m = statistics.mean(y2)
    y3m = statistics.mean(y3)
    y5m = statistics.mean(y5)
    y0d = statistics.stdev(y0)
    y2d = statistics.stdev(y2)
    y3d = statistics.stdev(y3)
    y5d = statistics.stdev(y5)
    c_code1 = c_code % 26
    gcolor = color_codes[c_code1]
#ここでループごとの色指定を行う。26色を用意しているので26で割ったあまりの数値で指定する。27番目で最初の色に戻るわけだ。
    x_ax = np.array([0, 1, 2, 3])
    y_ax = np.array([y0m,y2m,y3m,y5m], dtype=np.float16)
    y_err = np.array([y0d,y2d,y3d,y5d])
    ax.errorbar(x_ax, y_ax, y_err, capsize=3, linestyle=lineStyle, linewidth=lineWidth, color=gcolor, marker=mstyle, markersize=mSize, markeredgewidth=1, label = y.id)
    plt.xticks([0,1,2,3], ['control','1d','3d','7d'], rotation=90, fontsize=tSize, weight='bold')
    ax.legend(bbox_to_anchor=(1, 1), loc='upper left', borderaxespad=1, fontsize=lSize)
    plt.xlabel('samples', fontsize=xSize, weight='bold')
    plt.ylabel('FPKM', fontsize=ySize, weight='bold')
    plt.tick_params(axis='y', which='major', labelsize=ytSize)
    plt.tight_layout()
    plt.subplots_adjust(left=0.15, right=0.6)
    img_data1a = dir + "a_" + str(time.time()) + ".eps"
    plt.savefig(img_data1a)
    canvas = FigureCanvasAgg(fig)
    buf = io.BytesIO()
    canvas.print_png(buf)
    data = buf.getvalue()
    img_data = urllib.parse.quote(data)
#ここまででグラフを描写
    c_code = c_code + 1
#これで色指定をシフトさせる。

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

DataTablesの表にチェックボックスをつけて選択した内容から次のアクションを起こす

まず、DataTablesで表をweb上に表示できていることが前提。

いろいろなプラグインがあるのだけれど、ほぼ標準的に使われているButtonsに加え、Selectプラグインも入れておく。

<script type="text/javascript" src="{{ url_for('static', filename='js/datatables.min.js') }}"></script>
<script type="text/javascript" src="{{ url_for('static', filename='js/dataTables.select.min.js') }}"></script>
<script type="text/javascript" src="{{ url_for('static', filename='js/dataTables.buttons.min.js') }}"></script>

説明はコードを見ればわかるだろう。

<script type="text/javascript">
    $(document).ready(function() {
        var events = $('#events');
        var table = $('#result').DataTable({
            "lengthMenu": [[10, 25, 50, 100, -1], [10, 25, 50, 100, "All"]],
            scrollY: 600,
            scrollX: true,
            stateSave: true,
            stateDuration: 0,
            colReorder: true,
            columnDefs: [ {             #この辺がチェックボックスを表示させるコード
                    orderable: false,
                    className: 'select-checkbox',
                    targets:   0             
                } ],
                select: {
                    style:    'multi',           #複数行にチェックを入れる事ができる
                    selector: 'td:first-child',      #一列目にチェックボックスをつける
                    blurable: false,
                },
                order: [[ 1, 'asc' ]],
            buttons: [
                'selectAll',              #全行にまとめてチェック
                'selectNone',            #全行のチェックを外す
                {
                text: 'Go to selected data',      #チェックを入れた行から次のアクションを起こすボタンを作る
                action: function ( e, dt, type, indexes ) {
                    selData = table.rows({selected: true}, indexes).data().toArray();    #チェックを付けた行を集めてjavascriptの2次元配列を作る
                    const transpose = a => a[0].map((_, c) => a.map(r => r[c]));     #ある列の値だけを集めたいので行列入れ替えの関数を作る
                    selData2 = transpose(selData)[34]                 #この場合34列目(チェックボックス除く)の値を集める
                    window.open("/exp_diff/?id="+selData2, '_blanc')          #集めた文字列を引数にデータベースを再検索して表示させる
                    }
                },
                'copy',
                {
                text: 'Download TSV',
                extend: 'csvHtml5',
                fieldSeparator: '\t',
                extension: '.tsv'
                },
                {
                extend: 'colvis',                #カラムを選択するこの機能は廃止されたらしいのだが
                text: 'Column select'
                }
            ],
            dom: 'Biftlp'
          });
  } );
</script>

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