​Accurate quantification of circular RNAs identifies extensive circular isoform switching events

环状RNA的准确定量及识别广泛的环状RNA可变剪切事件

期刊：Nat. Commun.（IF=11.878）
发表单位：中国科学院北京生命科学研究院

导 读

本研究介绍了circRNA鉴定和定量的新算法CIRIquant。与其他同类算法相比，其在鉴定circRNA的过程中降低了假阳性率。在circRNA的定量和差异表达分析方面，CIRIquant也考虑到RNase R处理不均等因素对最终结果的影响，通过qPCR验证其结果相较于其他方法拥有更高的准确度。同时，CIRIquant还提供对circRNA与线性转录本比例的分析，为circRNA生物发生和调控的相关研究提供基础。最后利用肝癌癌与癌旁的样品进行分析，介绍了CIRIquant一般的分析流程。

摘 要

环状RNA（circRNA）的检测和定量面临若干重大挑战，包括较高的错误发现率，不均匀的rRNA消耗和RNase R处理效率，以及对反向剪接连接读数的低估。在这里，我们提出了一种新颖的算法CIRIquant，用于精确的circRNA定量和差异表达分析。我们进一步开发了实施两种独立措施的一站式差异表达分析流程，这有助于揭示circRNA及其线性对应物之间竞争性剪接的调控。我们将CIRIquant应用到肝细胞癌的RNA-seq数据集，并表征了线性循环切换和循环转录使用切换事件的两个重要组，这证明了探索肝脏肿瘤发生中广泛的转录组变化的有前途的能力。

前 言

circRNA具有在细胞代谢过程中发挥重要作用的巨大潜力，并且从高通量RNA-seq数据表征和定量circRNA已成为circRNA研究中的新问题。已经开发出许多用于表征circRNA的计算方法。这些方法大多数采用基于比对的策略来识别circRNA的反向剪接连接（BSJ）读取，并且在circRNA鉴定中灵敏度有限且显着的假阳性率高。因此，迫切需要可靠的计算工具来精确定量circRNA及其亲本线性转录本。circRNA在不同样品中的差异表达分析是circRNA研究中的常规分析。对于使用RiboMinus / RNase R处理过的RNA-seq文库进行的研究，RNase R处理步骤中富集系数的变化可能导致circRNA表达水平的估计偏差。

为了克服这些限制，我们建议使用CIRIquant对circRNA及其亲本转录本进行准确定量，并过滤BSJ读物的假阳性。CIRIquant可以使用当前广泛使用的工具（例如CIRI2，CIRCexplorer2，find_circ等）进行circRNA鉴定，然后为所鉴定的circRNA转录本生成伪参考序列，以重新排列推定的BSJ读段。利用针对参考基因组和伪circRNA转录物的序列比对结果，CIRIquant不仅可以更准确，更灵敏地鉴定BSJ序列，而且还可以可靠地定量circRNA的连接比率。而且，CIRIquant为一站式circRNA差异表达分析提供了便利的功能。我们应用CIRIquant来调查肝细胞癌（HCC）肿瘤样本和其相邻的正常组织之间的circRNA表达谱，这揭示了肝肿瘤发生中广泛的转录组学变化。我们相信CIRIquant提供了一种准确有效的定量方法来表征circRNA，并在校正实验和计算偏差后进行差异表达分析，这将大大提高我们对circRNA多样性和功能的理解。

结 果

1 定量circRNA表达的挑战

为了严格评估circRNA表达目前的挑战，我们从先前的六项研究中收集了63个转录组样本，涉及四种物种（人，小鼠，蝇和线虫）。使用HISAT2将所有这些RNA-seq数据集与其参考基因组和核糖体RNA（rRNA）序列进行比对，以评估定位速率和rRNA序列分数。来自四个物种的RNA-seq数据集显示rRNA序列消耗效率的差异非常高，这在很大程度上是由于目前RiboMinus转录组分离试剂盒的特异性和效率有限。对于每个样品，使用BWA-MEM算法将原始读数进一步与参考基因组比对，然后进行CIRI2进行circRNA检测和定量。我们选择CPM（每百万映射读取的计数）代表circRNA表达水平，以消除源自不同文库插入物大小和测序深度的偏倚。考虑到RNase R处理已广泛用于circRNA富集，我们比较了RNase R处理和未处理样品中基因和circRNA的表达水平，以研究RNase R处理是否会在转录定量中引入偏倚。考虑到RNase R处理有望提高circRNA检测的饱和度，因此鉴定出更多circRNA时，每个特定circRNA的环状读数越多，相对表达值越低。为了评估RNase R处理的效率，我们计算了不同RNA-seq数据集中circRNA的富集系数。富集系数在来自不同物种甚至同一物种内的样品之间表现出明显的差异。来自人类样品的数据集circases的RNase R富集效果最高，平均富集系数为0.95，其他物种的RNA-seq样品中RNase R系数的峰值通常低于0.9。为了进一步研究RNase R治疗中的随机效应，我们比较了实验重复中circRNA的富集系数。综上所述，这些发现表明，circRNA的定量受其丰度，rRNA去除率和RNase R处理效率的影响，这需要更有效的算法来解决这些问题**。（注：对包含4个物种的RNA-seq数据集进行常规的circRNA鉴定和定量分析发现，circRNA的定量受多个因素的影响，降低了可重复性）

