RNA-seq のde novo assembly

通常RNA-seqしたらreferenceのfastaファイルを使ってmappingして発現解析なりするわけだが、referenceが完備されていない種のseqはどうするのか？
近縁種のreferenceを使う、というのが簡単な手段なわけだが、今回mappingしてみるとmap rateが30％くらいしかなくて、ほとんどまともにマッピングできない、という事態となったので、次の手としてRNAのde novo assemblyを行って独自referenceを作ってしまおう、と考えた。
以前にも一度その近縁種のreferenceをガイドに使ってTrinityによるgenome guided assemblyを試していたのだが、これによって出来上がったreferenceではやはり30%くらいしかマップできなかった。要するにguideに使ったgenome referenceが全体の30％くらいしかカバーしていないので残りのread は結局mappingされないという訳のようだ。
なので、ガイドなしのde novo assemblyとなった。

Trinity - 井上 潤
とか
Trinity | de novo アセンブリー
とか
shortreadbrothers.blogspot.com
などを参考に

#!/bin/sh

Trinity --seqType fq \
        --single L001_trim.fastq.gz,L002_trim.fastq.gz,L003_trim.fastq.gz,L004_trim.fastq.gz \
        --CPU 8 \
        --max_memory 60G \
        --output ref/trinity \
        --monitoring \
        --no_cleanup 

こんな感じのスクリプトを作成し、実行。
メモリは100万リードあたり1Gくらいいるよ！という情報だったので60G盛ってるが、3000万リードを4分割してあるfastqを処理するとピークで15Gくらいメモリを消費している模様。まとめて3000万リード読み込むわけじゃないのかな。
というわけで16Gメモリのノードでもなんとかなりそうな感じ。
処理時間はかなり要することは間違いない。E3-1230v6の4コア8スレッド全部使っても丸１日？

出た。
単純なCPUスペック的にはE3-1230v6のほうが勝っているのだが、E5-2667の6コア12スレッド＊２CPUのほうが早く終わったよ。

3世代5CPU対決

RX1330M3 E3-1230 v6 3.5GHz 4C8T 224.0GFLOPS

18/11/06 13:06:26
[M::bwa_idx_load_from_disk] read 0 ALT contigs
[M::process] read 1815774 sequences (80000053 bp)...
[M::process] read 1833396 sequences (80000011 bp)...
[M::mem_process_seqs] Processed 1815774 reads in 69.270 CPU sec, 10.917 real sec
[M::process] read 1827330 sequences (80000067 bp)...
[M::mem_process_seqs] Processed 1833396 reads in 72.403 CPU sec, 11.981 real sec
[M::process] read 1837492 sequences (80000005 bp)...
[M::mem_process_seqs] Processed 1827330 reads in 72.093 CPU sec, 11.749 real sec
[M::process] read 419113 sequences (18201018 bp)...
[M::mem_process_seqs] Processed 1837492 reads in 70.634 CPU sec, 11.775 real sec
[M::mem_process_seqs] Processed 419113 reads in 16.366 CPU sec, 3.573 real sec
[main] Version: 0.7.17-r1188
[main] CMD: bwa mem -t8 Trinity-GG MNt_1_trim.fastq.gz
[main] Real time: 54.819 sec; CPU: 303.178 sec
18/11/06 13:07:21

55 sec

RX300S7 E5-2667 2.9GHz 6C12T*2 139.2GFLOPS/cpu

18/11/06 13:04:24
[M::bwa_idx_load_from_disk] read 0 ALT contigs
[M::process] read 5476498 sequences (240000039 bp)...
[M::mem_process_seqs] Processed 5476498 reads in 283.738 CPU sec, 42.144 real sec
[M::process] read 2256607 sequences (98201115 bp)...
[M::mem_process_seqs] Processed 2256607 reads in 119.757 CPU sec, 12.516 real sec
[main] Version: 0.7.17-r1188
[main] CMD: bwa mem -t24 Trinity-GG MNt_1_trim.fastq.gz
[main] Real time: 72.568 sec; CPU: 415.462 sec
18/11/06 13:05:37

73 sec

RX300S7 E5-2643 3.3GHz 4C8T*2 105.6GFLOPS/cpu

