2020-12-09

De novo transcriptomeのアセンブリ評価やフィルタリングのツール

2022/12/29　コマンド修正

2024/01/16 誤字修正

人の文章ばかり載せているので、たまには自分の言葉で書いてみる。

　ハイスループットシークエンシングのコストは年々下がっているが、それでも高品質なゲノム配列の解読と遺伝子予測、機能的注釈を得るまでには多くの人、時間、そして資金が必要になる。このゲノム解析の難易度の高さが、未知の生物の分子レベルの研究を行う上での大きな障害となっている。De novo transcriptome reconstructionとは、ゲノム配列情報なしで、De novo transcriptome assemblerによってRNA seqのシークエンシングリードから新規にRNA配列を再構成する方法である。得られたtranscriptome情報を使ってどのようなRNA分子が発現しているのかを調べれば、ゲノム配列がない生物でも生物学的な知見を得ることができるというわけである。この性質から、De novo transcriptome解析は、ゲノムが分からない生物でパイロット的に行われることが多い（実際、De novo transcriptome の論文ではタイトルやAbstractにファーストトライアル~かそれに類する表現が使われることが多い）。このような書き方をすると、transcriptomeの単なるスナップショットであり、ゲノム解析をするほどリソースが潤沢ではない非モデル生物の研究者の武器のように聞こえるかもしれないが、それは早計で、ゲノムが利用できない状態で素早く細胞のステータスを知るため（ゲノム解析と並行して）De novo transcriptome 解析を行う事もあるし、高品質なゲノムが利用できる生物においても、アノテーションがない新規の転写産物を探す目的で行うこともある。つまり、生物を問わず有用な研究アプローチである。有効性が高いため、このブログでも、De novo transcriptome 関連のツールは度々紹介してきた。

　De novo transcriptomeベースの発現解析では、遺伝子発現解析やGOエンリッチメント解析、パスウェイエンリッチメント解析、共発現解析などの解析手法の精度は、アドホックに調整されるトランスクリプトーム配列に完全に依存している。しかし、De novo transcriptome assemblerはそれほど完全なRNA配列の再構成手法ではなく、通常完全なトランスクリプトームは再構成できない。これは、入力するデータの性質を考えれば理解できる。まず、RNA seqのデータは特定の条件のRNAのスナップショットであり、各々の細胞内で発現しているRNA分子のコピー数は数桁のオーダーで幅がある。つまり、どのような条件でRNAを抽出しても大量に発現しているRNA分子と中程度に発現しているRNA分子、ほとんど発現していないRNA分子があり、低発現のRNA配列は、シークエンシングがほとんど行われないため、再構成する事が困難になる。これはgenomic DNAのシークエンシングとは明確に異なる。また、RNAは常に新陳代謝しており、合成途中のRNAと分解途中のRNA分子がある。合成や分解が活発に行われているRNA分子がシークエンシングされると、カバレッジは均等にならずRNA配列の一部しか再構成できない恐れがある（成熟前のRNAにはpolyAがないはずので、polyAから逆転写してライブラリ作成をしているなら合成中のRNAは読まれないはず）。これは短いプローブのみで定量していたマイクロアレイ時代には問題になりにくかった話だが、RNA seqが全ての配列をシークエンシングするため、RNA分子が部分的に分解された場合、対応するリード量の歪みが起こってアセンブルに影響を与えるということである。さらに、転写にはバリエーションが存在し、１遺伝子内に複数のRNAアイソフォームが存在することが多い。オルタナティブスプライシングによって特定のエキソンだけ異なり、あとは同じエキソンからなるアイソフォームのRNA配列が生じている場合、スプライシングバリアントを区別したアセンブルは難しく、特にショートリードシークエンシングのデータでは難易度が高い。さらに、ライブラリ調整法にもよるが、真の末端はタギングできなかったりする。このような解析上の困難の積み重ねで、多くのRNA配列は５’末端から３’末端までの全長は再構成されず、部分的な配列のみが再構成される。この影響は軽微ではなく、ゲノム配列にマッピングして進める発現解析と比較して、De novo transcriptome からの発現解析は精度面で劣る傾向がある（目当てのRNA配列が２つに分かれて再構成され、その各々が別々の配列としてアノテーション・定量されることを想像してみて下さい）。

