macでインフォマティクス

macでインフォマティクス

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

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

 

2、トリミング (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

 

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

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

 

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

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

 

5、パターン1、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

 

6、パターン2、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