Docker で Trinity を使って de novo assembly する

次世代シークエンスでは、100塩基くらいの短い配列が、数百万から数千万個得られます。しかし、この配列から遺伝子の発現量や変異といったデータを得るには、コンピュータを使った解析が必要になります。

本記事では、非モデル生物を想定して、次世代シークエンスで得られた cDNA の断片配列から、解析に用いるためのリファレンスとなる遺伝子の塩基配列を作成する方法を解説します。

なぜアセンブリが必要なのか
fastq ファイルの前処理
Trinity でアセンブリする

なぜアセンブリが必要なのか

次世代シークエンスは、生物から抽出した DNA を断片化し、断片の配列を並行して大量に決定する技術です。これを RNA から逆転写した cDNA から行うことで、特定の条件下で発現していた RNA 由来の配列情報を大量に入手することができます。

そして、それらの断片配列を元のゲノムや転写産物と比較して、その断片がどの遺伝子に由来するものかを調べることができれば、その遺伝子の発現量を計算したり、断片配列に含まれる塩基配列の変異から、その遺伝子にどれくらい SNP が含まれるかを計算したりすることができます。
このように、遺伝子の配列をリファレンスに対応付けることを、マッピングといいます。次世代シーケンスでは、複数の条件でサンプリングを行い、それぞれの条件で得られた断片配列をリファレンスにマッピングし、発現量などの違いを調べるのが基本的な流れになります。

しかし、マッピングには、リファレンスとなるゲノムや転写産物の配列が必要です。モデル生物であれば、公開されているゲノム配列を使用すれば問題なく解析に進めますが、非モデル生物では、遺伝子の情報がほとんどわかっていないこともあると思います。

そこで、次世代シークエンスで得られた断片から、元の転写産物の配列を組み立てることにします。この作業をアセンブリといい、非モデル生物のリファレンスを作成するために必要な作業となります。

fastq ファイルの前処理

それでは早速アセンブリを行います...と言いたいところですが、先にデータの前処理を行います。次世代シークエンスは、同時に大量の塩基配列を決定しますが、その中には、配列決定の精度が低い部分や、シークエンスのために必要なアダプターの配列がサンプルと一緒に読まれている部分が含まれている場合があります。それらの配列が含まれているままアセンブリを行ってしまうと、正確ではない遺伝子がリファレンスに含まれることになり、後の解析に悪影響を与える可能性があります。そのため、データの前処理によって、アダプターや低品質な部分を取り除いてからアセンブリを行うことにします。
このように、リードの品質を確かめる作業をクオリティチェック、低品質なリードを取り除く作業をクオリティコントロールと言います (どちらも QC と呼ばれていて、ちょっとわかりにくい...)。

QC にはいろいろなツールがありますが、今回は fastp を使用します。fastp はクオリティチェックとクオリティコントロールを同時に行うことができ、処理がかなり高速なのでおすすめのツールです。Anaconda (Miniforge) で簡単に導入できます。
(fastp について詳しく知りたい方は、Github や論文をご確認ください。)

それでは、以下のコマンドでインストールしましょう。Miniforge のインストールについては、こちらで解説しています。

# 仮想環境作成 (前回お試しで作ったのとは違う環境を使用しています。今後もこちらの環境を使用します)
mamba create -n bioinfo -y

# 仮想環境に入る
mamba activate bioinfo

# インストール
mamba install bioconda::fastp -y

インストールできたら、fastp を実行します。私の場合は、R1.fastq と R2.fastq のペアから QC を行い、fastp_1.fastq と fastp_2.fastq を作成しています。また、「-w」は使用する CPU の数なので、パソコンに合わせて数値を変えてください (確か最大16までだったはず)。

fastp -i R1.fastq -I R2.fastq -o fastp_1.fastq -O fastp_2.fastq -h fastp_report.html -w 8

実行が終わったら、fastp_1.fastq, fastp_2.fastq, fastp_report.html, fastp.json というファイルができていると思います。
fastp.json に関しては、私はあまり見たことがないので省略します。
fastp_report.html を開くと、web ブラウザからレポートを見ることができ、実行前後のリード数や塩基数、品質、アダプターなどの情報が載っています。これを全てのサンプルに対して行ったら、前処理は終了です。

Trinity でアセンブリする

前処理ができたら、次は Trinity を使用してアセンブリを行います。アセンブリには大きく分けて2つの種類があり、近縁種のデータを参考にする方法を Reference-guided assembly、何もないところからリファレンスを組み立てる方法を de novo assembly といいます。実は Trinity はどちらもできるのですが、今回はより汎用性が高い de novo assembly を行います。以下に、ChatGPT による初心者向けの解説を書きました。

-- ChatGPT --
Trinity は、RNA-Seq データからトランスクリプト（転写産物）を再構築するための代表的なツールです。特に「リファレンスゲノムがない生物（de novo 解析）」でも使えるのが大きな特徴です。

・処理の流れ
1. Inchworm：最も高いカバレッジのシーケンスから "seed" を作り、リニアなコンティグを構築
2. Chrysalis：コンティグをクラスタリングし、リード情報で構造的な関係性を解析
3. Butterfly：アイソフォームや分岐構造を復元し、最終的なトランスクリプトを出力
(以下略)
----

詳しく知りたい方は、Github や論文をご確認ください。

それでは、Trinity を使用していきます。Trinity は依存ツールが多く、インストールに手間がかかるので、Docker を使用します (Anaconda でもインストールできるみたいですが、私の環境では実行の途中でエラーになりました)。
Docker のインストールについては、こちらで解説しています。
前回は自分で Dockerfile を書きましたが、今回は Docker Hub のイメージを使用する方法を解説します。以下のコマンドを実行してください (Docker は Miniforge とは独立しているので、仮想環境はそのままでも大丈夫です)。

docker pull trinityrnaseq/trinityrnaseq

しばらく時間がかかりますが、エラーが出なければ、Trinity が使えるようになっています。Docker の Trinity を使用するには、以下のようにコマンドを実行してください。
データは /data/file_name_1.fastq のように入力して、メモリと CPU はパソコンに合わせて調整してください (メモリが少ない場合、CPU も少なめにしておいた方が止まりにくいです)。

docker run --rm -it -v $(pwd):/data trinityrnaseq/trinityrnaseq \
Trinity --seqType fq --left /data/fastp_1.fastq --right /data/fastp_2.fastq \
--max_memory 10G --CPU 8 --output /data/trinity_out

終わったら、trinity_out, trinity_out.Trinity.fasta, trinity_out.Trinity.fasta.gene_trans_map というファイルができているはずです。
trinity_out.Trinity.fasta がアセンブリした転写産物の配列になっています。

これをリファレンスに使用して、発現量解析などを行うことができます。
今回は以上です。

おまけ: 書籍紹介

fastq, fasta といった用語や、fastq の品質の確認方法などの勉強に役立つ本をいくつか紹介します。
リンクをクリックすると、Amazon のページに飛びます。

・バイオインフォマティクス入門
生物学、コンピュータ、バイオインフォマティクスについて、基礎的な知識を習得するのに役立つ本です。見開きのページごとで内容がまとめられているので、毎日少しずつ勉強したいという方にもおすすめです。

・RNA-Seqデータ解析 WETラボのための超鉄板レシピ
具体的な解析方法が載っている本です。