ゲノム情報はないが近縁種のゲノムが解読されているような生物でRNA seqを行うと決まったら、どんなワークフローで進めるべきだろうか？マイクロアレイと違い、RNA seqならde novoでも解析は不可能ではない。ゲノムがモデル生物とほぼ同じならば、深く考えなければ、なんだかできそうに思えてくる。
実際、SRAのデータベースを見ているとそういったことを狙ったであろう解析データが無数に登録されている。その中にはpublish情報がないものもある。これが何を意味するかは置いておいて、モデル生物のゲノム情報の充実に伴い、そのゲノム情報を利用した近縁種の解析はこれからもますます増えてくると予想される。
ただしモデル生物に近縁でも塩基置換が起こってるだけでアライメントは難しくなる。ではde novo assemblyして進めるか？、というとせっかくゲノムが決まっているのだからなんとか利用したいという気持ちもあり、ジレンマに陥る。場合によりけりだと思うが、どういったワークフローが妥当なのか１つの例を通して考えてみることにする。
今回はアラビドプシス属のA.kamchaticaのRNA seqデータを使う。シロイヌナズナとゲノムはかなり似ていると思われるが、ゲノムはまだ解読されていない。シロイヌナズナのecotypeでも塩基置換がかなり起きていることが知られており、A. kamchaticaのRNA seqデータがどれくらいナズナのリファレンスゲノムにアライメントされるかは未知数である。
今回はアセンブルしてリード数などの比較を行うが、それ以上は深入りせずプロトコルに抽象性を残しておく。その方が汎用性のあるプロトコルになると考える。
Arabidopsis. kamchatica (アラビドプシス属）
データ
Transcriptome analysis of arabidopsis kamchatica under normal and stress conditions
control1のデータ量
リード数 3600万 x 2
リード長 101-bp
ファイルサイズ 12.6GB x 2
万能なワークフローは存在しないが、今回はまずgenome guideのアセンブリと、de novoのアセンブリを実行し、それぞれの手法で作ったcDNA (contig)をTAIRのcDNAと比較する。比較するのは、生じたcDNAの数と長さ、core geneの存在数、TAIR cDNAとの相同性などになる。TAIR cDNAと比較してArabidopsis kamchatica固有の配列の割合を見積もることで、発現解析にはどのようなcDNAセットを使うのが妥当なのか分析する。（ここで言う優良なcDNA セットとは、単純に漏れが無くてダブりも無いセット）。
マッピングに使うアセンブルデータが決まったら、それをTAIR のcDNAとマージし、さらにcd-hit-estを使って95%以上同一の重複配列を廃除する。
１、contigの作成。
genomeガイドの方法１つと、de novo assemblyの手法２つをを心みる。
1、 genome-guided assembly
tophatでのmapping
bowtie2-build -f A_thaliana_TAIR10.fa chr
tophat2 -p 4 -I 50000 -i 20 -o tophat chr control1/DRR013381_1.fastq control1/DRR013381_2.fastq
マッピング率の分析
#アライメントされたリード
samtools flagstat accepted_hits.bam
15901037 + 0 in total (QC-passed reads + QC-failed reads)
1879156 + 0 secondary
0 + 0 supplementary
0 + 0 duplicates
15901037 + 0 mapped (100.00% : N/A)
14021881 + 0 paired in sequencing
7007096 + 0 read1
7014785 + 0 read2
6438422 + 0 properly paired (45.92% : N/A)
6950656 + 0 with itself and mate mapped
7071225 + 0 singletons (50.43% : N/A)
57762 + 0 with mate mapped to a different chr
22166 + 0 with mate mapped to a different chr (mapQ>=5)
#アライメントされなかったリード
samtools flagstat unmapped.bam
58401591 + 0 in total (QC-passed reads + QC-failed reads)
0 + 0 secondary
0 + 0 supplementary
0 + 0 duplicates
0 + 0 mapped (0.00% : N/A)
58401591 + 0 paired in sequencing
29204640 + 0 read1
29196951 + 0 read2
0 + 0 properly paired (0.00% : N/A)
0 + 0 with itself and mate mapped
0 + 0 singletons (0.00% : N/A)
0 + 0 with mate mapped to a different chr
0 + 0 with mate mapped to a different chr (mapQ>=5)
ナズナゲノムにアライメントされたリード1590万、アライメントされなかったリード5840万。マッピング率はかなり悪い。この数値だとリファレンスベースで進めるとかなりの取りこぼしが出ると予想される。
ただし今回はアライメントの感受性は変えずこのまま進める。デフォルトのパラメータでどこまで妥当なデータが出せるか考える。
Cufflinks
cufflinks --no-update-check --overlap-radius 1 -o cufflinks accepted_hits.bam
- –overlap-radius <int> Transfrags that are separated by less than this distance get merged together, and the gap is filled. Default: 50 (in bp).
- –no-update-check Turns off the automatic routine that contacts the Cufflinks server to check for a more recent version.
Cuffmerge
echo cufflinks/transcripts.gtf > assemblies.txt
cuffmerge -s A_thaliana_TAIR10.fa assemblies.txt
/merged_asm にmerged.gtfができる。featureの数は61958だった。
ここからbedに変換してbedtoolsでfasta配列を抽出し、de novoで作ったfastaと比較する。
cat merged.gtf | awk '{OFS = " "} {print $1,$4,$5,$3,$6,$7}' > merged.bed
#コマンド省略 perlのスクリプトでbedの各featureのサイズを計算、７カラム目に追加= > sized.bed として出力
#サイズ順にソート
sort -k 7 -n -r sized.bed > size_sorted.bed
#コマンド省略 100bp以下を除去して保存
#bedtoolsでfastaを抽出。
bedtools getfasta -fi ../A_thaliana_TAIR10.fa -bed size_sorted.bed |fold -w 60 > cufflinks.fasta
featureの数は30422だった。大きく数が減ったのは、cuffmergeで0bpの配列が大量にできており、また100bp以下の非常に短いfeatureも大量にあったためである。
2、 de-novo assembly - rnaspades
rnaspades.py -t 20 -k auto -o RNA_spades -1 R1.fastq -2 R2.fastq
合計68374 contigできた。そのうち 100bp以上のcontigは68176だった。diploidのspeciesに使ってしまってますが、本来は推奨されていません、ご注意ください。
3、 de-novo assembly - trinity
Trinity --seqType fq --left R1.fastq --right R2.fastq --max_memory 40G --CPU 20
合計78657 contigできた。 そのうち 100bp以上のcontigは78656だった（１つだけ21 bp）。
アセンブルでエラーを起こした。フォーマットが古いせいだと思われたので、ヘッダーを修正して３行目を+に修正してランした。
4、 genome-guided assembly - trinity
trinityでも、ゲノムを参照しながらアセンブルする方法があるのを忘れていました。ここでは結果なしになりますが、リンクを載せて起きます。
https://github.com/trinityrnaseq/trinityrnaseq/wiki/Genome-Guided-Trinity-Transcriptome-Assembly
２、buscoによるcore geneの探索。
BUSCO -g input.fasta -o OUTPUT -l embryophyta_odb9 -m trans -c 8
1、 genome-guided assembly
#BUSCO was run in mode: trans
Summarized benchmarks in BUSCO notation:
C:0%[D:0%],F:0%,M:0%,n:1440
13 Complete BUSCOs
13 Complete and single-copy BUSCOs
0 Complete and duplicated BUSCOs
113 Fragmented BUSCOs
1314 Missing BUSCOs
1440 Total BUSCO groups searched
13しかできていない。
2、 de-novo assembly - rnaspades
#BUSCO was run in mode: trans
Summarized benchmarks in BUSCO notation:
C:0%[D:0%],F:0%,M:0%,n:1440
937 Complete BUSCOs
668 Complete and single-copy BUSCOs
269 Complete and duplicated BUSCOs
164 Fragmented BUSCOs
339 Missing BUSCOs
1440 Total BUSCO groups searched
937見つかった。339 missingがある。
3、 de-novo assembly - trinity
#BUSCO was run in mode: trans
Summarized benchmarks in BUSCO notation:
C:0%[D:0%],F:0%,M:0%,n:1440
1078 Complete BUSCOs
797 Complete and single-copy BUSCOs
281 Complete and duplicated BUSCOs
105 Fragmented BUSCOs
257 Missing BUSCOs
1440 Total BUSCO groups searched
1078見つかった。257 missingがある。
参考; TAIR cDNA
Summarized benchmarks in BUSCO notation:
C:0%[D:0%],F:0%,M:0%,n:1440
1357 Complete BUSCOs
896 Complete and single-copy BUSCOs
461 Complete and duplicated BUSCOs
17 Fragmented BUSCOs
66 Missing BUSCOs
1440 Total BUSCO groups searched
1357見つかる。TAIRのcDNAにも66 missingがある。
まとめ
cuffmergeはかなり厳しい。このケースではゲノムガイドアセンブリは現実的ではないようだ。 trinityがベストのように思える。
#BUSCOは以前紹介しました。簡単な使い方をまとめています。
アセンブル結果をCore gene setの検出数で評価する BUSCO - macでインフォマティクス
３、contigの長さの比較
check_fasta_length.pl RNA_spades/transcripts.fasta |cut -f 2 - > spades.txt #spadesのfastaの長さを抽出
check_fasta_length.pl trinity/transcripts.fasta |cut -f 2 - > trinity.txt #trinityのfastaの長さを抽出
check_fasta_length.pl Arabidopsis_thaliana.TAIR10.cdna.all.fa |cut -f 2 - > TAIRcDNA.txt #TAIRcDNAのfastaの長さを抽出
Rにて
spades=scan("spades.txt")
trinity=scan("trinity.txt")
cuffmerge=scan("cufflinks_contig")
TAIRcDNA=scan("TAIRcDNA.txt")
hist(spades, col = "#ff00ff40", breaks = 150)
hist(cuffmerge, col = "#00A47450", breaks = 150, add = TRUE) #グラフに追加
hist(TAIRcDNA, col = "#2F343720", breaks = 300, add = TRUE) #グラフに追加
黒がTAIRのcDNAの分布で、ピンクがrnaspadesでアセンブルしたcontigの分布、緑がcuffmergeでゲノムベースにアセンブルしたcontigの分布。縦軸は頻度。横軸は長さ(bp)。
cuffmergeで作られた領域は短いものが非常に多い。
次にtrinityとrnaspadesの比較。
hist(spades, col = "#ff00ff40", breaks = 150)
hist(trinity, col = "#0000ff40", breaks = 150, add = TRUE) #グラフに追加
ピンクがrnaspadesでアセンブルしたcontigの分布、青がTrinityでセンブルしたcontigの分布。Trinityの方が長いcontigが多い。
４、contigのblast解析
３つの方法で作成したcontigをTAIRのcDNAと比較する。blastnを使う。
まずblastのデータベースを作成。
makeblastdb -dbtype nucl -in Arabidopsis_thaliana.TAIR10.cdna.all.fa -out TAIR10_cDNA_db
相同性検索を実行。出力指定は可能だがdefaultで行う（リンク）。
#RNA_spades
blastn -query RNA_spades/transcripts.fasta -evalue 1e-10 -outfmt 6 -max_target_seqs 1 -db TAIR10_cDNA_db > RNA_spades_out.txt
#trinity
blastn -query trinity/transcripts.fasta -evalue 1e-10 -outfmt 6 -max_target_seqs 1 -db TAIR10_cDNA_db > trinity_out.txt
- -outfmt 6 タブ区切り出力
出力されたファイルにはe-valueが1e-10以下を満たす配列が保存されている。中身を確認しておく。
| 1 | 2 | 3 | 4 | 5 | 6 | 7 | 8 | 9 | 10 | 11 | 12 |
| AT1G19090.1 | AT1G19090.1 | 99.88 | 2509 | 3 | 0 | 1 | 2509 | 1 | 2509 | 0 | 4617 |
１、２列目がクエリーとデータベースの名前である。さらに3列目が相同性(%)、４列目がhitした長さ、５列目がhitした領域中のミスマッチの数、６列目がギャップの数である。１２列目にスコア (Bits) も出ている。
blastの結果をまとめる。48359あるTAIRのcDNAとヒットしたのは (e-value <1e-10)、
cuffmergeの30422 contigのうち30011 contig
rnaspadesの68176 contigのうち50715 contig
Trinityの78657 contigのうち 63575 contig
hitしなかったcontigはシロイヌナズナ ゲノムにはない配列だったか、シロイヌナズナゲノムに存在する配列だが、TAIR cDNAと転写領域が異なるためcDNAカタログ化されていない配列と考えられる。
2−4のまとめ
ここまでの解析で、シロイヌナズナのゲノムベースでアセンブルするより、de novoでアセンブルするデータの方が漏れが少なく、cDNAのサイズも長くなることがわかった。
ちなみに今回はリードの長さを分析してblast解析するという泥臭い方法をとったが、TransRateというツールを使うと、リファレンスとの相同性やエラー率など複数の側面からアセンブリを評価してくれる。
TransRateはこちらで紹介しています。
アセンブル結果をリードのアライメントパターンから評価する TransRate - macでインフォマティクス
この後どうすべきかであるが、万能な方法はない。他の人はどうしているのだろう？例えばこの論文ではrudandantなcontigを除く方法を提案している。
最初書いた時はシロイヌナズナの情報を諦めきれず、contig配列とシロイヌナズナの配列を１つにしてcd-hit-estで選抜してより長い配列を使うことを考えていた。しかしよく考えてみると、cd-hitで選抜すると、完全長のcDNAになっているシロイヌナズナの配列が優先されてde novoで作ったcDNAは捨てられ、結局アライメントでうまく行かなくなってしまう可能性が高い。またそれこそ恣意的なデータセットを使っていると取られかねない。
ということで、みんながやっているようにdenovo assmeblyで作ったcontigのみ使うことにする。解析に供した全サンプルのリードを使ってアセンブルすれば、完全なcDNAコレクションにはなつていないとしても、少なくともそのデータに関しては定量漏れの遺伝子が出る可能性は低い。
たたし、denovo assmeblyのデータだけ使うなら、アセンブルの精度が実験の品質を完全に左右することになるので、アセンブル精度を客観的に評価できる手法でアセンブル精度を担保しておく必要が出てくる。
以下は最近アクセプトされた論文での流れ。
2016年 nature scientific report de nono transcriptomeのペーパー
１、trinityを使い、k-merが25、32でアセンブリ。
↓
２、できたcontigをBUSCOとTransRateで評価。
↓
３、アセンブリには適切なk-mer coverageがあるとの先行研究から、ランダムサンプリングでリード数を減らし再アセンブル。
↓
４、１のデータと３のデータをBUSCOとTransRateで再評価。
↓
５、ランダムサンプリングしてリード数を減らしk-mer=32でアセンブルのがベストだと判断。
↓
６、Trinotateでアノテーション。
↓
７、 blastで植物ゲノムとベストヒットしないものを排除。
このように、k-mer、カバレッジ、などを試行錯誤してアセンブルし、それから評価ツールで良好なアセンブル条件を比較する流れを取っている。
今回のケースで、近縁種のゲノムが出ていたとしても、直接鋳型としてgenome guide assemblyをするのは厳しいことを身をもって実感した。ゲノムガイドとde novoのどちらも検討してもよいが、最終的には、近縁種のモデル生物のゲノム情報はアノテーション情報として利用するのが精一杯と思われる。
追記１
下のリンク先の論文は、複数のアセンブルツールを併用し、cd-hit、TGI clustering tools、EvidentialGeneを使ってredundancyを下げることで、単一のアセンブルツールより高品質なRNAのアセンブルを行なっている。
こちらもそう
https://dl.sciencesocieties.org/publications/tpg/pdfs/8/1/plantgenome2014.10.0064
TGI clustering toolsのpaper
https://academic.oup.com/bioinformatics/article/19/5/651/239299
追記２
追記３
発現パラーンからクラスター化する。
追記４
アセンブルの評価
Transrate
BUSCO
追記５
Evaluation of de novo transcriptome assemblies from RNA-Seq data
https://genomebiology.biomedcentral.com/articles/10.1186/s13059-014-0553-5
追記６
Differential expression analysis of de novo assembled transcriptomes - Nadia Davidson from Australian Bioinformatics Network
以下は、cd-hitを使ってシロイヌナズナのcDNAと統合しようと考えていた時の流れである。論文投稿のレベルを満たした手法とは思えないが、残しておく。
５、マージ
ここから先は２つの手段を取る。
A、アセンブルして得られたcDNAからシロイヌナズナに見つからないcDNAだけ選り分け、それをTAIR cDNAとマージしてリファレンスとする。
B、アセンブルで得られたcDNAとTAIR cDNAをそのままマージし、CD-HIT-ESTで重複を除きリファレンスとする。
Aの方法を実行するには、シロイヌナズナのcDNAと相同性を示さないcDNAを集める必要がある。そこで、step４でblast解析した結果を相同性を示すcDNAのデータベースとして、no hitのcDNAをrnaspadesのアセンブルfastaから抽出する作業を行なった。そうゆうオプションやツールがありそうだが（絶対ある）、すぐには見つからなかったのでperlのスクリプトを書いて処理した。
得られたno-hitのcDNA.fastaをblastにかける。organismはシロイヌナズナに限定し、最初の100 contigをblast解析する。e-valueの閾値はデフォルトでランした。
no-hitのcDNAだけをクエリにしているので、 スクリプトが上手くワークしているなら、シロイヌナズナのゲノムと全体にわたってマッチする配列は0のはずである。
結果は大まかに３パターンに分かれた。
highly homologous with A.thaliana genome
partially homologous with A.thaliana genome
no siginificant holology with A.thaliana genome
予想通り多くのcDNAはナズナゲノムと高い相同性は示さなかった。このことは、Arabidopsis. kamchatica 固有の配列が抽出されていることを意味する。しかしながら、blastにかけた100 cDNAのうち１割ほどのcDNAはシロイヌナズナゲノムと配列全体に渡って高い相同性を示した。恐らくは、スプライシングバリアントというよりTAIRのcDNAカタログと全く異なる転写領域のため、blastで相同性を示さなかったのだと思われる。
（シロイヌナズナのcDNAの再解析の論文も報告されているが、そういったデータと比べるとどうなるだろうか？ http://www.biorxiv.org/content/biorxiv/early/2016/04/05/047308.full.pdf ）
得られたno-hit CDAとTAIR 10 cDNAをマージして、no-hit_and_TAIRcDNA.faとした。step6のマッピングに進む。
次に Bの方法について書く。
まずrnaspadesで得られたfastaとTAIR cDNAを１つにして"TAIR_and_rnaspades_merged.fa"を作った。このファイルを入力としてcd-hit-estでクラスター化する。
cd-hit-est -i TAIR_and_rnaspades_merged.fa -o TAIR_and_rnaspades_merged_cd-hit_clustered.fa -c 0.95 -T 12
- -i input filename in fasta format, required
- -o output filename, required
- -c sequence identity threshold, default 0.9 this is the default cd-hit's "global sequence identity" calculated as: number of identical amino acids in alignment
- -T number of threads, default 1; with 0, all CPUs will be used.
閾値は95%とした。これ以上の相同性を示すcDNAは１つにまとめられる。
解析が終わると、クラスター化されて95%以上似た配列が除去されたfastaと、どのfastaがクラスターになったかを示す~fa.clstrができる。中身を見てみる。
user$ head -10 cd-hit_est_clustered.fa.clstr
>Cluster 2
0 15464nt, >AT3G02260.4... at +/98.96%
1 15700nt, >AT3G02260.2... *
2 15527nt, >AT3G02260.3... at +/99.98%
3 943nt, >NODE_18797_length_9... at -/98.09%
4 412nt, >NODE_35931_length_4... at +/99.03%
5 284nt, >NODE_45670_length_2... at +/97.18%
6 221nt, >NODE_54874_length_2... at -/96.38%
7 120nt, >NODE_65872_length_1... at -/99.17%
8 109nt, >NODE_67550_length_1... at +/93.58%
このクラスターだと8つのcDNAからなる。それがアスタリスクの ">AT3G02260.2 "にまとめられている。ただしcDNAの長さや塩基が編集を受けるわけではない。AT3G02260.2はクラスター化する前も後も15700-bpで100%同一の配列である。
cd-hit-estにより、116692あった配列が69650に減った。
ヒストグラムを書く。
ピンクがcd-hit-estにかける前のcDNA、黒がcd-hit-estにかけた後のcDNAである。横軸はcDNAの長さ、縦軸が頻度である。特に短い配列が多く減っているように見える。
BUSCOを使い、 クラスター化でcore geneが減ってしまっていないか確認する。
BUSCO -g input.fasta -o OUTPUT -l embryophyta_odb9 -m trans -c 8
cd-hit-estにかける前のfasta
#BUSCO was run in mode: trans
Summarized benchmarks in BUSCO notation:
C:0%[D:0%],F:0%,M:0%,n:1440
1395 Complete BUSCOs
414 Complete and single-copy BUSCOs
981 Complete and duplicated BUSCOs
22 Fragmented BUSCOs
23 Missing BUSCOs
1440 Total BUSCO groups searched
TAIRのcDNAを含むので、core geneの保有数は1440のうち1395まで上がっている(rnaspades単独のデータだと937)。
cd-hit-estにかけた後のfasta
#BUSCO was run in mode: trans
Summarized benchmarks in BUSCO notation:
C:0%[D:0%],F:0%,M:0%,n:1440
1335 Complete BUSCOs
835 Complete and single-copy BUSCOs
500 Complete and duplicated BUSCOs
17 Fragmented BUSCOs
88 Missing BUSCOs
1440 Total BUSCO groups searched
core geneは60減って1335となった。ただしduplicateが大きく減り、その分、single copyが大きく増えた。これをどう捉えるか...
追記
この記事を書いた時点ではしりませんでしたが、非モデル生物のRNA seq解析なら、CORESTという１遺伝子由来のcontigをまとめてクラスター状態でカウント可能にするツールが使えそうです。
contigをまとめることの有用性は2016年のアライナーの論文でも指摘されています。"5 Application of quasi-mapping for clustering de novo assemblies" の段落を読んでみてください。
引用
SPAdes: A New Genome Assembly Algorithm and Its Applications to Single-Cell Sequencing
Anton Bankevich, Sergey Nurk, Dmitry Antipov, Alexey A. Gurevich, Mikhail Dvorkin, Alexander S. Kulikov, Valery M. Lesin, Sergey I. Nikolenko, Son Pham, Andrey D. Prjibelski, Alexey V. Pyshkin, Alexander V. Sirotkin, Nikolay Vyahhi, Glenn Tesler, Max A. Alekseyev, and Pavel A. Pevzner.
Journal of Computational Biology. May 2012, 19(5): 455-477.
TopHat: discovering splice junctions with RNA-Seq
Cole Trapnell Lior Pachter Steven L. Salzberg
Bioinformatics (2009) 25 (9): 1105-1111.
Transcript assembly and quantification by RNA-Seq reveals unannotated transcripts and isoform switching during cell differentiation
Cole Trapnell, Brian A Williams, Geo Pertea, Ali Mortazavi, Gordon Kwan, Marijke J van Baren, Steven L Salzberg, Barbara J Wold & Lior Pachter Nature
Biotechnology 28, 511–515 (2010) doi:10.1038/nbt.1621
The Sequence Alignment/Map format and SAMtools.
Li H1, Handsaker B, Wysoker A, Fennell T, Ruan J, Homer N, Marth G, Abecasis G, Durbin R; 1000 Genome Project Data Processing Subgroup.
Bioinformatics. 2009 Aug 15;25(16):2078-9. doi: 10.1093/bioinformatics/btp352. Epub 2009 Jun 8.
CD-HIT Suite: a web server for clustering and comparing biological sequences.
Ying Huang, Beifang Niu, Ying Gao, Limin Fu and Weizhong Li.
Bioinformatics, 2010(26): 680-682.full text
Cd-hit: a fast program for clustering and comparing large sets of protein or nucleotide sequences.
Weizhong Li and Adam Godzik.
Bioinformatics, 2006(22): 1658-1659. full text
Tolerating some redundancy significantly speeds up clustering of large protein databases.
Weizhong Li, Lukasz Jaroszewski and Adam Godzik.
Bioinformatics, 2002(18): 77-82. full text
Clustering of highly homologous sequences to reduce the size of large protein databases.
Weizhong Li, Lukasz Jaroszewski and Adam Godzik.
Bioinformatics, 2001(17): 282-283. full text
