本篇回顾一下转录组组装之前常用的一些方法，对您理解最新组装方法，如HISAT2+Stringtie也是有帮助的。分为两种主要方法：基于映射的组装+基于从头组装。

1、基于映射的组装

在这里，我们描述了两个软件包，Cufflinks（Stringtie为最新的类似于Cuufflinks软件）和Scripture,二者都是基于RNA-Seq测序数据映射构建全长转录本。这两种方法都可以用于ab initio构建转录本，即不需要额外的基因模型。这两个程序之间的主要区别在于解决异构体的方法：Scripture输出所有可能的异构体，而Cufflinks输出可以解释样本数据的最小可能异构体。输出BED或GTF格式文件，包含在参考序列上的转录本坐标。由于参考序列是已知，可以直接使用任何语言脚本利用转录序列坐标得到为FASTA文件，例如Python或Perl。映射由TopHat完成，本例使用的是2.0版。输入数据由18号染色体的FASTA 格式文件“chr18.fa”，双末端读段文件“chr18_1.fq” 和“chr18_2.fq”，由Burrows–Wheeler转化名为”chr18” 的索引文件，映射结果输出到“top2”目录中。因为读段为长度2 × 75 bp，文库片段插入长为200 bp，所以读段之间的距离为50 bp，故将TopHat参数“-r”值设定为50 。本例中，映射是用四个线程完成的。使用TopHat，必须调用SAMtools和Bowtie（在这里bowtie2），故必须将二者路径添加到全局变量PATH中。

$ bowtie2-build chr18.fa chr18

$ tophat2 -r 50 -p 4 -o top2 chr18 chr18_1.fq chr18_2.fq

1.1 Cufflinks

Cufflinks是利用C++语言编写的软件，已多次更新。最新版本可以在从

http://cuffflinks.cbcb.umd.edu下载，从网页可获取用户指南及更多的信息。通过构造一个图表对数据进行吝啬分析，也就是说，它得到一系列代表RNA-seq读段的最小转录本集，然后再通过得到的转录本计算表达丰度。

Cufflinks支持TopHat2（HISAT2）版本产生的BAM文件。为了能充分利用双末端测序信息，BAM文件的reads的名字应不包含reads两端后缀。TopHat可以准确读取BAM文件的双末端信息，即，双末端都映射上的序列会以 “=“标记。且即便没有被预定分隔符隔开，双末端的后缀也不会被移除。例如，“/1”和“/2”后缀会自动从读段名中移除，而“_1”和“_2”不会。为让Cufflinks利用双末端的信息，在BAM文件中，成对的读段标识符（文件的一列）应该一致。利用SAMtools 的view命令，以BAM文件作为输入，可以轻易检测。

为了使用Cufflinks,SAMtools的安装路径必须设定为全局变量，即添加到PATH变量中。

这里对几项参数进行了设定。为提高计算速度，在下面的命令中，采用4线程，其他参数采用默认值。 Cuffinks运行命令如下：

$ cufflinks –p 4 –o outdir top2/accepted _ hits.bam

基因模型被储存在输出目录的GTF文件中。有以下四个输出文件：

rw------- 1 somervuo 50K Jul 15 10:43 genes.fpkm_tracking

-rw------- 1 somervuo 67K Jul 15 10:43 isoforms.fpkm_tracking

-rw------- 1 somervuo 0 Jul 15 10:42 skipped.gtf

-rw------- 1 somervuo 898K Jul 15 10:43 transcripts.gtf

拥有外显子信息的转录本保存在“transcripts.gtf”文件中。本例中，得到634个基因产生的750条转录本，分别列举在 “genes.fpkm_tracking”和“isoforms.fpkm_tracking” 文件中。

如果有几个插入片段长度不同的文库，最好先利用Cufflinks分别进行组装，再将组装结果合并，而不是将其比对文件BAM合并后组装。Cuffmerge程序可以用来合并Cufflinks的运行结果。

合并独立的片段的数目可以通过参数-overlap-radius进行控制。默认值是50 bp。值越大，能合并距离越远的基因。

在上面的例子中，没有使用已有的基因模型。如果有这样信息存在，可以给Cufflinks的“-g”参数提供一个gtf文件作为导向。

为与将Cufflinks的输出结果和已有的基因模型进行比较，可以利用Cuffcompare程序。若参考基因模型文件名为“ref.gtf”,命令如下：

