一种DNA同源重组异常的检测方法及其应用

摘要

本申请提供了一种DNA同源重组异常的检测方法，包括：（1）SNP位点筛选；（2）为筛选到的SNP位点设计捕获探针；（3）基因组DNA提取和文库构建；（4）文库靶向富集；（5）高通量测序并分析测序数据，判断HRD状态时使用Kolmogorov‑Smirnov检验或者scarHRD。

说明书

技术领域

本发明属于医学分子生物学领域，具体涉及一种DNA同源重组异常的检测方法及其应用。

背景技术

DNA修复的方式主要有两种，其中一种是精确无误的同源重组修复（HomologousRecombination Repair，HRR），另一种则是容易产生错误的非同源染色体粘合修复（Non-Homologous End-Joining，NHEJ）。第一种修复方式常见于DNA单链断裂时，第二种修复常见于双链DNA断裂情况。HRR是一条涉及到多个步骤的复杂的信号通路，其中关键蛋白为BRCA1和BRCA2。如果BRCA基因出现突变导致BRCA1和BRCA2蛋白失去功能，就会引起HRR功能异常（Homologous Recombination Deficiency, HRD）。另外，其它HRR相关基因，如 PALB2,CDK12，RAD51, CHEK2, ATM等发生突变、或BRCA1基因启动子发生甲基化、以及其他暂未明确的原因，都会引起HRD，导致基因组不稳定。

HRD是多种肿瘤中常见特征，与肿瘤发生相关，也与PARP抑制剂等肿瘤治疗药物的敏感性相关；此外HRD的检测在基因组功能研究，疾病相关基因筛选也有广泛的用途。

发明内容

一方面，本申请提供一种DNA同源重组异常的检测方法，包括：

（1）SNP位点筛选；

（2）为筛选到的SNP位点设计捕获探针；

（3）基因组DNA提取和文库构建；

（4）文库靶向富集；

（5）高通量测序并分析测序数据。

进一步地，步骤（5）包括数据质控、数据比对和点突变识别。

进一步地，点突变识别包括（a）利用已有的SNP数据库，建立相关性模型，产生重校准表，输入已知的多态性位点数据库，屏蔽不需要重校准的部分；（b）使用模型对原始碱基进行调整，仅调整非已知SNP区域；功能识别SNP和InDel；评估HRD状态。

进一步地，评估HRD状态包括： Kolmogorov-Smirnov检验检测样本与对照样本位点突变频率分布差异。

进一步地，评估HRD状态使用scarHRD的R包。

进一步地，步骤（1）中的SNP位点筛选包括从一个或多个人群中基因组数据库中按照以下规则筛选SNP：

A．剔除Y染色体和线粒体SNP，过滤次等位基因频率（MAF）小于5%的SNP；

B. 过滤显著偏离Hardy-Weinberg equilibrium的SNP；

C. 过滤Insertion和Deletion；

D. 筛选位点上下游75bp范围内不包含重复区域的SNP；

E. 筛选位点上下游75bp范围内序列与人类基因组其它区域无同源性的SNP；

F. 间隔30KB，在间隔点2KB范围内筛选GC含量最接近0.5的SNP。

进一步地，步骤（2）中的SNP位点设计捕获探针包括：从SNP位置前后延伸75bp，提取每个位点的参考序列，去掉重复区域的序列，其中重复序列采用RepeatMask软件分析得到；从第一个碱基开始截取78bp的序列做探针，再一次往后移动 n 个碱基，截取78bp的序列做探针，直到最后一个78bp；每个区域根据外显子的GC含量不同变化n，GC含量太高或太低时n 较小。

进一步地，步骤（4）的文库靶向富集按照专利CN201811600116.3说明书实施例2的试剂和方法进行

进一步地，步骤（5）中的高通量测序为二代测序。

另一方面，本申请提供了用于检测HRD的SNP位点筛选方法，包括从一个或多个人群中基因组数据库中按照以下规则筛选SNP：