　より完全なDe novo transcriptome reconstructionを目指して、これまでに様々なDe novo transcriptome assemblerが開発されてきたが、現状では、De novo transcriptome assemblerによってそれぞれ部分的に異なる配列が再構成される。得られた配列セットからどれがより適切なトランスクリプトーム配列であるのか簡単に判断する方法もない。そのため、アセンブラが出したトランスクリプトーム配列をそのまま使って解析を進めると結果に思わぬバイアスが存在する可能性がある。このような事態を引き起こさないためにも、複数の評価手法やフィルタリングを使い分けて、妥当な配列セットを得る事が欠かせない。生物学的に重要な洞察を得るためには、サンプル調整やシークエンシング（ウエット）、発現解析（ドライ）に慎重を期するだけではなく、トランスクリプトーム配列の調整にも細心の注意を払う必要があるということになる。

　ここでは、常に情報が不足しがちな非モデル生物のDe novo transcriptome解析に焦点を当て、アセンブルした配列の品質管理とフィルタリング方法を中心に紹介していく。目的は、データ解析の落とし穴にはまらないための知恵を蓄え、共有することにある。間違いやアドバイスがあったら教えていただきたい。

ポイント

ゲノムアセンブラと違い、De novo transcriptome のアセンブラは低発現の転写物から高発現の転写物まで組み立てることができるように設計されている。また、strand specificにRNA seqが行われていれば、ストランド特異的にアセンブルできるのも特徴である（アセンブル時はパラメータに注意）。DNAと違って、センス側のRNA配列（タンパク質コーディング）とアンチセンスのRNA（ノンコーディング）の両方が転写されている場合があり、それを区別するためである。
De novo transcriptome assemblerは、低発現のRNAから高発現のRNAまで満遍なくアセンブルする必要があるため、アセンブル過程で多様な戦略を組み込んでいる。しかし、どのようなヒューリスティックであっても、発現量が低いRNAのアセンブルは難しく、本来細胞内には存在しない疑わしい配列がアセンブルされる可能性がある。このような配列がどれくらい存在するのか調べるために、客観的な評価手法を使って精査しながら進める必要がある。
特定の遺伝子に興味がある場合、de novo transcriptome assemblyによってアイソフォームレベルの配列再構成ができることはメリットが大きい。例えば、何らかの刺激によってスプライシングパターンに変化が起きているのかかどうか調べることができる。
現在のGO enrichmernt解析などは遺伝子単位でしか注釈が存在しないため、アイソフォームレベルの配列再構成ができても、下流の発現解析では扱いづらくなる(羊土社のRNA-Seq レシピ本p.97も参照)。対策として、メジャーなアイソフォームだけ選抜して冗長性を減らす方法が存在する（後述）。
ペアエンド情報を使いより長い配列を構築できるため、De novo transcriptome ではできる限りペアエンドシークエンシングが望ましい。配列長も75x2より150x2の方がベターと考えられる（比較した論文は探してませんので妄想かもしれない）。
rRNA除去を行なっていても、polyAから逆転写してライブラリを作成していても、rRNAのリードは一定数読んでしまう。サンプルによってはかなり読んでしまう可能性もある。アセンブル前に予めSortMeRNAなどを使って除いておくと、アセンブルの計算負荷を抑制できる。一方でrRNAは分類情報として活用できるため、サンプルQCの一環としてphyloflashでアセンブルし、SILVAで調べる手もある。

アダプター除去などの前処理、発現解析などのコアの解析部分については、既に多くの情報があるため言及しない。De novo transcriptome assemblerも既にたくさん紹介しているので、重複しての紹介は避ける。

De novo transcriptome assemblerのランについてはこちらを参照。

汚染配列

De novo transcriptomeはまだ純粋な培養や育種が確立していない系にも利用することが多く、シークエンシングされたリードの中には、潜在的な汚染配列や共生生物のリードが混ざっている可能性がある。異なる分類群の汚染（共生）があるのか調べるには、blobtoolsが有用である。blobtoolsは、uniprotのプロテオームなどをデータベースとして、入力配列（クエリ）をタンパク質レベルで相同性検索し、ベストマッチのtaxaからクエリに分類をアサインする方法である。各配列をカバレッジ（発現量）を縦軸にlogでとり、GCを横軸にとってプロットし、その時にtaxaをプロットの色で表現したグラフを出力する。これをblobplotと呼び、汚染があるのか、それはどの程度なのかを視覚化して調べることができる。

