macでインフォマティクス

macでインフォマティクス

NGS関連のインフォマティクス情報についてまとめています。

ゲノム情報はないが、モデル生物と近縁な生物のRNA seq 解析

 ゲノム情報はないが近縁種のゲノムが解読されているような生物でRNA seqを行うと決まったら、どんなワークフローで進めるべきだろうか?マイクロアレイと違い、RNA seqならde novoでも解析は不可能ではない。ゲノムがモデル生物とほぼ同じならば、深く考えなければ、なんだかできそうに思えてくる。

 実際、SRAのデータベースを見ているとそういったことを狙ったであろう解析データが無数に登録されている。その中にはpublish情報がないものもある。これが何を意味するかは置いておいて、モデル生物のゲノム情報の充実に伴い、そのゲノム情報を利用した近縁種の解析はこれからもますます増えてくると予想される。

ただしモデル生物に近縁でも塩基置換が起こってるだけでアライメントは難しくなる。ではde novo assemblyして進めるか?、というとせっかくゲノムが決まっているのだからなんとか利用したいという気持ちもあり、ジレンマに陥る。場合によりけりだと思うが、どういったワークフローが妥当なのか1つの例を通して考えてみることにする。

 

 今回はアラビドプシス属のA.kamchaticaのRNA seqデータを使う。シロイヌナズナとゲノムはかなり似ていると思われるが、ゲノムはまだ解読されていない。シロイヌナズナのecotypeでも塩基置換がかなり起きていることが知られており、A. kamchaticaのRNA seqデータがどれくらいナズナのリファレンスゲノムにアライメントされるかは未知数である。

 

今回はアセンブルしてリード数などの比較を行うが、それ以上は深入りせずプロトコルに抽象性を残しておく。その方が汎用性のあるプロトコルになると考える。

 

Arabidopsis. kamchatica (アラビドプシス属

f:id:kazumaxneo:20170730105553j:plain

 

 

データ

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%以上同一の重複配列を廃除する。

 

 

 

 

 

1、contigの作成。

genomeガイドの方法1つと、de novo assemblyの手法2つをを心みる。

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のサイズを計算、7カラム目に追加= > 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だった。

 

 

3、 de-novo assembly - trinity

Trinity --seqType fq --left R1.fastq --right R2.fastq --max_memory 40G --CPU 20

合計78657 contigできた。 そのうち 100bp以上のcontigは78656だった(1つだけ21 bp)。

アセンブルでエラーを起こした。フォーマットが古いせいだと思われたので、ヘッダーを修正して3行目を+に修正してランした。

 

 

 

2、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がある。

 

まとめ

f:id:kazumaxneo:20170804143209j:plain

cuffmergeはかなり厳しい。このケースではゲノムガイドアセンブリは現実的ではないようだ。 trinityがベストのように思える。

 

 

#BUSCOは以前紹介しました。簡単な使い方をまとめています。

アセンブル結果をCore gene setの検出数で評価する BUSCO - macでインフォマティクス

 

 

3、contigの長さの比較

 

fastaのサイズについてヒストグラムを書き比較する。

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) #グラフに追加 

f:id:kazumaxneo:20170727011204j:plain

 黒がTAIRのcDNAの分布で、ピンクがrnaspadesでアセンブルしたcontigの分布、緑がcuffmergeでゲノムベースにアセンブルしたcontigの分布。縦軸は頻度。横軸は長さ(bp)。

cuffmergeで作られた領域は短いものが非常に多い。 

 

次にtrinityとrnaspadesの比較。

hist(spades, col = "#ff00ff40", breaks = 150) 
hist(trinity, col = "#0000ff40", breaks = 150, add = TRUE) #グラフに追加

 

f:id:kazumaxneo:20170804135811j:plain

ピンクがrnaspadesでアセンブルしたcontigの分布、青がTrinityでセンブルしたcontigの分布。Trinityの方が長いcontigが多い。

 

 

 

 

4、contigのblast解析

3つの方法で作成した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