A．剔除Y染色体和线粒体SNP，过滤次等位基因频率（MAF）小于5%的SNP；

B. 过滤显著偏离Hardy-Weinberg equilibrium的SNP；

C. 过滤Insertion和Deletion；

D. 筛选位点上下游75bp范围内不包含重复区域的SNP；

E. 筛选位点上下游75bp范围内序列与人类基因组其它区域无同源性的SNP；

F. 间隔30KB，在间隔点2KB范围内筛选GC含量最接近0.5的SNP。

另一方面，本申请提供了用于检测HRD的SNP位点的捕获探针设计方法，包括：从SNP位置前后延伸75bp，提取每个位点的参考序列，去掉重复区域的序列，其中重复序列采用RepeatMask软件分析得到；从第一个碱基开始截取78bp的序列做探针，再一次往后移动n 个碱基，截取78bp的序列做探针，直到最后一个78bp；每个区域根据外显子的GC含量不同变化n，GC含量太高或太低时n 较小。

本发明的方法可用于检测PARP抑制剂类药物的敏感性，也可以用于科研等用途，如研究肿瘤分子特征图谱、RARP抑制剂耐药机理等。

本发明的有益效果

1. 本发明设计50000个snp捕获探针捕测序检测HRD，捕获效率高、捕获稳定性和均一性较好，减少测序成本。

2.通过中国人群HRD检测与分析，拓宽了PARP抑制剂敏感性的生物标志物，为PARP抑制剂敏感人群筛选及用药指导提供了更多参考标准；

具体实施方式

实施例1、DNA同源重组异常的检测方法建立

1、 探针设计

（1）SNP的筛选

从千人基因组数据库包括北京汉族人群（CHB）、南方汉族人群（CHS）、西双版纳的傣族人群（CDX）在内的301个中国人群中按照以下规则筛选50000个SNP位点（部分点见下表1）：

A．剔除Y染色体和线粒体SNP，过滤次等位基因频率（MAF）小于5%的SNP；

B.过滤显著偏离Hardy-Weinberg equilibrium的SNP；

C.过滤Insertion和Deletion；

D.筛选位点上下游75bp范围内不包含重复区域的SNP；

E. 筛选位点上下游75bp范围内序列与人类基因组其它区域无同源性的SNP；

F. 间隔30KB，在间隔点2KB范围内筛选GC含量最接近0.5的SNP；

表1 部分SNP位点

（2）探针设计

从UCSC数据库获取50000个SNP的位置，并前后延伸75bp，提取每个位点的参考序列（参考基因组版本hg38），去掉重复区域的序列，重复序列采用RepeatMask软件分析得到。从第一个碱基开始截取78bp的序列做探针，再一次往后移动 n 个碱基，截取78bp的序列做探针，直到最后一个78bp。每个区域根据外显子的GC含量不同，n 会有变化，GC含量太高或太低 n 越小，探针设计越密集以达到捕获的均一性提高。

2、 基因组DNA提取

新鲜组织基因组DNA提取：取25 mg新鲜组织按照通用型柱式基因组提取试剂盒（康为，CWY004）说明书提取基因组DNA后Qubit检测浓度。

石蜡组织或石蜡切片基因组DNA提取：取10μm厚的石蜡块或8-10片石蜡切片样本按照GeneRead DNA FFPE Kit（Qiagen，180134）说明书提取基因组DNA后

Qubit检测浓度

血液基因组DNA提取：取200μL新鲜或冷冻的抗凝血液样品按照通用型柱式基因组提取试剂盒（康为，CWY004）说明书提取基因组DNA后Qubit检测浓度。

3、 基因组DNA文库构建

将步骤2中组织或血液基因组DNA利用Covaris超声打断仪（Covaris，S220）按照参数Peak Incident Power 175W、Duty Factor 10%、Cycles per Burst 200、Treatment Time180s片段化至200bp左右，然后利用KAPA Hyper文库构建试剂盒制备待测样本基因组DNA文库。