　具体的には、Trinityなどのアセンブル配列をクエリとしてblobtoolsで分析する。blobtools plotで -r superkingdomを指定すれば、Eukaryota, Bacteria, Archaea, Virusの分類でプロットして、明らかな汚染（または共生生物）があるのかを定量的に判別できる。プロットする順番は変更できるので、otherのプロットやno-hit（ノンコーディング RNAやゲノムの汚染の可能性）が邪魔なら後ろに来るように変更する。データベース依存の方法なので、誤った分類群へのアサインも起きる可能性があることは知っておく。

　生のショットガンリードを直接分析するkrakenやcentrifugeなどの実装も存在する。これらはk-merの完全マッチに基づくアラインメントフリーの方法で（kaijuのようなタンパク質のアミノ酸の完全一致の実装もあり）、リソースが豊富な計算環境では非常に高速に汚染の推定ができるが、間違ったアサインも散見される。これは完全マッチなので複数の生物にアサインされることを防げないためである；例えばA genomeのk mer hash tableとマッチする場合、A と近いB genome のk mer hash tableにも確率的にマッチする事が高い。その場合、AとBが両方存在しているとレポートされる（どちらにより近いのか決定できない）。この延長で、多くの分類群に存在する配列があった場合、全てに”ベストマッチ”する（言い換えれば２値分類、二項分類の方法を分類群の推定に使っている）。その結果、krona等で可視化すると、多くの分類群のわずかな汚染があると誤解を与えてしまう。さらに公共のデータベースには種レベルでは間違った分類が含まれることも問題を複雑にしている。このため、De novo transcriptomeの解析には向いていない。使うにしても、判定は分類階級のある程度上のランクに留めた方が良いと思う。一方で、アセンブル配列をクエリにするblobtoolsのランでは、通常はdiamond blastxを使ってタンパク質レベルでベストマッチの配列を探す。結果は一意なヒットにはならないが、基本的にタンパク質相同性は近い分類群の方がより高くなるため、ベストマッチは一意となる（データベースのプロテオームに冗長さがない限り）。タンパク質レベルの検索なので、塩基配列の完全マッチとは異なり、データベースから距離がある非モデル生物でも機能する。

（下で紹介しているTRAPID 2.0も参照して下さい）

アセンブルの完全性

１、転写物の基本要約統計

上にも書いた通り、現状では用いるアセンブラによって異なる配列が出力される。それらのうちどれが不完全なトランスクリプトーム配列で、どれがより完全なトランスクリプトーム配列なのかを調べる必要がある。まず調べるべきは、どのくらいの数の転写物配列が出力できていて、N50や平均長はどれくらいである。この目的のために既存のツール（実装）が複数利用できる。

#例えば高速なseqkitを使う
seqkit stats -T -a assemlyA.fasta assemblyB.fasta > assembly_stats

#TransRate
transrate --assembly transcripts.fasta --threads 12

#Trinityのアセンブリならスクリプトが使える。一番長いisoformの長さが分かる。
TrinityStats.pl Trinity.fasta

#（Trinityのアセンブリ）trinity geneそれぞれのisoform数をカウント
 grep '^>' Trinity.fasta  | cut -f1,2,3,4 -d"_" |sort | uniq -c > count

例えば転写物の平均長が100bp程度しかなかった場合、平均33アミノ酸程度のタンパク質しかコードしていないことになる。一般的なBulk RNA seqなら、de novoアセンブルがうまくいっていないと判断できる。比較のために近いモデル生物のcDNAも一緒に分析する。

２、マッピング率

アセンブル配列数や長さを調べただけでは本当に失敗している場合を除いてあまり有用な情報を得られないことが多い。そこで次に調べるべきは、アセンブルに使ったリードのうちのどれくらいがアセンブルした配列にマッピングされるかである（バックマッピング）。マッピング率が低い場合、部分的な配列のアセンブルしか上手くいっていない可能性がある。アセンブルした配列にイントロンはないので、マッピング時にスプリットアラインメントは考慮しなくてOK。

