一种鉴定肿瘤细胞的染色体外环状dna组成基因的方法
阅读说明：本技术 一种鉴定肿瘤细胞的染色体外环状dna组成基因的方法 (Method for identifying extrachromosomal circular DNA (deoxyribonucleic acid) constitutive genes of tumor cells ) 是由 张莹 董科显 贾学渊 王冬 于景翠 白静 傅松滨 于 2021-07-29 设计创作，主要内容包括：本发明公开一种鉴定肿瘤细胞的染色体外环状DNA组成基因的方法,通过在二代测序中加入umi以修正PCR导致的误差,然后针对肿瘤样本的二代测序数据设计了过滤算法以及提出新的更适合分析染色体外DNA的质控标准,最后通过联合三代测序校正接头序列,以确保预测的准确性。(The invention discloses a method for identifying extrachromosomal circular DNA (deoxyribonucleic acid) constitutive genes of tumor cells, which is characterized in that umi is added in second-generation sequencing to correct errors caused by PCR (polymerase chain reaction), then a filtering algorithm is designed aiming at second-generation sequencing data of a tumor sample, a new quality control standard more suitable for analyzing extrachromosomal DNA is provided, and finally a linker sequence is corrected by combining third-generation sequencing to ensure the accuracy of prediction.)
技术领域
本发明涉及生物信息技术领域，具体而言，涉及一种鉴定肿瘤细胞的染色体外环状DNA组成基因的方法。
背景技术
染色体外环状DNA(extrachromosomal circular DNAs，eccDNAs)为产生于染色体上的序列，在起源上有较高的异质性，可以影响细胞的生命活动，促进肿瘤细胞演进和适应性进化，是肿瘤的重要基因组特征。
染色体外环状DNA，最早在1965年，由HOTTA Y通过对肿瘤细胞核型分析发现的，但是因为当时分子技术水平的局限，未能得到广泛的关注，后来由于FISH(Fluorescence insitu hybridization，荧光原位杂交技术)技术的兴起，人们进一步了解了染色体外环状DNA的分子结构。但是FISH技术存在如下缺点：1)FISH实验手段非常繁琐，造价高昂；2)其分辨率很低，只能检测到1M左右的区域；3)对实验样本的要求很高，如必须得到中期分裂的细胞，这对于临床上组织样本的检测是无法实现的。因此，由于分辨率的限制，目前的技术只能获得较大片段的染色体外环状DNA上的基因信息，急需更高分辨率，准确性更高的手段来探索染色体外环状DNA上基因，这对于肿瘤的诊断及治疗有着重要意义。
在2019年，Deshpande V发明了一款软件，AmpliconArchitect，其利用二代测序数据对染色体外DNA进行预测。但是该软件在使用时需要提供指定的bed(靶向测序中一个能够指示目标基因位置的文件类型)区域，而该软件目前没有任何可以预测bed区域的工具，另外，因为肿瘤细胞的基因扩增是异常的，而该软件并未就肿瘤细胞的拷贝数进行有效校正。
因此，现有技术对于染色体外环状DNA的结构、组成仍然缺乏较为广泛而深入的研究，需要一种能够对染色体的染色体外环状结构进行预测及验证的方法。
发明内容
为了解决上述问题，本发明提供一种鉴定肿瘤细胞的染色体外环状DNA组成基因的方法，通过在测序建库过程中，加入指定的umi(Unique Molecular Identifiers序列)，以修正PCR(Polymerase Chain Reaction，聚合酶链反应)导致的误差，并且针对肿瘤样本的二代测序数据设计了过滤算法以及提出新的更适合分析染色体外DNA的质控标准。
本发明的另一目的在于，提供了一种能够准确预测bed区域的算法，针对肿瘤细胞内异常扩增的基因组片段以及发生染色体结构变异片段进行定位，为后续利用图论进行环状染色体分析提供基础。
为达到上述目的，本发明提供了一种鉴定肿瘤细胞的染色体外环状DNA组成基因的方法，其包括以下步骤：
步骤1：通过illumina平台对待分析样本进行全基因组二代测序，得到过滤后满足条件的fastq格式文件的二代测序结果，并进行质控；
步骤2：通过CNVkit对步骤1中二代测序结果进行拷贝数变异统计，并对得到的cnr文件进行均一化处理；
步骤3：通过AmpliconArchitect软件进行染色体外DNA的结构预测，并根据hg19参考基因组得到染色体外环状DNA在基因坐标轴的范围以及对应范围上所包含基因；
步骤4：通过Delly软件进行SV分析，并对得到的vcf格式的结果文件进一步进行筛选；
步骤5：将步骤4得到的结果与步骤2得到的结果结合后，与步骤3的结果取交集，得到一区域交界处的接头序列，然后将环状结构进行组装得到一样本；
步骤6：对组装后得到的样本进行Pacbio长序列测序，并与hg19参考基因组比对后得到比对后的第二bam文件，再进行筛选校正；
步骤7：将步骤5中组装后的样本序列重命名并加入到hg19参考基因组，再次进行比对以确认是否为染色体外环状DNA。
在本发明一实施例中，其中，步骤1具体为：
步骤11：采用barcode+NGS的方式进行全基因组二代测序，具体为通过illumina平台对待分析样本进行全基因组二代测序，再通过脚本程序在得到的fastq文件中提取umi序列，对提取到的所有umi序列进行聚类和信息统计分析，将支持数的相同的umi序列进行合并；
步骤12：根据统计分析结果将测序得到的序列中相同的barcode的序列合并，将合并后的集合作为一致性序列的集合，并在集合内修正因PCR和测序误差对reads本身的影响；
步骤13：通过BWA软件对提取后的fastq文件与hg19参考基因组进行比对，得到比对后第一bam格式文件。
在本发明一实施例中，其中，步骤1中的质控标准具体为去除N的比例大于10％的reads，其中N为未成功识别到的碱基。
在本发明一实施例中，其中，步骤2中均一化处理的算法具体为：
步骤21：将所述第一bam文件通过pysam选取10kb为一个bin，计算所有bin的平均覆盖度，排除覆盖度为0或者覆盖度大于5倍整体平均覆盖度方差的区域；
步骤22：重新计算平均覆盖度、方差及标准差,并通过下式计算得到任一窗口区域的高斯核函数参数，
式中，j代表迭代次数，ri，rj代表窗口区域中的任意两个bin即bin i和bin j，Hb为bin的位序，Hr为bin的带宽，norm表示恒定系数；
步骤23：当ΔF为大于0时，认定对应区域为一个边界；当ΔF小于0时，将对应区域合并到下一个区域，继续根据上式进行迭代，以得到筛选区域，其中，ΔF为当前计算的高斯核函数参数的值与上一次计算的高斯核函数参数的值的差值；
步骤24：在cnr文件中对筛选出来的区域进行K-mean聚类得到cns文件，并选取Copynumber>5,Length>100的区域作为bed文件候选区域，其中，Copynumber为拷贝数，Length为长度。
在本发明一实施例中，其中，步骤3具体为：
步骤31：通过AmpliconArchitect软件对步骤13得到的第一bam文件进行染色体外DNA的结构预测；
步骤32：当测序深度>5X时通过向下采样将深度调整为5X，并使用hg19参考基因组得到染色体外环状DNA在基因坐标轴的范围，其中X为倍数；
步骤33：得到基因坐标轴的范围上所包含的基因。
在本发明一实施例中，其中，步骤4具体为：
步骤41：通过Delly软件对步骤13得到的第一bam文件进行SV分析，得到vcf格式的结果文件；
步骤42：对vcf格式的结果文件按照以下条件进行进一步的筛选：
选择SVTYPE为BND；
选择INFO为PASS；
选择RV比率>0.5,DV比率>0.8，其中，RV比率＝RV/(RV+RR)，DV比率＝DV/(DV+DR)，RR为split reads，RV为split reads的变异，DR为discordant reads，DV为discordantreads的变异；以及
属于步骤12具备一致性的序列。
在本发明一实施例中，其中，步骤5中得到一区域交界处的接头序列的过程具体为：将取交集的结果比对步骤12得到的一致性序列的集合，得到所述区域交界处的接头序列中的一致性序列。
在本发明一实施例中，其中，步骤6具体为：
步骤61：将步骤5得到的样本进行Pacbio长序列测序，其中，由于样本的序列长度满足正态分布，则测序长度峰值出现在2kb处；
步骤62：将测序结果与hg19参考基因组进行比对，得到比对后的第二bam文件；
步骤63：从得到的第二bam文件中筛选出步骤5中接头序列中的一致性序列；
步骤64：根据步骤61的长序列测序结果，对步骤63的结果进行校正。
在本发明一实施例中，其中，步骤7具体为：
步骤71：将步骤5中组装后的序列重命名为chrR，并将其加入到hg19参考基因组，得到新构建的参考基因组；
步骤72：将步骤6得到的第二bam文件与步骤71得到的新构建的参考基因组进行比对，若能够完全比对上的reads的数目>0，则确认对应reads来自于所预测的染色体外DNA。
本发明公开的鉴定肿瘤细胞的染色体外环状DNA组成基因的方法，与现有技术相比，至少具有以下优点：
1)在二代测序中加入umi可以修正PCR导致的误差，可以将干扰序列信号影响降到最低；
2)针对肿瘤样本的二代测序数据设计了过滤算法以及提出新的更适合分析染色体外DNA的质控标准；
3)联合三代测序校正接头序列，并制作环状染色体特异性参考基因组进行二次比对，可以确保预测的准确性。
附图说明
为了更清楚地说明本发明实施例或现有技术中的技术方案，下面将对实施例或现有技术描述中所需要使用的附图作简单地介绍，显而易见地，下面描述中的附图仅仅是本发明的一些实施例，对于本领域普通技术人员来讲，在不付出创造性劳动的前提下，还可以根据这些附图获得其他的附图。
图1为本发明一实施例的流程图；
图2为本发明一实施例执行步骤1得到的部分结果示意图；
图3为本发明一实施例执行步骤3得到的部分结果示意图；
图4为本发明一实施例执行步骤7得到的部分结果示意图；
图5为本发明一实施例执行步骤8得到的验证示意图。
具体实施方式
下面将结合本发明实施例中的附图，对本发明实施例中的技术方案进行清楚、完整地描述，显然，所描述的实施例仅仅是本发明一部分实施例，而不是全部的实施例。基于本发明中的实施例，本领域普通技术人员在没有付出创造性劳动前提下所获得的所有其他实施例，都属于本发明保护的范围。
实施例一
图1为本发明一实施例的流程图，如图1所示，本实施例提供了一种鉴定肿瘤细胞的染色体外环状DNA组成基因的方法，其包括以下步骤：
步骤1：通过illumina平台(illumina公司提供的二代测序平台)对待分析样本进行全基因组二代测序，得到过滤后满足条件的fastq格式(一种存储了生物序列以及相应的质量评价的文本格式)文件的二代测序结果，并进行质控；
在本实施例中，其中，步骤1具体为：
步骤11：采用barcode(DNA条形码)+NGS(Next Generation Sequencing，下一代测序技术，即高通量测序)的方式进行测序，具体为通过illumina平台对待分析样本进行全基因组二代测序，再通过脚本程序在得到的fastq文件(一种存储了生物序列以及相应的质量评价的文件格式)中提取umi(Unique Molecular Identifiers)序列，对提取到的所有umi序列进行聚类和信息统计分析，将支持数的相同的umi序列进行合并，以实现监控和移除因PCR(Polymerase Chain Reaction，聚合酶链反应)和测序错误而产生的umi序列的目的；
步骤12：根据统计分析结果将测序得到的序列中相同的barcode的序列合并，将合并后的集合作为一致性序列的集合，并在集合内修正因PCR和测序误差对reads(高通量测序产生的序列)本身的影响，从而将PCR和测序错误产生的干扰序列信号影响降到最低，实现高精度和超低频突变测序；由于进行PCR的过程可能会出现大量重复序列，在比对的过程中无法准确判定这些重复的reads应该比对到哪里，因此加入umi后，相当于给每条reads添加一个条形码(barcode)，这样可以知道这些reads是从哪里重复得来。
步骤13：通过BWA软件(一款将序列比对到参考基因组上的软件)对提取后的fastq文件与hg19参考基因组进行比对，得到比对后第一bam格式文件(一种序列比对文件的格式)。
在本实施例中，其中，步骤1中的质控标准具体为去除N的比例大于10％的reads，其中N为未成功识别到的碱基。通常情况下正常识别的用TCGA表示，未成功识别的用N表示。由于在人群中有一些区域本身的序列差异比较大，如免疫复合物基因区域，因此在分析的过程中需要排除这些区域。
步骤2：通过CNVkit(进行CNV分析的一种工具)对步骤1中二代测序结果进行拷贝数变异统计，在得到cnr文件后，对这个cnr文件进行均一化处理；
在本实施例中，其中，步骤2中均一化处理的算法具体为：
步骤21：将第一bam文件通过pysam(一种处理基因组数据的python模块)选取10kb为一个bin(可视为一个集合)，计算所有bin的平均覆盖度，排除覆盖度为0或者覆盖度大于5倍整体平均覆盖度方差的区域；
步骤22：重新计算平均覆盖度(μwindow size,μws)、方差(θwindow size,θws)及标准差(σwindow size),并利用公式(1)计算得到任一窗口区域的高斯核函数参数，
式中，j代表迭代次数，ri，rj代表窗口(window)区域中的任意两个bin即bin i和bin j，Hb为bin的位序，Hr为bin的带宽，norm表示恒定系数；
步骤23：当ΔF为大于0时，认定对应区域为一个边界；当ΔF小于0时，将对应区域合并到下一个区域，继续根据上式进行迭代，以得到筛选区域，其中，ΔF为当前计算的高斯核函数参数的值与上一次计算的高斯核函数参数的值的差值，例如ΔF可以是本次通过式(1)计算得到的值减去上一次通过式(1)计算得到的值的差值；
步骤24：在cnr文件中对筛选出来的区域进行K-mean聚类(k均值聚类算法)得到cns文件，并选取Copynumber>5,Length(kb)>100的区域作为bed文件(靶向测序中一个能够指示目标基因位置的文件类型)候选区域，其中，Copynumber为拷贝数，Length为长度。
步骤3：通过AmpliconArchitect软件(一种用于寻找基因组来源的eccDNA的软件)进行染色体外DNA的结构预测，并根据hg19参考基因组得到染色体外环状DNA在基因坐标轴的范围以及对应范围上所包含基因；
在本实施例中，其中，步骤3具体为：
步骤31：通过AmpliconArchitect软件对步骤13得到的第一bam文件进行染色体外DNA的结构预测；
步骤32：当测序深度>5X时通过向下采样将深度调整为5X，并使用hg19参考基因组得到染色体外环状DNA在基因坐标轴的范围，其中X为倍数；
步骤33：得到基因坐标轴的范围上所包含的基因。
其中，步骤3中可以得到一个cycle.txt文件以及一个edges.txt，从而得到具体的范围和其上的基因。
步骤4：通过Delly软件(一种变异检测软件)进行SV(Structure Variantions，基因组结构性变异)分析，并对得到的vcf格式的结果文件进一步进行筛选；
在本实施例中，其中，步骤4具体为：
步骤41：通过Delly软件对步骤13得到的第一bam文件进行SV分析，得到vcf格式的结果文件；
步骤42：对vcf格式的结果文件按照以下条件进行进一步的筛选：
选择SVTYPE(SV的类型)为BND(一种变异类型)；
选择INFO(附加信息列)为PASS；
选择RV ratio(比率)>0.5且DV ratio(比率)>0.8，其中，RV ratio(比率)＝RV/(RV+RR)，DV ratio(比率)＝DV/(DV+DR)，RR为split reads(一类特殊的read，其出现通常是由基因组中的结构变异造成的，这类read在映射中不再保持连续序列的形式，而是包含了一定长度的空位)，RV为split reads的变异(variant)，DR为discordant reads(不一致的reads，指的是R1和R2端比对之间的距离超过了期望的插入片段长度或者比对到了不同链的reads，这些reads比对情况可能是由于基因组结构变异引起的)，DV为discordantreads的变异(variant)；以及
属于步骤12具备一致性的序列(consensus)。
步骤5：将步骤4得到的结果与步骤2得到的结果结合后，与步骤3的结果取交集，得到一区域交界处的接头序列，然后将环状结构进行组装得到一样本；
在本实施例中，其中，步骤5中得到一区域交界处的接头序列的过程具体为：将取交集的结果比对步骤12得到的一致性序列的集合，得到所述区域交界处的接头序列中的一致性序列。
步骤6：对组装后得到的样本进行Pacbio(一种用于进行第三代测序的系统)长序列测序，并与hg19参考基因组比对后得到比对后的第二bam文件，再进行筛选校正；
在本实施例中，其中，步骤6具体为：
步骤61：将步骤5得到的样本进行Pacbio长序列测序，其中，由于样本的序列长度满足正态分布，则测序长度峰值出现在2kb处；
步骤62：将测序结果与hg19参考基因组进行比对，得到比对后的第二bam文件；
步骤63：从得到的第二bam文件中筛选出步骤5中接头序列中的一致性序列；
步骤64：根据步骤61的长序列测序结果，对步骤63的结果进行校正。
步骤7：将步骤5中组装后的样本序列重命名并加入到hg19参考基因组，再次进行比对以确认是否为染色体外环状DNA。
在本实施例中，其中，步骤7具体为：
步骤71：将步骤5中组装后的序列重命名为chrR，并将其加入到hg19参考基因组，得到新构建的参考基因组；
步骤72：将步骤6得到的第二bam文件与步骤71得到的新构建的参考基因组进行比对，若能够完全比对上的reads的数目>0，则确认对应reads来自于所预测的染色体外DNA。
实施例二
本实施例以人卵巢癌细胞系UACC-1598作为材料，鉴定肿瘤细胞的染色体外环状DNA组成基因的方法如下：
步骤1：经过二代测序得到fastq文件，在按照去除N的比例大于10％的reads的标准进行质控后，质量合格。图2为本发明一实施例执行步骤1得到的部分结果示意图，如图2所示的部分质控合格的reads。然后，通过BWA软件对提取后的fastq文件与hg19参考基因组进行比对，得到比对后第一bam格式文件。
步骤2：利用CNVkit筛选出Copynumber>3，Length(kb)>100的区域，如表1中展示的部分筛选出的实例。
表1
步骤3：利用AmpliconArchitect软件，通过输入步骤1的第一bam文件进行环状预测，得到其预测染色体外DNA上所容纳基因，图3为本发明一实施例执行步骤3得到的部分结果示意图，如图3所示，部分得到的预测染色体外DNA上所容纳的基因显示于图中。
步骤4：利用delly得出SV突变结果，经过过滤，得到符合条件的染色体变异信息，如下表2展示部分过滤后的结果。
表2
步骤5：综合步骤2、步骤3和步骤4的结果，如果步骤4中所得到的基因，在步骤2中有证实，则需要判断该基因是否在步骤3中存在合理的染色体变异信息，如果存在，则判断该基因位于染色体外环状DNA上。根据本实施例的材料，得出染色体外DNA上基因有ADAMTSL4、ECM1、ENSA、GOLPH3L、HORMAD1、MCL1、RPRD2、TARS2、MYCN、CLDN11、EIF5A2、RPL22L1、SLC7A14、SKIL、MYCNO、GACAT3、ADAMTSL4-AS1、RN7SL473P及FALEC，并且本实施例经常规实验证实，上述基因存在于UACC-1598的染色体外环状DNA上。
步骤6：利用Pacbio三代测序，得到比对后的第二bam文件，进一步得到接头处序列，选取其中质量值为60的序列，校正接头序列。
步骤7：将组装后的序列命名为chrR，加入到参考基因组得到新构建的参考基因组，将第二bam文件与新构建的参考基因组再次进行比对，在比对结果中，若chrR完全比对上的reads数目>0，则确认对应的序列为染色体外环状DNA。图4为本发明一实施例执行步骤7得到的部分结果示意图，如图4所示，由于该序列中完全比对上的reads数目>0，则确认该序列染色体外环状DNA。
步骤8：在完成上述步骤1～7后，利用免疫荧光原位杂交技术进行验证，图5为本发明一实施例执行步骤8得到的验证示意图，如图5所示，根据所预测染色体外环状DNA序列特异性探针设计，发现荧光信号出现在染色体外环状DNA上面，从而证明上述步骤1～7的预测方法是准确的。
本发明公开的鉴定肿瘤细胞的染色体外环状DNA组成基因的方法，通过在二代测序中加入umi以修正PCR导致的误差，然后针对肿瘤样本的二代测序数据设计了过滤算法以及提出新的更适合分析染色体外DNA的质控标准，最后通过联合三代测序校正接头序列，并制作环状染色体特异性参考基因组进行二次比对，以确保预测的准确性。
本领域普通技术人员可以理解：附图只是一个实施例的示意图，附图中的模块或流程并不一定是实施本发明所必须的。
本领域普通技术人员可以理解：实施例中的装置中的模块可以按照实施例描述分布于实施例的装置中，也可以进行相应变化位于不同于本实施例的一个或多个装置中。上述实施例的模块可以合并为一个模块，也可以进一步拆分成多个子模块。
最后应说明的是：以上实施例仅用以说明本发明的技术方案，而非对其限制；尽管参照前述实施例对本发明进行了详细的说明，本领域的普通技术人员应当理解：其依然可以对前述实施例所记载的技术方案进行修改，或者对其中部分技术特征进行等同替换；而这些修改或者替换，并不使相应技术方案的本质脱离本发明实施例技术方案的精神和范围。