一种植物线粒体基因组编码环状rna的识别算法
技术领域：
：1.本发明涉及生物信息
技术领域：
：，具体涉及一种植物线粒体基因组编码环状rna的识别算法。
背景技术：
：：2.核糖核酸(rna)是一种由脱氧核糖核酸(dna)转录而成的核糖核酸聚合物，包含信使rna(messagerna，mrna)、转运rna(transferrna，trna)、核糖体rna(ribosomalrna，rrna)和小核rna(smallnuclearrna，snrna)等种类。rna既有线性分子，又存在以共价键形式头尾相连的环状结构分子，即环状rna(circularrna，circrna)。circrna分子的封闭环状结构，使得其比线性rna更稳定，不易被核糖核酸外切酶降解。转录组测序(rna-sequencing，rna-seq)技术是高通量鉴定circrna的基础。rna-seq可以分为单端测序及双端测序两种类型，单端测序由一端读取rna序列，而双端测序则从rna的两端读取序列，并在后续通过两端的重叠序列将双端测序序列结合在一起。双端测序可以获得更长的读长(read)。3.从高通量测序结果中识别circrna是其功能研究的关键步骤，这一步需要适合的算法辅助。已有的研究表明，内含子反向剪接是产生circrna的主要途径。现在常用的circrna识别算法包括ciri21.、ciri2.、find_circ3.、circexplorer4.和uroborus5.等。这些算法首先通过不同原理识别反向剪接产物，然后根据细胞核内含子剪接位点的保守序列或者基因注释信息等进行过滤，从而获得可行度较高的候选circrna。4.植物线粒体基因组与细胞核基因组存在差异，主要包括以下几个部分：(1)植物线粒体基因组内含子多为自我剪接(self-splicing)的ii类内含子，尚无报道显示这类内含子存在反向剪接途径；(2)相对于细胞核基因组，植物线粒体基因组很小(比如：玉米b73的线粒体基因组大小约570kb，而其细胞核基因组大小约2.18gb)；(3)植物线粒体基因组的编码区存在大量rna编辑位点，而细胞核基因组不存在。这些差异使得植物线粒体基因组编码环状rna(circrna)具有与细胞核基因组编码circrna具有不同的特点。目前常用的算法多是针对细胞核基因组编码circrna设计，它们无法有效识别植物线粒体基因组编码circrna。5.以ciri2为例。ciri2是目前应用广泛的circrna识别算法之一。它使用基于多片段匹配(multipleseedmatching)的最大似然估计(maximumlikelihoodestimation)来识别含有反向剪接衔接点(back-splicedjunction，bsj)的读长，并通过计算错误发现率(falsediscoveryrate，fdr)降低由于重复序列和碱基错配所导致的假阳性。ciri2在识别植物线粒体基因组编码circrna时，存在多个问题，包括：(1)ciri2通过识别内含子反向剪接产物来筛选circrna，而植物线粒体基因组的ii类内含子未发现反向剪接现象存在；(2)植物线粒体基因组存在大量rna编辑位点，这导致ciri2预测植物线粒体基因组编码circrna时，fdr过高，进而导致大量真实存在的circrna被过滤；(3)植物线粒体基因组包含大量重复序列，ciri2会过滤源自这些重复区域的circrna。6.参考文献：7.[1]gaoy,zhangj,zhaof.circularrnaidentificationbasedonmultipleseedmatching[j].briefingsinbioinformatics.2018,19(5):803-810.[0008][2]gaoy,wangj,zhaof.ciri:anefficientandunbiasedalgorithmfordenovocircularrnaidentification[j].genomebiology.2015,16(1).[0009][3]memczaks,jensm,elefsiniotia,etal.circularrnasarealargeclassofanimalrnaswithregulatorypotency[j].nature.2013,495(7441):333-338.[0010][4]max,xuew,chenl,etal.circexplorerpipelinesforcircrnaannotationandquantificationfromnon-polyadenylatedrna-seqdatasets[j].methods.2021,196:3-10.[0011][5]songx,zhangn,hanp,etal.circularrnaprofileingliomasrevealedbyidentificationtooluroborus[j].nucleicacidsresearch.2016,44(9):e87.技术实现要素：[0012]本发明所要解决的技术问题是如何识别植物线粒体基因组编码环状rna(circrna)。[0013]为了解决上述技术问题，本发明首先提供了识别植物线粒体基因组编码circrna的方法。所述方法包括如下步骤：使用blast工具将rna-seq数据比对到参考基因组，得到比对结果文件；将所述比对结果文件中的数据进行序列特征筛选获得候选circrna文件；对所述候选circrna文件进行条件筛选获得植物线粒体基因组编码circrna。[0014]上文所述方法中，所述序列特征筛选可为去除所述比对结果文件中的只有一个序列片段比对到所述参考基因组上的读长，保留有且只有两个序列片段能方向相反比对到所述参考基因组的读长。[0015]上文所述方法中，所述序列特征筛选可包括如下步骤：[0016]a1)去除所述比对结果文件中的只有一个序列片段比对到所述参考基因组上的读长，保留两个以上序列片段能够方向相反比对到参考基因组的读长的比对结果文件。[0017]a2)提取所述比对结果文件中所述两个以上序列片段的位置信息数据和所述序列片段对应于参考基因组上的正负链信息数据。[0018]a3)对所述正负链信息数据进行处理，保留有且仅有两个序列片段方向相反的比对到参考基因组的读长，对所述序列片段按照所述位置信息进行排序整合，得到候选circrna。[0019]上文所述方法中，所述位置信息可为所述序列片段对应的读长的id信息。所述位置信息还可包括所述序列片段在所述读长上的位置信息、及其所述序列片段比对到参考基因组序列上的位置信息。所述方向可为核苷酸序列5’到3’的方向。[0020]上文所述方法中，所述序列片段的核苷酸个数大于等于11。[0021]上文所述方法中，a3)对所述正负链信息进行处理可包括如下步骤：将每个所述序列片段的所述正负链信息数据分别填入正链文件与负链文件两个文件；将所述正链文件中的所述序列片段进行整合得到正链整合文件，将所述负链文件中的所述序列片段的起始和终止位置信息进行交换，然后进行整合得到负链整合文件。[0022]上文所述方法中，所述blast工具的“e”参数设置范围可为10-5～5。[0023]所述“e”参数的设置具体可为2。[0024]上文所述方法中，所述筛选条件可为如下：[0025]b1)环化位点的重叠核苷酸与空缺核苷酸的核苷酸个数应小于等于3。[0026]b2)候选circrna长度小于等于10000个核苷酸，并且大于等于对应读长的长度。[0027]所述环化位点可为所述两段方向相反比对到参考基因组上的序列片段在读长上的衔接点。所述重叠核苷酸可为所述两个序列片段方向相反比对到参考基因组上的位置有重叠核苷酸。所述空缺核苷酸可为所述两个方向相反比对到参考基因组上的序列片段衔接处存在的不能比对到参考基因组上的核苷酸。[0028]上文所述方法中，所述候选circrna长度可为所述环化位点对应于比对到参考基因组上的起始核苷酸与终止核苷酸之间的核苷酸个数。所述起始核苷酸可对应于比对到参考基因组上的5’端核苷酸。所述终止核苷酸可对应于比对到参考基因组上3’端核苷酸。[0029]上文所述方法中，所述rna-seq数据可为fasta格式文件。所述fasta格式文件可为经过拆分的rna-seq数据文件。[0030]上文所述方法中，所述参考基因组的序列可为经过格式化的参考基因组的序列。所述格式化工具可为formatdb。[0031]为了解决上述技术问题，本发明还提供了识别植物线粒体基因组编码环状rna(circrna)的装置。所述装置可包括如下模块：[0032]c1)序列比对模块：用于使用blast工具将植物线粒体的转录组测序(rna-seq)数据比对到参考基因组，得到比对结果文件。[0033]c2)序列特征筛选模块：用于将所述比对结果文件中的数据进行序列特征筛选获得候选circrna文件。[0034]c3)条件筛选模块：用于对所述候选circrna文件进行条件筛选获得最终circrna。[0035]所述序列特征筛选可通过包括如下步骤的方法建立：为去除所述比对结果文件中的只有一个序列片段比对到所述参考基因组上的读长，只保留有两个序列片段能够方向相反比对到所述参考基因组的读长。[0036]所述blast工具的“e”参数设置范围可为10-5～5。[0037]所述“e”参数的设置具体可为2。[0038]上文所述装置中，所述条件筛选模块可通过包括如下步骤的方法建立：[0039]b1)环化位点的重叠核苷酸与空缺核苷酸的核苷酸个数应小于等于3。[0040]b2)候选circrna长度小于等于10000个核苷酸，并且大于等于对应读长的长度。[0041]所述环化位点可为所述两段方向相反比对到参考基因组上的序列片段在读长上的衔接点。所述重叠核苷酸可为所述两个序列片段方向相反比对到参考基因组上的位置重叠核苷酸。所述空缺核苷酸可为所述两个方向相反比对到参考基因组上的序列片段衔接处存在的不能比对到参考基因组上的核苷酸。[0042]上文所述装置中，所述候选circrna长度可为所述环化位点对应于比对到参考基因组上的起始核苷酸与终止核苷酸之间的核苷酸个数。所述起始核苷酸可对应于比对到参考基因组上的5’端核苷酸。所述终止核苷酸可对应于比对到参考基因组上3’端核苷酸。[0043]上文所述装置中，所述rna-seq数据可为fasta格式文件。所述fasta格式文件可为经过拆分的rna-seq数据文件。[0044]上文所述装置中，所述参考基因组的序列可为经过格式化的参考基因组的序列。所述格式化工具可为formatdb。[0045]为了解决上述技术问题，本发明还提供了一种存储有计算机程序的计算机可读存储介质。所述计算机程序使计算机建立如上文所述的方法的步骤或所述计算机程序使计算机建立如上文所述装置的模块。[0046]为了解决上述技术问题，本发明还提供了存储有计算机程序的计算机可读存储介质。所述计算机程序使计算机执行如上文所述的方法的步骤或所述计算机程序使计算机执行如上文所述装置的模块的步骤。[0047]针对已有方法的缺点，本发明旨在结合植物线粒体基因组特征，开发一款适合植物线粒体基因组编码circrna的识别算法，即mitochondrion-encodedcircrnaidentifier(简称meci)。该算法用于从植物rna-seq数据中鉴定高可信度的线粒体基因组编码circrna。[0048]ciri2通过识别内含子反向剪接产物来鉴定circrna，而植物线粒体基因组的内含子多数为ii类内含子，尚无报道显示这类内含子存在反向剪接途径。meci基于blast比对，检测rna-seq读长序列中的环化位点(circularizationjunction，即两段反向比对序列在读长上的衔接点；图1)。meci对circrna的识别不局限于内含子的反向剪接。植物线粒体基因组存在大量rna编辑位点，如果使用ciri2预测线粒体基因组编码circrna，会使得fdr过高，进而导致大量circrna被过滤。meci通过设置较大的blastne值(默认设置为2，可调节范围为10-5～5)，既保留长度较小(≥11nt)的比对片段，又保留存在错配的读长。植物线粒体基因组存在一定数量的重复序列，meci允许一个线粒体基因组编码circrna比对到基因组的不同位置，而ciri2会删除这些circrna。ciri2不允许环化位点(两段反向比对序列在读长上的衔接点)存在重叠或者空缺核苷酸的circrna，但许多植物线粒体基因组编码circrna的衔接处含有重叠核苷酸或者空缺核苷酸。meci允许circrna在环化位点处存在重叠与空缺核苷酸，并限制其数量不大于3。[0049]本发明所建立的植物线粒体基因组编码circrna的方法meci的创新点如下：[0050]1.meci基于blastn比对，并使用e值控制错误匹配和比对片段长度。[0051]2.meci对circrna的筛选依据植物线粒体基因组特征设置条件，并且不依赖于内含子反向剪接规则。具体条件如下：[0052]a.要求一个读长有且只含有两段反向片段比对到参考基因组；[0053]b.读长的两个反向比对片段衔接处出现不多于3个重叠或空缺核苷酸；[0054]c.e值默认设置为2，可调节范围为10-5～5，并且允许出现多个错配核苷酸；[0055]d.读长的两个反向比对片段的长度不小于11个核苷酸；[0056]e.读长两个反向比对片段之间的环化位点在参考基因组上的起始核苷酸与终止核苷酸之间的核苷酸个数大于该读长的核苷酸个数。[0057]与现有技术相比，本发明的有益效果：[0058]现在已有的circrna识别算法主要用于鉴定来自细胞核内含子反向剪接产物，而植物线粒体内含子多属于ii类内含子，不存在反向剪接途径。meci根据植物线粒体基因组的特点设计。首先，meci通过blast将rna-seq数据比对到参考基因组，并筛选只存在两段反向比对片段的读长；然后，通过一系列筛选条件，获得高可信度的读长。具体筛选条件如下：允许比对片段存在多个错配碱基，允许环化位点(两段反向比对序列在读长上的衔接点)处有不大于3个碱基的重叠或空缺，允许一个circrna可以比对到基因组多个位置，circrna长度(环化位点在参考基因组上的首位核苷酸与末位核苷酸之间的核苷酸个数)大于读长长度(读长的核苷酸个数)。[0059]本发明实施例中使用三种circrna识别算法，即meci、find_circ和ciri2，分别对ncbi下载的玉米和拟南芥rna-seq数据中鉴定线粒体基因组编码circrna。结果显示，在玉米3次重复数据中，meci分别鉴定到32686，23560和2638个线粒体基因组编码circrna，find_circ分别鉴定到1782，1197和1256个，而ciri2分别鉴定到263，219和184个。在拟南芥2次重复数据中，meci分别鉴定到49772和60289个线粒体基因组编码circrna，find_circ分别鉴定到2333和2782个，而ciri2分别鉴定到453和593个。通过设计发散引物进行rt-pcr扩增验证，结果显示，meci新检测的玉米和拟南芥的线粒体基因组编码circrna大部分是正确的，从而证明了meci预测植物线粒体基因组编码circrna具有较高的可信度。本发明建立的meci算法识别植物线粒体基因组编码circrna的能力优于已有的circrna识别方法，可以有效识别植物线粒体基因组编码circrna。附图说明[0060]图1为含有环化位点读长的确定和分类。如果1个读长包含两段方向相反比对到参考基因组的片段，则将其确定为候选circrna，两个片段在读长上的衔接点即为环化位点。根据环化位点处的区别，circrna可以分为三种类型。类型i：环化位点处含有重复比对到参考基因组的核苷酸(即重叠核苷酸)，类型ii：环化位点处含有不能比对到参考基因组的核苷酸(即空缺核苷酸)，类型iii：环化位点处不含重叠或者空缺核苷酸。[0061]图2为参数传入及检验流程图。[0062]图3为序列预处理流程图。[0063]图4为序列文件拆分流程图。[0064]图5为序列比对流程图。[0065]图6为blast比对数据解析流程图。[0066]图7为环状rna筛选流程图。[0067]图8为植物线粒体环状rna识别完整流程图。[0068]图9为meci，ciri2和find_circ共同或者独自从rna-seq数据中预测的高可信度线粒体基因组编码circrna的数量比较。[0069]图10为rt-pcr验证meci新预测的玉米线粒体基因组编码circrna。使用pcr验证源自玉米zmrrn26，zmcob，zmnad2t2和zmtrny的circrna。基因结构如图所示，粗线表示内含子(int)，成熟线性rna的5’和3’末端位置如箭头所指。黑色和蓝色直线分别代表meci和rt-pcr在该位点所检测circrna的位置。收敛引物和发散引物的位置如图所示。con：收敛引物，di：发散引物；gd：gdna，-和+：rnase–和rnase+cdna；箭头指向收敛引物所扩增目标片段；折线所指发散引物扩增条带被分别回收，克隆，转化和测序。m：dnamarker。[0070]图11为rt-pcr验证meci新预测的拟南芥线粒体基因组编码circrna。使用pcr验证源自拟南芥atatp1，atcox3，atnad4l，attrny和attrnd位点的circrna。基因结构如图所示，成熟线性rna的5’和3’末端位置如箭头所指。黑色和浅灰色直线分别代表meci和rt-pcr在该位点所检测circrna的位置。收敛引物和发散引物的位置如图所示。con：收敛引物，di：发散引物；gd：gdna，-和+：rnase–和rnase+cdna；箭头指向收敛引物所扩增目标片段；折线所指发散引物扩增条带被分别回收，克隆，转化和测序。m：dnamarker。具体实施方式[0071]下面结合具体实施方式对本发明进行进一步的详细描述，给出的实施例仅为了阐明本发明，而不是为了限制本发明的范围。以下提供的实施例可作为本
技术领域：
：普通技术人员进行进一步改进的指南，并不以任何方式构成对本发明的限制。[0072]下述实施例中的实验方法，如无特殊说明，均为常规方法，按照本领域内的文献所描述的技术或条件或者按照产品说明书进行。下述实施例中所用的材料、试剂等，如无特殊说明，均可从商业途径得到。[0073]实施例一、植物线粒体基因组编码circrna识别算法meci的建立[0074]1.数据输入及数据预处理[0075]1.1数据输入[0076]meci算法开始运算时，需要进行参数的传入及检验。第一步，检验植物线粒体rna-seq序列文件及参考基因组序列文件是否正确输入。第二步，检验各个选项参数(包括“‑in1”、“‑in2”以及“‑genome”等参数，表1)是否设定，未设定的参数填入默认设置参数。第三步，创建结果输出文件夹(流程如图2所示)。[0077]1.2数据预处理[0078]rna-seq分为单端测序和双端测序两种类型。因为两种类型测序文件的处理方式存在差异，所以需要对输入序列文件进行预处理。[0079]单端测序：将序列文件由fastq格式转化为fasta格式，进行后续处理。[0080]双端测序：首先，使用双端测序文件拼接程序flash(下载网址：https://jaist.dl.sourceforge.net/project/flashpage/flash-1.2.11.tar.gz)将双端测序的两个序列文件进行拼接，并得到一个拼接序列文件。然后，将拼接序列文件由fastq格式转化为fasta格式，用于后续处理(流程如图3所示)。[0081]表1.meci所需要的输入数据及可选选项[0082][0083]2.序列比对[0084]2.1序列文件拆分[0085]在进行blast比对之前，将fasta格式序列文件中的序列按照一定的大小拆分成多个文件，以提高后续序列比对效率。拆分文件的大小为测序文件读长数的0.1％～100％。[0086]计算步骤1.2fasta格式序列文件的读长数，并将读长数除以序列块大小参数(chunksize；一般设置为序列文件读长数的0.1％～100％)，从而获得需要拆分的序列块数量。创建拆分文件，设定文件名。编写并调用perl语言脚本，拆分fasta格式序列文件并填入之前准备好的拆分文件中，得到fasta格式的拆分文件(流程如图4所示)。[0087]2.2序列比对[0088]使用从ncbi获得的formatdb和blastall程序(两个程序的下载网址：https://ftp.ncbi.nlm.nih.gov/blast/executables/legacy.notsupported/2.2.26/blast-2.2.26-x64-linux.tar.gz)进行序列比对。使用formatdb程序将参考基因组序列格式化得到格式化的参考基因组序列，然后开始序列比对。[0089]创建多线程比对，计算fasta格式的拆分文件中每一条序列的碱基数，使用blastall将fasta格式的拆分文件中的序列与格式化的参考基因组序列进行批量序列比对，并输出比对结果文件(流程如图5所示)。其中，使用blastall进行比对过程中的e值默认设置为2，可选择范围为10-5～5。减小e值，则blastn比对片段(读长比对到参考基因组的序列片段)的错配率变小而长度增长；增大e值，则比对片段变短，错配率增加。读长比对片段的核苷酸个数大于等于11。[0090]表2.meci使用的程序[0091][0092][0093]3.序列特征筛选并获得候选环状rna[0094]此步骤对步骤2.2得到的比对结果文件进行处理，以筛选出符合circrna特征的比对结果。[0095]参考已有的报道，circrna识别算法主要通过读长的序列特征来鉴定circrna，即一个读长可以分为两段分别比对到参考基因组的序列，而且两段序列比对到参考基因组的位置方向相反。两段反向比对序列在读长上的衔接点即为circrna的环化位点(图1)。[0096]根据环化位点处的区别，circrna可以分为三种类型。类型i：环化位点处含有重复比对到参考基因组的核苷酸(即重叠核苷酸)，类型ii：环化位点处含有不能比对到参考基因组的核苷酸(即空缺核苷酸)，类型iii：环化位点处不含重叠核苷酸或者空缺核苷酸(图1)。[0097]在算法中实现，主要有以下几个步骤(流程如图6所示)：[0098]第一步，分析blast比对结果文件中的数据，删除只有一个序列片段能够比对到参考基因组的读长，即去除与参考基因组完全对齐的读长，仅保留候选读长(即有两个以上序列片段分别比对到参考基因组的读长)的比对结果文件。[0099]第二步，在包含候选读长的比对结果文件中，提取每个读长上能够比对到参考基因组的序列片段在读长上的位置信息及其比对到参考基因组序列上的位置信息，以及该读长的id信息数据；计算每个序列片段对应于参考基因组上的正负链信息，并将数据填入正链与负链两个文件。[0100]第三步，分别对正链文件与负链文件进行处理。对于正链文件，将其中的序列片段进行整合得到正链整合文件；对于负链文件，将其中的序列片段的起始和终止位置信息进行交换，然后进行整合得到负链整合文件。将正链整合文件和负链整合文件中的序列按读长的id排序(读长的id信息来源于原始rna-seq数据)，并保留有且仅有两个方向相反比对到参考基因组的序列片段的读长。将同一个读长上的两段序列片段整合在一起，获取候选circrna文件。[0101]4.对候选circrna进行条件筛选获得最终植物线粒体基因组编码circrna[0102]步骤3的候选circrna由序列特征筛选所获得，可能存在假阳性。为了获取高可信度的circrna，需要结合植物线粒体基因组特征，设置筛选条件对候选circrna进行过滤，从而获得最终植物线粒体基因组编码circrna(流程如图7所示)。筛选条件如下：[0103]第一，环化位点(两段方向相反比对到参考基因组上的序列片段在读长上的衔接点)的重叠核苷酸(即两个方向相反比对到参考基因组上的序列片段中间有重叠的几个核苷酸)与空缺核苷酸(即两个方向相反比对到参考基因组上的序列片段中间有不能比对到参考基因组上的核苷酸)的核苷酸个数应小于等于3(图1)。[0104]第二，circrna长度(环化位点对应于比对到参考基因组上的起始核苷酸与终止核苷酸之间的核苷酸个数)小于等于10000个核苷酸，并且大于等于对应读长的长度(序列核苷酸个数)。其中，起始核苷酸对应于比对到参考基因组上的5’端核苷酸；终止核苷酸对应于比对到参考基因组上3’端核苷酸。[0105]经上述条件筛选获得的植物线粒体基因组编码circrna填入输出表格，并生成两个结果文件：一个包含最终植物线粒体基因组编码circrna数据文件，命名为circrna_details.xls；另一个包含最终circrna对应测序读长数据的文件，命名为read_details.xls。[0106]本发明建立的植物线粒体基因组编码circrna识别方法meci的完整流程图如图8所示。[0107]实施例二、meci与现有技术的应用效果比较[0108]1.三种circrna识别算法预测植物线粒体基因组编码circrna的比较[0109]从ncbi下载玉米和拟南芥线粒体circrna的rna-seq数据(表3；两组数据的下载链接如下：www.ncbi.nlm.nih.gov/bioproject/prjna719584)。玉米和拟南芥rna-seq数据分别含有3次和2次生物学重复。[0110]表3.玉米和拟南芥线粒体circrna的转录组数据[0111][0112]使用三种circrna识别算法，即meci、find_circ和ciri2，分别从两组数据中鉴定玉米和拟南芥的线粒体基因组编码circrna。在玉米3次重复中，meci分别鉴定到32686，23560和2638个线粒体基因组编码circrna，find_circ分别鉴定到1782，1197和1256个，而ciri2分别鉴定到263，219和184个(表4)。在拟南芥2次重复中，meci分别鉴定到49772和60289个线粒体基因组编码circrna，find_circ分别鉴定到2333和2782个，而ciri2分别鉴定到453和593个。[0113]为了提高预测可靠性，取至少2次重复中都出现的circrna，作为高可信度线粒体基因组编码circrna。按照这个标准，meci，find_circ和ciri2在玉米中分别鉴定到7524，674和66个线粒体基因组编码circrna，而在拟南芥中分别鉴定到9819，685和105个。其中，37和42个线粒体基因组编码circrna为3种算法在玉米和拟南芥中分别共同预测(表4和图9)；7482和9749个由meci分别在玉米和拟南芥中新发现，而使用已有算法find_circ和ciri2未能识别获得的(图9)。[0114]表4.不同circrna识别方法在相同rna-seq数据中预测线粒体基因组编码circrna的数量比较[0115][0116]其中，find_circ算法的主要步骤如下：[0117](1)使用bowtie2，将rna-seq序列文件的每一个读长与参考基因组比对，丢弃与基因组完全匹配的读长；[0118](2)在剩余读长的5’及3’末端各取20碱基序列作为锚点(anchor)，分别与基因组比对。若这两个锚点能够分别比对到参考基因组，且在读长和参考基因组上的顺序相反，则扩展比对。[0119](3)将这个读长的全长序列与基因组比对。若比对结果包含步骤(2)中的两段比对序列，且环化位点(两段反向比对序列在读长上的衔接点)的侧翼序列为gt-ag保守位点，则认为这是一个候选circrna。[0120](4)候选circrna通过以下条件进一步筛选。第一，环化位点(两段反向比对序列在读长上的衔接点)应位置明确且只有一个；第二，每个circrna至少有两个独立的读长支持；第三，锚点应当仅能比对到基因组上的一个位置；第四，circrna长度(环化位点在参考基因组上的首位核苷酸与末位核苷酸之间的核苷酸个数)不大于100kb；第五，读长的两段序列与基因组比对时，最多允许两个错配碱基。[0121]ciri2算法的主要步骤如下：[0122](1)使用bwa-mem将rna-seq序列文件的每一个读长与参考基因组比对。当一个读长的部分片段(片段a)可以比对到基因组上，而侧翼的另一个片段(片段b)不能与参考基因组序列精确比对时，保留该读长用于下一步分析。[0123](2)将片段b分割为多个小片段，并将每个小片段比对到参考基因组上。当小片段能够比对到参考基因组上且与片段a在读长和参考基因组上的顺序相反时，即为反向剪接类型小片段。当小片段能够比对到参考基因组上且与片段a在读长和参考基因组上的顺序相同时，统计为正向剪接类型小片段。通过调节fdr阈值，可以控制小片段的数量、长度及其比对到参考基因组的错配率。[0124](3)分别统计正向剪接和反向剪接小片段数量，并通过最大似然估计法判断读长是否为反向剪接读长。[0125](4)计算反向剪接读长对应的circrna信息。在默认参数下，算法要求circrna长度(环化位点在参考基因组上的首位核苷酸与末位核苷酸之间的核苷酸个数)不大于200kb且环化位点(两段反向比对序列在读长上的衔接点)的侧翼序列存在gt-ag保守位点。[0126]2.结果验证[0127]2.1实验验证方法[0128]2.1.1主要试剂及其配方或者来源：[0129](1)dna提取液(配方见表5)[0130]表5.dna提取液配方[0131]组分终浓度尿素(urea)7.0m氯化钠(nacl)0.3mtris-hcl(ph8.0)50mmedta(ph8.0)24mm月桂酰肌氨酸钠(sarkosyl)1％[0132](2)线粒体提取缓冲液(配方见表6)[0133]表6.线粒体提取液配方[0134]组分终浓度蔗糖0.3m焦磷酸钠(na4p2o7)5mm磷酸二氢钾(kh2po4)10mm聚乙烯吡咯烷酮(pvp)1％(w/v)edta2mm牛血清蛋白(bsa)1％(w/v)半胱氨酸(cysteine)5mm抗坏血酸20mm[0135]使用h3po4调节ph至7.3。[0136](3)线粒体洗涤缓冲液(配方见表7)[0137]表7.线粒体洗涤液配方[0138]组分终浓度蔗糖0.3megta1mmmops10mm[0139]使用h3po4调节ph至7.2。[0140](4)主要试剂和耗材[0141]dna聚合酶2xtaqmastermix购买于诺唯赞生物科技股份有限公司公司(中国南京)，rna和dna纯化回收试剂盒购买于天根生化科技有限公司(中国北京)，trizol试剂购买于美国英杰生命技术有限公司(invitrogen，美国)，primescripttmii逆转录酶购买于宝日医生物技术有限公司(takara，中国北京)，rnaser购买于epicenter公司(美国)，miracloth购买于calbiochem公司(美国)。[0142]2.1.2植物材料和种植环境[0143]玉米(zeamays)材料的遗传背景为w22自交系，拟南芥(arabidopsisthaliana)材料的生态型为哥伦比亚(columbia)。玉米种植在华南农业大学试验田，而拟南芥种植在人工气候室(22℃，黑暗培养)。[0144]2.1.3基因组(genomicdna，gdna)提取[0145]取大约0.2g玉米15dap(daysafterself-pollination，自交授粉天数)籽粒或者拟南芥暗培养5dag(daysaftergermination，萌发后5天)幼苗，放入研钵中，使用dna提取液提取玉米和拟南芥gdna。[0146]2.1.4玉米和拟南芥线粒体的富集[0147](1)取10克玉米15dap籽粒或者拟南芥5dap暗培养幼苗，收集于冰预冷的50ml离心管。[0148](2)向预冷的研钵中加入3～4ml线粒体提取缓冲液，将离心管中的组织全部倒入研钵，在冰上充分研磨。[0149](3)利用miracloth过滤研磨匀浆，收集于冰浴冷的50ml离心管，并使用10ml左右线粒体提取缓冲液冲洗miracloth上的残留组织。[0150](4)8000g，4℃，离心10min。[0151](4)转移上清液至新的50ml离心管中，20000g，4℃，离心10min。[0152](5)小心移除上清，沉淀即为线粒体。[0153](6)向沉淀中加入3ml线粒体洗涤缓冲液，轻摇混匀，分装于3个1.5ml去rna酶离心管。20000g，4℃，离心10min。[0154](7)弃上清，使用液氮将沉淀速冻，保存于-80℃，用于后续实验。[0155]2.1.5线粒体rna提取[0156](1)取步骤2.1.4得到的-80℃保存的线粒体沉淀1管，加入1mltrizol试剂，充分振荡混匀，冰上放置10min。[0157](2)加入200μl氯仿，充分混匀，冰上放置5min。[0158](3)12000rpm，离心10min。转移上清至新的1.5ml去rna酶离心管。[0159](4)加入等体积异丙醇，充分混匀，冰上放置30min。[0160](5)14000rpm，离心15min，去上清，管中白色沉淀即为rna。[0161](6)使用80％乙醇清洗白色沉淀两次，最后使用移液枪把剩余液体除净。[0162](7)室温放置5min，加入30μl去rna酶水。[0163](8)使用液氮将rna速冻，保存于-80℃冰箱，供后续使用。线粒体富集：[0164]2.1.6去除gdna污染[0165]使用dnasei去除线粒体rna中的gdna污染，反应体系见表8。[0166]表8.dnasei处理线粒体rna[0167]组分终浓度dnasei3μl10xdnasei缓冲液3μl线粒体rna20μg去rna酶水补水至30μl[0168]反应条件如下：37℃孵育30min，加3μledta，75℃孵育10min，灭活dnasei。反应液使用rna纯化试剂盒进行回收，用于后续rnaser处理。[0169]2.1.7线粒体rna去线性化[0170]使用rnaser处理线粒体rna，去除其中线性rna，反应体系见表9。[0171]表9.rnaser去线性化[0172]组分终浓度rnaser2μl10xrnaser缓冲液2μl线粒体rna20μg去rna酶水补水至40μl[0173]反应条件：37℃孵育30min。反应液使用rna纯化试剂盒进行回收，用于后续cdna合成。[0174]2.1.8cdna第一链合成[0175]使用primescripttmii逆转录酶和随机引物，以rnaser处理和未处理的线粒体rnamastermix试剂进行pcr扩增。pcr反应体系和条件分别见表12和13。[0185]表12.pcr反应体系[0186][0187][0188]表13.pcr反应条件[0189][0190]2.1.11pcr产物分离，回收和克隆[0191]将pcr产物经胶回收纯化后连接至pmd18-t载体转化大肠杆菌dh5α，挑取转化后生长出的单克隆菌进行pcr鉴定，选择含有预期片段大小的单克隆菌送至公司测序。[0192]2.1.12pcr扩增产物的环化位点确定[0193]使用ncbi的blastn功能(https://blast.ncbi.nlm.nih.gov)，将测序序列与参考基因组进行比对。根据比对结果，判断单克隆的插入片段是否含有环化衔接点。[0194]2.2pcr扩增和测序结果分析[0195]为了验证步骤1中meci新预测的线粒体基因组编码circrna是否准确，随机选取来自玉米和拟南芥线粒体基因组的11个位点(玉米：zmrrn26，zmcob，zmnad2t2和zmtrny；拟南芥：atatp1，atcox3，atnad4l，atrpl5，attrny和attrnd)，使用rt-pcr扩增由这些位点编码的circrna(图10和11)。根据meci预测结果，线粒体基因组位点一般可以产生多个circrna异构体(图10和11，表14)。比如：来源于玉米zmnad2t2的circrna异构体数量为322个，而来自拟南芥atatp1的异构体数量多达982个(表17)。此外，大多数线粒体基因组编码circrna异构体环化衔接点的5’末端比较接近，它们的主要区别在于3’末端位置不同(图10和11)。[0196]表14.rt-pcr单克隆的测序结果分析[0197][0198]pcr扩增产物的琼脂糖凝胶电泳结果显示，每对收敛引物均可以从gdna，rnase-cdna和rnase+cdna中扩增单一片段，并且长度与预测大小一致。相对而言，发散引物只能从rnase-和rnase+cdna模板中扩增circrna，而且扩增片段成弥散状，这一特点与上述circrna异构体特征一致，即来自同一位点的circrna异构体具有相似5’末端和不同3’末端。不同发散引物扩增片段经过pcr产物回收，载体克隆，大肠杆菌感受态转化，阳性克隆鉴定等步骤，送至公司进行测序。单克隆测序和序列比对的结果显示，rt-pcr所鉴定的226个circrna中，有144个被meci算法所预测，验证比例为63.72％。相对而言，这些rt-pcr鉴定的circrna，被ciri2和find_circ算法预测的数量为0。这些结果说明，meci新检测的线粒体基因组编码circrna大部分是正确的，从而证明了meci预测植物线粒体基因组编码circrna具有较高的可信度。[0199]综上，本发明建立的meci算法识别植物线粒体基因组编码circrna的能力优于已有的circrna识别方法。[0200]以上对本发明进行了详述。对于本领域技术人员来说，在不脱离本发明的宗旨和范围，以及无需进行不必要的实验情况下，可在等同参数、浓度和条件下，在较宽范围内实施本发明。虽然本发明给出了特殊的实施例，应该理解为，可以对本发明作进一步的改进。总之，按本发明的原理，本技术欲包括任何变更、用途或对本发明的改进，包括脱离了本技术中已公开范围，而用本领域已知的常规技术进行的改变。按以下附带的权利要求的范围，可以进行一些基本特征的应用。当前第1页12当前第1页12
技术特征：1.识别植物线粒体基因组编码环状rna的方法，其特征在于：所述方法包括如下步骤：使用blast工具将转录组测序数据比对到参考基因组，得到比对结果文件；将所述比对结果文件中的数据进行序列特征筛选获得候选circrna文件；对所述候选circrna文件进行条件筛选获得植物线粒体基因组编码环状rna；所述序列特征筛选为去除所述比对结果文件中的只有一个序列片段比对到所述参考基因组上的读长，保留有且只有两个序列片段能方向相反比对到所述参考基因组的读长。2.根据权利要求1所述的方法，其特征在于：所述序列特征筛选包括如下步骤：a1)去除所述比对结果文件中的只有一个序列片段比对到所述参考基因组上的读长，保留两个以上序列片段能够方向相反比对到参考基因组的读长的比对结果文件；a2)提取所述比对结果文件中所述两个以上序列片段的位置信息数据和所述序列片段对应于参考基因组上的正负链信息数据；a3)对所述正负链信息数据进行处理，保留有且仅有两个序列片段方向相反的比对到参考基因组的读长，对所述序列片段按照所述位置信息进行排序整合，得到候选circrna。3.根据权利要求1或2所述的方法，其特征在于：所述blast工具的“e”参数设置范围为10-5
～5。4.根据权利要求1-3中任一权利要求所述的方法，其特征在于：所述筛选条件如下：b1)环化位点的重叠核苷酸与空缺核苷酸的核苷酸个数应小于等于3；b2)候选circrna长度小于等于10000个核苷酸，并且大于等于对应读长的长度；所述环化位点为所述两段方向相反比对到参考基因组上的序列片段在读长上的衔接点；所述重叠核苷酸为所述两个序列片段方向相反比对到参考基因组上的位置有重叠核苷酸；所述空缺核苷酸为所述两个方向相反比对到参考基因组上的序列片段衔接处存在的不能比对到参考基因组上的核苷酸。5.识别植物线粒体基因组编码环状rna的装置，其特征在于：所述装置包括如下模块：c1)序列比对模块：用于使用blast工具将植物线粒体的转录组测序数据比对到参考基因组，得到比对结果文件；c2)序列特征筛选模块：用于将所述比对结果文件中的数据进行序列特征筛选获得候选circrna文件；c3)条件筛选模块：用于对所述候选circrna文件进行条件筛选获得最终circrna；所述序列特征筛选通过包括如下步骤的方法建立：为去除所述比对结果文件中的只有一个序列片段比对到所述参考基因组上的读长，只保留有两个序列片段能够方向相反比对到所述参考基因组的读长。所述blast工具的“e”参数设置范围为10-5
～5。6.根据权利要求5所述的装置，其特征在于：所述条件筛选模块通过包括如下步骤的方法建立：b1)环化位点的重叠核苷酸与空缺核苷酸的核苷酸个数应小于等于3；b2)候选circrna长度小于等于10000个核苷酸，并且大于等于对应读长的长度；所述环化位点为所述两段方向相反比对到参考基因组上的序列片段在读长上的衔接点；所述重叠核苷酸为所述两个序列片段方向相反比对到参考基因组上的位置重叠核苷酸；所述空缺核苷酸为所述两个方向相反比对到参考基因组上的序列片段衔接处存在的不
能比对到参考基因组上的核苷酸。7.一种存储有计算机程序的计算机可读存储介质，所述计算机程序使计算机建立如权利要求1-4中任一权利要求所述的方法的步骤或所述计算机程序使计算机建立如权利要求5-6中任一权利要求所述装置的模块。8.存储有计算机程序的计算机可读存储介质，所述计算机程序使计算机执行如权利要求1-4中任一权利要求所述的方法的步骤或所述计算机程序使计算机执行如权利要求5-6中任一权利要求所述装置的模块的步骤。
技术总结本发明公开了一种植物线粒体基因组编码环状RNA的识别算法。本发明所建立的MeCi识别方法根据植物线粒体基因组的特点设计，使用blast工具将RNA-seq数据比对到参考基因组，得到比对结果文件；将比对结果文件中的数据进行序列特征筛选获得候选circRNA；对候选circRNA进行条件筛选获得植物线粒体基因组编码环状RNA。通过使用MeCi对玉米和拟南芥的转录组数据分析可以获得现有技术检测不到的新预测线粒体基因组编码circRNA。实验证据表明，新预测的线粒体基因组编码circRNA具有较高的可信度，MeCi优于现有方法，可以有效识别植物线粒体基因组编码circRNA。体基因组编码circRNA。
技术研发人员：张亚锋 廖珣 常冯瑞
受保护的技术使用者：华南农业大学
技术研发日：2022.03.21
技术公布日：2022/7/5