2 CircCRIM1在体外促进NPC细胞迁移和侵袭

目前关于circRNA识别的计算方法主要基于BSJ读码的检测。但是，这些策略在从RNA-seq数据定量circRNA表达方面显示出明显的缺点，因为它们使用的RNA-seq对齐器并非设计用于通过BSJ信号对读段进行映射，特别是对于跨越多个连接位点的序列。在这里，我们提出了一种新的高效方法，用于从转录组数据中准确识别和量化BSJ上的线性和环状转录本。首先，使用HISAT2将RNA-seq读数与参考基因组比对，然后将CIRI2或其他circRNA检测工具应用于鉴定假定的circRNA。为了准确量化circRNA的表达水平并过滤假阳性BSJ，我们通过串联BSJ区的两个全长序列生成了一个伪circRNA参考序列，将候选circRNA读数与该伪参考重新对齐，并确定BSJ读数是否可以线性且完全与BSJ区域对齐。此外，通过结合针对参考基因组和伪参考序列的比对结果，我们可以计算跨BSJ的环状剪接连接读数的百分比来确定每个circRNA的连接比率。最后，使用规范的RNA-seq数据分析管道来获取转录本水平的表达信息。

对于RNase R处理过的RNA-seq数据，由于在不同研究中RNase R处理效率不均，因此无法直接将circRNA BSJ表达值用于比较分析。因此，我们实施了一个高斯混合模型以拟合其效率分布，然后将拟合模型用作RNase R系数校正的后验分布。对于circRNA的差异表达分析，我们提出了两种策略来评估病例和对照样品中circRNA的差异表达（DE）和差异剪接（DS）。当没有可用的生物学复制品时，我们利用倍数变化和方差信息来利用广义倍数变化计算circRNA的DE和DS得分。对于生物复制样品，进行统计测试以评估circRNA表达值和连接比率变化的重要性。为了推断样品之间circRNA表达的真正差异，我们使用基因表达数据实施对数倍数变化（TMM）归一化的修整均值，以消除系统的批次效应。因此，应用edgeR中的广义线性模型来确定circRNA在整个实验条件下是否显着差异表达，并且精确的比率测试用作circRNA连接比率差异的显着性检验。（注：介绍CIRIquant的分析流程，鉴定circRNA时降低假阳性率，差异定量分析时降低RNase R处理效率的影响）

3 模拟研究

为了评估circRNA量化中现有算法的性能，我们生成了用于性能比较的模拟数据集。我们首先模拟了100-bp至250-bp不同读长的RNA-seq读段，然后应用CIRI2，CIRCexplorer2，DCC，find_circ和KNIFE评估它们对circRNA检测的敏感性，对于除CIRI2以外的大多数工具，检测灵敏度随着测序读取长度的增加而降低。我们选择了在经验表达分布下具有配对末端100bp读数和circRNA的模拟数据。我们应用了这五种方法从模拟数据中检测circRNA，然后将预测的circRNA坐标用作CIRIquant的输入，以过滤假阳性并定量circRNA表达。我们计算了鉴定出的circRNA的预测BSJ读数数与模拟BSJ读数数之间的Pearson相关系数，在这些方法中，CIRI2达到了最佳性能（0.97），而其他工具的相关系数则在0.87至0.92之间。通过CIRIquant调整后，所有五种方法的相关系数均显着提高。我们通过Sailfish-cir或CIRIquant计算了模拟覆盖率和预测TPM之间的皮尔逊相关系数，在模拟数据集上，CIRIquant（0.947）的性能比Sailfish-cir（0.618）好得多。综上所述，这些结果表明，CIRIquant在模拟和实验验证的数据集上均可在circRNA定量分析中获得更好的性能。（注：利用CIRIquant对不同circRNA鉴定方法的分析结果进行筛选，可以提高不同方法之间的相关性，并且CIRIquant的定量分析准确度高于其他方法）

