CN114067907B

CN114067907B - 一种准确鉴定rna病毒基因组变异的方法

Info

Publication number: CN114067907B
Application number: CN202010758480.3A
Authority: CN
Inventors: 彭鑫鑫; 季序我
Original assignee: Precision Scientific Technology Beijing Co ltd; Predatum Biomedicine Suzhou Co ltd
Current assignee: Precision Scientific Technology Beijing Co ltd; Predatum Biomedicine Suzhou Co ltd
Priority date: 2020-07-31
Filing date: 2020-07-31
Publication date: 2022-07-08
Anticipated expiration: 2040-07-31
Also published as: CN114067907A

Abstract

本发明提供了一种准确鉴定RNA病毒基因组变异的方法。所述方法基于二代测序技术，通过改造clean reads以及提供新的分析策略和流程，最大限度的减轻NGS测序引入的错误，实现RNA病毒基因组变异的准确检测。

Description

一种准确鉴定RNA病毒基因组变异的方法

技术领域

本发明涉及病毒检测技术领域，特别是涉及一种准确鉴定RNA病毒基因组变异的方法。

背景技术

2020年1月7日，中国疾病预防控制中心(China CDC)从患者的咽拭子样本中鉴定出一种新型冠状病毒，国际病毒分类委员会冠状病毒研究小组(CSG)将新冠病毒正式命名为“SARS-CoV-2”。大多数SARS-CoV-2肺炎患者的症状较轻，预后良好，然而也有一部分患者出现严重的肺炎，肺水肿，ARDS或多器官功能衰竭和死亡。并且在短短2个月后，中国科研团队最新发现显示：新冠病毒已于近期产生了149个突变点。理论上，RNA病毒的基因组变异的主要来源为：1)RNA病毒的自身突变：RNA病毒具有比较快的变异速度(10^-6到10^-4替换每核苷酸位点每次感染)，显著高于DNA病毒(10^-8到10^-6个替换每核苷酸位点每细胞感染)；2)RNA病毒的编辑事件：人体细胞存在ADAR1/ADAR2介导的A-to-I的碱基编辑以及APOBEC介导的C-to-U的碱基编辑两套系统，当RNA病毒进入人体细胞后，同样面临被碱基A被ADAR1/ADAR2酶改为I(碱基I会被翻译机器与选择性剪切机器识别为G),以及碱基C被编辑为U的可能性；3)RNA病毒不同毒株的基因组差异。准确鉴定RNA病毒变异，是理解RNA病毒演化的基础，对研究病毒宿主来源、传播机制、遗传多样性等具有重要意义。

第二代测序技术(Next Generation Sequencing,NGS)被广泛用于RNA病毒的检测与基因组研究。例如在新型冠状病毒SARS-CoV-2的检测中，使用SARS-CoV-2患者的肺泡灌洗液，抽提RNA样本，进行宏转录组(metatranscriptome)测序，可以10⁷左右的有效读长(reads)。传统方法是，将这些reads比对到病毒参考基因组序列，将与病毒参考基因组序列不同的点视为单核苷酸突变位点。而这种方法会产生大量假阳性，该结果主要由于如下因素导致：1)NGS方法本身导致：目前NGS的错误率高达10^-5(A>C/T>G,C>A/G>T以及C>G/G>C)到10^-4(A>G/T>C)。使用NGS进行RNA病毒变异检测时，通常碱基数达到10⁹个，这样序列中将混入至少10⁴个测序错误。2)reads上不同位置的碱基错误引入概率差异很大，尤其是末端碱基出错概率远高于中间部分碱基。这些错误在上述方法中将被视为RNA病毒的基因组变异，从而导致出现大量假阳性结果。在早期的研究结果中即可发现，NGS技术产生的reads上，末端的位置贡献的变异事件高于中间位置的碱基。因此，现有技术包含以下缺点：1)忽略NGS高错率对RNA病毒基因组变异鉴定的影响。2)忽略read末端碱基错误率对RNA病毒基因组变异鉴定的影响。

发明内容

为了解决现有技术中存在的问题，本发明提供了一种能够准确鉴定RNA病毒基因组变异的方法，所述方法为基于NGS产生的RNA病毒宏转录组测序数据鉴定RNA病毒基因组变异，通过改造clean reads以及提供新的分析策略和流程，最大限度的减轻NGS测序引入的错误，实现RNA病毒基因组变异的准确检测。

更为具体的，所述方法包括如下技术步骤：

1、收集病毒样本并提取RNA；

2、基因文库的建立；

3、使用Illumina HiSeq 2500/4000进行NGS测序；

4、测序数据分析：

1)将来自宏转录组测序的fastq原始文件进行质控。使用软件fastp(https:// github.com/OpenGene/fastp)去除低质量序列，去除reads的接头，产生clean reads。

2)使用软件fastv(https://github.com/OpenGene/fastv)从clean reads中提取病毒reads。

3)将clean reads的5’端去除1-10个碱基，对3’端去除1-4个碱基。