#例えばbowtie2のsensitive setting。
bowtie2-build -f assemlyA.fasta bowtie2_index
bowtie2 -p 20 --sensitive -x bowtie2_index -1 pair_1.fq.gz -2 pair_2.fq.gz |samtools sort -@ 8 -> assemlyA.bam
samtools index -@ 4 assemlyA.bam

#statistics
bamtools stats -in assemlyA.bam
bamtools stats -in assemlyB.bam

#samtools flagsatsならsecondary alignmentも見れる
samtools flagsts -@ 4 assemlyA.bam
samtools flagsts -@ 4 assemlyB.bam

The Tuxedo protocolだったかどうか忘れてしまったが、80%以上のマッピングがボーダーとして論文に記載されていたと記憶している。だがこれは低い閾値で、アセンブルした全配列にマッピングした場合、98%くらいのマッピング率になることもある。この簡単な評価方法は通常のRNA seq でも利用価値がある。イルミナのハイクオリティなRNA seqでマッピング率が低くなるなら、原因を探るべきである。

(モデル生物のハイクオリティなアノテーション情報を使って既知転写領域にマッピングすると、90％くらいのマッピング率になることもあったりする。これは、転写領域としてまだアノテーションに含まれていない領域がたくさんあるからと推測される。この仮定が本当かどうかは、ゲノムにそのままマッピングしてマッピング率を出せば分かる)。

３、BUSCO

　シングルかつ90%以上の種のゲノムに見つかる遺伝子（2倍範囲内のサイズ　サイズ分布の2 S.D範囲内）をシングルコピーオルソロガス遺伝子として定義し、BUSCOの研究チームはこれを様々な分類群と分類階級で行ない、データベースとして公開している。BUSCOは、このデータベースを使ったゲノム、RNA、タンパク質の評価ツールである。シングルコピーオルソロガス遺伝子がシングルコピーで完全長見つかるのかどうか、もしくは複数見つかるのか調べることで、DNA/RNA アセンブリやプロテオームの完全性、冗長性を評価できる。様々な分類群のゲノムで計算されており、真核生物全部のような大きなランクだけでなく、よく研究された分類群については、属（ジーナス）のような下のランクの階級でも利用できる場合がある。真核生物全部になるとリボソーマル関連タンパク質のような配列しかなくなりシングルコピーオルソロガス遺伝子は数百まで減るが、狭い分類クレードでは保存されている遺伝子が多いため、BUSCOのシングルコピーオルソロガス遺伝子セットは数千まで増え、解析精度が高くなる。よって、通常は狭い分類で計算されたデータセットを使う。とは言え、真核生物全部のような大きなランクも使えないわけではなく、ほとんどゲノムがわかっていない分類のアセンブルでは大きなランクを使わざるを得ないし、全ゲノム重複などのイベントを経験している分類群では、その上の方の分類群で定義されたBUSCOをあえて使うことで、かつてそのようなイベントが起きたかどうか推測できる可能性がある（duplicated BUSCOを見る）。BUSCOはゲノムで定義されているため、transcriptomeアセンブルに使うと、スプライシングバリアントの存在でduplicated BUSCOが増える可能性がある点に注意する：実際モデル生物のcDNAにはたくさんのアイソフォームも含まれるため、分析するとduplicated BUSCOが多く出てくる。

BUSCOは、v4までAugustusを使用して遺伝子予測を行なっていた。しかしこのツールは割とトラブルを起こしがちで、かつゲノムサイズが大きいとランタイムが長引く原因となっていた。そのためか、現在ベータのv5ではAugustusがデフォルト設定ではなくなっている（オプションで使用可能）。transcriptomeやproteomeではこの遺伝子予測が不要のため、ランタイムは大幅に短くなる。
論文では結果の要約、例えばC:97.2%[S:23.1%,D:74.1%],F:1.6%,M:1.2%,n:255だけが記載されることもある。
BUSCOで定義された遺伝子は種間での保存性が高いだけで、必ずしも構成的に発現しているとは限らない。
BUSCOの出力にはsummary.tsvが出力される、どのようなタンパク質配列がヒットして、どのようなタンパク質配列がヒットしなかったのか探ることができる。

f:id:kazumaxneo:20201210230533p:plain

４、Transrate

