一种区分RNA测序数据中基因表达差异与长拷贝数变异的方法

摘要

本发明公开了一种区分RNA测序数据中基因表达差异与长拷贝数变异的方法，包括以下步骤：先用已知不含CNV的样本和对照样本分别提取RNA，建库测序、质控、比对，得到比对结果BAM文件；统计原始reads count并进行标准化处理，得到矫正后reads count文件，过滤掉部分对照样本后合并得到原始reads count矩阵文件；使用差异表达检测软件检测测试样本中的差异表达基因，获得每个外显子的坐标区间，确定判断阈值，再用待测试样本替换已知不含CNV的样本重复操作，根据阈值来判断为拷贝数变异或表达差异。该方法判断准确，操作方便。

说明书

技术领域

本发明涉及生物信息学与精准医学全基因组变异检测技术领域，具体涉及一种区分RNA测序数据中基因表达差异与长拷贝数变异的方法。

背景技术

拷贝数变异(Copy number variation,CNV)是由基因组发生重排而导致的，一般指长度为1kb以上的基因组大片段的拷贝数增加或者减少，主要表现为亚显微水平的缺失和重复，是人类疾病的重要致病因素之一。异常的拷贝数变化(CNV)是许多人类疾病(如癌症、遗传性疾病、心血管疾病)的一种重要分子机制。作为疾病的一项生物标志，染色体水平的缺失、扩增等变化已成为许多疾病研究的热点，然而传统的方法(比如G显带，FISH，CGH等)存在操作繁琐、分辨率低等问题，难以提供变异区段的具体信息。

得益于高通量测序(又称下一代测序，Next Generation Sequencing,NGS)技术，在碱基层面分析拷贝数变异已经成为可能。基于NGS平台的变异检测方案众多，如全基因组测序(WGS)、全外显子测序(WES)、捕获测序(Targeted sequencing)、RNA测序等。目前各种测序方案各有优劣，比如WGS虽然能分析全基因组上的各种变异，但是市场价格比较昂贵，且对数据计算和存储要求较高，因此没有大规模在临床上使用。WES由于其捕获了基因组上的重要编码基因和一些与疾病相关的特殊区域，价格较低，数据量少，分析周期短，故使用场景比较普遍，但是WES只能在基因层面研究相关变异，一旦涉及到转录组水平就无计可施了。因此，越来越多的学术期刊开始力推RNA测序技术，该方案能够真正从转录组水平研究个体的变异情况。RNA测序可以全面快速地获取某一物种特定器官或组织在某一状态下的几乎所有转录本，这样就可以用来研究某一生理条件或时空条件下该物种转录水平和正常情况下的差异。这些测序结果相对于基因组来说，数据更少，更方便分析，此外转录组测序最大的好处是费用比较低。相对于高通量测序来说，费用要少很多。

目前基于RNA测序的研究方案更多是集中在可变剪接、融合基因、基因表达上面，几乎很少有学术期刊报道在RNA中检测CNV的案例。RNA数据检测CNV有其技术难点，比如检测到的CNV区域如何判定是表达差异导致的还是真实CNV，如果仅仅从RNA数据本身是很难做到有效区分的。

发明内容

本发明的目的在于，提供一种基于二代测序平台，使用RNA测序数据来检测长拷贝数变异(CNV)的方法，能够有效区分检测到的长CNV是由基因差异表达所导致，还是长拷贝数变异所致。

目前的检测方法无法区分真实CNV和表达差异导致的假CNV，其原因在于基因表达差异的检测结果和拷贝数变异的检测结果都会呈现出reads覆盖度上的异常。本发明可以在只利用RNA数据的基础上利用创新性的生物信息学算法区分拷贝数变异和基因表达差异，其原理是在合理的生物学规律下不同个体间发生异常表达的基因数量或者基因组区域是小范围的，如果较多的基因或者大片段的基因组区域发生表达量差异，则更有可能是拷贝数变异导致的。基于此本发明在对常规生物信息学软件检测到差异基因后通过注释更多信息去识别是否是长拷贝数变异(本发明中指6Mb以上的拷贝数变异)，充分利用RNA数据。

本发明技术方案详述如下：

