de novo assemblyで得たRNAのコンティグをクラスター化して、非モデル生物のDEG解析を可能にする Corset

　RNA seqデータをde novoでアセンブルすると、一般に同じ遺伝子のアイソフォームが区別され、それぞれを別々にアセンブルするため、似た複数のコンティグが生じてしまう（ SNPまたはindelだけが異なるコンティグを繰り返し報告する）。付け加えて、こうしてできたコンティグが異なるアイソフォームを真に反映している保証もない。こういったそっくりなコンティグは、アイソフォーム、二倍体生物や集団内のシークエンシングエラー、反復、カバレッジの変化または遺伝的変異などに起因して生じている。よってデノボアセンブリによって生成されるコンティグの数は、本来のトランスクリプトームの数より大きくなる。このようなデータを使ってDEG解析を行うと、より多くのコンティグの間でリードを割り当てなければならないので、コンティグあたりの平均カウントが減少する。また、統計的に有意な変動を示すコンティグが同定されたとしても、実際は複数のコンティグがリストに存在するので、解釈が難しくなる。

　Corsetは、コンティグに複数回マップされた（一連の読み取りごとに複数のアライメントが報告される）リードを対象に、発現パターンに基づいてコンティグを階層的にクラスター化する。

Github

https://github.com/Oshlack/Corset

wiki

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

ダウンロード

linuxなら以下からビルド済みバイナリファイルをダウンロードできる。

https://github.com/Oshlack/Corset/releases

ラン

1、テストデータのダウンロード。

prefetch SRR453569 & #ダウンロード1
prefetch SRR453570 & #ダウンロード2
prefetch SRR453571 #ダウンロード3

#sraからpiared-end.fastqに変換。
mkdir fastq_db #ここに収納
fastq-dump ~/ncbi/public/sra/SRR453569.sra -O fastq_db/ --split-files &
fastq-dump ~/ncbi/public/sra/SRR453570.sra -O fastq_db/ --split-files &
fastq-dump ~/ncbi/public/sra/SRR453571.sra -O fastq_db/ --split-files
 
seqkit stats fastq_db/*

$ seqkit stats fastq_db/*

file format type num_seqs sum_len min_len avg_len max_len

fastq_db/SRR453569_1.fastq FASTQ DNA 4,032,514 407,283,914 101 101 101

fastq_db/SRR453569_2.fastq FASTQ DNA 4,032,514 407,283,914 101 101 101

fastq_db/SRR453570_1.fastq FASTQ DNA 6,745,975 681,343,475 101 101 101

fastq_db/SRR453570_2.fastq FASTQ DNA 6,745,975 681,343,475 101 101 101

fastq_db/SRR453571_1.fastq FASTQ DNA 6,163,396 622,502,996 101 101 101

fastq_db/SRR453571_2.fastq FASTQ DNA 6,163,396 622,502,996 101 101 101

２、トリミング (trimmomaticはbrewで導入できます)。

trimmomatic PE -phred33 SRR453569_1.fastq SRR453569_2.fastq SRR453569_pair1.fastq U1.fastq SRR453569_pair2.fastq U2.fastq LEADING:20 TRAILING:20 MINLEN:50

trimmomatic PE -phred33 SRR453570_1.fastq SRR453570_2.fastq SRR453570_pair1.fastq U1.fastq SRR453570_pair2.fastq U2.fastq LEADING:20 TRAILING:20 MINLEN:50

trimmomatic PE -phred33 SRR453571_1.fastq SRR453571_2.fastq SRR453571_pair1.fastq U1.fastq SRR453571_pair2.fastq U2.fastq LEADING:20 TRAILING:20 MINLEN:50

３、fastqを１つにまとめ、Trinityでアセンブルする（プールが大きいほど様々なmRNAのアセンブルが期待できる）。今回は３データをまとめアセンブル。

cat *pair1.fastq | sed 's/ HWI/\/1 HWI/g' > all_1.fastq 
cat *pair2.fastq | sed 's/ HWI/\/2 HWI/g' > all_2.fastq

４、１つにまとめたfastqを使いtrinityでde novo アセンブル (Trinityもbrewで導入できます)。

Trinity --max_memory 30G --seqType fq --left all_1.fastq --right all_2.fastq

５、パターン１、salmonでリードをマッピング (salmonもbrewで導入可)

#index
cp trinity_out_dir/Trinity.fasta . #fastaをカレントに移動
salmon index --index Trinity --transcripts Trinity.fasta

#以下のシェルスクリプトsalomn_loop.shを実行する。
FILES=`ls SRR*_P1.fastq | sed 's/_P1.fastq//g'` 
for F in $FILES ; do 
  R1=${F}_P1.fastq 
  R2=${F}_P2.fastq 
  salmon quant --index Trinity --libType IU --dumpEq -1 $R1 -2 $R2 --output ${F}.out 
done

bash salmon_loop.sh #実行

salmonの出力ディレクトリのaux_info/eq_classes.txtを指定して、Corsetをラン。

corset -g 1,1,1,2,2,2 -n A1,A2,A3,B1,B2,B3 -i salmon_eq_classes SRR*/aux_info/eq_classes.txt

６、パターン２、bowiteでリードをマッピング

bowtieを使うなら、以下のように進める。

#index
cp cp trinity_out_dir/Trinity.fasta . #fastaをカレントに移動
bowtie-build Trinity.fasta Trinity

#以下のシェルスクリプトbowtie_loop.shを実行する。
FILES=`ls *_1.fastq | sed 's/_1.fastq//g'`
for F in $FILES ; do
　R1=${F}_1.fastq
　R2=${F}_2.fastq
　bowtie --all -S Trinity -p 20 -1 $R1 -2 $R2 > ${F}.sam #20スレッド
　samtools view -@ 10 -S -b ${F}.sam > ${F}.bam
done

bash bowtie_loop.sh

出力bamを指定してCorsetをラン。

corset -g 1,1,1,2,2,2 -n A1,A2,A3,B1,B2,B3 *.bam

修正途中のまま時間が経ってしまったのでまとめ直しました。

引用

Corset: enabling differential gene expression analysis for de novoassembled transcriptomes

Nadia M Davidson and Alicia Oshlack

Genome Biology201415:410

macでインフォマティクス

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

de novo assemblyで得たRNAのコンティグをクラスター化して、非モデル生物のDEG解析を可能にする Corset