4)使用BWA MEM(http://bio-bwa.sourceforge.net/)将病毒reads比对病毒基因组，产生bam文件。

5)使用samtools(http://www.htslib.org/)对上述bam文件产生mpileup文件。

6)对与参考基因组序列(来自NCBI公共数据库)相比，存在不同碱基的位点而言，如果n(A>G)>10×(n(A>C)+n(A>T))，则该位点被定义为：A>G的替换位点，其中所述n(A>G)为碱基A变为G的个数；n(A>C)为碱基A变为C的个数；n(A>T)为碱基A变为T的个数，其它类型依次类推。

7)对6)产生的位点清单，进行替换碱基类型read数目的过滤。要求有效位点需要有替换碱基类型read数目大于1，例如，含有变异的read数目可以为1，2，3，4等整数。

8)基于5)产生的mpileup文件计算碱基错配率。碱基错配率＝(与参考基因组不同的碱基数目)/(比对到参考基因组所有reads的碱基数据)。

9)对7)中明确定义的碱基替换位点，基于8)产生的整体碱基错配率，使用二项分布(binomial distribution)，计算观测该点碱基替换数目的概率。选取p<0.05的位点。例如，p值范围可以为<0.01，或者，p值范围可以为<0.001。此外，也可以使用FDR替代p值。

本发明采用上述技术方案，与现有技术中公开的传统方法相比，具有如下技术效果：显著提高变异鉴定的准确性，尤其是提高低频变异鉴定的准确性，以及显著降低假阳性率。这些优点能够确保准确甄别大量可靠的RNA编辑信号(包括A-to-I以及C-to-U的编辑)，以及RNA病毒基因组突变的信号。具有广阔的应用前景。

附图说明

附图用来提供对本发明的进一步理解，并且构成说明书的一部分，与本发明的实施例一起用于解释本发明，并不构成对本发明的限制。在附图中：

图1本发明分析流程；

图2实施例1检测结果：本发明所述方法相对于传统方法在鉴别nCov-RNA-2病毒RNA样品变异事件的差异对照(每组碱基柱状图中在前为传统方法，在后为本发明方法)；

图3实施例2检测结果：本发明所述方法相对于传统方法在鉴别nCov-RNA-2病毒RNA样品变异事件的差异对照(每组碱基柱状图中在前为传统方法，在后为本发明方法)。

具体实施方式

以下结合附图对本发明的优选实施例进行说明，应当理解，此处所描述的优选实施例仅用于说明和解释本发明，并不用于限定本发明。

实施例1

本实施例是采用第二代测序技术(Next Generation Sequencing,NGS)对从感染SARS-CoV-2患者1收集的支气管肺泡灌洗液(bronchoalveolar lavage fluid,BLAF)样本进行检测，具体操作步骤如下:

1、收集感染SARS-CoV-2患者(nCov-RNA-1)的全BALF样本，使用RNA提取试剂盒，提取RNA。

2、使用RNA-Seq文库准备试剂盒，进行逆转录与扩展，建立测序文库。

3、使用Illumina HiSeq 2500/4000进行NGS测序。

4、测序数据分析：

1)将来自宏转录组测序的fastq原始文件进行质控。使用软件fastp(https:// github.com/OpenGene/fastp)去除低质量序列，去除reads的接头，产生clean reads。测试数据来自https://trace.ncbi.nlm.nih.gov/Traces/sra？run＝SRR11059947。

2)使用软件fastv(https://github.com/OpenGene/fastv)从clean reads中提取病毒reads。测试数据获得病毒序列共57671380条。

3)将clean reads的5’端去除10个碱基，对3’端去除4个碱基。

4)使用BWA MEM(http://bio-bwa.sourceforge.net/)将病毒reads比对上病毒基因组。此处使用SARS-CoV-2参考基因组(NC_045512.2)，产生bam文件。

7)对6)产生的位点清单，进行替换碱基类型read数目的过滤。要求有效位点需要有替换碱基类型read数目大于1。

9)对7)中明确定义的碱基替换位点，基于8)产生的整体碱基错配率，使用二项分布(binomial distribution)，计算观测该点碱基替换数目的概率。选取p<0.05的位点。

检测结果如附图2所示。根据图2可明显看出，相对于传统方法(具体方法描述参见对比例1)，本发明所述方法在变异事件数目与变异位点数目检测上显著改善，大大降低变异鉴定的假阳性率。

实施例2

本实施例是采用第二代测序技术(Next Generation Sequencing,NGS)对从感染SARS-CoV-2患者2收集的支气管肺泡灌洗液(bronchoalveolar lavage fluid,BLAF)样本进行检测，具体操作步骤如下:

1、收集感染SARS-CoV-2患者(nCov-RNA-2)的全BALF样本，使用RNA提取试剂盒提取RNA。

