Hatena::ブログ(Diary)

kuroの覚え書き

[Vmax] [Macintosh/PC] [Network] [Make] [science] [Programing]

注:個人的覚え書きであり、いかなる内容も保証されるものではありません。
(ツッコミは歓迎しますが、全てに応えられるとも限りませんのであしからず。)

2018-05-21

[] CentOS7にNGS関連のソフトウェアをインストール

まずとにかくyumでwgetとnanoだけは入れておく。

そのうえでlinuxbrewを入れてしまえば、ほぼ問題なく殆どの環境が構築できてしまった。

http://linuxbrew.sh

Linuxbrewのインストール方法は上のサイトのままでOK

幾つかのソフトは

$ brew tap brewsci/science

$ brew tap brewsci/bio

でリポジトリを追加した上でinstallする必要があったが。

あとjavaSEは直接ファイルをダウンロードして、展開してやる必要があった。

HISAT2とFastQCはperlを要求するので、perlもbrewで入れた。その結果HISAT2は無事使えるようになったが、FastQCは

$ fastqc
-bash: /home/linuxbrew/.linuxbrew/bin/fastqc: /usr/bin/perl: 誤ったインタプリタです: そのようなファイルやディレクトリはありません

というエラーが出る。

/home/linuxbrew/.linuxbrew/bin/perlから/usr/binにリンクを張ってみたが

$ fastqc
perl: warning: Setting locale failed.
perl: warning: Please check that your locale settings:
	LANGUAGE = (unset),
	LC_ALL = (unset),
	LANG = "ja_JP.UTF-8"
    are supported and installed on your system.
perl: warning: Falling back to the standard locale ("C").
Exception in thread "main" java.awt.HeadlessException: 
No X11 DISPLAY variable was set, but this program performed an operation which requires it.
	at java.desktop/java.awt.GraphicsEnvironment.checkHeadless(GraphicsEnvironment.java:208)
	at java.desktop/java.awt.Window.<init>(Window.java:548)
	at java.desktop/java.awt.Frame.<init>(Frame.java:423)
	at java.desktop/java.awt.Frame.<init>(Frame.java:388)
	at java.desktop/javax.swing.JFrame.<init>(JFrame.java:180)
	at uk.ac.babraham.FastQC.FastQCApplication.<init>(FastQCApplication.java:63)
	at uk.ac.babraham.FastQC.FastQCApplication.main(FastQCApplication.java:332)

こんな表示が出るな。これじゃ駄目かな?

[][]テスト運転

さて、早速

見せてもらおうか、富士通の中古サーバーの性能とやらを

(注:ガ◯ダ◯って実は殆どみたことないんだが)

大体270MBくらいのfastq.gzファイルを4つ

trimmomatic SE \
        -threads $cpu -phred33 -trimlog log.trimlog \
        ./fastq_files/"$samp"_"$exp"_L001_R1_001.fastq.gz \
        ./trimmed/"$samp"_"$exp"_L001_R1_001_trim_np.fastq.gz \
        ILLUMINACLIP:"$adapter":2:40:15 \
        LEADING:3 \
        TRAILING:3 \
        SLIDINGWINDOW:4:15 \
        MINLEN:36

こんな感じでトリミングし、bwaでアラビcDNAにマッピング。出来上がったBAMファイルをマージまでやってやる。

DDBJスパコンで8CPU指定で2:03:57かかったところを自前鯖16スレッド使用で0:30:57で終了!

なんと赤い彗星の3倍速いを軽く上回り、4倍以上速いじゃないの。

内訳を見ると

DDBJ→富士通自前鯖

Trimmomatic 1:12:22→0:16:26

bwa 0:06:17→0:07:51

samtools merge 0:06:59→0:06:40

んーJAVAが4.5倍くらい速く、その他はいい勝負か、ちょっと遅い。

たしかbwaもsamtoolsもシングルコアで処理をする部分があって、そこで地のクロックスピードの差が出ちゃってるのかね。

結論:javaは速い

ただし1ノードしか無いから並列処理で負ける。多サンプル並列処理しないデータ解析ならそこそこいい感じだと思う。

2018-05-19

[] 自前鯖ゲット

というわけで出来心から中古ラックマウントサーバマシンを入手してしまった。