一种区分RNA测序数据中基因表达差异与长拷贝数变异的方法，包括以下步骤：

(1)选择已知不含CNV的样本作为阈值确定样本；

选择表型正常，即无测试样本致病表型、与测试样本组织类型一致的样本作为对照样本，对照样本数量在30个以上；

(2)阈值确定样本和对照样本分别提取RNA，建库测序，得到测序数据文件，即FASTQ文件；

对测序数据文件进行质控、比对，得到比对结果文件，即BAM文件；

利用比对结果文件和基因组注释文件，统计所有对照样本的每个exon区间内的原始reads count，并进行标准化处理，得到矫正后reads count文件，计算所有对照样本之间矫正后reads count文件的相关性系数，过滤掉相关性系数R＜0.7的对照样本；合并阈值确定样本和过滤后对照样本的原始reads count文件，得到原始reads count矩阵文件；

以原始reads count矩阵文件为输入，使用差异表达检测软件检测阈值确定样本中的差异表达基因，通过基因组注释文件获得该差异表达基因的每个外显子的坐标区间，然后以外显子为单位遍历整条染色体，把外显子差异表达状态相同、且基因与基因之间物理距离小于1Mb的基因合并为1个差异表达区，同时要求每个差异表达区内的95％以上的外显子的差异表达状态相同，找到阈值确定样本中最长的差异表达区间记为nMb，以及上调和下调的有意外显子最小占比值记为m％；

(3)将步骤(2)中的阈值确定样本替换为待测试样本，重复步骤(2)的操作，所得最长差异表达区间长度在nMb以上，并且上调和下调的有意外显子占比在m％以上，则确定待测试样本中该差异表达区域实际为拷贝数变异，其余的单个或多个相邻或不相邻的基因为表达差异；

所述差异表达状态，是指与对照样本相比，阈值确定样本和待测试样本该基因的表达量的上调或者下调的状态。

reads count：读段数量，是指利用二代测序技术对基因或者转录本进行测序，测序测到的每条序列为一个读段，即read，通过统计某一区域测到的reads数目即readscount。

可选或优选的，上述方法中，对原始reads count进行标准化处理的方法如下：

(1)利用基因组注释文件提取每个基因对应的全部外显子的坐标作为外显子窗口，计算每个外显子窗口内的原始reads count，并过滤掉reads count值为0的外显子窗口，对每个reads count值取Log2对数，得到对数矩阵matrixA；

(2)根据对数矩阵matrixA，计算每个基因的各个外显子在所有样本中的平均值，即几何平均数，获得相应的几何平均数矩阵matrixB；

(3)用对数矩阵matrixA中各个外显子窗口的每个reads count值取Log2对数所得值减去对应基因的外显子的几何平均数值，并得到矫正值矩阵matrixC；

(4)再基于矫正值矩阵matrixC，对每个样本的所有外显子窗口的原始readscount进行矫正，方法是用原始reads count除以matrixC对应的值，从而获得最终每个外显子窗口内的标准化reads count，得到矫正后reads count文件。

可选或优选的，上述方法中，步骤(2)的计算方法为：

公式中：

RC