2、使用RNA-Seq文库准备试剂盒进行逆转录与扩展，建立测序文库。

3、使用Illumina HiSeq 2500/4000进行NGS测序。

4、测序数据分析：

1)将来自宏转录组测序的fastq原始文件进行质控。使用软件fastp(https:// github.com/OpenGene/fastp)去除低质量序列，去除reads的接头，产生clean reads。测试数据来自https://trace.ncbi.nlm.nih.gov/Traces/sra？run＝SRR11059946。

2)使用软件fastv(https://github.com/OpenGene/fastv)从clean reads中提取病毒reads。测试数据获得病毒序列共4796888条。

3)将clean reads的5’端去除2个碱基，对3’端去除1个碱基。

4)使用BWA MEM(http://bio-bwa.sourceforge.net/)将病毒reads比对上病毒基因组。此处使用SARS-CoV-2参考基因组((NC_045512.2),产生bam文件。

6)基于上述mpileup文件，定义碱基替换类型。如果n(A>G)>10×(n(A>C)+n(A>T))，则该位点被定义为：A>G的替换位点，其中所述n(A>G)为碱基A变为G的个数；n(A>C)为碱基A变为C的个数；n(A>T)为碱基A变为T的个数，其它类型依次类推。

检测结果如附图3所示(作为对照采用的传统方法参见对比例2)。

对比例1

按照现有技术中传统检测方法检测实施例1中所述患者1样品中的SARS-CoV-2病毒基因组突变，具体步骤为：

步骤1-3同实施例1-2；

步骤4、测序数据分析：

2)将clean reads的使用软件fastv(https://github.com/OpenGene/fastv)提取病毒reads。测试数据获得病毒序列共57671380条。

3)使用BWA MEM(http://bio-bwa.sourceforge.net/)将病毒reads比对上病毒基因组。此处使用SARS-CoV-2参考基因组((NC_045512.2),产生bam文件。

4)使用samtools(http://www.htslib.org/)对上述bam文件产生mpileup文件。

5)统计四种碱基(A、G、T、C)分别发生替换事件的总数目与发生替换的位点数目。

对比例2

按照现有技术中传统检测方法检测实施例2中所述患者2样品中的SARS-CoV-2病毒基因组突变，具体方法同对比例1。

显然，本领域的技术人员可以对本发明进行各种改动和变型而不脱离本发明的精神和范围。这样，倘若本发明的这些修改和变型属于本发明权利要求及其等同技术的范围之内，则本发明也意图包含这些改动和变型在内。

Claims

1.一种准确鉴定RNA病毒基因组变异的方法，其特征在于，所述方法包括如下步骤：

1) 提取病毒样品中的总RNA；

2) 建立病毒总RNA文库；

3) 高通量测序；

4）将3）中得到的宏转录组测序fastq原始文件进行质控，使用软件fastp去除低质量序列，去除reads的接头，产生clean reads；

5) 使用软件fastv从clean reads中提取病毒reads；

6) 将clean reads的5’端去除1-10个碱基，对3’端去除1-4个碱基；

7) 使用BWA MEM将病毒reads比对病毒基因组，产生bam文件；

8) 使用samtools对步骤7）中所述bam文件产生mpileup文件；

9) 基于步骤8）中所述mpileup文件，定义碱基替换类型，对与参考基因组序列相比，存在不同碱基的位点而言，如果n(A>G) > 10×(n(A>C) + n(A>T))，则该位点被定义为：A>G的替换位点，其中所述n(A>G)为碱基A变为G的个数；n(A>C)为碱基A变为C的个数；n(A>T)为碱基A变为T的个数，

其它类型依次类推；

10）对步骤9）产生的位点清单，进行替换碱基类型read数目的过滤，要求有效位点需要有替换碱基类型read数目大于1且为整数；

11）基于步骤8）产生的mpileup文件计算碱基错配率，所述碱基错配率=与参考基因组不同的碱基数目/比对到参考基因组所有reads的碱基数目；

12）对步骤9）中明确定义的碱基替换位点，基于步骤11）产生的整体碱基错配率，使用二项分布法计算观测该点碱基替换数目的概率p，选取p<0.05的位点。

2.如权利要求1所述的方法，其特征在于，步骤3）所述的高通量测序为第二代高通量测序技术。

3.如权利要求1所述的方法，步骤6）为将clean reads的5’端去除10个碱基，对3’端去除4个碱基。

4.如权利要求1所述的方法，步骤10）中所述要求有效位点需要有替换碱基类型read数目为1、2、3或4。

5.如权利要求1所述的方法，步骤12）中选取的p值范围为小于0.01。

6.如权利要求1所述的方法，步骤12）中选取的p值范围为小于0.001。

7.如权利要求1所述的方法，步骤12）中使用FDR统计法替代二项分布法计算观测该点碱基替换数目的概率。