macでインフォマティクス

macでインフォマティクス

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

SuperTranscripts 其の2

 

Githubで紹介されている、Supertransctiptsを非モデル生物のRNA seq解析に適用する場合の流れを確認する。以下のようなフローをとる。

f:id:kazumaxneo:20180717215632p:plain

de novoアセンブリして得たcontigにマッピングして定量し、定量結果からクラスタリングする。そのクラスタリング情報を使ってSupertransctiptsを構築する。最後に、Supertransctiptsを使ってマッピングして定量し、発現量の異なるSupertransctiptsを統計的に絞り込む。

 

mac os 10.12環境でテストした。

 

1、ダウンロード

ダウンロードしてFastqに変換。

#SRAに登録されているRNA seqデータを使う。
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/

 

2、ペアエンドのfastqへの変換

forループで処理。pfastq-dump(紹介)を使えばより早く終わります。

mkdir fastq

for srr in $(cat sra_run_ids.txt)
do
fastq-dump --split-files sra_files/${srr}.sra -O fastq/
done

 

3、クオリティフィルタリング

高速なfastp(紹介)を使い、アダプターとクオリティトリミングしてfastq.gz出力。forループで処理。

for srr in $(cat sra_run_ids.txt)
do
fastp -i fastq/${srr}_1.fastq -I fastq/${srr}_2.fastq \
-o fastq/${srr}_clean1.fastq.gz -O fastq/${srr}_clean2.fastq.gz \
-h fastq/${srr}_report.html
done

 

4、de novo アセンブリ

コンカテネートしたfastqを作り、アセンブリを行う。

#concatenate。公式では、このstepでsedを使い、Trinityのランに必要な/1を付加している。今回はfastq作成時にシングルだが--split-filesをつけて付加済み
cat fastq/*fastq.gz > fastq/all.fastq.gz

#assembly
Trinity --max_memory 40G --CPU 16 -seqType fq --single fastq/all.fastq --full_cleanup
mv trinity_out_dir.Trinity.fasta Trinity.fasta

#Remove fastq file
rm all.fastq

 

5、contigへのマッピング

HISAT2を使う(紹介)。

#index
hisat2-build trinity.fasta genome_index

#mapping
for srr in $(cat sra_run_ids.txt)
do
hisat2 -q -x genome_index -U fastq/"$srr"_1.fastq -p 16 -S "$srr".sam
#bamにしておく。IGVでみるならソートも必要。
samtools sort -@ 12 -m 8G -O BAM "$srr".sam -o "$srr".bam
done

#samが要らなければ消す
rm *sam

 

6、クラスタリング

Corsetを使う(紹介)。

for srr in $(cat sra_run_ids.txt)
do
corset -r true-stop "$srr".bam &
done
wait

#最初の3つはグループ1、残り4つはグループ2。-nでサンプル名を指定。
corset -g 1,1,1,2,2,2 -n A1,A2,A3,B2,B2,B3 -i corset *.corset-reads

 

7、SuperTranscripts構築

Corsetで得たクラスター情報を使い、de novo assemblyで作ったcontigからSuperTranscriptsを構築する。本ツールを使う(紹介)。

python Lace.py trinity.fasta clusters.txt -t --cores 16 -o FlyLace/

 

8、SuperTranscriptsへのマッピング

今度はSuperTranscriptsへマッピングする。

#index
hisat2-build FlyLace/SuperDuper.fasta superTranscripts_index

#mapping
for srr in $(cat sra_run_ids.txt)
do
hisat2 -q -x superTranscripts_index -U fastq/"$srr"_1.fastq -p 16 -S "$srr"ST.sam
#bamにしておく。IGVでみるならソートも必要。
samtools sort -@ 12 -m 8G -O BAM "$srr"ST.sam -o "$srr"ST.bam
done

#samが要らなければ消す
rm *ST.sam

 

9、リードカウント。

featurecountを使う(リンク)。

mkdir Counts

featureCounts -T 20 -p -B -O --fraction -a FlyLace/SuperDuper.gff -f -o Counts/counts.txt \
SRR074548ST.bam SRR074547ST.bam SRR074532ST.bam \
SRR031714ST.bam SRR031715ST.bam SRR031728ST.bam SRR031729ST.bam
  • -B    Only count read pairs that have both ends aligned.
  • -O    Assign reads to all their overlapping meta-features (or features if -f is specified).

 

10、 DTU検出

これで標準的なDTU(DTU)検出パイプラインを実行できるようになった。GithubではDexSeq(pubmed)またはDiffSplice(pubmed)を使うことが勧められている。しかし、ここではexamleの記載通りedgeR(紹介)のdiffを使う(DexSeqDiffSpliceより高速かつより保守的)。

library(edgeR) 
counts <- read.table("Counts/counts.txt",header=TRUE,sep="\t")

#Define groups
treatment= c('1','1','1','2','2','2','2')
sample = c('1','2','3','4','5','6','7')

#Make DGElist and normalise
dx <- DGEList(counts[,c(7:12)])
dx<-calcNormFactors(dx,group=treatment)

#Make exon id
eid = paste0(counts$Chr,":",counts$Start)

#Define design matrix
design <- model.matrix(~treatment)

#Voom transform
vx <- voom(dx,design)
#Fit with limma
fx <- lmFit(vx,design)

ex <- diffSplice(fx,geneid=counts$Chr,exonid=eid)
res<-topSplice(ex,number="20")
write.table(res[,c("GeneID","P.Value","FDR")],"VoomSplice.txt",col.names = FALSE,row.names=FALSE,quote=FALSE,sep="\t")

 

5−6のstepはオーサーらによりもう1つ代替手段が提示されている。詳細はリンク先のStep3を参照。

  

引用

SuperTranscripts: a data driven reference for analysis and visualisation of transcriptomes

Davidson NM, Hawkins ADK, Oshlack A

Genome Biol. 2017 Aug 4;18(1):148.

 

Example: Differential Transcript Usage on a non model organism · Oshlack/Lace Wiki · GitHub