$ cuffcompare –r ref.gtf transcripts.gtf

输出文件包括总结和在两个文件基因模型中相似的基因信息。

2 Scripture

Scripture是利用Java编写的软件。可以从网页http:/ / www. broadinstitute. org/ software/ scripture/下载。Scripture基于剪切读段的信息来分选数据。基因组上由读段的剪切位点连接形成的岛域，可以用读段的双端的信息进一步连接起来。这些区域的异构体将被输出。

从构造连通图开始。以参考基因组序列的所有碱基做它的节点。如果在基因或转录本中有两个一致碱基相邻，就将这两个节点连接起来。含有剪切位点的读段给出了外显子和内含子的边界，且一个连接点必须被至少两个RNA-seq的Read支持。允许的供体/受体剪切位点包括典型的GT/AG和非典型的GC/AG 和AT/AC位经过数据统计。与作背景Reaad映射分布相比，对连接图中的通路的富集情况中进行数据统计。数据统计是通过以固定长度的窗口游览连接图和设定每个窗口的p值来实现。重要的窗口整合进创建的转录本图中，这些转录本图是利用双末端连接到之前的独立片段确定。

Scripture以经过排序的BAM文件和参考染色体的FASTA文件作输入。在写书阶段的Scripture最新版本2.0不能被获取，但初始版本可以通过作者获取。如法如下：

java –jar Scripture Version2.0.jar –task reconstruct \

-alignment top2/accepted_hits.bam –genome chr18.fa \

–out out –strand unstranded –chr 18

早期的Scripture版本输出两个文件，一个是包含基因模型的 BED文件，另外一个是包含转录本图的DOT文件，而2.0版本中输出四个文件。另外，新版中还创建了BAM文件同目录下的一个并列文件。四个输出文件如下：

-rw------- 1 somervuo 80K Jul 8 15:13 out.connected.bed

-rw------- 1 somervuo 250K Jul 8 14:09 out.paired Counts.txt

-rw------- 1 somervuo 229K Jul 8 14:09 out.paired Genes.bed

-rw------- 1 somervuo 104K Jul 8 14:09 out.scripture.paths.bed

“out.scripture.paths.bed” 文件中为只利用一个读段信息的初始转录本，“out.connected.bed”为采用了双末端读段信息的转录组。后面的文件中，包含504个基因产生的549个转录本。

2、 DE NOVO组装

在这,我们讲述了两种从头组装（de novo）构建全长转录本序列的软件，即无参考基因组序列辅助。两者都利用de Bruijn图。第一种软件由Velvet 和 Oases两个程序组成。Velvet是一个基因组中组装程序，其构建一个组装图将被第二个程序Oases利用，以查找能代表异构体的通路。另一软件Trinity有三个模块组成。首先，RNA-seq读段先被组装，再聚类，每个聚类代表一个基因组的位点。每个聚类构建一个de Bruijn图，抽提出线性转录本序列，以便同一位点可以产生多个异构体。两者都会在组装前将序列数据拷贝到同一个文件下，所以，如果处理大数据，请在组装前检查磁盘空间。

2.1 Velvet + Oases

Velvet由C语言编写，自定为基因组装配器[17]。后来，编写出以Velvet[9]输出结果进行转录本组装的另外一个程序Oases。Velvet可以从网页http:/ / www. ebi. ac. uk/ ~zerbino/velvet/下载，Oases 可从http:/ / www. ebi. ac. uk/ ~zerbino/oases/.下载。二者软件包中都有写好的使用手册。

Velvet 由两个程序组成：velveth 和 velvetg。前者计算数据的 k-mer数，后者寻找并抽取de Bruijn的叠连群。Oases分割de Bruijn图，并从每个位点抽取异构体。由 Oases得到的转录本序列通常比Velvet得到叠连群长很多。为在Velvet中充分利用双末端读段，必须对其进行隔行读取，也就是说，一个文件中的一对读段的两条序列的定位是相邻的。如果一对读段最初是被分别放在不同的文件中，Velvet包中的Perl 脚本会将其进行连接处理。该命令会创建一个新的文件 “chr18_12.fq”，再在该文件进行隔行读取。

$ shuffle Sequences_fastq.pl chr18_1.fq chr18_2.fq chr18_12.fq