18/11/06 12:54:28
[M::bwa_idx_load_from_disk] read 0 ALT contigs
[M::process] read 3649170 sequences (160000064 bp)...
[M::mem_process_seqs] Processed 3649170 reads in 185.107 CPU sec, 37.522 real sec
[M::process] read 3664822 sequences (160000072 bp)...
[M::process] read 419113 sequences (18201018 bp)...
[M::mem_process_seqs] Processed 3664822 reads in 189.660 CPU sec, 28.777 real sec
[M::mem_process_seqs] Processed 419113 reads in 21.859 CPU sec, 3.837 real sec
[main] Version: 0.7.17-r1188
[main] CMD: bwa mem -t16Trinity-GG MNt_1_trim.fastq.gz
[main] Real time: 84.943 sec; CPU: 406.781 sec
18/11/06 12:55:53

85 sec

RX200S7 E5-2620 2.00GHz 6C12T*2 96.0GFLOPS/cpu

18/11/06 13:18:06
[M::bwa_idx_load_from_disk] read 0 ALT contigs
[M::process] read 5476498 sequences (240000039 bp)...
[M::mem_process_seqs] Processed 5476498 reads in 395.009 CPU sec, 58.175 real sec
[M::process] read 2256607 sequences (98201115 bp)...
[M::mem_process_seqs] Processed 2256607 reads in 165.597 CPU sec, 17.457 real sec
[main] Version: 0.7.17-r1188
[main] CMD: bwa mem -t24 Trinity-GG MNt_1_trim.fastq.gz
[main] Real time: 98.595 sec; CPU: 576.243 sec
18/11/06 13:19:45

99 sec

RX200S6 E5630 2.54GHz 4C8T*2 no-data

18/11/06 13:01:24
[M::bwa_idx_load_from_disk] read 0 ALT contigs
[M::process] read 3649170 sequences (160000064 bp)...
[M::mem_process_seqs] Processed 3649170 reads in 247.487 CPU sec, 49.641 real sec
[M::process] read 3664822 sequences (160000072 bp)...
[M::process] read 419113 sequences (18201018 bp)...
[M::mem_process_seqs] Processed 3664822 reads in 241.258 CPU sec, 36.402 real sec
[M::mem_process_seqs] Processed 419113 reads in 27.730 CPU sec, 4.817 real sec
[main] Version: 0.7.17-r1188
[main] CMD: bwa mem -t16Trinity-GG MNt_1_trim.fastq.gz
[main] Real time: 109.888 sec; CPU: 529.980 sec
18/11/06 13:03:14

110 sec

ということでスペック通り
E3-1230 v6 > E5-2667 > E5-2643 > E5-2620 > E5630
multi cpuは単純にflops*cpu数にはやっぱりならんのな。

はてなダイアリーが終了するらしいので、はてなブログへ移行することにする

というわけで引っ越してきた。

2018年10月31日19:20現在の訪問者数　1192384

ゼロから再スタート。

はてなダイアリーでできていたことがこっちではできなかったりするのな。
結構リンク元参照とか楽しめたのに。思いがけない出会いがあったり。

カレンダーから記事をたどるのも結構使ってたけどね。

UPSを導入ーterastationとサーバの両方を連動してシャットダウンー

停電や不慮の電源トラブルに備えてオムロンのUPSを導入した。

全サーバー群の電源を確保できるような大容量のUPSは高価なため、WEBサービス用データを保持しているヘッドノードとNFSサーバとして使っているTeraStationだけでもつないでおけるようにと、750Wのものを手配した。

うちのTeraStationはUPS連動して電源異常を感知すると一定時間後にシャットダウンできる機能がついている。一方、UPSの方にも同様のユーティリティソフトが付いていて、サーバのシャットダウンができるようになっているらしい。

しかし、ここで困ったことが。UPSから異常検知シグナルを送るポートはUSB1ポートしかない。
付属ソフトでは複数台のサーバを連動してシャットダウンできる。
TeraStationは複数台ネットに接続していれば連動してシャットダウンできる。
ところがサーバとTeraStationは連動できないみたい。