4 RNase R处理校正

circRNA文库制备中不同样品的RNase R消化效率不一致可能导致circRNA表达定量偏倚。为了评估CIRIquant中RNase R效率校正的性能，我们将该方法应用于先前研究中生成的三个人类细胞系RNA-seq数据集。我们实施了高斯混合模型以适应RNase R消化系数的分布，然后使用其后验分布来估计RNase R处理之前每个circRNA的原始读取计数。在RiboMinus / RNase R数据中，由于RNase R对circRNA的富集，circRNA的表达水平往往被高估了，而RNase R处理效率的校正可以极大地降低RNase R处理样品与未处理样品之间的偏差。为了通过实验验证RNase R校正的可靠性，对四个RiboMinus文库中的五个随机选择的circRNA进行了20次RT-PCR实验，经过RNase R校正后，CIRIquant显着降低了RNase R处理引起的偏差（RMSE = 0.006），并且比所有其他工具都具有更好的性能。（注：使用qPCR验证CIRIquant降低了RNase R处理的偏差）

为了进一步研究RNase R处理对差异表达分析的影响，我们在HeLa细胞中进行了RNA干扰介导的TRA2B的敲除，其中使用了设计为最小化偏离的shRNA耗尽了TRA2B（一种重要的序列特异性丝氨酸/精氨酸剪接因子）。我们在敲除TRA2B之前和之后分离了总RNA，然后将每个样品分为两个部分，其中一个部分进行了核糖体RNA消耗和RNase R消化，而另一部分仅通过核糖体RNA处理。我们分别使用CIRI2和CIRIquant鉴定circRNA并定量其表达水平，然后分别计算了敲除TRA2B后RiboMinus和RiboMinus / RNase R文库的circRNA表达水平的对数倍变化。我们观察到对于在两个数据集中均可检测到的circRNA，RiboMinus和RiboMinus / RNase R组之间的log2倍变化显示出较低的相关性，表明RNase R处理应在随后的分析中引起明显的偏倚。因此，校正RNase R的治疗效率对于circRNA差异表达分析至关重要。（注：通过敲低实验的差异表达分析对RNase R处理前后样品分析发现鉴定的差异circRNA相关性较低，说明校正RNase R的重要性）

我们采用了两种不同的策略来处理RiboMinus和RiboMinus / RNase R数据集。首先，对于未经RNase R处理的RiboMinus数据，CIRIquant使用从GFOLD调整的广义倍数变化方法计算差异表达（DE）分数，其中差异表达的circRNA通过考虑倍数变化和后验方差进行排序log2倍数变化的分布。其次，对于RiboMinus / RNase R数据，使用RiboMinus和RiboMinus / RNase R数据中发现的circRNA估算RNase R消化效率的分布。然后，使用高斯混合模型（GMM）拟合分布。利用拟合模型作为先验分布，可以从两个数据集中的RiboMinus / RNase R数据中的BSJ读取数和正向拼接结（FSJ）读取来推断RiboMinus数据中的原始BSJ读取计数。最后，将病例和对照样本校正后的读数计数的后验分布用于DE得分计算。与仅从RiboMinus / RNase R数据得出的DE得分相比，RNase R校正的其他步骤可以滤除大多数倍数变化和P值相对较低的circRNA。总体而言，结果表明在CIRIquant中实施的RNase R校正可以有效过滤假阳性并生成更可靠的差异表达分析。（注：根据是否结果RNase R的处理，使用CIRIquant对应的两种分析方法，验证CIRIquant校正RNase R偏差的作用）

5 线性-circRNA同工型转换事件的识别

为了确定circRNA与其亲本线性转录本的比率，CIRIquant利用了两步读取比对策略，包括标准的RNA-seq读取比对和对伪参考序列的重新比对。因此，可以准确确定在连接位点的给定circRNA的BSJ读数与规范连接读数的比率，利用伪参考对齐策略的CIRIquant在结点比率估计方面获得了显着的高性能（PCC = 0.97）。（注：研究circRNA与mRNA的比例是研究circRNA发生调控的重要方法，CIRIquant增加了伪参考序列的比对，获得了比其他方法更高的可信度）