第一个步是创建一个哈希列表：这里我们设定k-mer的长度为 25 bp，设定输出目录为“vdir”，同时设定数据格式。本例中，双末端读段为FASTQ 格式。遍历de Bruijn 图和抽取重叠群在第二步完成。这里我们设定的插入片段长度为200 bp。Velvet中，插入序列长度即为文库片段长，包含了两端读长。将“–read_trkg” 设定 “yes” 非常重要，这样Oases才能利用读段跟踪信息。

$ velveth vdir 25 –fastq –short Paired chr18_12.fq

$ velvetg vdir –ins_length 200 –read_trkg yes

Oases用于产生de Bruijn 图。 Oases的输入文件为Velvet输出文件所在路径。此时，双末端读段的插入长度也必须进行设定。本例中，设定最小的转录本长度为200 bp。

$ oases vdir –ins_length 200 –min_trans_lgth 200

输出目录vdir包含文件如下。转录本序列储存在FASTA文件“transcripts.fa” 中 。每条FASTA序列的名字记录描述了基因位点和异构体信息。由Oases创建的另一文件是“contig-ordering.txt”。

-rw------- 1 somervuo 25M Jul 16 11:56 Graph2

-rw------- 1 somervuo 11M Jul 16 11:59 Last Graph

-rw------- 1 somervuo 1.2K Jul 16 11:59 Log

-rw------- 1 somervuo 5.5M Jul 16 11:56 Pre Graph

-rw------- 1 somervuo 34M Jul 16 11:55 Roadmaps

-rw------- 1 somervuo 84M Jul 16 11:55 Sequences

-rw------- 1 somervuo 1.3M Jul 16 11:59 contig-ordering.txt

-rw------- 1 somervuo 2.6M Jul 16 11:56 contigs.fa

-rw------- 1 somervuo 253K Jul 16 11:59 stats.txt

-rw------- 1 somervuo 1.6M Jul 16 11:59 transcripts.fa

例如一条FASTA序列的名字记录为“Locus_10_Transcript_1/3_Confidence_0.571_Length_3815”。表明在Locus 10位点产生了3个转录本，该序列为其第一条转录本。置信值在0到1（越高越好），转录本长度单位为bp。该例中，“transcripts.fa” 文件中转录本最小长度为200 bp，由862个基因位点产生1308条转录本。

在Oases0.2版中，能以不同k-mer长度进行组装再进行合并。在Oases包中，一个Python脚本执行该功能。这里我们取19到29之间的奇数作为k-mer的值。另外， 给定Velvet 和 Oases 参数“–d” 和 “–p” 。利用该Python脚本，不需要使用“–read_trkg”参数。

$ python oases_pipeline.py -m 19 -M 29 -o odir -d " -fastq \

-short Paired chr18_12.fq" -p "-ins_length 200 \

–min_trans_lgth 200"

结果产生对应不同k-mer值的输出文件，以及包含合并组装结果的目录”odirMerged” 。本例中，在odirMerged目录中的“transcripts.fa” 文件包含由827基因位点产生的4468条转录本序列。在创建不同k-mer值的输出目录后，可以无需从开头起始就能合并一些组装结果。这由Python脚本的”-r”参数完成。例如，仅对k >= 25的组装结果进行合并，命令如下：

$ python oases_pipeline.py -m 25 -M 29 –r -o odir

结果生成由783个基因位点产生的2159转录本序列。

在Velvet 中k-mer值默认的最大参数为31；但也可以使用更大的数值。例如，将k值设定到51，可以由以下命令实现：

$ make ‘MAXKMERLENGTH=51’

注：Oases也必须进行同样的设定

当利用不同的k值进行组装时，最好使用“–m” 和 “–M”参数，而非设定k值独立进行组装。因为Velvet拷贝读段数据到序列文件，但是设定“–m” 和 “–M”参数时，它仅拷贝一次，其他目录会创建第一次输出结果目录的字符连接。

欢迎关注生信人

TCGA | 小工具 | 数据库 |组装| 注释 |   基因家族  |  Pvalue

基因预测  |bestorf |  sci | NAR | 在线工具 | 生存分析 | 热图

 生信不死 | 初学者 | circRNA | 一箭画心| 十二生肖 | circos

 舞台|基因组 | 黄金测序 | 套路 | 杂谈组装 |  进化 | 测序简史