こりゃTeraStationの連動機能は使わずに、サーバのシャットダウンシークエンスにTeraStationのシャットダウンも組み込むしかないか？と思ったのだが、どうもTeraStationは普通に簡単なコマンド一発でシャットダウンできるようにはなっていないらしい。

どうすりゃいいのかというとwgetでwebベースのコントロールサイトにアクセスして、webからシャットダウンボタンを押したように見せかけないとならないみたい。
そのためにはBASIC認証でログインして〜と結構面倒そうだな。

http://www.argv.org/~chome/blog/noisefactory/2010/02/tera-stationremoteshutdown.html
古いTeraStationはそもそもUPS連動機能がなかったらしく、先人のスクリプトは発見した。
これをモディファイして、同じようなことができるようにしてやらないと。

RX300S7インストール覚書

RX300S7にCentOS7をインストール。
計算ノードとして。
300GB　SASドライブにシステムインストール。

一般ユーザーは作成しない。
NFSでヘッドノードの/home以下、/usr/local以下をマウント
NISでヘッドノードで作成したユーザ情報を共有

NFSのホスト：rx1330m3
NISのホスト：rx1330m3.nis
NAS（NFS)：192.168.1.10

今設定しようとしているドメイン：192.168.1.2 y0-rx200s7

# yum install ypbind rpcbind
# yum install nano
# nano /etc/hosts
# cat /etc/hosts
127.0.0.1   localhost localhost.localdomain localhost4 localhost4.localdomain4
::1         localhost localhost.localdomain localhost6 localhost6.localdomain6
192.168.1.1     rx1330m3 rx1330m3.nis
192.168.1.2    y0-rx200s7
192.168.1.3   y1-rx200s6
192.168.1.4   y2-rx200s6
192.168.1.5   y3-rx300s7
192.168.1.6     n1-rx300s7
192.168.1.7     n2-rx300s7

# authconfig --enablenis --nisdomain=nis --nisserver=rx1330m3 --enablemkhomedir --update
# systemctl start rpcbind ypbind
# systemctl enable rpcbind ypbind

NFSのauto mountを導入

# yum -y install nfs-utils
# nano /etc/idmapd.conf 
# cat /etc/idmapd.conf
[General]
#Verbosity = 0
# The following should be set to the local NFSv4 domain name
# The default is the host's DNS domain name.
Domain = y0-rx200s7
.........

# systemctl restart rpcbind
# nano /etc/fstab
# cat /etc/fstab
#
# /etc/fstab
# Created by anaconda on Tue Aug 13 13:07:11 2019
#
# Accessible filesystems, by reference, are maintained under '/dev/disk'
# See man pages fstab(5), findfs(8), mount(8) and/or blkid(8) for more info
#
/dev/mapper/centos-root /                       xfs     defaults	0 0
UUID=********-****-****-****-************ /boot                   xfs     defaults        0 0
/dev/mapper/centos-home /home                   xfs     defaults	0 0
/dev/mapper/centos-swap swap                    swap    defaults	0 0
rx1330m3:/usr/local  /usr/local              nfs     defaults	0 0
192.168.1.10:/mnt/array1/home /mnt/nfs  nfs     defaults  0 0

# yum -y install autofs
# nano /etc/auto.master
# cat /etc/auto.master
.........
/-    /etc/auto.mount

# nano /etc/auto.mount
# cat /etc/auto.mount
/home -fstype=nfs,rw rx1330m3:/home

# systemctl start autofs
# systemctl enable autofs

torque-client設定

# yum install -y epel-release
# yum install -y torque-client torque-mom
# nano /etc/torque/server_name
# cat /etc/torque/server_name
rx1330m3

# nano /var/lib/torque/mom_priv/config

# cat /var/lib/torque/mom_priv/config
# Configuration for pbs_mom.
$pbsserver rx1330m3
$usecp *:/home /home
$usecp *:/mnt/nfs /mnt/nfs
$log_file_suffix %h

# firewall-cmd --add-port=15001/tcp --zone=public --permanent
# firewall-cmd --add-port=15002/tcp --zone=public --permanent
# firewall-cmd --add-port=15003/tcp --zone=public --permanent
# firewall-cmd --add-port=15003/udp --zone=public --permanent
# firewall-cmd --add-port=15004/tcp --zone=public --permanent
# systemctl start pbs_mom
# systemctl enable pbs_mom