为了进一步证明CIRIquant在差异表达分析中的适用性，我们使用专门设计的shRNA敲除了HeLa细胞中三个众所周知的剪接因子（MBNL1，PTBP1和TRA2B），然后分离了三个敲除样品和一个对照样品的总RNA模拟转录组测序。为通过实验验证CIRIquant估算circRNA连接率的准确性，对5个随机选择的circRNA及其亲本线性转录本进行了60次定量实时RT-PCR分析。对于每个circRNA，使用向外引物扩增BSJ区域，同时设计了两对靶向5'和3'irexon的向内引物，以确定这些circRNA的连接比例。从qRT-PCR和CIRIquant获得的circRNA的连接率显示出很高的一致性，这表明CIRIquant在确定circRNA连接率方面的可靠性。（注：通过qPCR方法验证了CIRIquant在三个敲出样品中计算的circRNA比例的准确性）基于差异表达分析，我们在三个敲低样品中分别发现了162、346和97个DE-circRNA。然后我们计算了DS分数，发现DS-circRNA和DE-circRNA显示出很高的一致性。然而，我们观察到了仅通过使用DS或DE分数鉴定的一定比例的circRNA，我们分别将其称为DS特异性和DE特异性circRNA。在三个数据集中，DE特异性和DS特异性circRNA的比例有所不同，这表明这些剪接因子的敲除可能通过不同的机制影响circRNA的生物发生。在不同的细胞类型中或响应刺激，基因可以大大改变不同同种型的相对丰度，通常称为同种型转换，这可能会对生物学产生重大影响。我们将线性-****circ****同工型转换定义为在环形转录本及其亲本线性转录本之间发生的交换。设置了严格的阈值以检测高度可信的LC转换事件，其中连接比率在不同样品中应显示出显着变化（Δ接合比> 0.3），而LC转换circRNA应从主要剪接同工型改变（接合比> 0.5）到次要异构体（接合比<0.5）。我们在这三个敲低样品中发现了多个线性-circ同工型转换事件，表明宿主的竞争性调节这些剪接因子敲除后，基因的典型剪接和circRNA的产生。我们提出了一种新的方法来研究circRNA研究中的线性-circ同种型转换事件，这有助于我们了解circRNA生物发生的机制，并揭示circRNA与线性对应物之间竞争性剪接的调控。（注：通过对三个不同剪接因子敲低样品中计算的差异剪接circRNA，可以推断这些circRNA的生物发生与对应的剪接因子有关，同时结合功能富集分析可以分析各个剪接因子调控的下游基因通路）

6 肝细胞癌的转录组广泛变化

为了系统地研究circRNA在肝细胞癌（HCC）中的表达情况，我们将CIRI2和CIRIquant应用于来自肝癌患者的正常和肿瘤肝样品的40个RNA-seq数据集，并在这些样品中检测并定量了46,984个circRNA。使用mRNA表达水平，circRNA表达水平和circRNA的结合比率这三个特征可以清楚地区分肿瘤和邻近的正常样品，这表明mRNA，circRNA和连接比也可以用作将肿瘤与正常组织样品区分开的潜在生物标记。（注：利用40组HCC癌与癌旁样品的RNA-seq数据，鉴定其中的circRNA）

我们进一步进行了差异表达分析，并将承载circRNA的基因分为三类：差异表达基因，DE-circRNA的宿主基因和DS-circRNA的宿主基因。在正常和肿瘤样品之间，总共有1159个基因发生了显着变化，而只有459个基因承载着DE-circRNA，而284个基因承载着DS-circRNA。发现几个特征明确的与癌症相关的基因（例如ZKSCAN1，ABCB4和CYP2C8）产生的circRNA在肿瘤样品中发生了显着变化，其中，circZKSCAN1和circSMARCA5在HCC肿瘤样本中均被下调。（注：鉴定癌与癌旁样品中差异表达的circRNA）

为了从不同方面进一步探讨circRNA表达模式的变化，我们计算了每个circRNA的连接率和circRNA宿主基因的环状转录使用率（CTU）。对于每个circRNA宿主基因，仅考虑两个主要的环状转录物，并且使用具有最长连接位点距离的circRNA比例来确定环状转录物的使用情况。我们研究了所有LC / CTU转换circRNA的组织特异性表达模式，通过circAtlas43在每个组织中circRNA的相对表达水平进行了测量。我们注意到hsa_intergenic_006404，一种来自5号染色体的基因间circRNA，被认为是LC转换circRNA，在肿瘤样品中表现出明显的连接率降低。通过实时定量PCR（qRT-PCR）和Sanger测序，我们确认了hsa_intergenic_006404在14对肿瘤和邻近正常样品中的表达水平，并验证了其BSJ序列。总之，hsa_intergenic_006404作为肝脏特异性的circRNA，可能在HCC的发生中起重要作用。（注：将差异表达的circRNA的连接率变化分类，筛选出表达量和连接率都有变化的hsa_intergenic_006404，并通过qPCR验证了其特异性表达）

