macでインフォマティクス

macでインフォマティクス

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

de novo transcriptomeのcontigクラスタリングツール Corset

 

 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回実行で行えるため、使うのが簡単である。

f:id:kazumaxneo:20180614221350j:plain

非モデル生物でリードカウントベースの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 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データを使う。

 

1、ダウンロード

ダウンロードして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

目に見える問題は起きてない。

 

2、クオリティトリミング

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はアンペア。アンペアになったリードは使わない。

 

3、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をインストールする(なるべく手を抜く)。

 

4、マッピング

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を使う例は公式を参考にしてください。実際に動かして確認しましたが、流れは同じです。

 

5、クラスタリング

Corsetにbamの付属情報を与える。今回の6サンプルは、SRR453566-68の3つがbatch growth、SRR453569-71の3つがchmostat 条件になっている。条件が同じか違うかどうかは-gフラグで与える。次のコマンドでは、12は違うグループと見なされる。

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、2列目がアサインされた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で、それ以降の列がそのクラスターにアサインされたリードのカウント。

 

6、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

 

 

7、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が出力される。

f:id:kazumaxneo:20180616093103j:plain

 

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