- 我这里只是用了4个数据进行测试,分成两组(随意写的,仅供测试,在ballgown可以看到);
- hisat2比对后,samtools sort生成sort.bam,用于stringtie的输入;
- stringtie的推荐版本使用bam生成了gtf格式的文件,后利用gtf生成非冗余转录本,用于下一步分析;
- stringtie的fast版本则直接生成了最终结果;
1.Hisat2
- 上游利用了hisat2进行比对(网上教程颇多,这里不再赘述),整个采用的是hisat2+StringTie+ballgown;
2.StringTie
- 差异分析应用,StringTie有两种模式,官方推荐和fast模式;
- StringTie --merge 功能可生成在所有样本中均出现的非冗余转录本;输入文件是上一步生成的组装文件(GTF格式)和参考注释文件(我这里使用的是gencodeV29);
2.1.StringTie-recommend
image2.1.1生成每个样本的组装文件(gtf格式)
R=/teach/database/gtf/gencode.v29.annotation.gtf
for file in $(ls ~/airway/04align/*.sort.bam);
do
out_file1=${file##*/}
out_file=${out_file1%.sort.bam}
#echo $out_file1 $out_file
stringtie ${file} -l ${out_file} -p 5 -o ~/airway/stringtie_recommend/${out_file}_hg38.gtf -G $R 2>~/airway/log/stringtie_${out_file}.log2
done
2.1.2 生成mergelist.txt(即之前生成的gtf文件)
ls ~/airway/stringtie_recommend/*gtf > mergelist.txt
##基于组装结果生成非冗余转录本
R=/teach/database/gtf/gencode.v29.annotation.gtf
stringtie --merge -p 8 -G $R -o stringtie_merged.gtf ~/airway/stringtie_recommend/mergelist.txt
2.1.3生成ballgown需要的结果
- -B/-b 估计转录本丰度并生成Ballgown所需的文件
- -e(只采纳给定的参考转录本)
R=/teach/database/gtf/gencode.v29.annotation.gtf
for file in $(ls ~/airway/04align/SRR1039513.sort.bam);
do
out_file1=${file##*/}
out_file=${out_file1%.sort.bam}
stringtie -e -B -p 8 -G stringtie_merged.gtf -o ~/airway/stringtie_recommend/${out_file}/${out_file}_hg38.gtf $file 2>~/airway/log/stringtie_re_${out_file}.log2
done
image
2.2.StringTie-fast
imageR=/teach/database/gtf/gencode.v29.annotation.gtf
for file in $(ls ~/airway/04align/SRR1039513.sort.bam);
do
out_file1=${file##*/}
out_file=${out_file1%.sort.bam}
mkdir ~/airway1/08stringtie/${out_file}
#echo $out_file1 $out_file
stringtie ${file} -B ~/airway1/08stringtie/${out_file} -o ~/airway1/08stringtie/${out_file}/${out_file}.gtf -G $R 2>~/airway1/log/stringtie_${out_file}.log2
done
image
3.ballgown
因为我使用的数据是抽取了部分reads往下进行的,所以,这里未展示差异分析结果;
library(ballgown)
bg = ballgown(dataDir='~/ballgown/08stringtie', samplePattern='SRR', meas='all')
pData(bg) = data.frame(id=sampleNames(bg), group=rep(c(1,0), c(2,2)))
stat_results = stattest(bg, feature='gene', meas='FPKM', covariate='group',getFC = T)
head(stat_results)
参考内容:
1.stringtie_manual
网友评论