SuperTranscripts 其の２ - macでインフォマティクス

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

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

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

１、ダウンロード

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

２、ペアエンドの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

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

高速な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

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

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

６、クラスタリング

Corsetを使う（紹介）。

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

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

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

９、リードカウント。

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).

１０、 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")