Stringtie的使用

Posted 2023-04-30

参考技术A

本文用作个人学习， 全部摘抄于 Stringtie说明书中文翻译版
参考链接： Stringtie说明书
参考链接： Stringtie说明书中文翻译版

StringTie的基本用法 ： stringtie <aligned_reads.bam> [options]*

其中， aligned_reads.bam 是输入文件，该输入文件要求必须按其基因组位置排序， HISAT2的输出文件则需经过samtools sort生成的bam文件才可当做输入文件。

其他可选参数：

其中， aligned_reads.bam 是输入文件，该输入文件要求必须按其基因组位置排序，如TopHat的输出文件 accepted_hits.bam 可直接当做输入文件，而 HISAT2的输出文件则需经过samtools sort生成的bam文件才可当做输入文件。

输入BAM文件中的每个 spliced read 比对（即跨越至少一个连接点的比对）必须包含标签XS，用以指示测序产生的read是来源于基因组序列上的哪条链产生的RNA。由TopHat和 HISAT2 (需参数 --dta，该参数用于发现剪接位点) 产生的比对结果中已经包含标签XS。但是，有的mapping程序(read mapper)未必含有标签XS，所以，用户在进行下一步分析时需要进行检查。

注意：一定要使用-dta选项来运行HISAT2，否则结果将会受到影响。

作为选项，可以向StringTie提供GTF / GFF3格式的参考注释基因组文件。在这种情况下，StringTie更喜欢使用注释文件中的这些“已知”基因，对于那些被表达的基因，它将计算coverage，TPM和FPKM值。它还会产生额外的转录本，而注释文件中并没有这些转录本。请注意，如果不使用选项-e，那么参考转录本就需要被reads 完全覆盖，以便包含在StringTie的输出中。在这种情况下，其他通过StringTie从数据中组装的转录本，且不在注释文件中的转录本也会输出。

注意：如果用户正在分析注释较好的基因组，例如人类，小鼠或其他模型生物，则强烈建议您提供注释文件。

主要输出文件有：

1、 GTF文件： 记录组装的转录本信息

2、 Tab文件： 记录基因丰度信息

3、 GTF文件：完全覆盖与参考注释基因组文件所匹配的转录本信息

4、 *.ctab文件：用于下游Ballgown软件做差异表达分析的输入文件

5、 GTF文件：在合并模式下，生成一个合并的GTF文件

如果StringTie使用-A <gene_abund.tab>选项运行，则返回包含基因丰度的文件。

如果StringTie与 -C <cov_refs.gtf> 选项一起运行（需要选项-G

如果StringTie与 -B 选项一起运行，它将返回Ballgown输入文件，包含以下文件：(1) e2t.ctab, (2) e_data.ctab, (3) i2t.ctab, (4) i_data.ctab, and (5) t_data.ctab。

如果StringTie使用 --merge 选项运行，它将多个GTF / GFF文件作为输入，并将这些转录本合并和组装成非冗余转录本集合。

HISAT2+StringTie+Ballgown安装及使用流程

HISAT2+StringTie+Ballgown安装及使用流程

2015年Nature Methods上面发表了一款快速比对工具hisat，作为接替tophat和bowtie的比对工具，它具有更快的比对速度和更高的比对率，最近把这个流程走完一遍，感觉优势还是很明显的。 
一、HISAT2： 
1、下载安装： 
hisat2下载地址：ftp://ftp.ccb.jhu.edu/pub/infphilo/hisat2/downloads/hisat2-2.1.0-Linux_x86_64.zip 
hisat2官方手册：http://ccb.jhu.edu/software/hisat2/manual.shtml 
下载完成后解压缩： 
unzip hisat2-2.0.5-Linux_x86_64.zip 
进入hisat2-2.0.5文件夹： 

这里面的绿色文件都是可执行文件，所以只需要把目录添加到环境变量中即可： 
vim进入编辑bashrc文件，在文本中输入红色方框内的内容，保存退出，然后source ~/.bashrc 即可 

此时我们就可以直接调用hisat2命令了。 
2、建立索引： 
如同tophat一样，比对之前需要利用bowtie建立index，hisat2同样需要建立索引： 
首先提取gtf文件中的剪切位点和外显子位置： 
extract_splice_sites.py gencode.vM4.annotation.gtf >gencode.vM4.annotation.for.hisat2.ss 
extract_exons.py gencode.vM4.annotation.gtf >gencode.vM4.annotation.for.hisat2.exon 
建立索引： 
hisat2-build -p 30 --ss gencode.vM4.annotation.for.hisat2.ss --exon gencode.vM4.annotation.for.hisat2.exon GRCm38.p3.genome.fa mouseGencodeIndex 
##如果电脑内存<200G，那么可以不用--ss/--exon参数，但是在比对的时候需要加 
--known-splicesite-infile参数。3、比对： 
我的数据是双段的无链特异性数据，此处需要把sam文件转化为bam文件，所以需要提前安装samtools： 
        hisat2 --known-splicesite-infile gencode.vM4.annotation.for.hisat2.ss --dta -t -p 24 -x mouseGencodeIndex -1 samp_1.fq.gz -2 samp_2.fq.gz -S accepted_hits.sam &> alignment_summary.txt 
       samtools view -bS accepted_hits.sam > accepted_hits.bam 
       samtools sort accepted_hits.bam -o accepted_hits_sorted.bam 
       rm accepted_hits.bam 
       rm accepted_hits.sam 

