よく似たゲノム情報を使い不完全なゲノム情報しか持たない種のRNA seq解析の精度を上げる自動化されたツール Necklace

2018 10/31、11/2 タイトル、コード等修正、docker追加

2021 3/9 出力例追記

　シーケンシングされた種の数が増加しているが、ゲノムの大部分は不完全である。それらにはギャップが含まれていても、配置されていない領域が残っていてもよく、アノテーションが不十分な場合もある。ゲノムを持つ種のRNAシーケンシング（RNA-seq）はモデル生物と同じ手順に従う。つまりリードをゲノムにアライメントし、注釈付き遺伝子とオーバーラップするリードをカウントし、リードカウントに基づいて発現変動があるかをテストする[論文より ref. 1]。しかし、このアプローチは、多くの生物にとって重要なbiologyを見逃す可能性がある。リファレンス配列のギャップまたはアノテーションの欠落のために、遺伝子のセグメントが見逃される可能性がある。下流のdifferential expression 解析は、遺伝子数が過小評価されているため、統計分析の能力を低下させる可能性が高い。著者らは、遺伝子が複数のアセンブリscaffoldsにまたがる場合に、異なるアノテーションを持つ遺伝子としてアノテートされる例を観察した。しかし、これは遺伝子が単一scaffolds内に収まる場合でも起こり得る。最悪の場合、遺伝子全体が見逃される可能性がある。

　RNA-seq解析は、発現された遺伝子に関する情報をデータ自体から抽出することで、利用可能な遺伝子モデルを修復するゲノムガイドの方法かde novo assemblyの方法を使うのが理想的である。しかし、アセンブリを含む解析は、複数のアセンブリを統合する必要がある場合に複雑になる（一部略）。

　Davidson et alの論文[ef. 4]（ツール紹介）で、著者らはsuperTranscriptomeの概念を導入した。そこでは、各遺伝子はすべてのエクソンを含む1つの配列によって表される。スーパートランスクリプトは、アセンブリやアノテーションの異なるソースからのトランスクリプトームをコンパクトで統一されたリファレンスに結合できる便利な手段を提供した。チキンRNA seqに適用すると、チキンリファレンスゲノムには存在しない何百もの遺伝子セグメントを回収できることが示された。

　ここでは、不完全なリファレンスゲノムを持つ種について、[ref.4]で説明したプロセスを自動化するNecklaceというソフトウェアを紹介する。 NecklaceはRNA-seqリードへのパス、リファレンスゲノム、および1つ以上のリファレンスゲノムアノテーションを含む設定ファイルを入力として受け取る。デノボアセンブリはエラーが起こりやすいため、関連する（よくアノテーションけされた）種のコード配列の中に新たにアセンブリされた遺伝子が発見される必要がある。したがって、関連する種のゲノムおよび注釈もNecklaceに提供されなければならない。Necklaceは、ゲノムガイドおよびデノボアセンブリに含まれるステップを実行し、アセンブリされたトランスクリプトームを、ユーザーが関心ある種のアノテーションと組み合わせる。 superTranscriptomeを構築した後、Necklaceは、リードをsuperTranscriptomeにアライメントしてリードをカウントし、edgeR [ref.5]、DEseq [ref.6]、またはDEXseq [ref.7]などのよく確立されたツールによる試験の準備を行う。

　新しいデータセットにNecklaceを適用することを実証するために、著者らはsheep milkの公開RNA-seqデータを分析した。ヒツジリファレンスゲノムのみを使用する場合と比較して、Necklaceを使った分析では、遺伝子にアサインされたリードが18％多く、発現変動遺伝子は19％多く検出された。

necklaceパイプライン。論文より転載。