富士通のPRIMERGYシリーズで、5年保守期間が切れて市場に出たものの、大量に売れ残って投げ売りされている2010年ごろのモデル

RX200 S6がそれだ。

http://jp.fujitsu.com/platform/server/primergy/product-navi/html/rx200s6-201005-pgr2062e63.html

http://jp.fujitsu.com/platform/server/primergy/pdf/20100531/rx200_s6.pdf

標準状態からかなりカスタマイズされたそれは

E5503 (2GHz)→E5630 (2.53GHz)×2

2GB DDR3 1333 UDIMM →4GB

SAS HDD:146.8GB(10krpm)×1→SAS HDD:146.8GB(15krpm)×4

となっている。

新品で納入されたときには100万を超えていそうな構成だけど、今回入手価格は6800円w

メモリ、HDDの容量はちょっとあれだけど、4コア(8スレッド)x2ということでゲノム解析に使うには現状の4コア4スレッドのiMacでやるより「3倍は速い」と思う。

さて、電源を入れてみると予想外にWindows server 2008が立ち上がってきた。ただしAdministratorのパスワードがわからないのでどうしようもないし、そもそもWindowsを使う意味もないので、さくっとLinuxを入れることにする。

事前調査にて、Fujitsuが非サポートながら動作確認を取っているらしいリストが公開されているらしいことを知った。

http://ottoserver.com/pcserver1/wp/archives/2264

それと、SASのRAIDコントローラのドライバが要りそうなこともそこには書かれていたが。

http://jp.fujitsu.com/platform/server/primergy/software/linux/products/distribution/free-os.html

まずはリストを見てみると、流石に古い機種なので最新のバージョンまではテストされていず、CentOSだと6.2まで動作実績があるようだ。多分枯れた構成となっているから、7でも問題なくインストールできそうな気もするが、とりあえずそこまで最新を追っかけることも必要としていないので、6.2で行ってみることにする。なお、なぜCentOSなのかといえば、これまでにそれなりに触っていて敷居が低くなっているからと言うしかない。

インストールするためにisoイメージをダウンロードしてくる。

http://mirror.nsc.liu.se/centos-store/6.2/isos/x86_64/

一応DVD2もダウンロードしておいたが、最小構成でインストールする分には必要はなさそうだ。

isoイメージをDVDーRに焼き焼きしておもむろにDVDドライブにセット、電源オン。とくに問題なくインストーラが起動し、メディアのチェックをするかと聞かれるので、一応やってもらうことにする。メディアは正常ということでインストールに進む。富士通のサイトの動作情報によると、インストールに際し、特にドライバを当てたりする必要もなさそうだったので、そのままストレージをフォーマットしてインストールに進めばいいようだ。

インストールはデフォルト構成(minimal install)でやっておく。今時のLinuxは必要な機能をあとから追加するのもそう苦ではないので、HDD容量も厳しいゆえ、できるだけいらないものは入れない方向でやっていく。基本クライアントからsshで利用する前提なのでGUIも入れない。

さて、インストールも拍子抜けするぐらい何事もなく終わり、再起動となる。

BIOSが終わり、CentOSの画面らしきグラフィックが表示され、プログレスバーが進んでいくが・・・・何故かいきなりkernel panic

https://mynotebook.h2np.net/post/130

https://qiita.com/murabiton/items/d44fc01b9e1ac2fdcfab

https://rfs.jp/server/security/selinux01.html

どうもSELinuxというのが問題のような感じがする。

ということで起動時にGRUBの画面で適当なキーを押して止めて起動オプションを入れてやる。


中略(あとで書く)

で、起動を再開させると、今度はプログレスバーが最後まで行って、無事に起動できた。


ntpの設定

https://web-dev.hatenablog.com/entry/linux/centos/time-jp


https://qiita.com/maruware/items/eb659266a45021cf486c

[] CentOS7にすることにした。

CentOS6.2にて環境構築を試みるも、やっぱりちょっとあちこち古くて、構築に難儀しそうだったので、最新の7を入れることにした。そうしたところ自分の知っているコマンドと色々変わっているらしく、それはそれで色々苦労しそうな予感。

まずはネットワークの設定から。

https://qiita.com/taro0219/items/8d3be39c5cb882d1ba5d