リファレンスフリーのトランスクリプトームアセンブリ評価ツールとしてRSEM-evalがある(Li et al）。RSEM-evalは、リードデータを与えられたアセンブリ尤度を提供し、同じ入力データから生成されたアセンブリの比較を可能にする。TransRateはRSEM-evalから発展して、TransRateコンティグスコアとTransRateアセンブリスコアという新しいリファレンスフリーの統計を用いて、これらの特徴を評価する。

繰り返しになるが、de novo のtranscriptomeアセンブリではよく似た転写物のアセンブリが難しいため、スプライシングバリアントや複数の遺伝子ファミリーのアセンブリが成功しない恐れがある。厄介なケースだと、ハイブリッドな配列（キメラ）が作成されてしまう。このエラーは、転写物配列内のリードのマッピング状況を見ることで調べることができる。例えばよく似たファミリー遺伝子が同時に発現している場合でも、よく似たファミリー遺伝子も組織や時間によって発現量は大きく異なる場合が多くあり、特定のスナップショットのRNAseqリードをハイブリッドな転写物配列にマッピングすると、コンティグに沿ってリードカバレッジに変化が生じる点が観察される可能性がある。リードがコンティグの端から外にずれている場合、配列が不完全であることが示唆される。ペアエンドのリードが2つの転写物にブリッジしてアラインメントされる場合、転写物が２つに分かれていることが示唆される。冗長なマッピングは、1つの転写物が複数のコンティグによって繰り返し表現されている場合に起きる。またユニークにマッピングされたリードにhomoのindel変異があれば、アセンブルにエラーがある証拠となる。このような統計からアセンブルの精度を評価する。

５、TIN value

2016年、Liguo Wangらによって、RNA分解を測定するためのトランスクリプトインテグリティナンバー（TIN）に関する論文が発表された。RNAの分解速度は転写産物に特異的であり、それはいくつかの内因性および外因性因子のほか、AU-rich配列、転写産物の長さ、GC含有量、二次構造、RNA-タンパク質複合体などの他の因子によっても調節される。RNAの分解速度は、機能グループ間で、転写物間では最大で10倍の差があることがわかっている（論文より）。

https://bmcbioinformatics.biomedcentral.com/articles/10.1186/s12859-016-0922-z

TIN値を求める機能はRSeQCにも実装されている。

RSeQCにはstranded RNA seqかどうかを判定するinfer_experiment.py等のサポートスクリプトもある。

６、DETONATEスコアを使用する。これはトランスクリプトームアセンブリの品質を厳密に計算し評価するもので、異なるパラメータ設定や全く異なるツールを使用して複数のアセンブリを実行する場合に有用。DETONATEスコアが最も高いアセンブリが最良とみなせる。

７、（Trinityのアセンブリ向け）blastヒットから完全長転写産物の組み立てに成功している配列の割合を推定する。これはTrinityのマニュアルで紹介されている方法であり、組み立てた転写産物をその近縁のトランスクリプトームと比較して、全長率を推定する。

wget ftp://ftp.uniprot.org/pub/databases/uniprot/current_release/knowledgebase/complete/uniprot_sprot.fasta.gz
gzip -dv uniprot_sprot.fasta.gz
makeblastdb -in uniprot_sprot.fasta -dbtype prot

#blastx
blastx -query Trinity.fasta -db uniprot_sprot.fasta -out blastx.outfmt6 -evalue 1e-20 -num_threads 20 -max_target_seqs 1 -outfmt 6

-max_target_seqs 1はDNAでは期待したように動作しないが、Trinityのマニュアルによると、タンパク質データベース検索では直感通りに動作するように見える（らしい）。

組み立てられた転写物をすべての既知のタンパク質に対して整列させ、その長さのX%以上にわたって整列するユニークなトップマッチタンパク質の数を調べる。そのためのスクリプト"analyze_blastPlus_topHit_coverage.pl"が付属している。

analyze_blastPlus_topHit_coverage.pl blastx.outfmt6 Trinity.fasta uniprot_sprot.fasta

blastx出力ファイルに、トップヒットの長さ、Trinity転写産物へのアライメントに含まれるヒットの長さのパーセント、そのデータベースエントリーのヘッダー記述を含むフィールドが追加される。さらにトップマッチのデータベースエントリーの長さカバー率の分布が提供される。Trinityのマニュアルでは、さらに１つのタンパク質配列に不連続に並んでアラインされる場合、それらをグループ化してカウントするスクリプトの使用方法について説明されている。

フィルタリング

１、Transrate

上で紹介したTransrateには、de novo transcriptomeアセンブリを評価し、品質の良い配列だけを別に出力する実装がある。

transrate --assembly transcripts.fasta --left pair_R1.fq.gz --right pair_R2.fq.g --threads 20

#dockerを使う（こちらがオススメ）
sudo docker run --rm -itv $PWD:/data/ -w /data arnaudmeng/transrate:1.0.3 \
 transrate --threads 20 --assembly /data/Trinity.fasta \
 --left /data/pair_R1.fastq.gz \
 --right /data/pair_R2.fastq.gz

ランが終わると、アセンブリ評価スコアとともにgood.transcripts.fastaとbad.transcripts.fastaが出力される。de novo transcriptomeの配列を使うと、数十%の配列がフィルタリング対象になることもある。Transrateのラン後にBUSCOやマッピング率などを改めて分析する。

２、Trinityのfilter_low_expr_transcripts.pl

De novo transcriptome assemblerは一定の冗長性を許与してアセンブルを行うため、出力にはほぼ同一の転写物が複数出力されている事がある。一般的にこれらを丸め込んだほうが下流解析の精度が上がる。これは、GO termなどの語彙は遺伝子レベルでしか付かないためである（アイソフォームが複数あると、過剰に表現された語彙が生じやすいということ）。

Trinityに付属するスクリプトを使うと、Trinityのアセンブルからメジャーなアイソフォームを選抜できる。ランするにはまずalign_and_estimate_abundance.plにより、RSEMによるリードカウントを実行する（RSEMはde novo transcriptomeの定量が可能な初の実装）。確率的リードカウント値の取得後、メジャーなアイソフォームを出力するスクリプト:filter_low_expr_transcripts.plを実行する。２つのスクリプトを実行するだけで低発現アイソフォームをフィルタリングできる。

align_and_estimate_abundance.pl --transcripts Trinity.fasta \
--seqType fq \
--left all-1.fq.gz --right all-2.fq.gz \
 --est_method RSEM \
 --aln_method bowtie2 \
 --trinity_mode \
 --prep_reference \
 --coordsort_bam \
 --thread_count 20 \
 --output_dir RSEM_outdir

filter_low_expr_transcripts.pl --transcripts Trinity.fasta \
 --highest_iso_only \
 --trinity_mode \
 --matrix RSEM.genes.results \
 > output.fasta

Trinityはこのようなヘルパースクリプトが充実しているのも素晴らしい。

３、CD-HIT-EST

CD-HIT-ESTは似た配列を丸めこむ目的でよく使われる。発現量を重みとして加えず長さだけで残す配列を選抜するため、先に1と2のランを行ってから実行した方が良いと思われる。

cd-hit-est -i transcripts.fasta -c 0.98 -T 20 -M 20000 -o output.fasta

先に１、２の処理を実行している場合は、通常CD-HIT-ESTでほとんど配列は減らない、または全く変化しない。この処理で配列数が大きく減るなら、何かのコマンドが上手くいっていない可能性がある。

リードカウント

妥当な配列が選抜できたら、リードをトランスクリプトーム配列にマッピングして定量する。カウントにはDe novo transcriptomeに対応したRSEMを使うのが望ましい。要は上のコマンドを、それぞれのサンプルのペアエンドfastqを使って再び走らせる。

フィルタリングして得られた転写物配列を指定する。

#sample1
align_and_estimate_abundance.pl --transcripts filtered.fasta \
--seqType fq --left sample1-R1.fq.gz --right sample1-R2.fq.gz \
 --est_method RSEM \
 --aln_method bowtie2 \
 --trinity_mode \
 --prep_reference \
 --coordsort_bam \
 --thread_count 20
 --output_dir RSEM_outdir_sample1

#sample2
align_and_estimate_abundance.pl --transcripts filtered.fasta \
--seqType fq --left sample1-R1.fq.gz --right sample2-R2.fq.gz \
--est_method RSEM \
--aln_method bowtie2 \
--trinity_mode \
--prep_reference \
--coordsort_bam \
--thread_count 20
--output_dir RSEM_outdir_sample2

--seqType <string> fq|fa
--est_method <string> alignment_based: RSEM. alignment_free: kallisto|salmon
--trinity_mode Setting this mode will automatically generate the gene_trans_map and use it.
--prep_reference prep reference (builds target index)
--coordsort_bam provide coord-sorted bam in addition to the default (unsorted) bam.
--thread_count number of threads to use (default = 4)

if alignment_based est_method:

--aln_method <string> bowtie|bowtie2 (note: RSEM requires either bowtie or bowtie2)

例えばsamtools idxstats ででカウントする場合、contigごとのカウントはidxstatsコマンドで取得できるが、確率的なカウントアサイン方法ではないので、よく似た配列のマルチマッピングリードのカウントでバイアスが強まると予想される。まずはRSEMを試すことを勧める。

アノテーション

TransDecoderを使ってコード領域を予測し、推定タンパク質コード配列を得る。

得られたprotein.fastaにeggNOG mapperやtrinotate等で機能的注釈をつける。コード領域が見つからなかった配列はnon-coding RNAの可能性がある。これらは、本当であれば別に機能的注釈をつけることが望ましいが、タンパク質コードの配列より難易度は高い。

追記

2020/12

TranSuiteというツールも発表されている（現在はまだコードが整備されていない）。

TranSuite: a software suite for accurate translation and characterization of transcripts | bioRxiv

その後、紹介しました。

https://kazumaxneo.hatenablog.com/entry/2021/11/21/202254

リフトオーバー

2020年12月にBioinfomaticsジャーナルにアクセプトされたLiftoffというツールは、ヒトのタンパク質をコードする遺伝子の98.3%以上を、配列同一性が98.2%のチンパンジーゲノムアセンブリにリフティングできたとの記載がある。Liftoffを使うことで近縁なモデル生物のアノテーションを高い精度で移せるかもしれない。LiftoffはPreprintの時に本ブログで紹介した。

DRAP（紹介）にかけて、メガトランスクリプトームを統合するという進め方。

データセットには断片的な配列や汚染された配列が含まれていることが多く、汚染配列があると解析が難しくなる。これらの課題を軽減するために、TRAPID 2.0が開発された。TRAPID 2.0は、アセンブルされたトランスクリプトームデータを高速かつ効率的に処理するためのウェブアプリケーションで、各トランスクリプトに構造情報、機能情報、分類学的情報からなる数層のアノテーションを提供する。

予測されたタンパク質の機能的アノテーション

Interproscanデータベースを使う事でGO termをアサインできる。ここでは、GOやKEGG、その他の情報を一度にアサインできるeggNOG mapperを勧める。

DeepNOGも興味深い。

2021 12/24 追記

TrinityとTrinotateには、Trinotate（Trinityの開発者のBrian Haasさんが開発）のアノテーション結果を読み込んでGO enrichment解析を行うスクリプトが整備されている（リンク）。

git clone https://github.com/Trinotate/Trinotate.git
perl Trinotate/util/Trinotate_get_feature_name_encoding_attributes.pl Trinotate_report.xls  > annot_feature_map.txt

また、TrinotateにはTrinotateWebが含まれており、トランスクリプトームアノテーションをナビゲートし、Trinity/RSEM/Bioconductor解析フレームワークを用いてトランスクリプト発現および差分発現を解析するためのローカル駆動のWebベースのグラフィカルインターフェースが提供されている。

GOがアサインされたら、DEGsのリストについてGO enrichment解析を行う。

GO enrichment解析後、GOterm<tab>p-valueの順にソートしたリストを作成する。

GO:0009409 4.20E-10
GO:0009266 1.10E-08
GO:0009415 1.10E-08
GO:0009414 4.70E-08
GO:0009628 8.80E-07
GO:0006950 8.00E-06
GO:0009631 2.50E-05

REVIGOで視覚化する。

STRINGdbにある生物、もしくはそれに非常に近い生物のDe novo transcriptome解析なら、ランクベースの機能的エンリッチメント解析をSTRING上で簡単に実行できる。ブログで簡単に説明している。

NCBIのblastnやblastpではTranscriptome Shotgun Assembly Sequence Database（TSA）を選ぶことができる。

f:id:kazumaxneo:20210113201618p:plain