Necklaceは、bpipeフレームワーク[ref.8]を使用して構築されたパイプラインである。 Necklaceはいくつかのサブステージで構成されている。初期ステージのゲノガイドアセンブリおよびde novoアセンブリ、遺伝子群への転写産物のクラスタリング、スーパートランスクリプトームの構築のための再構築、およびリードをマッピングしてカウントし発現変動遺伝子とスプライシングバリアントの分析に備える。アライナーやアセンブラなどの外部ソフトウェアと、C / C ++で書かれた独自のユーティリティセットを管理する。入力として、Necklaceはraw RNA-seqリードと該当する種のリファレンスゲノムおよびアノテーションファイルをとる。それに加えて、Necklaceは、よく研究されている種であるヒト、ショウジョウバエ、酵母などからリファレンスゲノムとアノテーションを取り込む。

リードカウントまでやってくれます。正規化は自分で実行する必要があります。

wiki

https://github.com/Oshlack/necklace/wiki

関連ツイート

インストール

ubuntu16.04を使ってテストした（docker使用。ホストOS: mac os10.12）。

依存

bash command line
Java 1.8 - You can also use earlier versions of Java, but need to perform de novo transcriptome assembly outside Necklace. See Options.
wget - Only required for installing Necklace and not for running it. If the required tools are installed manually, this is not needed.
g++

本体 Github

下記インストールスクリプトを実行するとcondaコマンドを途中実行する。他の依存とのコンフリクトが起きないようにも、仮想環境で実行するのが無難。javaの最新版は2018年現在version10だが、本ツールは8を要求する。以下の方法で1.8を導入した（確認は"java version"）

#JDK1.8導入
sudo apt-get install software-properties-common
sudo apt-get update
sudo add-apt-repository ppa:webupd8team/java
sudo apt-get update
sudo apt-get install oracle-java8-installer
java -version

#本体
git clone https://github.com/Oshlack/necklace.git
cd necklace/

#linuxの場合
./install_linux64.sh

#macの場合 (30分ほどかかる)
./install_mac.sh

#パスが表示されるかcheck
cat tools.groovy

$ cat tools.groovy

// Path to tools used by the pipeline

bpipe="/home/kazu/necklace/tools/bin/bpipe"

hisat2="/home/kazu/necklace/tools/bin/hisat2"

stringtie="/home/kazu/necklace/tools/bin/stringtie"

gffread="/home/kazu/necklace/tools/bin/gffread"

blat="/home/kazu/necklace/tools/bin/blat"

cluster="/home/kazu/necklace/tools/bin/cluster"

python="/home/kazu/necklace/tools/bin/python"

lace="/home/kazu/necklace/tools/bin/lace"

gtf2flatgtf="/home/kazu/necklace/tools/bin/gtf2flatgtf"

samtools="/home/kazu/.pyenv/shims/samtools"

Trinity="/home/kazu/necklace/tools/bin/Trinity"

bowtie2="/home/kazu/necklace/tools/bin/bowtie2"

make_blocks="/home/kazu/necklace/tools/bin/make_blocks"

featureCounts="/home/kazu/necklace/tools/bin/featureCounts"

上では直しているが、samtoolsは導入に失敗していたので、別途導入して、tools.groovyのbowtie2="" 行にパスを記載した (*1)。

> tools/bin/bpipe -h

# tools/bin/bpipe

Bpipe Version 0.9.9.2 Built on Mon Jul 18 03:11:39 UTC 2016

usage: bpipe [run|test|debug|execute] [options] <pipeline> <in1> <in2>...

retry [test]

stop

history

log

jobs

checks

status

cleanup

query

preserve

diagram <pipeline> <in1> <in2>...

diagrameditor <pipeline> <in1> <in2>...

-b,--branch <arg> Comma separated list of branches to

limit execution to

-d,--dir <arg> output directory

-f,--filename <arg> output file name of report

-h,--help usage information

-l,--resource <resource=value> place limit on named resource

-L,--interval <arg> the default genomic interval to execute

pipeline for (samtools format)

-m,--memory <arg> maximum memory in MB, or specified as

<n>GB or <n>MB

-n,--threads <arg> maximum threads

-p,--param <param=value> defines a pipeline parameter, or file of

paramaters via @<file>

-r,--report generate an HTML report / documentation

for pipeline

-R,--report <arg> generate report using named template

-t,--test test mode

-u,--until <arg> run until stage given