4、 文库靶向富集

将待测样本组织和血液的基因为文库按照专利201811600116.3实施例2的试剂和方法进行靶向捕获，得到靶向捕获文库，在Illumina平台测序。

5、 生信分析

将步骤4得到的靶序列捕获文库通过Nextseq500、X Ten、NovaSeq等二代测序平台进行高通量测序，得到测序原始数据进行以下分析。

基本分析：

（1）碱基识别

使用Illumina官方软件bcf2fastq（version 2.15.0.4），根据样本index序列，将Illumina测序仪下机二进制BCF格式文件转化并拆分为单个样本可读文件fastq格式。

（2）数据质控

使用cutadapt（version 1.16）去除测序接头，删除低质量碱基，生成clean reads。其中 cutadapt（version 1.16）的参数为（-q10,10 --nextseq-trim=10 -aATCTCGTATGCCGTCTTCTGCTTG-A AGATCGGAAGAGCGTCGTGTAGGGAAAGAGTGTAGATCTCGGTGGTCGCCGTATCATT），序列长度小于80。

（3）数据比对

使用序列比对软件BWA（version 0.7.12-r1044）mem功能将clean reads比对至人基因组hg19，其中BWA的软件参数为（mem -M -t）。

使用samtools sort（version 1.2-99-ge2bb18f）功能，根据序列比对位置对序列进行排序。其中参数为（samtools fixmate -O bam sample.sam - | samtools sort -\@5 -m 1G - sample.sort）

对产生的sample.sort.bam文件建立索引，其中参数为（samtools indexsample.sort.bam > sample.sort.bam.bai）

使用GATK的AddOrReplaceReadGroups对提取出来的sample.sort.bam添加readgroup，在后续分析中GATK要求bam文件的header必须包含@RG，参数为（java GATKAddOrReplaceReadGroups -I sample.sort.bam -O sample.sort.header.bam -LBgenome.fa -PL ILLUMINA -SM GP1 -PU GRP1 -SO coordinate --CREATE_INDEX true）

使用bamtools对sample.sort.header.bam文件进行过滤，将一些不合格的序列删除，参数为（bamtools filter -isMapped true -isPaired true -isProperPair true -insample.sort.header.bam -out sample.sort.header.flt.bam）

使用软件bammarkduplicates2过滤测序实验过程中产生的重复序列，参数为（/disk1/software/biobambam2/bin/bammarkduplicates2 I= sample.sort.header.flt.bam O=sample.rmdup.sorted.bam M= sample.duplication-report2.txt markthreads=threadindex=1 rmdup=1）

（4）点突变识别

利用已有的snp数据库，建立相关性模型，产生重校准表，输入已知的多态性位点数据库，用于屏蔽那些不需要重校准的部分。使用软件GATK（version 4.1.4.0）BaseRecalibrator功能建立校准表，参数为（java GATK BaseRecalibrator -R ref_fa -Isample.rmdup.sorted.bam -known-sites dbsnp_150.hg19.vcf -O recal.table -OBItrue -L bed.bed）

根据这个模型对原始碱基进行调整，只会调整非已知SNP区域。使用软件GATK（version4.1.4.0）ApplyBQSR模块对bam文件进行调整，参数为（java GATK ApplyBQSR -R ref_fa -I sample.rmdup.sorted.bam -bqsr recal.table -OBI true –O sample.recal.bam）

使用软件GATK（version 4.1.4.0）Mutect2功能识别SNP和InDel，其参数为（java GATKMutect2 -R ref_fa -I sample.recal.bam -L bed.bed --af-of-alleles-not-in-resource 0.00003125 -O sample.Mutect2.raw.vcf）

HRD状态评估：

结合组织样本和血液样本Kolmogorov-Smirnov检验与scarHRD 的计算结果，对HRD状态进行评估。