しかしネットワーク設定さえしてしまえば、あとはほとんど不都合なく使えそな感じはある。

とにかく必要そうなプログラムを列記しておこう。

pigz

java

python3

MARIADB

linuxbrew

samtools

tophat

bowtie

cufflinks

FastQC

fastx-tools

HISAT2

BWA

SRA Toolkit

R

Trimmomatic

2018-05-04

[] [RaspberryPi] Raspberry PiにDAC

Raspberry Piのヘッドフォンジャックはなかなかにひどい音しか出ない。HDMIから出せばそこそこまともな音が出るのにどうしてなんだろうね。

ってことで、DACを積んで出力をしてみることにする。USBにつなぐDACというのがかんたんなんだろうけど、折角ラズパイなのでピンヘッダから直接デジタル信号を取り出してアナログオーディオに変換してやりたい。

I2Sインターフェイスでかんたんにできるようなので、アマゾンで基板を買ってみた。

配線とソフトウェアの設定は

http://www.neko.ne.jp/~freewing/raspberry_pi/raspberry_pi_i2s_dac_audio_ti_pcm5102/

こちらの通りにやってみた。

$ omxplayer -o alsa:hw:0,0 hoge.mp3

Audio codec mp3 channels 2 samplerate 44100 bitspersample 16

Subtitle count: 0, state: off, index: 1, delay: 0

こんな感じで音声を出力できるところまではできた。

さて問題はここから。

Raspbianの標準ブラウザChromiumの音をこのインターフェイスから出すにはどうしたらいいんだろう。

YouTubeの音をもうちょっとマシな音で聞きたいんだけどなあ。

[] [RaspberryPi] Raspberry PiにDAC その2

YouTubeの音声をもうちょっとキレイに聴きたいという動機で着手したRaspberry PiのDAC搭載だが、Chromiumからの音を出せないでいたわけだが、どうにか音をだすことに成功した。

結論から言うと元々の内蔵音源をオフにしてやる必要があったようだ。

/boot/config.txtの

dtparm-audio=on

をコメントアウト

dtoverlay=hifiberry-dac

を追記

~/.asoundrc

pcm.!default{

type hw

card 0

}

ctl.!default{

type hw

card 0

}

これでreboot

難点はUIの音量調節が効かなくなること。

そして最も困ったことがタッチパネルの電源をラズパイのUSBから取ると電流不足となってしまうらしいこと。

つまりタッチパネルを別電源で取らないと使えないということ。

ますますタブレットの気軽さとはかけ離れていくな。


追記

モニタに供給する電流の設定をデフォルトから下げてやるとラズパイからのUSB電源供給でもフリーズすることがないことがわかった。

config.txt

#KUMAN HDMI Display setting

max_usb_current=0.5

のように通常1だったところを0.5に下げてやるだけでOKだ。

2018-04-19

[] fastaから配列を取り出す

$ samtools faidx genome.fasta chr1:15,000-50,000

こんな感じでfastaファイルから取り出したい配列を指定して取り出すことが可能。

ある遺伝子を特定したとして、その上流プロモーターの配列をfastaから取り出したいなどの用途で使える。

2018-04-03

[][] SQLとBLASTの連携

まず、SQLで抽出された遺伝子のIDを元にtranscriptのFASTAファイルから配列を抽出。

抽出した配列をfasta形式で一時保存し、Arabidopsisの遺伝子データベースに対してtblastxでサーチ。

結果をwebに表示。

こんな流れを構築した。

https://pypi.python.org/pypi/biopython/1.70

https://qiita.com/MaedaTaro_Umiushi/items/06f9d847ba025c175547

2018-03-29

[] DDBJスパコン運用〜ファイル転送2

大量のデータをやり取りしたいときはASPERAで行うのがいいが、ちょっとしたファイルをアップロードダウンロードするならsshfsでディスクをマウントしてしまうのが手っ取り早い。

コマンドラインでマウントポイント ~/ddbjにマウントするには

$ sshfs username@gw2.ddbj.nig.ac.jp:/target/directory/here ~/ddbj

てな感じで良い。

MacのFinderから使う方法はDDBJのFAQに出ていた。はずだがどこに出ていたかわからなくなった。

あったあった。

