扫码联系客服
6月24日的《自然·通讯》在线两篇文章，同时把注意力放在Pacbio的单分子测序技术Iso-Seq，分别用来做玉米和高粱的转录组研究。标题分别为“Unveiling the complexity of the maize transcriptome by single-molecule long-read sequencing”和”A survey of the sorghum transcriptome using single-molecule longreads”.今天笔者就来解读一下这两篇文章。
关键词：single-module long-read; Pacbio; Iso-Seq;Transcriptome;isoform
先来看第一篇文章：“Unveiling the complexity of the maize transcriptome by single-molecule long-read sequencing”.
【研究背景】
玉米是世界上最重要的农作物之一，也是一个重要的研究转录组调控网络的模式生物，对它的研究有助于人们从时空，发育，环境等多维度增加对转录组的认知。2009年首个玉米基因组参考基因组公布，并构建了遗传图谱和物理图谱。在基因组注释方面，采用了证据支持（cDNA，EST，RNA-seq等多种数据）和从头预测的方法，构建了一个非冗余的基因集。
随着NGS技术的发展，积累了越来越多的转录组数据，但是问题在于，NGS数据的读长短，·难以为每条转录本提供全长序列，限制了其在可变剪接方面研究应用。而且在某些案例中，短读长序列常常组装出低质量的转录本，导致基因组的注释不完整或者错误。所以目前的痛点是：二代测序数据的读长短，在转录本的可变剪接方面研究受限，准确度低；在基因组注释方面，注释不完整，可靠度差。正是这些限制。越来越多的研究者把注意力集中在长读长的第三代测序技术上，尤其是Pacbio公司的测序技术。
【研究现状】
在真核生物中，转录本的可变剪接是很普遍的。可变剪接大幅度地增加了蛋白，调控分子，细胞，生理和发育过程/通路的指令。在动物基因组，通过多外显子的可变剪接增加转录本的亚型（isoform）而不是大幅增加基因数目来完成生物功能复杂化的趋向。比如在果蝇中，DSCAM基因有超过38,000个isoform，两倍于整个基因集的数目。在人和拟南芥中，分别有~95%和61%的多外显子基因是可变剪接的。已有的研究中表明：可变剪接是生物体一种重要的调节机制。在果蝇中，Bcl-x 基因有两个isoform，一个抑制细胞凋亡，而另一个激活凋亡。在拟南芥中，ZIFL1基因有两个异构转录本，ZIFL1.1 与生长素的转运相关，而ZIFL1.3调节气孔的运动。可变剪接主要存在以下几种模式（如Figure 1所示）：
Figure1 可变剪接的几种模式
根据目前的研究表明：在动物中Exon skipping占主导地位，而在植物中Intron retention占主要地位。
【研究思路】
【样本选取】
选取玉米 Inbred line B73不同发育阶段的6个组织。分别是root,pollen,embryo,endosperm,immature ear及immature tassel组织。每个组织做三个生物学重复。
之所以选择多个组织做分析，是因为每个组织都会特异性表达转录本，只占全部基因的一部分。多组织，多发育时间段，可以尽可能获得更多基因的转录本，更全面的去了解一个物种的transcript landscape。
值得一提的是，作者在文章中也提到为什么没有选择叶片做分析。因为叶片的细胞中含有大量的叶绿素，用以光合作用，因此和光合作用相关的基因表达丰度很高。需要更高的测序深度，才能捕获到转录本的多样性。
【建库测序】
为了和二代short-read sequencing 做一个很好的对比。同时采用二代和三代建库测序。二代和三代样本选取一致。
二代测序采用Illumina 公司的Hiseq 2000， 采用PE101测序策略。
三代测序用Pacbio Iso-Seq策略，用多级大小片段建库测序。这是为了避免pacbio平台对短转录本的偏好性。
总共产生3,716,604条reads，过滤掉低质量的reads，总共获得643,330 条全长的，高质量的转录本序列(FL)。其中606,145 (94.2%)的FL用GMAP软件比对到玉米参考基因组上，获得111,151个非冗余的isoform。
【结果展示】
1. 在111,151个isoform中鉴别出829个转座元件。大部分是LTR序列，Gypsy (59%) ,Copia（33%）.在6个组织中，胚乳有最多的TE元件，雄穗最少。
2. 将剩下的110,322个isoform比对到Ref_v3基因集（FGS）上，确定出这些isoform的基因座，占注释基因座的~70%。为了进一步和参考基因集做比较，把这些isoform归为8类，详见Figure 4。
Figure4 Comparison of PacBio and RefGen_v3 isoforms
a) 其中3%的转录本为新发现的基因座的转录本；
b) 57%为新的isoform，即在参考基因组中可以找到，至少有一个共同的剪接位点（splice site）；
c) 17%的isoform和参考基因集为同一个基因模型；
d) 5.6%的isoform和参考基因有外显子有交集，但没有共同的剪接位点；
e) 0.5%的isoform 位于参考基因的内含子区；
f) 4%的isoform在参考基因的相反链上有外仙子区域的交集；
g) 3.9%的isoform是参考基因转录本的的部分区域；
h) 2.9%的isoform是包含好几个参考基因座位置。这些可能是当时基因集注释出现的错误，把一个基因注释成多个基因。
注意: 此处为了验证新的isoform，从Gramene数据库下载了玉米844条公认的split gene model序列。结果证明了新的isoform是没有问题的。同时也说明了V3基因集存在很多注释问题。
3. 长度统计表明，Iso-Seq预测出来的转录本整体上比V3基因集的要长，如图Figure 5。
Figure5 Comparison of isoform length between RefGen_v3 and PacBio data
在目前的V3基因集中一个基因平均有2.84个isoform，而Iso-Seq数据显示，一个基因平均有6.56个isoform，是前者的两倍多。同时数据也表明：基因的isoform数目和其intron数目呈现正相关，而isoform数目和基因表达量没有显著关联。
4. 数据饱和度分析
640+k 条FL已经接近可检测基因数目的饱和值。但是Pacbio的数据只检测到参考基因集中~70%的基因。目前很多表达数据显示，V3基因集中>6%的基因没有支持证据或者是假阳性预测。这说明，70%的基因检测率很可能是被低估了。为了进一步探究测序深度和转录本发现率之间的关系，对6个组织的数据做了Rarefaction分析，结果显示花粉比较容易达到检测饱和，而根组织需要更多的数据量才能达到饱和。大片段文库更容易触及检测饱和，文库片段越小，转录本多样性越高。
5. 组织特异性isoform和可变剪接模式分析
Isoform的组织特异性分析显示：花粉有更高的组织特异性，而根的特异性最低。5种可变剪接模式分析显示：IR（intron retention）在植物中占绝对优势。并随机选择了IR事件相关的基因做qPCR，结果都得到验证。而ES（Exonskipping）在植物中比较少见，和之前的研究比较吻合。剪接模式在6种组织中表现并不完全一致，其中胚乳组织中Mutually exclusive exons（外显子互斥，见Figure 1）占主导比例。
6. 转录因子预测
在玉米的V3参考基因集中，转录因子数目为2,624，分为57个家族。Iso-Seq数据在53个家族中发现了新的isoform，将转录因子的数据增加到5,423个。根据表达模式分为20个cluster。
7. LncRNA鉴定
鉴定出878个lncRNA,其中11个是以前报道过的，867个是新的lncRNA。这些lncRNA的平均长度为1.1kb，比之前的认知的lncRNA要长很多(平均400+kb)。基于在参考基因集上位置，这些lncRNA分为三类。A）58%在基因间区；B）21%来自反义链；C）16%来自正义链；D）5%来自内含子区域。另外，大部分lncRNA 是单外显子的。lncRNA也展示出组织特异性表达，较non-lncRNA更少数目的高表达。Lnc-RNA和编码基因在染色体上的分布比较相似，富集在着丝粒外围区，和低重复区及低CG/CHG甲基化区域相关。
8. 融合转录本鉴定
从Pacbio数据中鉴定出1,430个融合转录本。其中143个被illumina数据支持。结果表明，融合事件多发生在染色体间。随机挑选若干转录本并用RT-PCR和sanger 测序验证。
9. Pacbio isoform的验证和定量
主要利用illumina 数据来验证。用TopHat2.0.8b软件比对到参考基因组上，和Iso-Seq数据做比较。发现Iso-Seq数据中有86%以上的剪接点得到支持。考虑到Tophat2可检测的splice motif比较少，又用全能的START软件做比对，然后发现，有90%以上的剪接点得到illumina数据的支持。此时也表现出组织差异性，花粉和胚乳支持率较低,GC/AG模式在这两个组织中占比最高。用Cufflinks软件对转录本定量，计算出FPKM值。同时用6大组织及SRA数据库的转录组数据（叶，芽，顶端分生组织及幼苗）做表达量聚类分析，发现花粉和其他组织的表达模式差异较大。这是因为花粉高度特异化，仅包含三个细胞的缘故。
10. Pacbio isoform和short-read 组装比较
Short-read组装软件选择有参组装的Cufflinks和无参组装的Trinity。第一轮比较，考察Pacbio和short-read组装软件都能获得的基因座，评估转录本仅仅基于其剪接点。结果发现两者能组装出来的比例很低（Cufflinks: 22%; Trinity: 8%）。第二轮比较，考察这两个软件对Pacbio isoform 的还原能力。结果发现，随着isoform的复杂度增加，两个软件的敏感性急剧下降。
11. Isoform 甲基化分析
之前的报道认为：剪接位点的甲基化抑制了转录本的可变剪接。利用已发表的甲基化数据做分析，结果表明：CHG甲基化确实抑制可变剪接；CpG甲基化水平反而和isoform的数据正相关，甲基化促进可变剪接；CHH甲基化和可变剪接没有显著相关性。改变之前有限的认知。
进一步分析发现：lnc-RNA 和non-lncRNA基因的起始和终止位点甲基化水平较低。Non-lncRNA在基因主体部分表现出更高的CpG甲基化模式，lnc-RNA则是更高的CHG甲基化模式。两者甲基化模式的不同，也反映出两者基因表达水平的差异。
【重要结论】
1. 无论是NGS数据还是EST数据又或者FL cDNAs数据，都在一定程度上限制了基因注释，转录本精确重构。而long-read 技术则在规模，取样深度，转录本完整性和价格等需求上得到了多重满足，成为转录组研究，基因注释的利器；
2. Pacbio数据可以在很大程度上改善目前的玉米（或其他物种）基因组注释，校正之前的错误注释；
3. Pacbio isoform分析表明，长转录本是很普遍存在的，比之前认知的比例要高。LncRNA的预测分析也表明lncRNA的长度高于先前的认知。
4. Pacbio 的数据改变了此前对转录本甲基化的认知。阐明了甲基化在基因不同位置的特征；
5. 基因融合现象在玉米中比之前的认知要普遍的多。融合转录本的形成，增加了玉米转录组的复杂性。融合事件中，反式-剪接的机制出现的比例多。
6. 转录组在不同组织中有较大的差异，无论是转录本表达数目和表达量，还是可变剪接的模式，或者融合事件，甚至建库测序策略也会有所不同。
再来看第二篇文章 “A survey of the sorghum transcriptome using single-molecule long reads”.
【研究背景】
1. 转录后的mRNA前体，经常会进行可变剪接（AS）或者可变多聚腺苷酸化(APA)，增加了转录组的多样性。转录组的复杂性在决定基因的编码潜能和多机制调节基因表达方面起着重要的作用。目前在人和植物中分别发现~95%和60%以上的多外显子基因存在可变剪接。研究表明植物的非生物胁迫能显著影响转录本前体的可变剪接，可变剪接的基因会出现过量表达。
2. 目前转录本研究技术和工具的受限，对AS和APA知之甚少。而Pacbio 的单分子长读长技术为精确研究AS和APA事件提供了好的工具。
3. 高粱是世界上非常重要的C4作物，也是重要的耐胁迫的谷类，重要的研究非生物胁迫的模式生物。本研究以高粱为研究材料，采用Pacbio的Iso-Seq策略，并开发出一款好用的转录组异构体测序分析流程（TAPIS）,调研高粱转录组特征，改进高粱基因集注释。
【研究思路】
【样本选择】
选取高粱BTx623 品系的幼苗作为研究，设置干旱胁迫处理组和对照组，没有生物学重复。
【测序建库】
Pacbio 平台： PacBio RS II instrument. Iso-Seq 文库，1–2 and 2–6kb；
Illumina平台： SE88，处理和对照50M以上的reads数目，总计110M reads。
【结果展示】
1. 数据产出
Pacbio测序获得1,838,330个ROIs（reads of insert），其中884,638是全长的。RIO长度从20-3,886不等，平均长度为1,042。 由于Iso-Seq是环状cDNA模板多次测序，测序的错误率还是比较低的。和高粱参考基因比对，确定单碱基的错误率为2.34%。
2. TAPIS流程
TAPIS流程使用Python语言写成，从Iso-Seq数据校正，比对基因组到isoform鉴定做了一个相对完整的开发，如Figure 7所示。
Figure7 Schematic workflow of the transcriptome assembly and analysispipeline for Pacific Biosciences Isoform Sequencing reads.
TAPIS流程中对Pacbio数据的校正效果比较好。如Table1所示。
以illumina 数据为基准，无论用LoRDEC或者Proovread方法，最终能获取的reads 在~80%左右。而TAPIS-ref以参考基因组为基准来比对和校正，能获得95%左右的reads。然后用混合的方法测试，显示用illumina的数据做校正，然后用TAPIS-ref的方法校正，无论用LoRDEC或者Proovread，都可以获得~96%的reads，和单纯的TAPIS-ref相比提高不大。这表明TAPIS-ref方法校正效果很好，可以不需要用illumina数据来校正了。
3. 数据特征分析
从867,089条含polyA位点的reads中，总共鉴定出27,860个转录本，这些转录本分布如下：18%是没有内含子的；85.3%和参考基因集有交集；11.3%可能为新检测到的基因；3.4%的转录本是跨越两个基因模型的。可能是目前基因集的注释错误或者通读的转录本。总之，用Pacbio数据检测到14,550个基因的转录本，而用110M reads的illumina数据，用Cufflinks可以检测到17,750个转录本表达（FPKM>2.2），这意味着我们可以预测到表达基因的82%以上了。
4. AS事件分析
根据之前的报道，在植物中80%以上的编码基因和很多非编码RNA都存在内含子保留的剪接（玉米文章也有类似的结论）。在高粱中大约1,500个mRNA前体经历了AS事件。我们的研究发现，有10,053isoform经历了AS事件，其中只有2,950个isoform和现有的基因模型吻合。如Figure 8所示。Pacbio数据表明在高粱转录组AS事件比之前的认知要普遍很多。其中有超过7,000 个之前未报道的AS事件。在27,860个转录本中，11,342（40.7%）是新发现的，7,065 (25.4%)是全长isoform。统计表明大部分基因（69.9%）只有一个isoform，只有1%左右的基因有5个以上的isoform。为了验证本研究可靠性，特异挑选了6个在参考基因集中是单isoform，而在Pacbio数据中是多isoform的基因做RT-PCR验证，如Figure 9所示，Pacbio结果和电泳检测基本吻合。另外也发现在非生物胁迫下，isoform呈现出差异性（表达量或者isoform数目的变化）
Figure 8 Alternative splicing envets
5. 用illumina数据佐证Iso-Seq
为了验证Pacbio结果的可靠性，同时也和illumina的数据（用Cufflinks组装）做了比较分析，发现虽然有93.1%以上的Pacbio确定的可变剪接在cufflinks组装中找到，但是真正的只有59%的isoform被重构出来。由此可见，Cufflinks能检测到的isoform数目不到全部的2/3。Cufflinks再次被吊打。
6. APA事件研究
APA事件对转录本的修饰改造也是非常重要的。PolyA尾对mRNA的转运，定位，稳定，翻译都有重要影响。近来的高通量数据研究显示，APA事件通过在基因的编码区或者3-UTR区产生isoform，增加了转录组的复杂度。APA也影响了植物的生长发育，尤其是开花事件。由于Iso-Seq文库构建时采用Oligo dT 引物做cDNA的合成，因此所有的poly A都被呈现在reads中。所以Iso-Seq也是研究APA的利器。利用TAPIS流程，在可检测的14,450个基因中发现11,013个至少有一个支持的ploy A位点。平均有5个poly A reads比对到注释基因上，最终确定出38,736 polyA 位点。其中2,301(20.9%)的检测到ploy A 位点基因注释出起始密码子和终止密码子，有5,641个poly A 位点。针对这些位点做分布分析发现，96.4%的位点是在3-UTR区，3.5%的位点在编码区，还有0.1%的5-UTR区，与之前拟南芥的研究结果比较类似（Figure 10）。同时用EST数据评估了这些结果。发现EST中71.3%的poly A被Pacbio数据支持。大约7700个基因有2个或2个以上为poly A 位点，占表达基因数目的50%以上。用3-RACE方法对这些结果做随机验证，发现有较好的支持。同时，APA事件在处理组和对照组也表现出差异化，说明APA事件是一个受调控的现象。3-UTR区的长度在miRNA依赖的基因表达调控中具有重要的角色。虽然APA事件不能显著影响转录本的长度，但是影响到3-UTR的长度，而3-UTR长度影响了mRNA的稳定性。对poly A 位点的上下游50nt的碱基做偏好性分析，发现在poly A切割位点之后，显现出A突增， GC含量低的特征，与拟南芥的研究结果一致。在切割点上游25nt处有一个AAUAAA信号序列，在35nt处有一个UGUA信号序列。
Figure 10 拟南芥PolyA cluster 在基因体上的分布
7. 差异表达分析
虽然Iso-Seq产出的reads量只有常规的RNA-seq的1/20,但是长读长弥补了reads数目的劣势。把Iso-Seq的数据用GFOLD软件做差异表达分析，确定出186个差异表达基因（85上调，101下调），并用RT-qPCR做验证。（可惜没有深入分析差异表达的基因）。
8. 新基因发现
原来的参考基因集有~34,500个基因，Iso-Seq检测到其中的14,550个，占比42% (玉米的检测率大约为70%，大部分原因是玉米文章采用了多组织测序分析，基因表达的多样性会高很多，所以取样策略改进，可以检测/预测到更多的isoform及新基因)。并且鉴定出2,171个新基因。
9. Non-coding RNA 鉴定
总共鉴定出49个miRNA，20个是新发现的，14个有内含子，11个表现出可变剪接。用RT–PCR验证其中的3个miRNA，处理组样本的isoform均表现出较高的丰度。
【重要结论】
1. 采用Pabio Iso-Seq 方法，对高粱的转录组认知达到空前的高度。鉴定出27,860个转录本，其中1,342个是新发现的转录本isoform；10,053个转录本经历AS事件；7,700多个基因经历APA事件；鉴定出2,171个新基因；
2. Pacbio Iso-Seq 方法在转录本AS和APA鉴定方面表现出巨大的优势，本研究的数据为高粱基因集的改进注释提供资源；
3. TAPIS流程是一个好用的Iso-Seq分析工具，在没有illumina数据的情况下也能对三代数据做有效的校正。
4. 以后利用更多的组织样本及不同条件下的转录组研究，将会发现更多的AS及APA转录本，有助于高粱基因集的更好注释。
5. 当然，文库构建还是有改进空间的，可以更好的检测到低丰度的AS及APA事件，对研究植物的多种胁迫适应性至关重要。
【小结】
1. Pacbio Iso-Seq 策略为转录组研究提供了新的方法和视野。全长转录组的深入研究将会改写目前的很多认知。对基因集的注释也会有很大的帮助。
2. 第一篇文章的研究思路相对更好，分析更到位，即验证了之前的正确认知，也否定了以前的错误认知。更好的契合了标题中的Unveiling，而第二篇的可借鉴之处是，开发了较完整的流程，助力Iso-Seq的深入挖掘，对转录组本身的研究从非生物胁迫入手，更像是玩票，是一个抛砖引玉的Survey性质的工作。
3. 从这两篇文章的工作看，转录组远比目前的认知要复杂得多。未来转录组研究可以更加深入；
4. 目测不久就会有更多的物种加入全长转录组研究的俱乐部；当然，我相信基因组重新注释的工作也会纷纷上马了。当笔者刚看完这两篇文章，PNAS就在线了一篇用Pacbio平台研究玉米串联基因拷贝的文章“Analysis of tandem genecopies in maize chromosomal regions reconstructed from long sequence reads”。
5. “工欲善其事，必先利其器”。随着Pacbio为代表的三代单分子测序技术（已经开始预热四代了哈）的发展和改进，会把基因组、转录组等组学研究推向一个新的境界。
【参考文献】
1 Wang, B. et al. Unveiling the complexity of themaize transcriptome by single-molecule long-read sequencing. Nature communications 7, 11708, doi:10.1038/ncomms11708(2016).
2 Abdel-Ghany, S. E. et al. A survey of the sorghumtranscriptome using single-molecule long reads. Nature communications 7,11706, doi:10.1038/ncomms11706 (2016).
3 Wu, X. et al. Genome-wide landscape ofpolyadenylation in Arabidopsis provides evidence for extensive alternativepolyadenylation. Proceedings of theNational Academy of Sciences of the United States of America 108, 12533-12538,doi:10.1073/pnas.1019732108 (2011).
2016/6/27
原文转载自薛猫_柳叶刀的新浪博客http://blog.sina.com.cn/s/blog_4d43fb020102wjtn.html
欢迎关注生信人
