2019 6/11 インストール追記
RNA-seqは、トランスクリプトームの様々な側面を研究するための強力な技術である。それは、遺伝子発見、選択的スプライシングイベントの検出、継時的発現分析、融合の検出、SNPおよび転写後エディティングなどの変異の同定を含む広範囲の用途を有する[ref.1]、[ref.2]。マイクロアレイのような従来技術に比べてRNA-seqの利点の1つは、データの生成と分析にリファレンスゲノムとアノテーションが必要ないため、非モデル生物のトランスクリプトーム全体の分析が可能になることである。リファレンスゲノムが利用できない場合、トランスクリプトームはRNA-seqのリードから直接de novo assemblyされる[ ref.3]。新たなトランスクリプトームアセンブリのためのいくつかのプログラムが存在する：Velvet [ref.6]およびAbyss [ref.7]ゲノムアセンブラを拡張するOases [ref.4]およびTrans-abyss [ref.5]ならびにTrinity [ref.8] ]。これらのプログラムは、コンティグ（contig）と呼ばれる転写配列へ数百万のショートリードをアセンブリすることができる。
RNA-seqの一般的かつ生物学的に重要な用途の1つは、2つ以上の条件の間で発現が異なる遺伝子を同定することである[ref.9]。しかし、デノボアセンブリされたトランスクリプトームに対するdifferential expression解析を行うことは、遺伝子あたり複数のコンティグが報告されているため困難である。トランスクリプトームアセンブラーは、同じ遺伝子のアイソフォームを区別し、それぞれを別々に出力するため、共通の配列を持つ複数のコンティグが生じる。さらに、しばしば、異なるアイソフォームを真に代表してはいないコンティグを作るが、これは、二倍体またはプールされた集団内のシークエンシングエラー、リピートのカバレッジ変動や変異などのアーティファクトから生じている。その結果、トランスクリプトームのアセンブラは、転写物の断片化されたバージョン、またはSNPかindelだけ異なるコンティグを繰り返し報告することがよくある。驚くべきことに、シークエンシングエラー、SNP、または選択的スプライシングを伴わないデータのアセンブリでさえ、遺伝子ごとに複数のコンティグを生成することがシミュレーションによって示されている[ref.10]。したがって、デノボアセンブリによって生成されるコンティグの数は、典型的には大きい。例えば、8,000万リードのアセンブリから数十万のコンティグが生じる[ref.11]。
de novoトランスクリプトームアセンブリによって生成される必然的に長いコンティグのリストは、differential expression解析のためにいくつかの問題を引き起こす：i）リードは重複配列に対して明確にアラインできないため、リードのアライメント起源の決定には誤りが多くなる。 ii）より多くのコンティグ間でリードが割り当てられなければならないので、コンティグあたりの平均カウントが減少しdifferential expressionの統計試験の能力が低下する、; iii）複数テストの調整がより厳しくなる。そしてiv）differential expressionのコンティグが同定されても、多数の遺伝子がリストに複数回出現するので、解釈が困難となる。Differential expression解析をコンティグよりも遺伝子を使い行えば、これらの困難は克服されるであろう。しかしながら、de novo アセンブリされたコンティグセットから遺伝子レベルの発現を評価するための手順は、これまでの文献では十分に検討されていない。
De novoアセンブリされたトランスクリプトームからdifferential expression genes (DEG) を同定するには、いくつかの段階がある（論文 図1（下の図））。最初にアセンブリされ、リードはコンティグへマッピングされる。コンティグはクラスター化され遺伝子にし、その後、各遺伝子の発現レベルを解析する（図1）。DEGを検出するために統計試験を行う。
いくつかの研究では、この分析パイプラインの個々のステップを比較している。例えば、様々なアセンブラの検討や、アセンブリ前の品質管理のメリットが検討されている[ref.12]〜[ref.15]。同様に、リードカウントベースの統計的な検定を行う方法も比較評価されている[ref.16]、[ref.17]。しかし、トランスクリプトームアセンブリから遺伝子レベルのカウントを得るための系を比較したり、提案したりする研究はほとんどない[ref.18]、[ref.19]。（一部略）
この論文では、任意のde novoトランスクリプトームアセンブリから遺伝子レベルのカウントを得るための方法とソフトウェアであるCorsetを紹介する。Corsetは、新しくアセンブリされたトランスクリプトーム配列にマルチマッピングされたリードのセットを受け取り、共通するリードの比率および発現パターンに基づいてコンティグを階層的にクラスタリングする。発現パターンが異なる場合、リードの配列がシェアされているパラログなどの遺伝子間を区別可能である。マップされたリードを使用して、Corsetは遺伝子レベルのカウントを出力する。次いで、遺伝子レベルのカウントは、edgeRおよびDESeqなどのカウントベースのフレームワークを使用して、differential expressionの試験を容易に実行できる。著者らは、Corsetが、複数の評価基準に従って、ほかのクラスタリング手法と比較して一貫して優れた性能を発揮することを示している。Corsetはアセンブラに依存しないメソッドであるため、さまざまなソースからコンティグや転写産物を組み合わせることができる。また、Corsetはクラスタリングとカウントのステップが1回実行で行えるため、使うのが簡単である。
非モデル生物でリードカウントベースのDEG検出を行う流れ。論文より転載。
簡単に言うと、同じ転写産物には、contigが複数に分かれていてもリードが冗長にマッピングされうることを利用してクラスタリングする手法です。複数mappingが全て報告されているbamを使います。
インストール
cent os6とmac os10.12でテストした。pythonのバージョンはは3.5.2 、Anaconda custom (64-bit)。
本体 Github
https://github.com/Oshlack/Corset
git clone https://github.com/Oshlack/Corset.git
cd Corset/
./configure
make
make install
#condaでも導入できる（リンク）。
conda install -c bioconda -y corset
> ./corset
$ ./corset
Running Corset Version 1.06
No input files specified
Corset clusters contigs and counts reads from de novo assembled transcriptomes.
Usage: corset [options] <input bam files>
Input bam files:
The input files should be multi-mapped bam files. They can be single, paired-end or mixed
and do not need to be indexed. A space separated list should be given.
e.g. corset sample1.bam sample2.bam sample3.bam
or just: corset sample*.bam
If you want to combine the results from different transcriptomes. i.e. the same reads have
been mapped twice or more, you can used a comma separated list like below:
corset sample1_Trinity.bam,sample1_Oases.bam sample2_Trinity.bam,sample2_Oases.bam ...
Options are:
-d <double list> A comma separated list of distance thresholds. The range must be
between 0 and 1. e.g -d 0.4,0.5. If more than one distance threshold
is supplied, the output filenames will be of the form:
counts-<threshold>.txt and clusters-<threshold>.txt
Default: 0.3
-D <double> The value used for thresholding the log likelihood ratio. The default
value will depend on the number of degrees of freedom (which is the
number of groups -1). By default D = 17.5 + 2.5 * ndf, which corresponds
approximately to a p-value threshold of 10^-5, when there are fewer than
10 groups.
-m <int> Filter out any transcripts with fewer than this many reads aligning.
Default: 10
-g <list> Specifies the grouping. i.e. which samples belong to which experimental
groups. The parameter must be a comma separated list (no spaces), with the
groupings given in the same order as the bam filename. For example:
-g Group1,Group1,Group2,Group2 etc. If this option is not used, each sample
is treated as an independent experimental group.
-p <string> Prefix for the output filenames. The output files will be of the form
<prefix>-counts.txt and <prefix>-clusters.txt. Default filenames are:
counts.txt and clusters.txt
-f <true/false> Specifies whether the output files should be overwritten if they already exist.
Default: false
-n <string list> Specifies the sample names to be used in the header of the output count file.
This should be a comma separated list without spaces.
e.g. -n Group1-ReplicateA,Group1-ReplicateB,Group2-ReplicateA etc.
Default: the input filenames will be used.
-r <true/true-stop/false>
Output a file summarising the read alignments. This may be used if you
would like to read the bam files and run the clustering in seperate runs
of corset. e.g. to read input bam files in parallel. The output will be the
bam filename appended with .corset-reads.
Default: false
-i <bam/corset> The input file type. Use -i corset, if you previously ran
corset with the -r option and would like to restart using those
read summary files. Running with -i corset will switch off the -r option.
Default: bam
-l <int> If running with -i corset, this will filter out a link between contigs
if the link is supported by less than this many reads. Default: 1 (no filtering)
-x <int> If running with -i corset, this option will filter out reads that align to more
than x contigs. Default: no filtering
Citation: Nadia M. Davidson and Alicia Oshlack, Corset: enabling differential gene expression
analysis for de novo assembled transcriptomes, Genome Biology 2014, 15:410
Please see https://github.com/Oshlack/Corset/wiki for more information
ラン
最低bamだけ指定すればランできる。
corset input.bam
cluster.txtが出力される。
ここでは公式の例に従いテストランを行う。データのダウンロードからアセンブリ、マッピング、クラスタリング、DEG検出まで行う本格的なものになっている。yeastのRNA seqデータを使う。
１、ダウンロード
ダウンロードしてFastqに変換。
cat >sra_run_ids.txt <<EOF
SRR453566
SRR453567
SRR453568
SRR453569
SRR453570
SRR453571
EOF
#ダウンロード
prefetch --option-file sra_run_ids.txt
mkdir sra_files
mv ~/ncbi/public/sra/SRR4535* sra_files/
forループでfastqに変換。
#fastqディレクトリを作成
mkdir fastq
for srr in $(cat sra_run_ids.txt)
do
fastq-dump --split-3 --skip-technical --readids
--defline-seq '@$sn/$ri' --defline-qual '+' sra_files/${srr}.sra
-O fastq/
done
- --split-3 Legacy 3-file splitting for mate-pairs: First biological reads satisfying dumping conditions are placed in files *_1.fastq and *_2.fastq If only one biological read is present it is placed in *.fastq Biological reads and above are ignored.
> ls -alh fastq/
$ ls -alh fastq/
total 16G
drwxr-xr-x 2 uesaka user 4.0K Jun 15 12:16 .
drwxr-xr-x 5 uesaka user 4.0K Jun 15 12:10 ..
-rw-r--r-- 1 uesaka user 1.3G Jun 15 12:13 SRR453566_1.fastq
-rw-r--r-- 1 uesaka user 1.3G Jun 15 12:13 SRR453566_2.fastq
-rw-r--r-- 1 uesaka user 1.7G Jun 15 12:14 SRR453567_1.fastq
-rw-r--r-- 1 uesaka user 1.7G Jun 15 12:14 SRR453567_2.fastq
-rw-r--r-- 1 uesaka user 1.3G Jun 15 12:15 SRR453568_1.fastq
-rw-r--r-- 1 uesaka user 1.3G Jun 15 12:15 SRR453568_2.fastq
-rw-r--r-- 1 uesaka user 916M Jun 15 12:15 SRR453569_1.fastq
-rw-r--r-- 1 uesaka user 916M Jun 15 12:15 SRR453569_2.fastq
-rw-r--r-- 1 uesaka user 1.5G Jun 15 12:16 SRR453570_1.fastq
-rw-r--r-- 1 uesaka user 1.5G Jun 15 12:16 SRR453570_2.fastq
-rw-r--r-- 1 uesaka user 1.4G Jun 15 12:17 SRR453571_1.fastq
-rw-r--r-- 1 uesaka user 1.4G Jun 15 12:17 SRR453571_2.fastq
目に見える問題は起きてない。
２、クオリティトリミング
trimmomatic（brewで導入可）を使っているが、ペアエンドを扱えるなら他のツールでもOK。
for srr in $(cat sra_run_ids.txt)
do
trimmomatic PE -phred33 fastq/${srr}_1.fastq fastq/${srr}_2.fastq fastq/${srr}_P1.fastq fastq/${srr}_U1.fastq
fastq/${srr}_P2.fastq fastq/${srr}_U2.fastq
LEADING:20 TRAILING:20 MINLEN:50
done
Pはペアエンド、Uはアンペア。アンペアになったリードは使わない。
３、de novoアセンブリ
全データをマージし、Trinityを使いde novoアセンブリを実行する。
#マージ。公式では、このstepでsedを使い、Trinityのランに必要な/1と/2を付加している。今回はダウンロード時に付加済み
cat fastq/*_P1.fastq > fastq/all_1.fastq
cat fastq/*_P2.fastq > fastq/all_2.fastq
#アセンブリを実行
Trinity --seqType fq --left fastq/all_1.fastq --right fastq/all_2.fastq --max_memory 100G --CPU 20 --full_cleanup
- --full_cleanup only retain the Trinity fasta file, rename as ${output_dir}.Trinity.fasta
- --CPU number of CPUs to use, default: 2
- --max_memory suggested max memory to use by Trinity where limiting can be enabled. (jellyfish, sorting, etc)
メモリは100GB指定している（CPU指定が増えると比例してメモリ使用量も増加）。カレントにtrinity_out_dir.Trinity.fastaが出力される。
参考: Trinityがエラーを起こし解決が難しければ、Docker環境で動かすのが楽かもしれません。Dockerならオーバーヘッドが比較的少なくランできます。例えばcentosのコンテナ入れて、そこにpyenv入れてcondaを導入し、conda installでTrinityをインストールする（なるべく手を抜く）。
４、マッピング
salmon（紹介）やbowtie2を使う。全てのアライメントをレポートするフラグを立てる必要がある。
#index
bowtie-build trinity_out_dir.Trinity.fasta Trinity
#mapping
for srr in $(cat sra_run_ids.txt)
do
bowtie -p 20 --all -S Trinity -1 fastq/${srr}_P1.fastq -2 fastq/${srr}_P2.fastq |samtools sort -@ 20 -O BAM -o fastq/${srr}.bam -
done
- --allをつけ、全てのアライメントを出力するようにしている。fastq/にbamが出力される。salmonを使う例は公式を参考にしてください。実際に動かして確認しましたが、流れは同じです。
５、クラスタリング
Corsetにbamの付属情報を与える。今回の６サンプルは、SRR453566-68の３つがbatch growth、SRR453569-71の３つがchmostat 条件になっている。条件が同じか違うかどうかは-gフラグで与える。次のコマンドでは、1と2は違うグループと見なされる。
cd fastq/
#サンプルごとにクラスタリング。6サンプルを並列実行。子プロセスが終わるのを待つwaitをつけている。
for FILE in `ls *.bam` ; do
corset -r true-stop $FILE &
done
wait
#統合
corset -g 1,1,1,2,2,2 -n A1,A2,A3,B1,B2,B3 -i corset *.corset-reads
clusters.txt、count.txtが出力される。
> head clusters.txt
$ head clusters.txt
TRINITY_DN1053_c0_g1_i1 Cluster-0.0
TRINITY_DN1053_c0_g1_i2 Cluster-0.0
TRINITY_DN2782_c0_g1_i1 Cluster-1.0
TRINITY_DN1713_c0_g9_i1 Cluster-2.0
TRINITY_DN1710_c0_g11_i2 Cluster-3.0
TRINITY_DN1836_c0_g1_i1 Cluster-4.0
TRINITY_DN180_c0_g1_i1 Cluster-5.0
TRINITY_DN947_c0_g2_i1 Cluster-6.0
TRINITY_DN865_c0_g1_i1 Cluster-7.0
1列目が転写産物のID、２列目がアサインされたcluster ID。
>head count.txt
$ head counts.txt
A1 A2 A3 B1 B2 B3
Cluster-0.0 361 541 397 361 458 562
Cluster-1.0 0 1 2 2 3 8
Cluster-2.0 0 3 3 3 4 4
Cluster-3.0 0 5 1 5 0 0
Cluster-4.0 10 5 4 9 15 25
Cluster-5.0 35 25 19 867 1663 1949
Cluster-6.0 7 6 2 3 3 10
Cluster-7.0 1650 2430 1856 608 702 829
Cluster-8.0 1041 1575 1214 734 849 1132
1列目がクラスターIDで、それ以降の列がそのクラスターにアサインされたリードのカウント。
６、edgeRを使ったDEG判定
> R #Rを立ち上げる
library(edgeR)
count_data<-read.delim("counts.txt",row.names=1)
group<-factor(c(1,1,1,2,2,2))
dge<-DGEList(counts=count_data,group=group)
dge<-DGEList(counts=count_data,group=group)
keep <- rowSums(cpm(dge)>1) >= 3
dge <- dge[keep,,keep.lib.sizes=FALSE]
dge<-calcNormFactors(dge)
dge<-estimateCommonDisp(dge)
dge<-estimateTagwiseDisp(dge)
results<-exactTest(dge)
topTags(results)
Comparison of groups: 2-1
logFC logCPM PValue FDR
Cluster-3549.0 7.781056 8.703321 0.000000e+00 0.000000e+00
Cluster-1221.0 6.381804 6.825066 0.000000e+00 0.000000e+00
Cluster-3976.0 4.793833 7.872592 0.000000e+00 0.000000e+00
Cluster-5221.0 4.487102 7.319362 7.386281e-321 9.764664e-318
Cluster-34.2 7.148952 10.719010 5.103542e-317 5.397506e-314
Cluster-2899.0 4.225118 7.564345 2.930408e-310 2.582666e-307
Cluster-1599.0 5.453314 6.612733 2.553171e-309 1.928738e-306
Cluster-4134.0 8.640848 8.879755 4.277488e-309 2.827419e-306
Cluster-502.0 8.604537 9.306601 2.350392e-296 1.380986e-293
Cluster-2771.0 5.183932 7.621110 2.961453e-291 1.566016e-288
続き。
allTags<-topTags(results,n=dim(results$table)[1])1
#ここではFDR0.05で切る。
deClusters=rownames(allTags)[allTags$FDR<0.05]
#保存
write.table(deClusters,"clusters_of_interest.txt",quote=F,row.names=F,col.names=F)
> head clusters_of_interest.txt
$ head clusters_of_interest.txt
Cluster-3549.0
Cluster-1221.0
Cluster-3976.0
Cluster-5221.0
Cluster-34.2
Cluster-2899.0
Cluster-1599.0
Cluster-4134.0
Cluster-502.0
Cluster-2771.0
７、FASTAへの変換
用意されているスクリプト（Corset-tools）を使えば、上のcluster IDをfastaに変換できる。
#Corset-toolsをクローンする。
git clone https://github.com/Adamtaranto/Corset-tools
python Corset-tools/fetchClusterSeqs.py -i trinity_out_dir.Trinity.fasta -t clusters_of_interest.txt -o contigs_of_interest.fasta -c clusters.txt
（python2.7環境で実行）
cluster IDに該当するMulti-FASTAが出力される。
8、アノテーション
Multi-FASTAが得られたら、blast2goやTrinnotate（紹介）を使ってアノテーションをかける。
以前一度紹介しましたが、中途半端だったのでまとめ直しました。
引用
Corset: enabling differential gene expression analysis for de novo assembled transcriptomes
Davidson NM, Oshlack A
Genome Biol. 2014 Jul 26;15(7):410.
corset-project - Example.wiki
https://code.google.com/archive/p/corset-project/wikis/Example.wiki