https://sc.ddbj.nig.ac.jp/ja/system-guide/faq/02mghj

とにかく

http://macfusionapp.org

をクライアントとして使えと。

[]BLASTを自前で

バイオインフォマティックスをやってるならBLASTくらい自分で構築しないとな、

というわけでもないけど、やはりアノテーションデータベースを構築したいので、blastをMacに入れることにする。

さらっとぐぐってみたらHomebrewで行けそうなので

$ brew install blast
Updating Homebrew...
==> Auto-updated Homebrew!
-------
$ blastn -version
blastn: 2.7.1+
 Package: blast 2.7.1, build Feb 16 2018 14:27:28

素晴らしい。ちゃんとインストールできた。

データベースを作らないと。マッピングにはゲノムコンティグを使ったけど、ついでにダウンロードしておいたtranscriptのfastaファイルを使って、

$ mkdir blast
$ cd blast/
$ makeblastdb -in /Users/kuro/rnaseq/ref/transcript.fa -dbtype nucl -out mydb

Building a new DB, current time: 03/29/2018 17:43:14
New DB name:   /Users/kuro/local/bin/blast/mydb
New DB title:  /Users/kuro/rnaseq/ref/transcript.fa
Sequence type: Nucleotide
Keep MBits: T
Maximum file size: 1000000000B
Adding sequences from FASTA; added 30364 sequences in 1.38843 seconds.

$ ls
mydb.nhr  mydb.nin	mydb.nsq

こんな感じでデータベースが一瞬で出来上がる。

さあ、blast!

適当な遺伝子配列をfastaで保存しておいて、

$ blastn -query test.fasta -db mydb
BLASTN 2.7.1+


Reference: Zheng Zhang, Scott Schwartz, Lukas Wagner, and Webb
Miller (2000), "A greedy algorithm for aligning DNA sequences", J
Comput Biol 2000; 7(1-2):203-14.



Database: /Users/kuro/rnaseq/ref/transcript.fa
           30,364 sequences; 48,838,386 total letters



Query= xxxxx1Scf01234 xxxxx1Scf01234:68915..70156 (- strand) class=mRNA
length=1242

Length=1242
                                                                      Score     E
Sequences producing significant alignments:                          (Bits)  Value

  yyyyy.106780.1 pacid=16959505 locus=yyyyy.106780 annot-version=...  196     3e-49
  yyyyy.349830.1 pacid=16979855 locus=yyyyy.349830 annot-version=...  187     2e-46


> yyyyy106780.1 pacid=16959505 locus=yyyyy.106780 annot-version=v1.0
Length=1716

 Score = 196 bits (106),  Expect = 3e-49
 Identities = 224/281 (80%), Gaps = 7/281 (2%)
 Strand=Plus/Plus
=====中略=====
Lambda      K        H
    1.33    0.621     1.12 

Gapped
Lambda      K        H
    1.28    0.460    0.850 

Effective search space used: 58427489152


  Database: /Users/kuro/rnaseq/ref/transcript.fa
    Posted date:  Mar 29, 2018  5:43 PM
  Number of letters in database: 48,838,386
  Number of sequences in database:  30,364



Matrix: blastn matrix 1 -2
Gap Penalties: Existence: 0, Extension: 2.5

できた。


参考ページ

http://bioinfo-dojo.net/2016/10/25/blast_makeblastdb/

[]BLASTでトップヒットだけを抽出する

アノテーションデータベースが作りたいので、A種のgene aについてB種の遺伝子をBLASTで調べて、一番上に来るものを1:1対応させたい。

http://bioinfo-dojo.net/2016/03/25/blast_besthit_outfmt7/

ここにあるように、blastの結果から1行目だけをawkで取り出す。

その上でcutで先頭の2カラムを抽出すればいけるよね。

ワンライナーなら

$ blastn -query test.fasta -db mydb -outfmt 7 | awk '/hits found/{getline;print}' | grep -v "#" | cut -f 1,2 > anno.txt

これでばっちり。

あとはA種のgeneのfastaからすべての配列を順にblastにかけてannotationを作るスクリプトを考えないと。

と思ったけど、考えるまでもなく、普通に複数の遺伝子の並んだfastaを上のワンライナーにかけるだけで、目的のリストが完成してしまった。かかった時間も数十秒。ほんとにちゃんとできているのか心配になるな。

