通常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のほうが早く終わったよ。
18/11/06 19:24:10
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").
______ ____ ____ ____ ____ ______ __ __
| || \ | || \ | || || | |
| || D ) | | | _ | | | | || | |
|_| |_|| / | | | | | | | |_| |_|| ~ |
| | | \ | | | | | | | | | |___, |
| | | . \ | | | | | | | | | | |
|__| |__|\_||____||__|__||____| |__| |____/
Single read files: $VAR1 = [
'./fastq_trimmed/L001_R1_001_trim.fastq.gz',
'./fastq_trimmed/L002_R1_001_trim.fastq.gz',
'./fastq_trimmed/L003_R1_001_trim.fastq.gz',
'./fastq_trimmed/L004_R1_001_trim.fastq.gz'
];
Trinity version: Trinity-v2.5.1
** NOTE: Latest version of Trinity is Trinity-v2.8.4, and can be obtained at:
https://github.com/trinityrnaseq/trinityrnaseq/releases
Tuesday, November 6, 2018: 19:24:12 CMD: java -Xmx64m -XX:ParallelGCThreads=2 -jar /home/linuxbrew/.linuxbrew/Cellar/trinity/2.5.1/util/support_scripts/ExitTester.jar 0
Tuesday, November 6, 2018: 19:24:12 CMD: java -Xmx64m -XX:ParallelGCThreads=2 -jar /home/linuxbrew/.linuxbrew/Cellar/trinity/2.5.1/util/support_scripts/ExitTester.jar 1
Tuesday, November 6, 2018: 19:24:12 CMD: mkdir -p /mnt/nfs/rnaseq/ref/trinity
Tuesday, November 6, 2018: 19:24:13 CMD: mkdir -p /mnt/nfs/rnaseq/ref/trinity/chrysalis
STARTING COLLECTL
Tuesday, November 6, 2018: 19:24:13 CMD: rm -rf /mnt/nfs/rnaseq/collectl
I'M THE PARENT, COLLECTL_PID=17752
I'M THE CHILD RUNNING TRINITY
Running CMD: bash -c "set -eou pipefail && cd /mnt/nfs/rnaseq/collectl && exec nohup /home/linuxbrew/.linuxbrew/Cellar/trinity/2.5.1/trinity-plugins/COLLECTL/collectl/collectl --procopts w --procfilt u1001 --interval :60 -sZ -oD | egrep 'Trinity|trinityrnaseq|Inchworm|Chrysalis|Butterfly|jellyfish|samtools|sort|bowtie' | egrep -v 'egrep' > /mnt/nfs/rnaseq/collectl/collectl.dat "
----------------------------------------------------------------------------------
-------------- Trinity Phase 1: Clustering of RNA-Seq Reads ---------------------
----------------------------------------------------------------------------------
---------------------------------------------------------------
------------ In silico Read Normalization ---------------------
-- (Removing Excess Reads Beyond 50 Coverage --
---------------------------------------------------------------
中略
/Trinity.fasta.tmp to /mnt/nfs/rnaseq/ref/trinity/read_partitions/Fb_0/CBin_270/c27066.trinity.reads.fa.out/Trinity.fasta
succeeded(90209) 100% completed.
All commands completed successfully. :-)
** Harvesting all assembled transcripts into a single multi-fasta file...
Wednesday, November 7, 2018: 08:58:51 CMD: find /mnt/nfs/rnaseq/ref/trinity/read_partitions/ -name '*inity.fasta' | /home/linuxbrew/.linuxbrew/Cellar/trinity/2.5.1/util/support_scripts/partitioned_trinity_aggregator.pl TRINITY_DN > /mnt/nfs/rnaseq/ref/trinity/Trinity.fasta.tmp
-relocating /mnt/nfs/rnaseq/ref/trinity/Trinity.fasta.tmp to /mnt/nfs/rnaseq/ref/trinity/Trinity.fasta
Wednesday, November 7, 2018: 13:53:19 CMD: mv /mnt/nfs/rnaseq/ref/trinity/Trinity.fasta.tmp /mnt/nfs/rnaseq/ref/trinity/Trinity.fasta
Wednesday, November 7, 2018: 13:53:19 CMD: /home/linuxbrew/.linuxbrew/Cellar/trinity/2.5.1/util/support_scripts/get_Trinity_gene_to_trans_map.pl /mnt/nfs/rnaseq/ref/trinity/Trinity.fasta > /mnt/nfs/rnaseq/ref/trinity/Trinity.fasta.gene_trans_map
###################################################################
Trinity assemblies are written to /mnt/nfs/rnaseq/ref/trinity/Trinity.fasta
###################################################################
TERMINATING COLLECTL, PID = 17752
finish
18/11/07 13:53:36417.1MB+422.5MB+422.4MB+412.1MB=1673.1MB
start:18/11/06 19:24:10
end: 18/11/07 13:53:36
1.67GBのファイルで18時間29分26秒
出来上がったfastaファイルは144.5MB
logファイルだけで43.3MBのテキストという超大仕事だな。