（1）组织样本和血液样本Kolmogorov-Smirnov检验

组织样本和血液样本K-S检验是由于对两样本的经验分布函数的位置和形状参数的差异都敏感而成为比较两样本的最有用且常规的非参数方法之一。检验统计量为:

累积分布函数：

其中

采用该方法进行检测样本与对照样本位点突变频率分布是否有差异。HRD状态评估过程中，将bed区域按照每15M进行划分，分别对检测样本和对照样本VCF文件的突变位点进行突变频率、深度、突变区域等信息进行提取，由于关注的是低频等位基因频率分布情况，所以大于0.5的突变频率按照1-0.5进行计算。因此

（a）提出假设H0：F

（b）计算样本累积频率与理论分布累积概率的绝对差，令最大的绝对差为

（c）如果

根据K-S检验结果，统计检测样本与异常样本差异区域的数量，计算差异区域占比，进行LOH/TAI/LST三指标的评分，根据评分，计算HRD的结果。

注：1.“基因组稳定/基因组不稳定”判断主要基于与同源重组修复缺陷（HRD）密切相关的杂合性缺失（LOH）、端粒等位不平衡（TAI）和大片端迁移（LST）等基因组疤痕的综合情况；

2.“LOH/TAI/LST三指标评分”通过综合分析LOH/TAI/LST三种基因组疤痕事件进行评分；

3.“LOH评分”通过分析LOH事件影响的基因组范围进行评分；

4.“LOH/TAI/LST三指标评分”和“LOH评分”是国际主流的基因组稳定性评价标准，只要一个以上指标表现为基因组不稳定，则HRD状态评估结果为“基因组不稳定”。

（2）scarHRD

scarHRD是一个基于二代测序数据的R包，可以通过杂合性缺失、端粒等位基因不平衡以及大片端迁移的数量进行HRD的评估。该方法是本专利评估HRD状态的第二种方法。

计算成对的正常样本与检测样本等位基因的拷贝数

使用软件sequenza-utils 的bam2seqz功能计算成对的正常样本与肿瘤样本的相对拷贝数，调用参数为（sequenza-utils bam2seqz -gc hg19.gc50Base.wig.gz --fastagenome.fa -n control.rmdup.sorted.bam --tumor case.rmdup.sorted.bam -C chr1chr2 chr3 chr4 chr5 chr6 chr7 chr8 chr9 chr10 chr11 chr12 chr13 chr14 chr15chr16 chr17 chr18 chr19 chr20 chr21 chr22 chr23 chr24 chrX | sequenza-utilsseqz_binning -w 50 -s - | gzip > sample.small.seqz.gz）

根据拷贝结果导入scarHRD包，计算HRD状态评分，调用参数为（scar_score("sample.small.seqz.gz",reference = "grch37", seqz=TRUE)）

通过HRDresults.txt文件判断评估样本HRD状态。

实施例2、使用实施例1的方法对11例卵巢癌患者的HRD检测

利用实施例1的方法对11例卵巢癌患者(均经患者知情同意)的HRD情况进行检测。检测结果显示：实施例1的试剂及方法对对目标区域具有高捕获率，目标区域的平均有效测序READ量达到20Mb，目标区域的平均测序深度为3000X以上(见表2)。

表2、利用本发明方法对11例卵巢癌患者检测的质控数据结果

通过分析分析11例卵巢癌患者的BRCA基因突变情况以及HRD特征性基因组疤痕事件，如杂合性缺失（Loss of Heterozygosity）、端粒等位基因不平衡（Telomeric AllelicImbalance）、大片端迁移（Large-scale State Transitions）等，分别根据LOH在基因组中的影响范围以及三种事件的发生频次，综合评估基因组不稳定性，并判断受检者HRD情况（见表3）。

表3、利用本发明方法对11例宫颈癌患者检测的HRD结果