Log2(RC

E

N

可选或优选的，上述方法中，步骤(3)的计算方法如下：

公式中：

RC

Log2(RC

E

N

可选或优选的，上述方法中，步骤(4)中矫正后的reads count矩阵如下表所示：

表格中，NRC

与现有技术相比，本发明具有如下有益效果：

本发明方法利用RNA测序数据先检测出异常表达的基因区域，再获得这些区域的坐标区间，根据最长差异表达区判断是否为基因表达差异导致的假CNV，从而发现真实CNV，该方法在仅利用RNA数据和特殊的生物信息学算法的情况下，就可以有效区分出RNA数据中的长CNV，不需要对样本进行全外显子测序(WES)或者全基因组测序(WGS)，节约成本，提高了RNA数据利用率。

附图说明

图1为Case1已知结果基因区间内差异表达基因的reads count信号，黑色曲线代表Case1，灰色曲线代表所有对照样本的平均值Ct(average)；横坐标是所有基因的所有exon，纵坐标是各个exon的矫正后reads count值，对照样本使用所有对照的矫正覆盖度的平均值。

图2为Case2已知结果基因区间内差异表达基因的reads count信号，黑色曲线代表Case2，灰色曲线代表所有对照样本的平均值Ct(average)；横坐标是所有基因的所有exon，纵坐标是各个exon的矫正后reads count值，对照样本使用所有对照的矫正覆盖度的平均值。

图3为实施例3上调和下调的有意外显子占比的比值分布图。

图4为本发明方法整体流程图。

图5为实施例2中reads count覆盖度矫正流程图。

图6为实施例3有意外显子阈值统计模块流程图。

具体实施方式

下面结合具体实施例对本发明作进一步说明，以使本领域的技术人员可以更好的理解本发明并能予以实施，但所举实施例不作为对本发明的限定。

实施例中的测试样本：

第一批测试样本：用在实施例1中的Case1样本和Case2样本，是已知具体CNV区间的阳性CNV样本；

第二批测试样本：用在实施例2中的Case1样本和Case2样本，也就是第一批测试样本；

第三批测试样本：用在实施例3中，是已经验证不含CNV的一批样本，相当于阈值确定样本。

实施例中的对照样本：

选择表型正常(无测试样本致病表型)、与测试样本组织类型一致的样本作为对照样本，对照样本数量在30个以上；对照样本用在所有实施例。

实施例1常规方法检测基因差异表达

本实施例以经过其他方法验证的已知真实CNV具体区间的两个样本，即阳性CNV样本作为第一批测试样本，分别标记为Case 1和Case2，已经确定CNV区间为拷贝数变异，而非差异化表达引起。测试样本的CNV区间和坐标信息如表1所示。

表1：Case 1和Case2样本的CNV区间

注：①Chr,Start,End表示染色体坐标；

②Size(～100kb)表示区间长度除以100kb。

上述两个Case 1和Case2样本使用差异表达检测软件DROP进行差异表达基因检测所得的差异表达基因情况如表2所示。

表2：DROP软件检测Case 1和Case2样本CNV区间内差异表达基因

根据检测结果可知，DROP软件检测到的差异表达基因与阳性CNV区间重叠率(overlap率)很高，对于DEL(拷贝数缺失)和DUP(拷贝数增加)，DROP分别检测到了DownExp(表达量下降)和UpExp(表达量上升)的结果。

通过以上验证可知，常规的差异表达检测软件所获得的检测结果，并不能区别出所得结果是基因表达差异导致还是拷贝数变异导致，其原因在于基因表达差异的检测结果和拷贝数变异的检测结果都会呈现出reads覆盖度(reads count)上的异常，覆盖度信号见实施例2。从表1-2中看出，表达量下降的表现与真实CNV拷贝数缺失的表现一致，而表达量上升的表现与真实CNV拷贝数增加的表现一致。

实施例2使用本发明方法检测基因差异表达与CNV的覆盖度信号

本实施例2中使用的测试样本记为第二批测试样本，是使用实施例1的已知CNV区间的阳性CNV样本即Case 1和Case 2样本，对已知真实CNV区间内的基因(下文标记为“已知结果基因”)进行分析。

第一部分：第二批测试样本和对照样本进行RNA测序，得到fastq文件后，计算基因的reads count，过滤对照样本

1、筛选对照样本，对照样本数量要求在30个以上，本实施例选择30个对照样本，对照样本需要满足以下条件：

(1)对照样本应该不表现出测试样本的致病表型；

(2)对照样本尽可能与测试样本同一批次建库、测序；

(3)若是公共对照，则对照样本的组织类型应该与测试样本一致。

2、测试样本和对照样本分别提取RNA，建库测序，得到测序数据文件，即FASTQ文件，包括测试样本的FASTQ文件和所有对照样本的FASTQ文件。

使用常规数据质控软件对FASTQ文件进行质控，如cutadapter，包括去接头序列、去除低质量碱基、去除含有N碱基比例高的reads(读段)等。

用STAR软件比对，主要使用不去冗余reads模式，得到比对结果文件，即bam文件。

利用BAM文件和基因组注释文件作为输入，通过htseq-count软件，统计所有对照样本的每个外显子(exon)区间内的原始reads count。

本实施例中，Case 1和Case 2作为第二批测试样本均已知CNV区间，为消除数据波动，在本实施例中对第二批测试样本的已知结果基因reads count也进行了标准化处理。

以Case1为例展示，Case2的标准化矫正计算公式同Case1。

标准化处理算法如下：

2.1获得已知结果基因的坐标后，利用基因组注释文件(GTF文件)提取每个基因的对应的外显子(exon)坐标作为外显子窗口(exon窗口)；

2.2分别从对照样本和Case 1的BAM文件中提取已知结果基因区域内的reads，并计算每个exon窗口内的原始reads count；同时计算所有样本每个exon窗口内的原始readscount。

如表3，是统计Case 1和对照样本已知结果基因相应exon窗口内原始reads count的矩阵，表中的RC

表3：对照样本和Case1样本的基因A、B在各自exon窗口内的原始reads count

2.3对原始reads count进行标准化处理，主要是为了进一步降低批次效应和过度扩增导致的数据异常，具体步骤如下：

(1)对表3中的每个reads count值并过滤掉reads count值为0的exon窗口后取Log2对数，得到对数矩阵matrixA，如表4所示。

表4：对数矩阵matrixA

(2)根据对数矩阵matrixA，计算每个基因各个exon在所有样本中的平均值，即相应的几何平均数矩阵matrixB，然后用对数矩阵matrixA中各个exon窗口的Log2值减去对应基因的exon的几何平均数Eij，并得到矫正值矩阵matrixC，如表5，表5中对照样本i的j基因的exonE窗口reads count的矫正值表示为N

样本i的基因j的第E号exon的Eij值的计算公式为：

表5：标准化矩阵matrixC

根据几何平均数矩阵对样本i的基因j的第E号exon的reads count对数值进行标准化的计算公式为：

(3)再基于矫正值矩阵matrixC对每个样本的所有已知结果基因的所有exon窗口的原始reads count进行矫正，从而得到最终每个exon窗口内的标准化reads count，从而降低其他过度异常表达基因reads count过高或过低的影响。

利用原始reads count矩阵和matrixC，对原始reads count(raw reads count)进行矫正得到标准化原始reads count(normalized reads count)的计算公式如下：

利用上述标准化算法对Case1和Case2样本的原始reads count进行矫正后，得到更能真实反应reads count变异的矩阵数值并绘制reads count差异图，如图1和图2所示。

由于Case1样本的已知结果基因区域范围chr1:6460786-12512047(参见表2的Case1Start及End)内并非所有基因都表达，且部分基因表达量过低或过高，在此我们过滤掉了低质量基因reads count信号，过滤指标为：①Case1样本中表达量为0的基因；②所有对照样本中表达量为0的基因；③对照样本中某个基因表达量高于或低于其他对照样本50％的基因。如果是不知道结果的测试样本，是不需要进行此过滤步骤的。

Case1样本经过上述过滤后该已知结果基因区域内的总基因数量是87个，全部基因(Case1样本的全部过滤后基因，包括已知结果基因区域)的exon总数是570个，exon的reads count变化情况表6。

表6：Case1差异表达基因区域内exon覆盖度变化

注：有35个exon的表达状态不明确，故不计算在内。

对Case2进行与Case1相同条件的标准化处理矫正和过滤后，Case2的已知结果基因区域内有98个有明确reads count变化信号的基因，总共1098个exon，详细结果见表7。

表7：Case2差异表达基因区域内exon覆盖度变化

注：有71个exon的表达状态不明确，故不计算在内。

根据表6和表7可知，在由拷贝数变异(CNV)导致的已知结果基因区域内，基因确实会发生差异表达，且表达状态与CNV状态是相符的，但并非所有基因的reads count信号都能体现CNV状态，如Case1样本中有95.68％(510/533)的exon reads count信号与CNV状态保持一致；在Case2样本中也有着类似的结果，有92.11％(946/1027)的exon覆盖度信号与CNV状态保持一致。

因此，如果只根据reads count信号变化就认为是拷贝数变异的方法是不严谨的，这会受到表达差异的长度的影响，因为从理论上一个生物个体如果在没有真实拷贝数变异的情况下，在较大基因组范围内的基因表达量与其他个体相比是相对固定的，即便存在基因差异表达，那也是个别基因。

至此，仍然不能区分差异表达和拷贝数变异。也就是说，仅使用测试样本和对照样本的reads count数据，即使经过标准化矫正，也无法区分差异表达和拷贝数变异。

继续执行本发明的方法步骤如下：

利用上述标准化方法分别对Case1、Case2和所有对照样本的reads count进行标准化矫正以后，得到矫正后reads count文件，计算所有对照样本之间矫正后reads count文件的相关性系数，过滤掉相关性系数R＜0.7的对照样本；合并测试样本和过滤后对照样本的原始reads count文件，得到原始reads count矩阵文件。本步骤的目的是过滤掉不符合筛选要求的对照样本，以免影响最终结果的判定。

以原始reads count矩阵文件为输入，使用差异表达检测软件DROP检测Case1、Case2中的差异表达基因，通过基因组注释文件获得该差异表达基因的每个外显子的坐标区间，然后以外显子为单位遍历整条染色体，把外显子差异表达状态相同、且基因与基因之间小于1Mb的基因合并为1个差异表达区，同时要求每个差异表达区内的95％以上的外显子的差异表达状态相同，找到Case1、Case2中最长的差异表达区间，如表8。

表8：Case1和Case2的最长差异表达区

判断方法：如果最长的差异表达区间长度在5.6Mb以上，并且上调和下调的有意外显子占比在75％以上，则确定该差异表达区域实际为拷贝数变异，其余则的单个或多个相邻或不相邻的基因为表达差异。根据表8的结果可知，Case1和Case2的CNV区域使用本发明的方法是可以进行判断的。

实施例3确定差异表达与拷贝数变异判断阈值

使用已经验证不含CNV的一批样本作为本实施例3中的测试样本，即第三批测试样本，通过把真实检测到的差异表达区合并找到最长的差异表达区数值和有意外显子最大比例。本实施例3对测试样本和对照样本进行类似实施例2中的方法提取原始reads count，矫正，DROP软件检测差异表达区间。然后以外显子为单位历遍整条染色体，把外显子差异表达状态相同、且基因与基因之间小于1Mb的基因合并为1个差异表达区，同时要求每个差异表达区内的95％以上的外显子的差异表达状态相同，找到测试样本中最长的差异表达区间。

参见表9、表10，分别为第三批测试样本中表达下调和表达上调的合并后的差异表达区。

表格中，ExpType表示差异表达区内的绝大部分基因的差异表达状态，DownExp表示下调，UpExp表示上调；Chr、DiffExp_Start、DiffExp_End、DiffExp_Len(～Mb)分别表示差异表达区的坐标和长度；ExpGene_Num、TotalExon_Num、DiffExon_Num分别表示差异表达区内的有表达信号的总基因数、总共的外显子数量、有表达信号的总外显子数量；UpExp、DownExp分别表示该划定基因区内发生上调和下调的基因数量比例。

表9：表达下调的测试样本的差异表达区

表10：表达上调的测试样本的差异表达区

表9和表10分别统计了这批次测试样本中能够合并的差异表达区间，以及区间内能够提供有效reads count差异信号的外显子数和比值。表中只展示了合并后区间长度在1Mb以上的差异表达区。

根据表9和表10展示的特定区域内不含有CNV的第三批测试样本检测结果可知，最长的差异表达区间长度没有超过5.6Mb，并且上调和下调的有意外显子占比都在75％以下。从而确定了将5.6Mb、75％作为判断差异表达和CNV的阈值是可行的，如果在全基因组范围内能够找到的最长差异表达区长度在5.6Mb以上，并且上调或下调的有意外显子占比在75％以上，则确定该差异表达区域实际为拷贝数变异，其余的单个或多个相邻或不相邻的基因为表达差异。

以上所述实施例仅是为充分说明本发明而所举的较佳的实施例，本发明的保护范围不限于此。本技术领域的技术人员在本发明基础上所作的等同替代或变换，均在本发明的保护范围之内。本发明的保护范围以权利要求书为准。