基于Stringtie的转录组数据分析与后续处理指南
一.StringTie介绍
StringTie 是用于 RNA-seq 的转录本组装和定量软件,StringTie 可以看做是cufflinks软件的升级版本,其功能和Cufflinks是一样的,包括下面两个主要功能:转录本组装和定量;相比Cuffinks, 其运行速度更快。该软件的官网:https://ccb.jhu.edu/software/stringtie/index.shtml。
1、Stringtie通过使用genome指导的组装的方法与从头组装的概念结合的新方法来改善转录组组装。
2、Stringtie的输入不仅可以是经过比对的结果,也可以是Stringtie通过从头组装read所得到的contig,当这两种输入都用到的时候,我们下面称之为stringtie+SR。
3、对于很多使用参考基因组辅助组装的方法,组装的的策略都是先对read进行cluter,然后建立一个graph model来推测每个基因所有可能的isoform,最终通过不同的graph的解析方法得到对转录本的组装结果。
4、有名的cufflinks用的是overlap graph,该模型中nodes代表fragment,如果两个fragment存在overlap并存在兼容的剪切模式,则对应的node连接起来。其解析方法为一种保守的算法,可以产生能够解释所有read的最少的转录本,尽管这种方法很吸引人,但是它没有考虑到转录本的丰度并且对于某些isoform来说该方法没有办法组装!
5、stringtie采用了组装转录本和估计表达量同步进行的方法,这不同于cufflinks的先组装后定量的策略。
6、首先,stringtie将read聚成cluster,然后采用了splice graph,其中node代表外显子或外显子的一部分,path将graph中可能 的剪切现象都连起来,最终对每个转录本通过创建一个网络流的方法,利用最大流算法(maximum flow algorithm)估计转录本的表达量。
7、最大流的问题是最优理论中的经典问题,但是目前还没有应用到转录本定量中。
8、与其它组装软件相比,stringtie具有很高的准确性和新型isoform的发现能力,其优势在于使用了网络流算法,同时stringtie也支持将read从头组装成更长的片段,这进一步提高了其组装的正确性。
9、其另一个优势在于它的最优化策略,它平衡了每次组装中每条转录本的覆盖度,这样可以对组装算法产生一定的限制,因为在组装基因组时,覆盖度是很重要的一个参数因为它需要被用来限制算法,否则组装器可能将重复的片段错误地堆叠到一起,相似地转录组装也是如此,在isoform中的每一个外显子需要有相似的覆盖度,如果忽略这个参数可能会产生一些保守但是错误的转录本,其中含有大量剪切位点的基因组装起来尤其困难。
常用的参数及描述:
-o [<path/>]<out.gtf> #设置StringTie组装转录本的输出GTF文件的路径和文件名。此处可指定完整路径,在这种情况下,将根据需要创建目录。默认情况下,StringTie将GTF写入标准输出。
-p <int> #指定组装转录本的线程数(CPU)。默认值是1
-G <ref_ann.gff> #使用参考注释基因文件指导组装过程,格式GTF/GFF3。输出文件中既包含已知表达的转录本,也包含新的转录本。选项-B,-b,-e,-C需要此选项(详情如下)
-l <label> #将<label>设置为输出转录本名称的前缀。默认:STRG
-A <gene_abund.tab> #输出基因丰度的文件(制表符分隔格式)
-C <cov_refs.gtf> #输出所有转录本对应的reads覆盖度的文件,此处的转录本是指参考注释基因文件中提供的转录本。(需要参数 -G).
-B #应用该选项,则会输出Ballgown输入表文件(* .ctab),其中包含用-G选项给出的参考转录本的覆盖率数据。(有关这些文件的说明,请参阅Ballgown文档。)如果选项-o 给出输出转录文件的完整路径,则* .ctab文件与输出GTF文件在相同的目录下。
-b <path> #指定 *.ctab 文件的输出路径, 而非由-o选项指定的目录。注意: 建议在使用-B/-b选项中同时使用-e选项,除非StringTie GTF输出文件中仍需要新的转录本,-B和-b选一个使用就行。
-e #限制reads比对的处理,仅估计和输出与用-G选项给出的参考转录本匹配的组装转录本。使用该选项,则会跳过处理与参考转录本不匹配的组装转录本,这将大大的提升了处理速度。
--merge #转录本合并模式。在合并模式下,StringTie将所有样品的GTF/GFF文件列表作为输入,并将这些转录本合并/组装成非冗余的转录本集合。这种模式被用于新的差异分析流程中,用以生成一个跨多个RNA-Seq样品的全局的、统一的转录本。如果提供了-G选项(参考注释基因组文件),则StringTie将从输入的GTF文件中将参考转录本组装到transfrags中。(个人理解:transfrags可能指的是拼接成更大的转录本片段,tanscript fragments)。在此模式下可以使用以下附加选项:
-G <guide_gff> #参考注释基因组文件(GTF/GFF3)
-o <out_gtf> #指定输出合并的GTF文件的路径和名称 (默认值:标准输出)
-m <min_len> #合并文件中,指定允许最小输入转录本的长度 (默认值: 50)
-c <min_cov> #合并文件中,指定允许最低输入转录本的覆盖度(默认值: 0)
-F <min_fpkm> #合并文件中,指定允许最低输入转录本的FPKM值 (默认值: 0)
-T <min_tpm> #合并文件中,指定允许最低输入转录本的TPM值 (默认值: 0)
-f <min_iso> #minimum isoform fraction (默认值: 0.01)
-i #合并后,保留含retained introns的转录本 (默认值: 除非有强有力的证据,否则不予保留)
-l <label> #输出转录本的名称前缀 (默认值: MSTRG)
二.StringTie应用
1.单样本转录本组装
输入文档是BAM格式的比对结果文档,该文档必需经过排序。这些文档可以是来源于Tophat比对的结果文档,也可以是hisat2的结果文档经过转换和排序的文档(使用samtools)。我们前面使用的就是hisat2比对后用samtools排序后的bam文件。除此以外,我们还需要gtf注释文件。关于gtf注释文件格式参考文章:生信中常见的数据文件格式。
gtf注释文件可以去genecode下载你需要的gtf文件,我这里下载的是小鼠的。
genecode官网:https://www.gencodegenes.org/
我的gtf文件在/data/mouse_annotation/ 目录下:
单样本组装:
stringtie -p 8 -G /data/mouse_annotation/gencode.mouse.annotation.gtf -o cleandata/stringtiedata/CK-4.gtf cleandata/samtools_bam/CK-4_sort.bam
批量处理:
for i in CK-7 CK-8 HGJ-10 HGJ-6 HGJ-9;do stringtie -p 8 -G /data/mouse_annotation/gencode.mouse.annotation.gtf -o cleandata/stringtiedata/${i}.gtf cleandata/samtools_bam/${i}_sort.bam;done
查看输出结果:
ll -h cleandata/stringtiedata/
2.多样本转录本整合
stringtie --merge [options] gtf.list :转录组merge模式,在该模式下,Stringtie可以利用输入的一个gtf list并将他们中的转录本进行非冗余的整理。可以在处理多个RNA-seq样本的时候,由于转录组存在时空特异性,可以将每个样本各自的转录组进行非冗余的整合,如果-G提供了参考gtf文件,可以将其一起整合到一个文件中,最终输出成一个完整的gtf文件。
将文件名的完整路径输入到一个文件中:
find /data/RNAseq/cleandata/stringtiedata/ -name *.gtf > /data/RNAseq/cleandata/stringtiedata/merglist.txt
查看一下merglist.txt文件中的内容:
cat /data/RNAseq/cleandata/stringtiedata/merglist.txt
stringtie --merge -p 8 -G /data/mouse_annotation/gencode.mouse.annotation.gtf -o cleandata/stringtiedata/stringtie_merged.gtf cleandata/stringtiedata/merglist.txt
3. 使用gffcompare检验数据比对到基因组上的情况(可选)
程序gffcompare可用于比较、合并、注释和估计一个或多个GFF文件(“查询”文件)的准确性。
gffcompare [options]* {-i <input_gtf_list> | <input1.gtf> [<input2.gtf> .. <inputN.gtf>]}
参数:
-i provide a text file with a list of (query) GTF files to process instead of expecting them as command-line arguments (useful when a large number of GTF files should be processed).
-h/-help Prints the help message and exits
-v/--version Prints the version number and exits
-o <outprefix> All output files created by gffcompare will have this prefix (e.g. .loci, .tracking, etc.). If this option is not provided the default output prefix being used is: “gffcmp”
-r An optional “reference” annotation GFF file. Each sample is matched against this file, and sample isoforms are tagged as overlapping, matching, or novel where appropriate. See the .refmap and .tmap output file descriptions below.
-R If -r was specified, this option causes gffcompare to ignore reference transcripts that are not overlapped by any transcript in one of input1.gtf,…,inputN.gtf. Useful for ignoring annotated transcripts that are not present in your RNA-Seq samples and thus adjusting the “sensitivity” calculation in the accuracy report written in the file.
-Q If -r was specified, this option causes gffcompare to ignore input transcripts that are not overlapped by any transcript in the reference. Useful for adjusting the “precision” calculation in the accuracy report written in the file.
-M discard (ignore) single-exon transfrags and reference transcripts (i.e. consider only multi-exon transcripts)
-N discard (ignore) single-exon reference transcripts; single-exon transfrags are still considered, but they will never find an exact match
-D discard "duplicate" (redundant) query transfrags (i.e. those with the same intron chain) within a single sample (and thus disable "annotation" mode)
-s <genome_path> path to genome sequences (optional); this will enable the "repeat" ('r') classcode assessment; <genome_path> should be a full path to a multi-FASTA file, preferrably indexed with samtools faidx; repeats must be soft-masked (lower case) in the genomic sequence
-e <dist> Maximum distance (range) allowed from free ends of terminal exons of reference transcripts when assessing exon accuracy. By default, this is 100.
-d <dist> Maximum distance (range) for grouping transcript start sites, by default 100.
-p <tprefix> The name prefix to use for consensus/combined transcripts in the <outprefix>.combined.gtf file (default: 'TCONS')
-C Discard the “contained” transfrags from the .combined.gtf output. By default, without this option, gffcompare writes in that file isoforms that were found to be fully contained/covered (with the same compatible intron structure) by other transfrags in the same locus, with the attribute “contained_in” showing the first container transfrag found. (Note: this behavior is the opposite of Cuffcompare's -C option).
-A Like -C but will not discard intron-redundant transfrags if they start on a different 5' exon (keep alternate transcript start sites)
-X Like -C but also discard contained transfrags if transfrag ends stick out within the container's introns
-K For -C/-A/-X, do NOT discard any redundant transfrag matching a reference
-T Do not generate .tmap and .refmap files for each input file
-V Gffcompare is a little more verbose about what it's doing, printing messages to stderr, and it will also show warning messages about any inconsistencies or potential issues found while reading the given GFF file(s).
--debug Enables debug mode, which enables -V and generates additional files: <outprefix>.Q_discarded.lst, <outprefix>.missed_introns.gtf and <outprefix>.R_missed.lst
参考:http://ccb.jhu.edu/software/stringtie/gffcompare.shtml
我们将前面我们整合的转录本文件检验数据比对到基因组上的情况
gffcompare -r /data/mouse_annotation/gencode.mouse.annotation.gtf -G cleandata/stringtiedata/stringtie_merged.gtf
会在输入的文件所在文件夹下产生几个文件:
其中gffcompare.annotated.gtf存储的是StringTie组装的转录本与注释文件内的转录本的差别信息,通过class_code来表示:
head cleandata/stringtiedata/gffcompare.annotated.gtf |grep class_code | cut -d ";" -f 5-8
(RNAseq1) [root@localhost RNAseq]# head cleandata/stringtiedata/gffcompare.annotated.gtf |grep class_code | cut -d ";" -f 5-8
ref_gene_id "ENSMUSG00000116111.1"; cmp_ref "ENSMUST00000229069.1"; class_code "="; tss_id "TSS1"
ref_gene_id "ENSMUSG00000115879.1"; cmp_ref "ENSMUST00000229227.1"; class_code "="; tss_id "TSS2"
ref_gene_id "ENSMUSG00000116214.1"; contained_in "ENSMUST00000229087.1"; cmp_ref "ENSMUST00000229989.1"; class_code "="
对照表如下:
gffcmp.stats文件存储比对结果的准确性和预测率。
cat cleandata/stringtiedata/gffcompare.stats
# gffcompare v0.11.2 | Command line was:
#gffcompare -r /data/mouse_annotation/gencode.mouse.annotation.gtf -G cleandata/stringtiedata/stringtie_merged.gtf -o cleandata/stringtiedata/gffcompare
#
#= Summary for dataset: cleandata/stringtiedata/stringtie_merged.gtf
# Query mRNAs : 155156 in 54744 loci (127176 multi-exon transcripts)
# (20889 multi-transcript loci, ~2.8 transcripts per locus)
# Reference mRNAs : 143543 in 54276 loci (116302 multi-exon)
# Super-loci w/ reference transcripts: 53940
#-----------------| Sensitivity | Precision |
Base level: 100.0 | 97.3 |
Exon level: 94.0 | 96.1 |
Intron level: 100.0 | 97.7 |
Intron chain level: 100.0 | 91.4 |
Transcript level: 99.5 | 92.1 |
Locus level: 99.7 | 98.3 |
Matching intron chains: 116302
Matching transcripts: 142855
Matching loci: 54117
Missed exons: 0/451928 ( 0.0%)
Novel exons: 3037/433435 ( 0.7%)
Missed introns: 6/288208 ( 0.0%)
Novel introns: 2186/295137 ( 0.7%)
Missed loci: 0/54276 ( 0.0%)
Novel loci: 804/54744 ( 1.5%)
Total union super-loci across all input datasets: 54744
155156 out of 155156 consensus transcripts written in cleandata/stringtiedata/gffcompare.annotated.gtf (0 discarded as redundant)
4.重新组装转录本并估算基因表达丰度
利用merge得到的gtf重新对各个样本做定量,并创建ballgown可读取文件。
mkdir cleandata/ballgown
for i in CK-4 CK-7 CK-8 HGJ-10 HGJ-6 HGJ-9; do stringtie -e -B -p 8 -G cleandata/stringtiedata/stringtie_merged.gtf -o cleandata/ballgown/${i}/${i}.gtf cleandata/samtools_bam/${i}_sort.bam; done
查看ballgown目录下生成的文件:
每一个样本对应一个文件夹。
查看其中一个:
5.read count数据输出
这里需要prepDE.py这个脚本。
prepDE.py -i cleandata/ballgown/
会在当前文件夹下产生2个csv文件。
分别查看一下:
推荐阅读
-
基于Stringtie的转录组数据分析与后续处理指南
-
windows下进程间通信的(13种方法)-摘 要 本文讨论了进程间通信与应用程序间通信的含义及相应的实现技术,并对这些技术的原理、特性等进行了深入的分析和比较。 ---- 关键词 信号 管道 消息队列 共享存储段 信号灯 远程过程调用 Socket套接字 MQSeries 1 引言 ---- 进程间通信的主要目的是实现同一计算机系统内部的相互协作的进程之间的数据共享与信息交换,由于这些进程处于同一软件和硬件环境下,利用操作系统提供的的编程接口,用户可以方便地在程序中实现这种通信;应用程序间通信的主要目的是实现不同计算机系统中的相互协作的应用程序之间的数据共享与信息交换,由于应用程序分别运行在不同计算机系统中,它们之间要通过网络之间的协议才能实现数据共享与信息交换。进程间通信和应用程序间通信及相应的实现技术有许多相同之处,也各有自己的特色。即使是同一类型的通信也有多种的实现方法,以适应不同情况的需要。 ---- 为了充分认识和掌握这两种通信及相应的实现技术,本文将就以下几个方面对这两种通信进行深入的讨论:问题的由来、解决问题的策略和方法、每种方法的工作原理和实现、每种实现方法的特点和适用的范围等。 2 进程间的通信及其实现技术 ---- 用户提交给计算机的任务最终都是通过一个个的进程来完成的。在一组并发进程中的任何两个进程之间,如果都不存在公共变量,则称该组进程为不相交的。在不相交的进程组中,每个进程都独立于其它进程,它的运行环境与顺序程序一样,而且它的运行环境也不为别的进程所改变。运行的结果是确定的,不会发生与时间相关的错误。 ---- 但是,在实际中,并发进程的各个进程之间并不是完全互相独立的,它们之间往往存在着相互制约的关系。进程之间的相互制约关系表现为两种方式: ---- (1) 间接相互制约:共享CPU ---- (2) 直接相互制约:竞争和协作 ---- 竞争——进程对共享资源的竞争。为保证进程互斥地访问共享资源,各进程必须互斥地进入各自的临界段。 ---- 协作——进程之间交换数据。为完成一个共同任务而同时运行的一组进程称为同组进程,它们之间必须交换数据,以达到协作完成任务的目的,交换数据可以通知对方可以做某事或者委托对方做某事。 ---- 共享CPU问题由操作系统的进程调度来实现,进程间的竞争和协作由进程间的通信来完成。进程间的通信一般由操作系统提供编程接口,由程序员在程序中实现。UNIX在这个方面可以说最具特色,它提供了一整套进程间的数据共享与信息交换的处理方法——进程通信机制(IPC)。因此,我们就以UNIX为例来分析进程间通信的各种实现技术。 ---- 在UNIX中,文件(File)、信号(Signal)、无名管道(Unnamed Pipes)、有名管道(FIFOs)是传统IPC功能;新的IPC功能包括消息队列(Message queues)、共享存储段(Shared memory segment)和信号灯(Semapores)。 ---- (1) 信号 ---- 信号机制是UNIX为进程中断处理而设置的。它只是一组预定义的值,因此不能用于信息交换,仅用于进程中断控制。例如在发生浮点错、非法内存访问、执行无效指令、某些按键(如ctrl-c、del等)等都会产生一个信号,操作系统就会调用有关的系统调用或用户定义的处理过程来处理。 ---- 信号处理的系统调用是signal,调用形式是: ---- signal(signalno,action) ---- 其中,signalno是规定信号编号的值,action指明当特定的信号发生时所执行的动作。 ---- (2) 无名管道和有名管道 ---- 无名管道实际上是内存中的一个临时存储区,它由系统安全控制,并且独立于创建它的进程的内存区。管道对数据采用先进先出方式管理,并严格按顺序操作,例如不能对管道进行搜索,管道中的信息只能读一次。 ---- 无名管道只能用于两个相互协作的进程之间的通信,并且访问无名管道的进程必须有共同的祖先。 ---- 系统提供了许多标准管道库函数,如: pipe——打开一个可以读写的管道; close——关闭相应的管道; read——从管道中读取字符; write——向管道中写入字符; ---- 有名管道的操作和无名管道类似,不同的地方在于使用有名管道的进程不需要具有共同的祖先,其它进程,只要知道该管道的名字,就可以访问它。管道非常适合进程之间快速交换信息。 ---- (3) 消息队列(MQ) ---- 消息队列是内存中独立于生成它的进程的一段存储区,一旦创建消息队列,任何进程,只要具有正确的的访问权限,都可以访问消息队列,消息队列非常适合于在进程间交换短信息。 ---- 消息队列的每条消息由类型编号来分类,这样接收进程可以选择读取特定的消息类型——这一点与管道不同。消息队列在创建后将一直存在,直到使用msgctl系统调用或iqcrm -q命令删除它为止。 ---- 系统提供了许多有关创建、使用和管理消息队列的系统调用,如: ---- int msgget(key,flag)——创建一个具有flag权限的MQ及其相应的结构,并返回一个唯一的正整数msqid(MQ的标识符); ---- int msgsnd(msqid,msgp,msgsz,msgtyp,flag)——向队列中发送信息; ---- int msgrcv(msqid,cmd,buf)——从队列中接收信息; ---- int msgctl(msqid,cmd,buf)——对MQ的控制操作; ---- (4) 共享存储段(SM) ---- 共享存储段是主存的一部分,它由一个或多个独立的进程共享。各进程的数据段与共享存储段相关联,对每个进程来说,共享存储段有不同的虚拟地址。系统提供的有关SM的系统调用有: ---- int shmget(key,size,flag)——创建大小为size的SM段,其相应的数据结构名为key,并返回共享内存区的标识符shmid; ---- char shmat(shmid,address,flag)——将当前进程数据段的地址赋给shmget所返回的名为shmid的SM段; ---- int shmdr(address)——从进程地址空间删除SM段; ---- int shmctl (shmid,cmd,buf)——对SM的控制操作; ---- SM的大小只受主存限制,SM段的访问及进程间的信息交换可以通过同步读写来完成。同步通常由信号灯来实现。SM非常适合进程之间大量数据的共享。 ---- (5) 信号灯 ---- 在UNIX中,信号灯是一组进程共享的数据结构,当几个进程竞争同一资源时(文件、共享内存或消息队列等),它们的操作便由信号灯来同步,以防止互相干扰。 ---- 信号灯保证了某一时刻只有一个进程访问某一临界资源,所有请求该资源的其它进程都将被挂起,一旦该资源得到释放,系统才允许其它进程访问该资源。信号灯通常配对使用,以便实现资源的加锁和解锁。 ---- 进程间通信的实现技术的特点是:操作系统提供实现机制和编程接口,由用户在程序中实现,保证进程间可以进行快速的信息交换和大量数据的共享。但是,上述方式主要适合在同一台计算机系统内部的进程之间的通信。 3 应用程序间的通信及其实现技术 ---- 同进程之间的相互制约一样,不同的应用程序之间也存在竞争和协作的关系。UNIX操作系统也提供一些可用于应用程序之间实现数据共享与信息交换的编程接口,程序员可以通过自己编程来实现。如远程过程调用和基于TCP/IP协议的套接字(Socket)编程。但是,相对普通程序员来说,它们涉及的技术比较深,编程也比较复杂,实现起来困难较大。 ---- 于是,一种新的技术应运而生——通过将有关通信的细节完全掩盖在某个独立软件内部,即底层的通讯工作和相应的维护管理工作由该软件内部来实现,用户只需要将通信任务提交给该软件去完成,而不必理会它的具体工作过程——这就是所谓的中间件技术。 ---- 我们在这里分别讨论这三种常用的应用程序间通信的实现技术——远程过程调用、会话编程技术和MQSeries消息队列技术。其中远程过程调用和会话编程属于比较低级的方式,程序员参与的程度较深,而MQSeries消息队列则属于比较高级的方式,即中间件方式,程序员参与的程度较浅。 ---- 4.1 远程过程调用(RPC)
-
实例解析:使用Hisat2, StringTie与Ballgown进行拟南芥转录组数据的分析步骤