2018-03-27

[][] BLASTとの連携

生物種間で発現を比較しようとすると、アノテーションデータベースで橋渡ししてデータを見るということになるわけだが、アノテーションにはどうしても限界がある。

例えば生物Aで発現変動が見られた遺伝子について、生物Bでの発現を見る分にはアノテーションされた遺伝子を頼りに検索すれば良いのだが、更にC種での発現変動を調べたいとして、AからCへのアノテーションが構築されてなくて、CからBのアノテーションしかなかった場合、Bを経由した検索ではCでの発現データにはたどり着けないことが往々にして起こる。

Aの遺伝子aのアノテーションがBの遺伝子aだったとして、一次配列を元にCの配列にBLAST検索をかけてhitするCの遺伝子は遺伝子aなのに、Cの遺伝子aのアノテーションはBの遺伝子a'となっていて、一致しないわけだ。

確かにCの遺伝子aの一次配列でBの配列にBLAST検索をかけるとa'がトップヒットし、2番めにaがヒットしていたりするので、アノテーションが間違っていた、というわけではないのだ。

ややこしい。

やはりデータベースからBのアノテーションではなく直接BLAST検索結果で繋いでやらないと、この作業は正確に実施できないのだな。

データベースを拡張してBLASTサーバも自前で立ち上げる必要があるのだろうか。

それともC→A間のアノテーションを独自に構築するほうが簡単だろうか。

[] sshの活用

https://qiita.com/ik-fib/items/0af46a19e3083bb0325e

まだまだ使いこなせてなかった技が色々。

手始めにDDBJのスパコンのUGEのモニタリングを行うqmonを試す。

これはsshでアクセスする際に-Yオプションを付けておくだけでよい。

$ ssh -y gw2.ddbj.nig.ac.jp

これだけでhostで実行したXアプリケーションをローカルに表示してくれる。

$ xeyes

f:id:k-kuro:20180327190121p:image

なんか懐かしい感じ。

[]DDBJのスパコンで解析したデータをGCEから利用

DDBJの利用規定を見る限り問題ないように思うのだが、DDBJスパコンで解析したデータをGCEにおいたWEBアプリから参照してデータベース解析を出先からでも行えるようにしてみた。nfsでマウントできればそれに越したことはないのだが、DDBJのほうが想定していない使い方だと思われる。一方sshfsでマウントしてデータを読み書きする方法はDDBJのサイト内にちゃんと書かれていたので、特に問題ないはずだ。早速sshfsクライアントをGCEでインストールして、所定の手続きによって無事にGCEにDDBJのディスクがマウントできることを確認した。

ただ、結局無料枠では実行速度が実用的とは言えないのがなあ。

2018-03-19

[][] DDBJスパコンのbcftoolsが古い!

パイプラインが走らないとウンウンと苦悩したのだが、なんとbcftoolsのバージョンが0.1.17だった。ローカルにインストールしているのが1.7なので、古すぎる。どうりでbcftools callがunknown commandと言われるわけだ。

2時間返せ〜

2018-03-16

[] DDBJスパコン運用〜ファイル転送

DDBJのスパコンのアカウントを取ったのだけど、その環境整備をやっている時間が取れなかったので、そのまま放置していた。ようやくちょっと時間が取れたのでアクセスしてみようと思ったら、メンテナンスで停止していて入れないし。

ようやくメンテが終わったらしいので、ちょっといじってみよう。

データ解析関連のソフトウェアは大方入っているので、まずは自分のデータをアップロードするところから。

ASPERAという高速ファイル転送サーバ、クライアントシステムが使えるらしいのでクライアントをインストールしておく。

DDBJのサイトからリンクが貼られているのは

http://downloads.asperasoft.com/download_connect//

ここで、コマンドラインからファイルを転送するもの。scpコマンドと同じような感じでファイルを転送するのだな。

しかし、こっちはMacだしGUIからもファイルをやり取りできたほうが便利ではあるので、自分でASPERAのサイトを探索してみると、使えそうなソフトが2つ。Aspera DriveというものとAspera desktop clientというもの。

http://asperasoft.com/software/transfer-clients/drive/#overview-4250