-v,--verbose print internal logging to standard error

-y,--yes answer yes to any prompts or questions

#導入が少し手間だったので、dockerhubにイメージを上げておきます。mac osでも使えるはずです。イメージにはNecklaceの下流解析のツールは入れていません。下流解析もコンテナ中で行うなら、必要なものを入れてcommitし直して使って下さい、

docker pull kazumax/necklace

#ホストのカレントディレクトリとイメージの/dataをシェアして起動
docker run -itv $PWD:/data/ kazumax/necklace
cd /root/necklace
tools/bin/bpipe -h

上のdockerイメージを使ってもテストランは正常に終わりますが、ラン中の進捗メッセージが非常に多く、何か見過ごしがあるかもしれません。注意して使って下さい。

テストラン

本ツールはBpipeというワークフローマネジメントシステムを使って構築されている。詳細はbpipeのペーパーとdocument参照。ランは、bin/のbpipeを叩くことで実行する。このテストランにはmac pro 2012で1時間ほどかかる。

#テストデータのダウンロードと解凍（羊のRNA seq。自身のゲノム以外に、ヒトゲノムをrelated speciesのゲノムとして使っている）
wget https://github.com/Oshlack/necklace/wiki/asserts/sheep_small_demo_data.tar.gz
tar -zxvf sheep_small_demo_data.tar.gz

#ラン
tools/bin/bpipe run necklace.groovy data/data.txt

情報はconfigファイルに記載する。

> cat data/data.txt

$ cat data/data.txt

// sequencing data

reads_R1="data/SRR2932539_small_1.fastq.gz,data/SRR2932540_small_1.fastq.gz,data/SRR2932541_small_1.fastq.gz,data/SRR2932542_small_1.fastq.gz,data/SRR2932561_small_1.fastq.gz,data/SRR2932562_small_1.fastq.gz,data/SRR2932563_small_1.fastq.gz,data/SRR2932564_small_1.fastq.gz"

reads_R2="data/SRR2932539_small_2.fastq.gz,data/SRR2932540_small_2.fastq.gz,data/SRR2932541_small_2.fastq.gz,data/SRR2932542_small_2.fastq.gz,data/SRR2932561_small_2.fastq.gz,data/SRR2932562_small_2.fastq.gz,data/SRR2932563_small_2.fastq.gz,data/SRR2932564_small_2.fastq.gz"

//The genome and annotation

annotation="data/Ovis_aries.Oar_v3.1.90.chr14.gtf"

genome="data/Ovis_aries.Oar_v3.1.dna.toplevel.chr14.fa"

//The genome and annotation of a related species

annotation_related_species="data/Homo_sapiens.GRCh38.90.CDS.chr19.gtf"

genome_related_species="data/Homo_sapiens.GRCh38.dna.toplevel.chr19.fa"

例えばfastqディレクトリに２サンプルのペアエンドfastq、refディレクトリにリファレンスとアノテーション情報があるなら以下のように記載する。

#===========================================================

// sequencing data

reads_R1="fastq/sample1_1.fq,fastq/sample2_1.fq"

reads_R2="fastq/sample1_2.fq,fastq/sample2_2.fq"

//The genome and annotation
annotation="ref/Ovis_aries.Oar_v3.1.90.chr14.gtf"
genome="ref/Ovis_aries.Oar_v3.1.dna.toplevel.chr14.fa"

//The genome and annotation of a related species
annotation_related_species="ref/Homo_sapiens.GRCh38.90.CDS.chr19.gtf"
genome_related_species="ref/Homo_sapiens.GRCh38.dna.toplevel.chr19.fa"

#===========================================================

ジョブの最後には以下のように表示される。

|| Running time : 0.00 minutes ||

|| ||

|| Read assignment finished. ||

|| ||

|| Summary of counting results can be found in file "counts/block.counts" ||

|| ||

===================== http://subread.sourceforge.net/ ======================//

======================================= Stage get_gene_info ========================================

======================================== Stage get_bp_info =========================================

====================================== Stage get_mapping_info ======================================