あとはbrewではなくyumでそれぞれのノードにインストールしているものを入れる。

# yum -y install java
# yum -y install perl

bedtoolsでread coverage

bedtoolsでbamファイルのread coverageを求める方法は以前にも書いた。
現在のbedtoolsのバージョン(2.27.1)ではbedファイルではなくbamファイルから直接read coverageを出すことが可能になっている

$ bedtools coverage -a referense.bed -b sample.bam > sample_cov.bed

というフォーマットで問題なく解析できるのだが。
ただしこの操作は非常にメモリを食うらしく、4GB超のbamファイルを解析しようとすると64GBのメモリを使い切ってスワップも15GBくらい使って、どうにもならなくなって落ちる。

ところが、

$ bedtools bamtobed -i sample.bam > sample.bed

と一旦bedに変換してから、先程のコマンドにつなげると、bedへの変換にはほとんどメモリを消費せず、coverageの計算も30GB程度のメモリで十分に収まることがわかった。

というわけで、bam直の解析はやめといたほうが良い。

IGV.jsの使い方は大体つかんだが、時々思うようにいかない

ということでjavascriptベースの他の選択肢もあたってみたい。

pileup.js
http://www.hammerlab.org/pileup/
これはシンプルでいいかもしれないな。
しかしnpmというパッケージシステムでインストールする？Node.jsってなに？また新しいことをいろいろ覚えないとならないのかね。

genomicviz
https://github.com/akmorrow13/genomicviz
pileup.jsをpythonに埋め込んでる？

biodalliance
http://www.biodalliance.org
これいいんじゃない？

そもそもIGV.jsでうまく表示されないのはBAMファイルからbedtools intersectでデータを処理したあとなんだが、bedtoolsの処理に問題があってBAMファイルがおかしくなっているんじゃないかな。

その点を調査してみよう。

Safari12でプラグインが使えなくなったときWEB認証のVPNをどうやって使うか

どうやらSafariがバージョンアップし、adobe flash player以外のプラグインがことごとく使用不能になったらしい。
そこで困ったのが大学のサーバへのVPN接続。
webからVPN接続のための認証を行っていたのに、mac_sslvpnが使えなくなって接続できなくなった。

そういえば接続が成立するとBIG-IPなんちゃら〜という表示が出ていたことを思い出し、WEBを検索。
スタンドアロンのVPNクライアントソフトウェア（接続ツール）というものがあるらしく、
一応App storeにも登録されたソフトウェアがあった。
https://itunes.apple.com/us/app/f5access/id1243219105?mt=12
これを起動し、Manage VPN Configurationsを開く。
いつもメールで送られてくる接続情報のリンクをサーバとして登録
ID・パスワードは都度変わるので空欄のままでよい
Connect to [作成した接続名]というところをクリックしてID・パスワードを入力。
無事接続。

reference-guided de novo assembly

ゲノムデータをreferenceにRNA-seqデータをマッピングしているのだが、ゲノムデータで使われているものとはことなる品種のRNA-seqデータをマッピングしてみたところ、かなりたくさんのSNPが含まれていることがわかった。
普通の発現解析なら、SNPがあろうがマッピングさえできれば良いわけだけど、いまやろうとしている仕事では100％マッチでマッピングされていることが重要なので、これではまずい。

ということでreferenceを使用品種のRNAseqデータで再構築してやりたい。本来的にはゲノムseqデータでde novo assemblyすれば良いわけだけど、これは大規模すぎるため、まずはRNA-seqデータでexon領域だけ、それも発現している遺伝子の部分だけ再構築してやることを目指す。

調べてみるとreference-based assemblyを行えば、referenceを再構築できそうなのだが、
https://bi.biopapyrus.jp/rnaseq/assembly/

肝心の方法はなかなか情報が得られない。
上のページでもできる、とは書いておきながら方法は書かれていないわけで。

http://seqanswers.com/forums/showthread.php?t=73620
ここでようやくそれらしい記述が。

