2018, Vol. 16
植物生长发育过程基因表达调控研究近年来随着二代测序的迅速发展得到广泛深入开展，RNA-Seq高通量测序给获得植物组织器官甚至单细胞全基因组转录信息带来了革命性突破[1-2]，随着测序成本的下降，测序数据量剧增，截至2017年8月，NCBI收录的SRA/RNA-Seq测序数量达到了65万条以上，其中人(Homo sapiens)与小鼠(Mus musculus)均已接近20万条，拟南芥1.35万条，玉米6.7千条，水稻3.8千条，很多数据目前可能并没有递交到NCBI SRA数据库。从高通量RNA测序结果中进行基因表达等分析近年来开发了很多工具[3-6]，包括没有参考基因组的序列从头组装工具[7-9]。很多工具算法复杂，完成多项分析任务，运行时间较长，研究者可能无法完全得知返回数据的运算细节，这对于从结果分析中发现基因表达的更多信息有一定难度[10-13]。
高通量测序序列较短而量大，利用序列本身直接从测序序列群体进行匹配会更有效找到比对序列，由于测序错误率引起匹配误差则可考虑不使用全长查询序列而使用局部序列进行匹配。为了加快RNA高通量测序获得基因表达量及进行从头组装，本文探讨利用开放源代码操作系统的文本过滤命令组合脚本进行RNA-Seq测序序列比对及拼接分析，分析过程清晰、简捷，不需要复杂的代码，可广泛用于进行快速RNA-Seq基因从头组装及表达量等分析。
1 材料与方法 1.1 植物材料山茶花品种克瑞默大牡丹和银凯旋的开花期花枝成熟全展叶片和未打开花苞的花瓣，采于佛山市植物园，由佛山市林业科学研究所提供。样品采集后马上液氮速冻并通过干冰运输。
1.2 RNA-Seq测序由北京奥维森基因科技有限公司提供测序服务，包括总RNA提取及建库, 双端测序(Paired-End, Illumina HiSeq 4 000)。序列格式fastq，提交6G高质量数据(Clean data)及7G未处理原始测序数据(Raw data)，每条序列长度150 mer，合并双端得到每个样品测序序列约5 000万条。
1.3 操作系统FreeBSD 11.0(amd64) (http://www.freebsd.org)，计算机内存8G, 硬盘1T，CPU为Intel(R) Xeon(R) CPU E3-1230 v5 @ 3.40GHz。
1.4 Unix文本处理命令包括awk, sed, cat, cut, comm, split, uniq, sort, tr等。所有序列处理与分析由FreeBSD操作系统的文本过滤命令及运行命令组合脚本完成。
1.5 序列处理提取序列(选用6G clean data), 仅保留序列，每条序列加上编号，然后逐级按照每10→100→5万条(“split”命令)进行随机排序并按随机方式合并序列。其中，每个100万分成若干个目录单独进行每1万条的随机排序，再随机组合合并所有序列。
1.6 表达量分析随机选取上述准随机排序后的100万条进行比对后统计表达量。1)比对：按照每10万条分割，选取其中任何一个10万条序列，每隔5个核苷酸取20个连续的核苷酸短序列(20 mer，共27个)，随机选取9个20 mer短序列进行与100万条的序列匹配比对，统计序列编号去除重复编号获得每条序列与100万条匹配后的数量。同时，获得9个20 mer序列的互补链亦进行与100万条序列匹配比对。2)计算表达量：分别计算两次(正链与互补链)比对的统计数即可初步获得每个查询序列与100万条序列比对的匹配序列，即匹配该条序列的正链与互补链的表达量(每100万条序列转录物数量)。
1.7 重叠群(基因)序列拼接将10万条序列与100万条序列比对得到的序列进行重叠群拼接，利用150 mer查询序列的5′和3′两端20 mer序列进行匹配，切除原始序列后找出剩下最长的序列进行连接。然后，将近10万条重叠群(部分可能没有获得)相互比对即可以去除重复的重叠群，即将重叠群序列再次进行9x20 mer随机片段选取并与其它重叠群序列比对，按照第一个序列号合并重叠群。并再次进行首尾20 mer拼接。
1.8 重叠群序列注释将组装拼接得到的序列与NCBI远程服务器核酸(nr/nt)数据库进行比对(https://blast.ncbi.nlm.nih.gov/Blast.cgi?PROGRAM=blastn&PAGE_TYPE=BlastSearch&LINK_LOC=blasthome)获得该序列的比对及注释信息。每次上传500条进行比对，保存比对结果进行注释提取。
2 结果 2.1 RNA-Seq测序基因表达分析流程RNA-Seq测序数据量大，选取较合适的序列数据大小来反映基因表达情况很必要，由于计算机运行内存的局限大多数程序都将序列数据进行切割，利用100万条进行比对及组装，获得重叠群拼接序列的表达量。RNA-Seq测序分析流程如图 1。
|
图 1 开源系统文本比对分析高通量RNA-Seq测序流程图 Figure 1 Flowchart for high-throughput RNA-Seq sequencing by text-filtering in open-source system |
山茶花叶片与花瓣RNA测序返回数据每个样品近5 000万条(双端测序合并，6G高质量数据, 即clean data)，进行每1万条准随机排序(见方法介绍)后如何选择适合的数据大小进行表达分析是很必要的，为此，分析过程进行了多个100万条及一个300万条的比对比较分析。另外，如果利用全部选取的100万条与100万条比对将耗尽计算机资源运行会非常慢甚至停止。因此，分析中从其中的100万条序列中选取了10万条叶片(银凯旋)、10万条花瓣(克瑞默)分别与多个不同的100万条序列进行比对，每个比对得到10万条比对结果，每条找到与其序列匹配的序列首先合并其序列号，后面分析再提取其序列进行拼接得到重叠序列。考虑到测序过程是cDNA正负链都可能被测序，所以也将10万条序列的互补链与100万条序列进行了比对。两次比对结果可以合并用于计算每个重叠序列的表达量。由于本文主要介绍文本过滤分析得到基因表达量，测序的质控分析结果在此不列出，下面分别介绍各个比对结果。
2.2.1 10万条序列与多个100万条序列比对选取多长的序列进行比对合适也是很关键的，如果选取全长150 mer序列进行比对计算机运行会很艰难，因此，考虑用20 mer的长度，要尽可能从150 mer长度的序列找到其他匹配序列，故考虑将150 mer隔5 mer分割将产生27个20 mer短序列，下面比较了24条序列(银凯旋叶片测序)用其全部的27个20 mer及分开的3个9个20 mer、随机的9个20 mer的比对结果，见表 1。
|
表 1 银凯旋叶片24条序列的20 mer与1个叶片100万条序列比对的匹配序列数量 Table 1 Matchig sequence amount of 20 mers from 24 sequences of Yangkaixuan blades aligned to 1 million sequences from 1 blade |
从表 1可以看出，用9个随机的20 mer短序列与27个所有的20 mer序列比对得到的比对匹配数量很一致(A, E~I)，而不进行随机选取的3个9个20 mer的匹配数量与27个20 mer的差异较明显(A, B~D)。因此，为了加快比对进程，后面的比对结果全部选取了150 mer序列中随机选取9个20 mer进行比对。
表 2用同样的银凯旋叶片测序的24条序列与银凯旋叶片测序的3个100万条(A/B/C)、1个300万条(D)和克瑞默的1个叶片100万条(G), 及2个银凯旋花瓣的100万条(F-双端读1/G-读2)进行比对，每个100万条是随机选取的，300万条序列里包含了A和B及另外的100万条(不在列表)。从表 2可以看到，3个银凯旋叶片的匹配数量也是很一致的，较高表达的24条序列与300万条的匹配数量也成倍增加，也说明基因表达并不是简单的线性增加。初步也可以看到两个品种的叶片表达数量是相似的，有些序列与花瓣则有明显的差异。从表 1和表 2结果可以看出，用100万条序列进行比对是可以反映高通量测序基因表达水平的。
|
|
表 2 24条序列9个随机20 mer与多个100万条及一个300万条序列比对的匹配序列数量 Table 2 Matching sequence amount of 9 sets of randomly-selected 20 mers from 24 sequences aligned to multiple sets of 1 million sequences and 1 set of 3 million sequences |
利用10万条序列与多个100万条序列比对得到的匹配结果能否反映基因的表达情况，需要进一步对其进行拼接后才能确定其表达量。另外，10万条序列是否包含了所有表达的基因序列亦需要进一步分析。利用银凯旋叶片10万条序列分别与多个100万条(两个品种及其叶片、花瓣)序列进行了比对，下面主要列出了其中部分结果。
银凯旋叶片10万条序列与包含自身序列的叶片100万条进行比对获得的9个比对结果如图 2，比对过程大约一周(每天约1.5万条，可同时运行6个相同任务计算机运行其他程序不受影响)，对比对结果进行排序去重复(剩下92 025条), 及去除没有匹配的(9 701条)，得到90 298条匹配(见图 3)，最高匹配序列6 615个, 大于1 000个的2 371个，小于1 000的87 927个。
|
图 2 银凯旋叶片10万条序列与1个100万条序列(包含自身10万条)部分序列比对匹配数 Figure 2 Partial matching amount of 100 000 sequences aligned to 1 set of 1 million sequences(including its own 100 000 sequences) of Yinkaixuan blade |
|
图 3 银凯旋叶片10万条序列与1个100万条序列比对结果排序去重复后匹配序列 Figure 3 Matching sequence after sorting and removing duplicated ones of 100 000 sequences aligned to one set of 1 million sequences of Yinkaixuan blade |
利用得到的90 298条比对结果的序列号从序列文件找出各自150 mer的序列，利用这些序列将其首尾20 mer序列作为匹配序列与找到的每一行内各自的匹配序列进行拼接，连接后得到加长的拼接序列，得到90 298条拼接序列，拼接过程运行时间较长约10天(按照每个窗口1万条分别进行拼接时间只需要3 d)。为了进一步合并拼接序列，去除了含有10个连续polyA/C/G/T的拼接序列，其中没有得到重叠群的序列有7 912条，序列超过150 mer的77 954条(重叠群序列见图 4), 拼接长度最长410 mer。
|
图 4 叶片10万条序列与100万条序列比对匹配序列第一次拼接部分序列 Figure 4 Partial sequences in first-time assembly of the blade's aligned matching sequences of 100 000 sequences and 1 million sequences |
将获得的重叠群序列(77 954条)进行再次相互比对(按照第一次比对的方法)，没有得到匹配的序列有16 478条，得到超过1条以上重叠群的有37 256条，合并第一个序列号后为16 754条，将这些序列选取第一个序列头尾20 mer进行拼接得到新的重叠群，长度最长为1 174 mer(部分序列截图见图 5), 长度超过500 mer的有9 729条。因为是按照第一个序列号组合的，新的16 574条重叠群仍然有匹配的序列群，所以需要进一步拼接，按照初次拼接的方法用第一个序列的150 mer进行20 mer分割随机选取9个去找到匹配的重叠群，再次拼接后得到6 125个具有一个以上的匹配群序列，10 629个不再出现匹配，而且最多匹配序列号只有218个。进一步可以对6 125条序列进行clustalw多重比对后可以获得单一的重叠群序列(本文没有列出)。
|
图 5 叶片10万条序列与100万条序列比对匹配序列第二次拼接部分序列 Figure 5 Partial sequences in second-time assembly of the blade's aligned matching sequences of 100 000 sequences and 1 million sequences |
上面的拼接中，得到的单一拼接序列包括再次拼接得到的16 478条及将含有2个以上的16 754条合并第一个序列号再比对匹配拼接后得到的10 629条，仍有6 125条序列含有2条及以上的匹配序列，如果进一步合并可以去除重复的匹配序列。因此，从10万条序列与100万条序列的比对匹配后拼接可以大致获得3万条左右的拼接序列。利用这些拼接序列的匹配序列按照拼接过程合并其表达数量，就得到它们的表达量计数。由于测序不排除正负链，所以计算表达量应该将正负链的表达进行计数。下面介绍得到的所有拼接序列的表达量。
第二次拼接得到的16 478条单一拼接序列正负链的表达量相当，最高3 369条，大于100的有53条，而大于1 000的只有16条(即千分之一)，小于等于10的最多有14 762条。限于篇幅，图 6列出了大于100的所有序列正负链表达量。
|
图 6 第二次拼接得到的单一序列表达量大于等于100的正负链表达量柱状图 Figure 6 Histogram of plus and minus strands with single sequence expression greater than or equal to 100 in second-time assembly |
第二次拼接得到的含有2条以上匹配序列的16 754条序列作为拼接序列合并其表达量，正负链也是相当的，而表达量普遍较高，大于1 000的有547条，大于等于10 000有20条，大于100的有3 040，小于等于100的有13 739条，小于等于10的4 204条。图 7也列出了50条序列的表达量。表达量较高是因为合并了多个匹配序列组，如图 7的第一个序列(正链表达量为25 520)拼接合并了86组首序列号相同的序列组，其中单组表达序列数大于500的就有58组(最高为758个)。
|
图 7 第二次拼接得到含有2个几以上匹配序列组合并后最高表达的50个序列的正负链表达量柱状图 Figure 7 Histogram of plus and minus strands expression of 50 highest expression sequences with more than 2 matching sequences after combination in second time assembly |
RNA-Seq测序过程首先对RNA进行了片段化，所以测序得到的序列是随机选取任何RNA片段获得的。由于一些序列比对没有找到匹配，部分序列进行第二次拼接没有找到匹配，因此，表达量的表示直接用第一次比对的150mer查询序列得到的匹配比对序列数作为该序列的表达量。
2.3 重叠群拼接序列注释初步分析单一拼接序列与含有2个以上拼接序列组的拼接序列其是否真实反映了某个基因，需要进一步利用该序列与NCBI核酸数据库进行blast比对，从返回比对序列及其注释可以初步推知。
每次用500条序列(fasta格式)与NCBI的blastn服务器进行比对，保存text格式的比对结果，提取得到了所有拼接序列的blast比对及其注释，完成所有3万条拼接序列的远程blalstn比对只需要半天。Blastn程序只是返回比对打分大于70%以上的比对，前面的两组拼接序列中再次拼接得到的16 478条单一序列blastn比对返回8 889条相似打分(Identities)高于70%的结果，而再次拼接大于多于两条拼接序列的16 754条序列进行blastn比对返回12 255条，说明较长的序列可以得到更多匹配的blastn比对。对从16 478条单一拼接序列进行blastn比对返回的8 889条注释进行第一条注释提取去重复后得到7 511条不同注释，其中一组53条为叶绿体全基因组序列注释，其他的均小于10条，出现2条及以上相同注释的有1 034条，只有1条注释的6 476条，注释中出现山茶属(Camellia)植物的有906条。多于两条拼接序列的16 754条序列进行blastn比对返回12 255条进行注释提取去除重复后得到8 418条注释，其中两组(162+36)为两个不同植物叶绿体全基因组序列注释，只有1条注释的6 462条，2条及以上的1 954条，出现山茶属(Camellia)植物的有1 962条。因此，两组拼接序列合并的有用基因注释只有15 678条，即是用10万条叶片测序得到的序列与100万条序列比对后拼接得到的3.3万条拼接序列得到的基因注释，没有得到有效注释的拼接序列为17 554条。表达量超过10 000的有20条，其中16条得到有效注释，结果见表 3，注释中仍然有相同基因名的。本文的注释分析是初步的(只是对blastn返回注释的第一个进行了提取)，而且其中出现Camellia的注释仍有重复，需要进一步分析。
|
|
表 3 最高表达的16条重叠拼接序列Genbank注释、blastn比对打分及表达量(每100万) Table 3 Genbank annotations, blastn aligned scores, and expression (per million) of 16 contig-joining sequences with highest expression |
本文利用开源系统文本过滤命令组合脚本对高通量RNA测序(RNA-Seq)的单次测序结果(6G高质量序列，约2 500万条)进行基因序列拼接及表达量分析，用10万条对任意100万条(按照每1万条进行了随机重排序)序列进行比对(利用随机9组20 mer长度寻找匹配序列)得到的匹配序列中拼接得到了重叠拼接序列，并获得了正负链的表达量。序列拼接长度最长达到1 174 bp, 10万条序列只有8千条没有得到匹配序列，拼接后得到单一拼接序列为16 478条，多于1条的有16 754条，为此这些序列再次进行了拼接，去除重复后得到10 629条单一拼接序列，而仍然有6 125条出现多于1条的匹配序列。
利用随机9组20 mer长度寻找匹配序列(测序长度150 mer，进行每隔5mer切割可以得到27组20mer长度序列)，得到的表达量与用27组20 mer序列进行比对得到的表达量很一致(见表 1)，因而为了加快运行速度只是选用随机9组20 mer进行比对就够了。利用叶片10万条随机9组20 mer对多个叶片或花瓣的100万条及300万条进行的比对也可以看到表达量很一致(见表 2)。从叶片10万条比对结果进行拼接得到的单一(16 478条)和多于1条(16 754条)的拼接序列进行了表达量分析，正负链的表达量一致(见图 6、图 7)。因此，从以上结果分析可以看到，对高通量测序进行拼接和表达量只随机选取100万条及分析正链就可以了。表示表达量的方法很多，如使用TPM (Tran-scripts per million)来表示高通量测序基因表达量[3]。本文对RNA-Seq测序的数据进行了随机打散(如图 1), 用每个150 mer查询序列对100万条序列进行比对并对比对匹配序列进行了重叠群序列拼接，第一次拼接可以获得最长410 mer序列，其表达量直接用查询序列进行比对获得的匹配序列数表示(见图 6)。我们进行了克瑞默的花瓣和叶片的各自10万条序列分别比对花瓣和叶片的3组100万条序列，对所有序列的匹配数进行了方差分析，差异不显著，比对匹配数高达5千的其标准差不超过125(结果未附)。因此，考虑到测序时RNA片段选取随机性及对测序序列进行了随机打散，表达量的表示用每个150 mer查询序列与单个100万条序列比对获得的比对匹配数是合适的。
3.2 开放源代码操作系统用于分析高通量基因表达量及从头拼接的优势开源系统(本文使用FreeeBSD系统)的文本过滤功能非常强大，组合命令行脚本可以完成很多分析工作。本文将测序序列作为文本进行处理，从上面结果可以看出使用文本过滤命令用于分析高通量测序序列完全可行的，不少序列分析平台也借助开源系统完成特定任务[6, 11]。
高通量测序由于数据量大，使用开源系统分析就具有明显优势。本文首先对测序序列按照每1万条进行随机排序合并(主要为split和sort命令)，从表达量分析结果可以看到随机排序是很有效的。比对选取全长测序序列(150 mer)显然会减慢比对进程，而较短序列也会得到太多非特异性匹配比对，而简单选取某些区段又会减少比对的匹配数量，为此，本文进行了每隔开5个碱基选取20 mer序列作为查询序列，150 mer测序序列可以得到27组20 mer序列组，全部用于比对对运算显然也不利，从表 1中看到，用27组20 mer与用随机9组20 mer进行比对得到的表达量是一致的，对运算速度影响不大，可以每次同时分析6组比对计算机运行速度不受影响，平均速度每天完成1.5万条比对(单组)，一周就可以同时得到6组10万条序列与100万条的比对。如果将每10万条比对100万条分成10个窗口只需要2天就可以。
本文也利用文本过滤命令组合脚本进行了从比对匹配序列中分析获得重叠序列并进行拼接，银凯旋叶片10万条比对匹配序列约9万多条很有多于1条匹配序列的，利用头尾20 mer进行对接也实现了重叠群序列的拼接，拼接长度最长1 174 bp(见图 4、图 5)。拼接过程比比对较耗时，约每天每个窗口完成7 000~1万条，完成9万条比对的拼接约需要10 d左右，得到77 954条拼接序列(先行去除含有连续10个以上的polyA/C/G/T)，7 912条没有得到拼接。第二次拼接用得到的拼接序列再次进行首尾20 mer拼接，得到单一的序列16 478条，多于1条的按照第一个序列号相同的进行合并后得到16 754条，16 754条序列组的序列利用150 mer测序序列按照随机9组20 mer再次进行比对及首尾拼接后得到10 629条单一序列，6 125条仍然含有多于1条以上的序列。初步对得到的拼接序列进行blastn比对(NCBI)获得了其注释(最高表达的16条序列的注释见表4)，超过半数拼接序列可以得到有效注释。
本地比对用blast2(本文使用的系统亦安装了NCBI的blast-2.2.26 -ia32-freebsd)也可以进行，利用RNA-Seq测序全长150 mer进行其比对速度仍然很快(结果未附)。由于blast对双链进行查找比对，不能对正负链的表达量进行分开，与本文利用文本比对得到的匹配转录本数量因而不能直接相比较，也不利于利用本文的拼接方法进行拼接。另外，部分序列利用blast无法比对而相同序列用文本比对可以得到匹配序列，且得到的拼接序列进行远程blast比对可以获得基因注释(结果未附)。进一步对两者得到的匹配序列是否一致，用银凯旋叶片10万条序列查找100万条叶片序列中9 701条未得到任何正链匹配的序列(结果2.2.2)，用blast去进行比对，发现其中6 247条亦找不到匹配序列，找到一个匹配序列的1 547条，只有2个序列多于100条(174和197)，其他匹配介于2~60条(结果未附)，这些匹配很可能来自同时blast查找互补链，符合正常情况下正链匹配少负链也是相对少的情形。上面的比较分析也进一步说明本文的文本比对得到的比对匹配转录数能够确切反映组织内基因的表达水平。
4 结论综合以上分析表明，本文建立的利用开源系统文本过滤命令组合脚本对高通量RNA-Seq测序进行序列比对及拼接，可以快速获得基因的表达信息。分析方法简捷、高效，用于进行RNA-Seq测序分析是完全可行的。
1) 本文进行的序列比对前提是序列在RNA-Seq测序中被重复测序，直接用序列本身进行文本匹配即可；
2) 本文通过将短序列进行每隔5个核苷酸取连续20个核苷酸序列并进行随机取其中9个20 mer与100万条目标序列比对，较好反映匹配目标相似序列数；
3) 本文方法可同时进行正链与其互补链分别与目标序列匹配查找，获得两者的表达量相当，表达量用每个查询序列与100万条序列比对匹配的序列数表示；
4) 本文利用查询序列两端20 mer进行重叠群组内拼接，进行两次拼接可获得大于1 Kb重叠序列，因此，进行转录组从头组装是完全可行的；
5) 根据从NCBI远程BLAST比对分析可提取重叠序列的注释，结合重叠群表达量可得到组织转录组信息。
| [1] |
GRIFFITH M, WALKER J R, SPIES N C, et al. Informatics for RNA Sequencing: A web resource for analysis on the cloud[J]. PLoS Computational Biology, 2015, 11(8): e1004393. DOI:10.1371/journal.pcbi.1004393 (  0) |
| [2] |
TSIRIGOS A, HAIMININ N, BILAL E, et al. GenomicTools: a computational platform for developing high-throughput analytics in genomics[J]. Bioinformatics, 2012, 28(2): 282-3. DOI:10.1093/bioinformatics/btr646 ( 0) |
| [3] |
CONESA A, MADRIGAL P, TARAZONA S, et al. A survey of best practices for RNA-Seq data analysis[J]. Genome Biology, 2016, 17(1): 181. DOI:10.1186/s13059-016-0881-8 ( 0) |
| [4] |
POPLAWSKI A, MARINI F, HESS M, et al. Systematically evaluating interfaces for RNA-Seq analysis from a life scientist perspective[J]. Briefings in Bioinformatics, 2016, 17(2): 213-223. DOI:10.1093/bib/bbv036 ( 0) |
| [5] |
TRAPNELL C, ROBERTS A, GOFF L, et al. Differential gene and transcript expression analysis of RNA-Seq experiments with TopHat and Cufflinks[J]. Nature Protocol, 2012, 7(3): 562-78. DOI:10.1038/nprot.2012.016 ( 0) |
| [6] |
GRNING BA, FALLMANN J, YUSUF D, et al. The RNA workbench: best practices for RNA and high-throughput sequencing bioinformatics in Galaxy[J]. Nucleic Acids Research, 2017, 45(Web Server issue): W560-W566. DOI:10.1093/nar/gkx409 ( 0) |
| [7] |
WANG S, GRIBSKOV M. Comprehensive evaluation of de novo transcriptome assembly programs and their effects on differential gene expression analysis[J]. Bioinformatics, 2017, 33(3): 327-333. DOI:10.1093/bioinformatics/btw625 ( 0) |
| [8] |
HONAAS L A, WAFULA E K, WICKETT N J, et al. Selecting superior de novo transcriptome assemblies: lessons learned by leveraging the best plant genome[J]. PLoS One, 2016, 11(1): e0146062. DOI:10.1371/journal.pone.0146062 ( 0) |
| [9] |
CABAU C, ESCUDIÉ F, DJARI A, et al. Compacting and correcting Trinity and Oases RNA-Seq de novo assemblies[J]. PeerJ, 2017, 5(7): e2988. DOI:10.7717/peerj.2988 ( 0) |
| [10] |
DAVIDSON N M, HAWKINS A D K, OSHLACK A, et al. SuperTranscripts: a data driven reference for analysis and visualisation of transcriptomes[J]. Genome Biology, 2017, 18(1): 148. DOI:10.1186/s13059-017-1284-1 ( 0) |
| [11] |
SREEDHARAN V T, SCHULTHEISS S J, JEAN G, et al. Oqtans: the RNA-seq workbench in the cloud for complete and reproducible quantitative transcriptome analysis[J]. Bioinformatics, 2014, 30(9): 1300-1. DOI:10.1093/bioinformatics/btt731 ( 0) |
| [12] |
ICAY K, CHEN P, CERVERA A, et al. SePIA: RNA and small RNA sequence processing, integration, and analysis[J]. BioData Mining, 2016, 9(1): 1-18. DOI:10.1186/s13040-016-0099-z ( 0) |
| [13] |
PERTEA M, KIM D, PERTEA G M, et al. Transcript-level expression analysis of RNA-Seq experiments with HISAT, StringTie and Ballgown[J]. Nature Protocol, 2016, 11(9): 1650-67. DOI:10.1038/nprot.2016.095 ( 0) |