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

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

【技术实现步骤摘要】
一种区分RNA测序数据中基因表达差异与长拷贝数变异的方法
本专利技术涉及生物信息学与精准医学全基因组变异检测
，具体涉及一种区分RNA测序数据中基因表达差异与长拷贝数变异的方法。
技术介绍
拷贝数变异(Copynumbervariation,CNV)是由基因组发生重排而导致的，一般指长度为1kb以上的基因组大片段的拷贝数增加或者减少，主要表现为亚显微水平的缺失和重复，是人类疾病的重要致病因素之一。异常的拷贝数变化(CNV)是许多人类疾病(如癌症、遗传性疾病、心血管疾病)的一种重要分子机制。作为疾病的一项生物标志，染色体水平的缺失、扩增等变化已成为许多疾病研究的热点，然而传统的方法(比如G显带，FISH，CGH等)存在操作繁琐、分辨率低等问题，难以提供变异区段的具体信息。得益于高通量测序(又称下一代测序，NextGenerationSequencing,NGS)技术，在碱基层面分析拷贝数变异已经成为可能。基于NGS平台的变异检测方案众多，如全基因组测序(WGS)、全外显子测序(WES)、捕获测序(Targetedsequencing)、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区间内的原始readscount，并进行标准化处理，得到矫正后readscount文件，计算所有对照样本之间矫正后readscount文件的相关性系数，过滤掉相关性系数R＜0.7的对照样本；合并阈值确定样本和过滤后对照样本的原始readscount文件，得到原始readscount矩阵文件；以原始readscount矩阵文件为输入，使用差异表达检测软件检测阈值确定样本中的差异表达基因，通过基因组注释文件获得该差异表达基因的每个外显子的坐标区间，然后以外显子为单位遍历整条染色体，把外显子差异表达状态相同、且基因与基因之间物理距离小于1Mb的基因合并为1个差异表达区，同时要求每个差异表达区内的95％以上的外显子的差异表达状态相同，找到阈值确定样本中最长的差异表达区间记为nMb，以及上调和下调的有意外显子最小占比值记为m％；(3)将步骤(2)中的阈值确定样本替换为待测试样本，重复步骤(2)的操作，所得最长差异表达区间长度在nMb以上，并且上调和下调的有意外显子占比在m％以上，则确定待测试样本中该差异表达区域实际为拷贝数变异，其余的单个或多个相邻或不相邻的基因为表达差异；所述差异表达状态，是指与对照样本相比，阈值确定样本和待测试样本该基因的表达量的上调或者下调的状态。readscount：读段数量，是指利用二代测序技术对基因或者转录本进行测序，测序测到的每条序列为一个读段，即read，通过统计某一区域测到的reads数目即readscount。可选或优选的，上述方法中，对原始readscount进行标准化处理的方法如下：(1)利用基因组注释文件提取每个基因对应的全部外显子的坐标作为外显子窗口，计算每个外显子窗口内的原始readscount，并过滤掉readscount值为0的外显子窗口，对每个readscount值取Log2对数，得到对数矩阵matrixA；(2)根据对数矩阵matrixA，计算每个基因的各个外显子在所有样本中的平均值，即几何平均数，获得相应的几何平均数矩阵matrixB；(3)用对数矩阵matrixA中各个外显子窗口的每个readscount值取Log2对数所得值减去对应基因的外显子的几何平均数值，并得到矫正值矩阵matrixC；(4)再基于矫正值矩阵matrixC，对每个样本的所有外显子窗口的原始readscount进行矫正，方法是用原始readscount除以matrixC对应的值，从而获得最终每个外显子窗口内的标准化readscount，得到矫正后readscount文件。可选或优选的，上述方法中，步骤(2)的计算方法为：公式中：RCcti-j-eE表示对照样本i的j基因的exonE上的原始readscount数；Log2(RCcti-j-eE)表示对对照样本i的j基因的exonE上的原始readscount数取Log2对数值；Ecti-j-eE表示RCcti-j-eE的几何平均数；Ncti-j-eE对照样本i的j基因的exonE上的原始readscount数的矫正值。可选或优选的，上述方法中，步骤(3本文档来自技高网...

【技术保护点】
1.一种区分RNA测序数据中基因表达差异与长拷贝数变异的方法，其特征在于，包括以下步骤：/n(1)选择已知不含CNV的样本作为阈值确定样本；/n选择表型正常，即无测试样本致病表型、与测试样本组织类型一致的样本作为对照样本，对照样本数量在30个以上；/n(2)阈值确定样本和对照样本分别提取RNA，建库测序，得到测序数据文件，即FASTQ文件；/n对测序数据文件进行质控、比对，得到比对结果文件，即BAM文件；/n利用比对结果文件和基因组注释文件，统计所有对照样本的每个exon区间内的原始reads count，并进行标准化处理，得到矫正后reads count文件，计算所有对照样本之间矫正后reads count文件的相关性系数，过滤掉相关性系数R＜0.7的对照样本；合并阈值确定样本和过滤后对照样本的原始reads count文件，得到原始reads count矩阵文件；/n以原始reads count矩阵文件为输入，使用差异表达检测软件检测阈值确定样本中的差异表达基因，通过基因组注释文件获得该差异表达基因的每个外显子的坐标区间，然后以外显子为单位遍历整条染色体，把外显子差异表达状态相同、且基因与基因之间物理距离小于1Mb的基因合并为1个差异表达区，同时要求每个差异表达区内的95％以上的外显子的差异表达状态相同，找到阈值确定样本中最长的差异表达区间记为nMb，以及上调和下调的有意外显子最小占比值记为m％；/n(3)将步骤(2)中的阈值确定样本替换为待测试样本，重复步骤(2)的操作，所得最长差异表达区间长度在nMb以上，并且上调和下调的有意外显子占比在m％以上，则确定待测试样本中该差异表达区域实际为拷贝数变异，其余的单个或多个相邻或不相邻的基因为表达差异；/n所述差异表达状态，是指与对照样本相比，阈值确定样本和待测试样本该基因的表达量的上调或者下调的状态。/n...

【技术特征摘要】
1.一种区分RNA测序数据中基因表达差异与长拷贝数变异的方法，其特征在于，包括以下步骤：
(1)选择已知不含CNV的样本作为阈值确定样本；
选择表型正常，即无测试样本致病表型、与测试样本组织类型一致的样本作为对照样本，对照样本数量在30个以上；
(2)阈值确定样本和对照样本分别提取RNA，建库测序，得到测序数据文件，即FASTQ文件；
对测序数据文件进行质控、比对，得到比对结果文件，即BAM文件；
利用比对结果文件和基因组注释文件，统计所有对照样本的每个exon区间内的原始readscount，并进行标准化处理，得到矫正后readscount文件，计算所有对照样本之间矫正后readscount文件的相关性系数，过滤掉相关性系数R＜0.7的对照样本；合并阈值确定样本和过滤后对照样本的原始readscount文件，得到原始readscount矩阵文件；
以原始readscount矩阵文件为输入，使用差异表达检测软件检测阈值确定样本中的差异表达基因，通过基因组注释文件获得该差异表达基因的每个外显子的坐标区间，然后以外显子为单位遍历整条染色体，把外显子差异表达状态相同、且基因与基因之间物理距离小于1Mb的基因合并为1个差异表达区，同时要求每个差异表达区内的95％以上的外显子的差异表达状态相同，找到阈值确定样本中最长的差异表达区间记为nMb，以及上调和下调的有意外显子最小占比值记为m％；
(3)将步骤(2)中的阈值确定样本替换为待测试样本，重复步骤(2)的操作，所得最长差异表达区间长度在nMb以上，并且上调和下调的有意外显子占比在m％以上，则确定待测试样本中该差异表达区域实际为拷贝数变异，其余的单个或多个相邻或不相邻的基因为表达差异；
所述差异表达状态，是指与对照样本相比，阈值确定样本和待测试样本该基因的表达量的上调或者下调的状态。


2.根据权利要求1所述的方法，其特征在于，对原始readscount进行标准化处理的方法如下：<...

【专利技术属性】
技术研发人员：鲍远亮，王义亭，王佳，
申请(专利权)人：赛福解码北京基因科技有限公司，
类型：发明
国别省市：北京;11