http://asperasoft.com/software/client-options/desktop-client/

Driveのほうはディスクをマウントして使えるようにする感じなんだろうか?これは起動してみたがうまく接続してくれず、駄目っぽい。

一方のclientのほうはHDDのバックアップツールとかFTPクライアントのようなインターフェースでファイルをサーバとコピーしあうような感じのソフトだな。こちらはサーバ名とアカウント名を入れて、公開鍵認証するとすんなりと接続ができた。

早速ファイルを転送しているが、ネットの帯域が狭く、40Mbps程度しか速度が出てないな。500GBのデータ転送で丸一日以上かかる計算になる。

と思ったんだけど、なんか速度のリミットを制限できるようになっているっぽい。これ、どこまであげられるんだろう?最初45Mbpsに設定されていたのを100Mbpsに上げたら、転送速度が95Mbpsまで上がったので、調子に乗って1Gbpsまで上げてみたけど、100Mbpsは超えられないようだ。ローカルの環境がその程度なんだろうかね。f:id:k-kuro:20180316111950p:image

2018-03-15

[][] python3を整備

python3をlinuxbrewでインストールしたが、どうもすっきりしないのでここは一旦brew uninstallして王道どおりyumでインストールし直してみる。

https://weblabo.oscasierra.net/python3-centos7-yum-install/

基本このページの指示通りにインストールするのだがやはりtkinterでつまづく。

tk-develをインストールするとかいろいろ試してみたけどそれだけでは解決せず、結局

$ yum search python36u
Loaded plugins: fastestmirror
Loading mirror speeds from cached hostfile
 * base: mirrors.gigenet.com
 * epel: mirror.steadfast.net
 * extras: centos.mirrors.tds.net
 * ius: ord.mirror.rackspace.com
 * updates: mirror.team-cymru.org
============================================= N/S matched: python36u ==============================================
python36u-debuginfo.x86_64 : Debug information for package python36u
python36u-lxml-debuginfo.x86_64 : Debug information for package python36u-lxml
python36u-mod_wsgi-debuginfo.x86_64 : Debug information for package python36u-mod_wsgi
python36u-psycopg2-debuginfo.x86_64 : Debug information for package python36u-psycopg2
python36u-setproctitle-debuginfo.x86_64 : Debug information for package python36u-setproctitle
python36u-test.x86_64 : The self-test suite for the main python36u package
uwsgi-plugin-python36u-debuginfo.x86_64 : Debug information for package uwsgi-plugin-python36u
python36u.x86_64 : Interpreter of the Python programming language
python36u-debug.x86_64 : Debug version of the Python runtime
python36u-devel.x86_64 : Libraries and header files needed for Python development
python36u-gunicorn.noarch : Python WSGI application server
python36u-libs.x86_64 : Python runtime libraries
python36u-lxml.x86_64 : XML processing library combining libxml2/libxslt with the ElementTree API
python36u-mod_wsgi.x86_64 : A WSGI interface for Python web applications in Apache
python36u-pip.noarch : A tool for installing and managing Python packages
python36u-psycopg2.x86_64 : A PostgreSQL database adapter for Python
python36u-redis.noarch : Python interface to the Redis key-value store
python36u-setproctitle.x86_64 : Python module to customize a process title
python36u-setuptools.noarch : Easily build and distribute Python packages
python36u-tkinter.x86_64 : A GUI toolkit for Python
python36u-tools.x86_64 : A collection of tools included with Python including 2to3 and idle
uwsgi-plugin-python36u.x86_64 : uWSGI - Plugin for Python support
  Name and summary matches only, use "search all" for everything.

ここをよく見ると、python36u-tkinterというのがあることに気が付き、これをインストールに含めてやったところエラーが出なくなった。

$ sudo ln -s /usr/bin/python3.6 /usr/bin/python3

とpython3のリンクを作って

$ sudo python3 manage.py runserver

で無事に起動することができた。


速度であるが、ローカルにあるiMac 2.7 GHz Intel Core i5 16G mem (4core)で0.125secかかるSQLの処理に0.170secかかっており、webへのレンダリングにはかなり時間を要している。やっぱりメモリが0.6Gしかないのはきついのか?ラズパイといい勝負な気がする。