图 肝细胞癌中差异表达的circRNAs。a 利用基因表达水平、circRNA表达值和连接比进行多维分析。样本之间的距离根据对数2倍变化计算。肿瘤和正常标本分别呈红色和灰色。b DE-circRNAs和DS-circRNAs基因与宿主基因的重叠。c 曼哈顿DE-和DScircRNAs的p-值图。水平虚线表示前25个最重要的DE-和DS-circRNA。DE-specific和DS-specific circRNAs分别以红色和蓝色显示。DE和DS组共享的圆圈显示为黑色。d 线性-circ（LC）切换和circ转录使用（CTU）切换循环的层次聚类。分别绘制了20对肝癌组织中这些环状rna的连接比和CTU图。e LC切换和CTU的组织特异性表达水平。相对表达水平代表特定组织中表达值与全部20个组织的比例。f LC切换电路实例，hsa-基因间。条形图显示了这种环状rna的肝脏特异性表达模式。g 蜂群图显示其在14对肿瘤组织和正常组织中的表达水平。绘制了上下游300nt区的H3K27Ac标志，揭示了hsa-间生区的潜在生烃机制。与大多数rbp相关的前两个CTU8。h CTU的RBP表达水平与CTU的相关性。热图中的颜色梯度表示spearman相关测试中的-log10 p值。

综上所述，通过同时使用circRNA连接率和环状转录本，CIRIquant可以在circRNA分析中提供新的视角，从而增强了我们探索肿瘤发生过程中广泛转录组变化的能力**。

讨 论

circRNA的准确定量是circRNA相关研究中一个具有挑战性的方面。尽管已经开发了几种用于circRNA检测的计算方法，但是，circRNA定量的精确度仍无法很好地表征。为了准确定量circRNA的表达，我们开发了CIRIquant，这是一种用于circRNA定量分析的多功能工具包。通过利用RNA序列比对仪将接头读段的拆分图谱转换为线性读段的剪接比对，CIRIquant在反向剪接的接头读段检测中显示出显着提高的灵敏度。在CIRIquant中，实现了BSJ读值计数和连接率的一般倍数变化，用于差异表达分析，这被证明是对有/无重复样品中差异表达的circRNA进行排序的有效方法。通过对模拟数据集和真实数据集的全面评估以及qRT-PCR验证，与以前的方法相比，CIRIquant在circRNA检测和定量方面显示出高效率并降低了假阳性率。

材料方法

连接读数的检测和定量

CIRIquant需要YAML格式的配置文件，其中包含FASTA格式的参考序列和GTF格式的注释，以及参考序列的BWA和HISAT2索引以进行读取比对。第一步，CIRIquant筛选出具有剪接比对的证据（最小映射片段长度> 5 bp），并将未映射的读数用作circRNA检测的候选物。默认情况下，CIRIquant使用CIRI2进行circRNA检测。但是，还支持通过其他任何检测工具生成的BED格式的手动拼接后的接合位点。通过将反向剪接区域的全长重复两次来生成循环序列的伪参考。然后，使用HISAT2将所有候选读段与伪参考序列比对，其中在连接位点的10 bp区域上一致地映射的读对被视为循环读段。此外，跨反向剪接连接位点排列的非圆形读段用于计算circRNA连接率。

circRNA的差异表达分析

对于没有生物学重复的样品，Fisher精确测试是在circ的BSJ读数的2×2表和总映射的读数上进行的，以确定在两种不同条件下circRNA的表达水平之间是否存在显着差异。R的比率测试包用于估算circRNA连接率的变化。对于具有生物学重复的RNA-seq数据集，在CIRIquant中实现了edgeR的改进版本，用于差异表达分析。简而言之，通过使用M值的均值（TMM）修整后的平均值来计算归一化因子，从而最大程度地减少样品之间对数倍变化，从而消除样品的特异性效应。这些因素用于circRNA丰度归一化。接下来，执行默认线性广义模型协议（GLM），以估计两组样本之间的差异的色散率和似然比测试。

长按识别二维码关注我们

image