1、2列目がクエリーとデータベースの名前である。さらに3列目が相同性(%)、4列目がhitした長さ、5列目がhitした領域中のミスマッチの数、6列目がギャップの数である。12列目にスコア (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配列とシロイヌナズナの配列を1つにしてcd-hit-estで選抜してより長い配列を使うことを考えていた。しかしよく考えてみると、cd-hitで選抜すると、完全長のcDNAになっているシロイヌナズナの配列が優先されてde novoで作ったcDNAは捨てられ、結局アライメントでうまく行かなくなってしまう可能性が高い。またそれこそ恣意的なデータセットを使っていると取られかねない。

 

 ということで、みんながやっているようにdenovo assmeblyで作ったcontigのみ使うことにする。解析に供した全サンプルのリードを使ってアセンブルすれば、完全なcDNAコレクションにはなつていないとしても、少なくともそのデータに関しては定量漏れの遺伝子が出る可能性は低い。

 

たたし、denovo assmeblyのデータだけ使うなら、アセンブルの精度が実験の品質を完全に左右することになるので、アセンブル精度を客観的に評価できる手法でアセンブル精度を担保しておく必要が出てくる。

 

以下は最近アクセプトされた論文での流れ。

 2016年 nature scientific report de nono transcriptomeのペーパー

1、trinityを使い、k-merが25、32でアセンブリ

2、できたcontigをBUSCOとTransRateで評価。

3、アセンブリには適切なk-mer coverageがあるとの先行研究から、ランダムサンプリングでリード数を減らし再アセンブル

4、1のデータと3のデータをBUSCOとTransRateで再評価

5、ランダムサンプリングしてリード数を減らしk-mer=32でアセンブルのがベストだと判断。

6、Trinotateでアノテーション

7、 blastで植物ゲノムとベストヒットしないものを排除。

 

 このように、k-mer、カバレッジ、などを試行錯誤してアセンブルし、それから評価ツールで良好なアセンブル条件を比較する流れを取っている。 

 

 

 

今回のケースで、近縁種のゲノムが出ていたとしても、直接鋳型としてgenome guide assemblyをするのは厳しいことを身をもって実感した。ゲノムガイドとde novoのどちらも検討してもよいが、最終的には、近縁種のモデル生物のゲノム情報はアノテーション情報として利用するのが精一杯と思われる。

 

 

 

 

 

 

 

 

 

 

 

以下は、cd-hitを使ってシロイヌナズナのcDNAと統合しようと考えていた時の流れである。論文投稿のレベルを満たした手法とは思えないが、残しておく。 

 

5、マージ

ここから先は2つの手段を取る。

 

A、アセンブルして得られたcDNAからシロイヌナズナに見つからないcDNAだけ選り分け、それをTAIR cDNAとマージしてリファレンスとする。

B、アセンブルで得られたcDNAとTAIR cDNAをそのままマージし、CD-HIT-ESTで重複を除きリファレンスとする。

 

Aの方法を実行するには、シロイヌナズナのcDNAと相同性を示さないcDNAを集める必要がある。そこで、step4でblast解析した結果を相同性を示すcDNAのデータベースとして、no hitのcDNAをrnaspadesのアセンブルfastaから抽出する作業を行なった。そうゆうオプションやツールがありそうだが(絶対ある)、すぐには見つからなかったのでperlスクリプトを書いて処理した。

得られたno-hitのcDNA.fastaをblastにかける。organismはシロイヌナズナに限定し、最初の100 contigをblast解析する。e-valueの閾値はデフォルトでランした。

f:id:kazumaxneo:20170727144655j:plain

 

no-hitのcDNAだけをクエリにしているので、 スクリプトが上手くワークしているなら、シロイヌナズナのゲノムと全体にわたってマッチする配列は0のはずである。

 

結果は大まかに3パターンに分かれた。

 

highly homologous with A.thaliana genome

f:id:kazumaxneo:20170727150550j:plain

partially homologous with A.thaliana genome

f:id:kazumaxneo:20170727150611j:plain

no siginificant holology with A.thaliana genome

f:id:kazumaxneo:20170727150656j:plain

 予想通り多くのcDNAはナズナゲノムと高い相同性は示さなかった。このことは、Arabidopsis. kamchatica 固有の配列が抽出されていることを意味する。しかしながら、blastにかけた100 cDNAのうち1割ほどの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を1つにして"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は1つにまとめられる。

解析が終わると、クラスター化されて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に減った。

ヒストグラムを書く。 

f:id:kazumaxneo:20170729112752j:plain

ピンクが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が大きく増えた。これをどう捉えるか...

 

 

 

 

 

 

 

引用

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