velvetというとさっきのページでも紹介されていたゲノムのde novo assemblyをde Bruijn graph（ド・ブラングラフと読むらしい）で行うソフトウェアらしいのだが、これに付随するColumbus extensionを使うとreference-guided assemblyもできるらしい。

てことで
https://www.ebi.ac.uk/~zerbino/velvet/
からダウンロードして環境構築してみる。

普通に適当なところに解凍し、makeするだけで使えるようになる。
あとはとにかくマニュアルに従うだけなのだが、まずはColumbus extensionのマニュアルを紐解く。

1 For impatient people
> head myRegions.fa
>chr1:123456789-123457789
ATGTGTGTACTAGCTAGCGCGCTAGCTAGTCATGTGTGTACTAGCTAGCGCGCTAGCTAGTC
[etc ...]
> sort myReads.sam > mySortedReads.sam
> velveth my_dir 21 -reference myRegions.fa \
-shortPaired -sam mySortedReads.sam
> velvetg my_dir [etc ...]

せっかち向けって・・・
これだけ見てもまあどうすりゃいいかは案外わかるな。

3世代のサーバの能力を検証してみた

ベンチマークソフトを使ってもいいけど、実際に仕事に使うスクリプトを処理するのにかかる時間を計測したほうが意味があるだろう。
ということで、試しにBAMファイルをcufflinksにかけて遺伝子発現量を算出させてみた。
処理してみたBAMファイルは
1, 4.4GBリード数74139744
2, 416KBリード数8653
3, 11KBリード数116
の3ファイル
サーバのスペックは
A, Intel(R) Xeon(R) CPU E3-1230 v6 @ 3.50GHz 4core 8thread mem:64GB
B, Intel(R) Xeon(R) CPU E5-2620 0 @ 2.00GHz 6core 12thread x 2 mem:16GB
C, Intel(R) Xeon(R) CPU E5630 @ 2.53GHz 4core 8thread x 2 mem:26GB

さて、

$ cat arrayjob180901.sh.o379-15 
18/09/01 21:58:26
[21:58:27] Loading reference annotation.
[21:58:29] Inspecting reads and determining fragment length distribution.
Processed 31122 loci.                       
> Map Properties:
>	Normalized Map Mass: 29205450.00
>	Raw Map Mass: 29205450.00
>	Fragment Length Distribution: Truncated Gaussian (default)
>	              Default Mean: 200
>	           Default Std Dev: 80
[22:00:12] Estimating transcript abundances.
Processed 31122 loci.                       
finish cufflinks
 
[22:02:38] Loading reference annotation.
[22:02:39] Inspecting reads and determining fragment length distribution.
Processed 31122 loci.                       
> Map Properties:
>	Normalized Map Mass: 8924.00
>	Raw Map Mass: 8924.00
>	Fragment Length Distribution: Truncated Gaussian (default)
>	              Default Mean: 200
>	           Default Std Dev: 80
[22:02:44] Estimating transcript abundances.
Processed 31122 loci.                       
finish cufflinks2

[22:03:03] Loading reference annotation.
[22:03:05] Inspecting reads and determining fragment length distribution.
Processed 31122 loci.                       
> Map Properties:
>	Normalized Map Mass: 2074.00
>	Raw Map Mass: 2074.00
>	Fragment Length Distribution: Truncated Gaussian (default)
>	              Default Mean: 200
>	           Default Std Dev: 80
[22:03:08] Estimating transcript abundances.
Processed 31122 loci.                       
finish cufflinks3
18/09/01 22:03:27
finish
18/09/01 22:03:27

こんな感じに処理は行われた。
A,
1, 4min11sec
2, 25sec
3, 24sec
total 5min1sec

B,
1, 5min59sec
2, 36sec
3, 33sec
total 8min10sec

C,
1, 6min41sec
2, 45sec
3, 42sec
total 8min10sec

データを書き込む処理などで実際に計算にかかった時間より余分がついてはいるがざっとこんな感じかな。
Bの24スレッドがあんまり伸びないのはメモリがボトルネックになっているのかもしれない。クロックも控えめだしな。
なんだかんだいって最新CPUが一番速いという当たり前のような結果となった。
まあそうじゃなかったら新しい機種がどんどん開発されている意味がないわな。