二、StringTie： 
比对完生成了sam文件，我们利用samtools将它转化为了排好序的bam文件，下一步就需要量化和确定表达值了，这里用到的StringTie相比之前的cufflinks来说功能强大了好多。 
1、下载安装： 
stringtie下载地址：http://ccb.jhu.edu/software/stringtie/dl/stringtie-1.3.3b.Linux_x86_64.tar.gz 
stringtie官方手册：http://ccb.jhu.edu/software/stringtie/index.shtml?t=manual 
直接下载解压就可以用了，它是可执行文件，也可以按上述方法将路径添加到环境变量中方便调用。 
2、运行： 
stringtie accepted_hits_sorted.bam -o outRes.gtf -p 28 -G gencode.vM4.annotation.gtf -A gene_abund.tab -B -e 
运行后每个样本文件夹下结果如下： 

这里我生成了结果gtf文件outRes.gtf和ballgown需要的.ctab文件，还有基因的表达量文件gene_abund.tab，该文件包括基因的表达量FPKM以及TPM等。当然如果你想要转录本的表达量，直接打开t_data.ctab这个文件，这里面有转录本的FPKM值。 
当然如果我们想利用DESeq2或者edgeR等计算差异表达，那我们就需要得到原始counts值矩阵来作为输入，此时我们需要利用StringTie自带的脚本prepDE.py来计算counts值，它可以同时对多个样本做： 
prepDE.py -i stringtieRes/ -g countsRes/gene_count_matrix.csv -t countsRes/transcript_count_matrix.csv 
stringtieRes/文件夹下面是我所有的样本的文件夹。 

*这里我们能得到所有样本的count matrix，但是只能拿到每个样本对应的FPKM值，又有什么方法能得到合并在一起的FPKM matrix呢？这就需要借助ballgown了。 
三、Ballgown： 
1、安装： 
首先你需要下载安装R，我的是3.4.0版本。 
source("https://bioconductor.org/biocLite.R") biocLite("ballgown") 
这里可能提示你安装XML包的时候会出现错误提示：Cannot find xml2-config 
这就需要你在自己电脑上安装相应的模块了，我的是centos7，于是安装相应的模块：
yum install libxml2-devel 
顺利安装上ballgown包。 
2、使用： 
读取所有样本到ballgown对象中：de> 
bg = ballgown(dataDir=de>de>de>YSde>, samplePattern=‘YT1‘, meas=‘all‘); 
#其中de>de>YS是我的所有样本的父目录，每个样本文件夹名字都包含YT1。 
#计算转录本和基因的FPKM值 
de>de>transcript_fpkm = texpr(bg, ‘FPKM‘) 
row.names(de>de>de>transcript_fpkmde>) = transcriptNames(bg) 
write.csv(de>de>transcript_fpkm,"de>de>de>transcript_fpkm_matrix.csvde>")
de>de>gene_expression = gexpr(bg) 
de>de>write.csv(de>de>de>gene_expressionde>,"de>de>de>de>gene_fpkmde>_matrix.csvde>") 
任务完成。 
3、差异表达分析： 
ballgown可以做case/control两两比较的差异表达，也可以做多组比较的差异表达（此时不能计算Fold Change值）， 
当然也可以做时间序列的差异。 
de>de>de>pData(bg) = data.frame(id=sampleNames(bg), group=rep(c(1,0), each=10)) 
#这里是条件矩阵，每行是一个样本，第二列是条件，如果是case/control那么就是0/1. 
de>de>de>de>stat_results = stattest(bg, feature=‘transcript‘, meas=‘FPKM‘, getFC=TRUE, covariate=‘group‘) 
#注意getFC在多组比较时候不能用，feature参数可以对基因‘gene‘或者转录本‘transcript‘或者外显子‘exon‘做 
差异表达分析。 
de>de>de>de>de>Data(bg) = data.frame(pData(bg), time=rep(1:10, 2)) #dummy time covariate timecourse_results = stattest(bg, feature=‘transcript‘, meas=‘FPKM‘, covariate=‘time‘, timecourse=TRUE)de> 
de> 
但是我个人不太推荐使用ballgown，喜欢使用DESeq2和edgeR来计算。 
de>