CN106715711B - 确定探针序列的方法和基因组结构变异的检测方法 - Google Patents

确定探针序列的方法和基因组结构变异的检测方法 Download PDF

Info

Publication number
CN106715711B
CN106715711B CN201480080426.0A CN201480080426A CN106715711B CN 106715711 B CN106715711 B CN 106715711B CN 201480080426 A CN201480080426 A CN 201480080426A CN 106715711 B CN106715711 B CN 106715711B
Authority
CN
China
Prior art keywords
region
candidate
probe
target sample
snp
Prior art date
Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
Active
Application number
CN201480080426.0A
Other languages
English (en)
Other versions
CN106715711A (zh
Inventor
李剑
王煜
李尉
李金良
赵霞
陈仕平
张现东
刘赛军
Current Assignee (The listed assignees may be inaccurate. Google has not performed a legal analysis and makes no representation or warranty as to the accuracy of the list.)
BGI Shenzhen Co Ltd
Original Assignee
BGI Shenzhen Co Ltd
Priority date (The priority date is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the date listed.)
Filing date
Publication date
Application filed by BGI Shenzhen Co Ltd filed Critical BGI Shenzhen Co Ltd
Publication of CN106715711A publication Critical patent/CN106715711A/zh
Application granted granted Critical
Publication of CN106715711B publication Critical patent/CN106715711B/zh
Active legal-status Critical Current
Anticipated expiration legal-status Critical

Links

Images

Classifications

    • CCHEMISTRY; METALLURGY
    • C12BIOCHEMISTRY; BEER; SPIRITS; WINE; VINEGAR; MICROBIOLOGY; ENZYMOLOGY; MUTATION OR GENETIC ENGINEERING
    • C12QMEASURING OR TESTING PROCESSES INVOLVING ENZYMES, NUCLEIC ACIDS OR MICROORGANISMS; COMPOSITIONS OR TEST PAPERS THEREFOR; PROCESSES OF PREPARING SUCH COMPOSITIONS; CONDITION-RESPONSIVE CONTROL IN MICROBIOLOGICAL OR ENZYMOLOGICAL PROCESSES
    • C12Q1/00Measuring or testing processes involving enzymes, nucleic acids or microorganisms; Compositions therefor; Processes of preparing such compositions
    • C12Q1/68Measuring or testing processes involving enzymes, nucleic acids or microorganisms; Compositions therefor; Processes of preparing such compositions involving nucleic acids

Landscapes

  • Chemical & Material Sciences (AREA)
  • Organic Chemistry (AREA)
  • Life Sciences & Earth Sciences (AREA)
  • Zoology (AREA)
  • Wood Science & Technology (AREA)
  • Proteomics, Peptides & Aminoacids (AREA)
  • Health & Medical Sciences (AREA)
  • Engineering & Computer Science (AREA)
  • Microbiology (AREA)
  • Immunology (AREA)
  • Physics & Mathematics (AREA)
  • Molecular Biology (AREA)
  • Biotechnology (AREA)
  • Biophysics (AREA)
  • Analytical Chemistry (AREA)
  • Biochemistry (AREA)
  • Bioinformatics & Cheminformatics (AREA)
  • General Engineering & Computer Science (AREA)
  • General Health & Medical Sciences (AREA)
  • Genetics & Genomics (AREA)
  • Measuring Or Testing Involving Enzymes Or Micro-Organisms (AREA)

Abstract

本发明提供了基于参考序列确定探针序列的方法和基因组结构变异的检测方法。其中,基于参考序列确定探针序列的方法包括:基于多个离散高频SNP位点,构建第一候选探针集,第一候选探针集由多个候选探针构成,并且多个候选探针中的每一个均含有至少一个离散高频SNP;将第一候选探针集中的多个候选探针与参考序列进行比对,以便获得比对结果;基于比对结果,对第一候选探针集进行第一筛选,获得第二候选探针集;将参考序列划分为多个具有预定长度的窗口,分别将第二候选探针集中的多个候选探针分配至各自匹配的窗口,以确定多个候选探针各自的位置信息;基于所说的位置信息以及离散高频SNP的等位基因频率,对第二候选探针集进行第二筛选,以便确定所述探针序列。

Description

确定探针序列的方法和基因组结构变异的检测方法
优先权信息
技术领域
本发明涉及基因组学及生物信息学技术领域,具体涉及确定探针序列的方法和基因组结构变异的检测方法。
背景技术
DNA拷贝数变异(copy number variation,CNV)和杂合性丢失(Loss ofheterozygosity,LOH)是不同类型的基因组变异。CNV是一种常见基因组结构变异,片段从1kb到几Mb不等,主要表现为亚显微水平的缺失和重复。LOH是指一对染色体上某一个染色体上基因缺失,与之配对的染色体上仍然存在,表现为在DNA很长一段区域只出现纯合子SNP。当LOH没有发生拷贝数的变化,即只从一个亲本遗传两个副本,被称单亲二倍体(uniparental disomy,UPD)。CNV,LOH,和UPD与许多常见的遗传性疾病,癌症和其他复杂疾病相关。建立一种准确、全面、高效、快速、简单、经济的检测CNV、LOH和UPD的方法,对于研究染色体变异事件,明确相关疾病的病因和采取相应的治疗方案,都具有重要的价值。
目前已有一些检查技术,比如PCR技术,包括实时荧光定量PCR技术和多重连接扩增技术(Multiplex Ligation-dependent Probe Amplification,MLPA),实时荧光PCR技术每次检测分析一个或数个靶点,MLPA一次能够分析40多个序列,灵敏度高,检测范围受限于探针所针对的染色体和区域;FISH技术,FISH技术一般用于检测特定的几条染色体,无法检测未知区域;基于芯片的技术,包括基于芯片的比较基因组杂交技术(array-basedComparative Genomic Hybridization,aCGH)和基于SNP芯片的技术(SNP-array),aCGH可检测全基因组范围内的CNV,不能检测出多倍体,小片段的丢失的漏检率高;以及测序技术,基于全基因组测序(whole genome sequnecing,WGS)检测全基因组范围的结构变异和基于目标区域测序检测目标区域的变异,主要有四种方法分析CNV,包括:配对末端映射(paired-end mapping),读长深度分析(read-depth analysis),分开读长策略(split-read strategies)和序列组装比较(sequence assembly comparisons)。
随着测序技术的发展,有必要研究基于测序结果特别是局部区域测序结果来发现基因组结构异常的手段,包括发现染色体非整倍性、CNV、***缺失(insertion-deletion,indel)、LOH、UPD以及SNP的手段。
发明内容
本发明的一方面提供一种基于参考序列确定探针序列的方法,包括以下步骤:基于多个离散高频SNP位点,构建第一候选探针集,第一候选探针集由多个候选探针构成,并且多个候选探针中的每一个均含有至少一个离散高频SNP;将第一候选探针集中的多个候选探针与参考序列进行比对,以便获得比对结果;基于比对结果,对第一候选探针集进行第一筛选,获得第二候选探针集;将参考序列划分为多个具有预定长度的窗口,分别将第二候选探针集中的多个候选探针分配至各自匹配的窗口,以确定多个候选探针各自的位置信息;基于所说的位置信息以及离散高频SNP的等位基因频率,对第二候选探针集进行第二筛选,以便确定所述探针序列。其中,离散高频SNP位点为等位基因频率大于10%,优选的不大于90%,并且与任意另外一个离散高频SNP位点在参考基因组上的物理距离不小于候选探针长度,候选探针长度为50-250mer。
利用本发明的确定探针序列的方法获得的探针,用于杂交捕获基因组获得多个基因组局部区域,捕获得的多个局部区域能够代表全基因组、能够反映全基因组变异信息,用于发现全基因范围的结构变异的发生。
本发明的另一方面提供了一种检测基因组结构变异的方法,适用于检测染色体非整倍性、拷贝数变异和***缺失,包括以下步骤:对目标样本基因组核酸进行测序,以获得基因组测序结果,所说的基因组测序结果由多个读段构成,可选地,所说的测序包括采用探针进行筛选,其中,探针是通过本发明一方面提供的基于参考序列确定探针序列的方法获得的。基因组测序结果,可以通过提取基因组DNA,依据现有高通量平台指导手册进行文库构建及上机测序获得;基因组测序结果也可以通过探针捕获目标样本的基因组并进行测序获得的,探针可以通过本发明一方面提供的基于参考序列确定探针序列的方法获得;将参考基因组分为m个区域,利用基因组测序结果中落入区域i的读段计算目标样本基因组区域i的覆盖深度TDi,其中,m和i为自然数,1≤i≤m,10<m;基于目标样本基因组区域i的覆盖深度与k个参照样本的区域i的覆盖深度的差异程度,判断目标样本区域i结构变异的发生,其中,k为自然数,k≥2,各参照样本的区域i的覆盖深度的得来方法可参照目标样本区域i的覆盖深度的获得方法。通过合并邻近发生结构变异的区域,进一步检测合并后的区域是否发生大的结构变异,或者说进一步检测发生在区域i的结构变异是否横跨几个区域。
本发明的再一方面提供了适用于检测另一种基因组结构变异——杂合性丢失的方法,包括以下步骤:获取目标样本的基因组测序结果,可选地,所说的基因组测序结果是通过探针捕获目标样本的基因组并进行测序获得的,探针是按照本发明一方面提供的基于参考序列确定探针序列的方法获得的;将基因组分成m’个区域,基于基因组测序结果中落在区域i中的读段和群体区域i数据,获得目标样本基因组区域i和群体区域i共有的SNP集,分别计算目标样本和群体的共有SNP集中的各个SNP位点所在片段的杂合度,获得目标样本基因组区域i的杂合度集Ui,和群体区域i的杂合度集U0i,比较目标样本Ui和群体U0i以确定目标样本区域i杂合性丢失是否发生;其中,共有SNP集中的每个SNP的等位基因频率都大于0.1,共有SNP集中的一个SNP位点所在片段是以与该SNP相邻的上下游两个SNP为边界点的,m’和i为自然数,m’≥i≥1,m’≥6。抽取多少样本能够真实反映群体,可根据检测所需的精确度、统计方法、样本数据分布情况等确定,群体数据由多个同物种的样本数据构成,可通过全基因组测序、或者依据获得目标样本数据的方法、或者从已完成已公开的数据库或网站获得,比如千人基因组数据。
本发明的再一方面提供一种计算机可读存储介质,用于存储供计算机执行的程序,本领域普通技术人员可以理解,在执行该程序时,通过指令相关硬件可完成上述检测基因组结构变异的各种方法的全部或部分步骤。所称存储介质可以包括:只读存储器、随机存储器、磁盘或光盘等。
根据本发明的最后一方面提供检测基因组结构变异的装置,包括:数据输入单元,用于输入数据;数据输出单元,用于输出数据;存储单元,用于存储数据,其中包括可执行的程序;处理器,与上述数据输入单元、数据输出单元及存储单元数据连接,用于执行存储单元中存储的可执行的程序,程序的执行包括完成上述检测基因组结构变异的各种方法的全部或部分步骤。
利用本发明的基于参考序列确定探针序列的方法获得的探针,利用探针或者包含这些探针的固相/液相芯片进行目标区域捕获测序,能够低测序成本的实现在全基因组范围内检测结构变异,包括覆盖人的23对染色体检测CNV、LOH和UPD,而且检测分辨率能根据需求通过调整探针的平均间距分布即增加/减少SNP位点进行调整。利用本发明的目标区域捕获测序结合生物信息分析方法实现了在全基因组范围内进行高分辨率、高准确性、高通量、低成本的CNV、LOH和UPD检测,同时本发明的基因组结构变异检测方法也适用于染色体非整倍性变异、SNP和Indel的检测,适用于基于全基因测序数据的结构变异分析检测。
本发明的附加方面和优点将在下面的描述中部分给出,部分将从下面的描述中变得明显,或通过本发明的实践了解到。
附图说明
本发明的上述和/或附加的方面和优点,结合下面附图对实施方式的描述将变得明显和容易理解,其中:
图1是本发明一个实施方式中的SeTR探针在全基因组上的特性的示意图,(A)SeTR探针序列的长度分布图;(B)SETR探针中两两探针的物理距离分布图。
图2是本发明一个实施方式中的SeTR探针的测试结果图,(A)目标区域的覆盖深度分布图(B)支持ref碱基型和非ref碱基型的reads分布图。
图3是本发明的一个实施方式中的CNV、LOH和UPD的检测流程示意图。
图4是本发明的一个实施方式中的Ri基准线示图。
图5是本发明一个实施方式中的检测到的一个样本(GM50275)的基因组结构变异的示意图,圆环由外到里,依次为I)染色体信息,II)ri值的变化(波浪线);III)Rhet对应的P值变化,IV)Rhet值变化(点)。
发明详细描述
下面详细描述本发明的实施例。下面通过参考附图描述的实施例是示例性的,仅用于解释本发明,而不能理解为对本发明的限制。
需要说明的是,术语“第一”、“第二”仅用于描述目的,而不能理解为指示或暗示相对重要性或者隐含指明所指示的技术特征的数量。由此,限定有“第一”、“第二”的特征可以明示或者隐含地包括一个或者更多个该特征。进一步地,在本发明的描述中,除非另有说明,“多个”的含义是两个或两个以上。
根据本发明的一种实施方式,提供一种基于参考序列确定探针序列的方法,包括以下步骤:
步骤一:构建第一候选探针集
利用分布于基因组的离散高频SNP位点构建第一候选探针集,第一候选探针集中的每一条候选探针包含至少一个离散高频SNP位点,离散高频SNP位点为等位基因频率大于10%、并且与任意另外一个离散高频SNP位点在参考基因组上的物理距离不小于候选探针长度,候选探针长度为50-250mer。
在本发明的一个具体实施方式中,离散高频SNP是通过千人基因组数据获得的,也可以从其它已公开的基因组数据或者获得的进一步选择等位基因频率小于90%的离散高频SNP位点,确定候选探针长度为100mer。
在本发明的一个具体实施方式中,每个候选探针包含一个离散高频SNP位点,并且离散高频SNP位点位于所说的候选探针的中段。这样每条候选探针只包含一个高频SNP位点,相邻候选探针之间可能有重叠也可能没有重叠。这里的“中段”,是相对于“前段”和“后段”来说的,可以按常规理解,比如一条序列,其上、下游1/3分别定为“前段”和“后段”,中间的1/3为“中段”;进一步的,离散高频SNP位点位于所说的候选探针的中点,这里的“中点”位置,比如一条序列包含2n+1个核苷酸,中点即为第n+1核苷酸的位置,而当一条序列含有2n个核苷酸,序列的中点为第n或第n+1个核苷酸的位置,这样可以增强探针对目标离散高频SNP位点的捕获效率。
在本发明的一个具体实施方式中,基于第一候选探针集中的候选探针序列的GC含量和/或单碱基重复对第一候选探针集进行预筛选,保留了第一候选探针集中的GC含量为35%-65%和/或单碱基重度小于7的候选探针。单碱基重复度是指在一段序列中一个碱基类型连续出现的次数,比如TGAAAAAAAAGC中,其中的A连续出现8次,该序列的A碱基重复度为8。序列GC含量偏高或偏低、高杂合度容易影响该序列的PCR或者杂交捕获过程,带来GC偏向性(GC bias)等,使捕获特异性降低,经此预筛选保留的第一候选探针集将不会与这些序列杂交,从而免除GC bias或低特异性捕获对结果产生的影响。
步骤二:将第一候选探针集与参考序列进行比对,以便获得比对结果
将第一候选探针集与参考序列进行比对,获得比对结果,获得第一候选探针集在参考序列上的位置信息。所使用的参考序列是已知序列,可以是预先获得的目标样本所属生物类别中的任意的参考模板。比如,目标样本是人类的,参考序列可选择美国国家生物技术信息中心(NCBI)提供的HG18或者HG19,进一步的可以预先配置包含更多参考序列的资源库,在进行序列比对前,先依据目标样本的性别、人种、地域等因素选择更接近的参考序列,有利于获得更有针对性的探针序列。
步骤三:对第一候选探针集进行第一筛选,以便获得第二候选探针集
在本发明的一个具体实施方式中,经过第一筛选保留的候选探针需满足以下两个条件中的任一个:1)第一候选探针集中的比对到参考基因组唯一位置的候选探针;2)第一候选探针集中的比对到参考序列多个位置、并且与参考序列多个位置中的至少两个位置的错配比例都小于10%;比如候选探针长度100mer,10个碱基错配即错配比例10%,错配率低用于杂交时能与目标区接近完全互补配对,捕获效果佳,特异性高。
步骤四:将参考序列划分为多个窗口,将第二候选探针集分配至各自匹配的窗口
将参考序列划分成多个具有预定长度的窗口,利用比对,将第二候选探针集中的多个候选探针分配到匹配上的窗口,获得各个候选探针在各自窗口上的位置信息。
多个预定长度的窗口的长度可以一致可以不一致,可以重叠可以不重叠,在本发明的一个具体实施方式中,参考序列为参考基因组,将参考基因组划分为多个一致长度的窗口,窗口长度为10Kb,且相邻两个窗口连接但不重叠。
步骤五:基于所说的位置信息以及离散高频SNP的等位基因频率,对第二候选探针 集进行第二筛选,确定探针序列
在本发明的一个具体实施方式中,进行第二筛选包括两个步骤,(a)如果存在多个候选探针位于同一个窗口,则确定离散高频SNP的等位基因频率最高的候选探针;(b)如果仅存在一个离散高频SNP的等位基因频率最高的候选探针,则选择该离散高频SNP的等位基因频率最高的候选探针作为探针,如果存在多个离散高频SNP的等位基因频率最高的候选探针,则选择多个离散高频SNP的等位基因频率最高的候选探针中距离窗口中心最近的候选探针作为所述探针。候选探针与窗口中心的距离可以是候选探针的中点与该窗口中心的距离。目标位置尽可能处于探针序列的中心位置,有利于提高捕获效率。
在本发明的一个具体实施方式中,对第二候选探针集进行第二筛选之后,当第二候选探针集中的分别落入相邻两个窗口的相邻两条候选探针在参考基因组上的距离大于相邻两窗口中任一窗口的长度时,可选择地,进一步将参考基因组上的位于相邻两条候选探针之间的短串连重复序列或者短串联重复序列的一部分添加到经第二筛选后的第二候选探针集中,一起构成探针序列。这样,利用这些设计获得的探针序列捕获全基因组时,能使捕获得的区域的间距呈现相对均匀的分布,能使捕获确定的区域组合更好的全面反映整个基因组信息。
根据本发明的另一种实施方式,提供一种基因组结构变异的检测方法,所说的基因组结构变异包括染色体非整倍性、拷贝数变异和***缺失的至少之一,包括以下步骤:
(一)对目标样本基因组核酸进行测序,以便获得基因组测序结果,所说的基因组测序结果由多个读段构成,基因组测序结果可以通过全基因测序获得,比如通过提取基因组DNA,依据现有高通量平台的指导手册,比如利用Illumina Hiseq2000/2500、Roche 454、Life technologies Ion Torrent、单分子或纳米孔测序平台等进行文库构建及上机测序获得读段(reads);或者通过探针捕获所述目标样本的基因组并进行测序获得,探针可以通过本发明一方面提供的探针的确定方法进行设计确定,接着按照现有的方法合成或制备而得的。
(二)将参考基因组分为m个区域,利用测序结果中的读段中落入区域i的读段计算目标样本基因组区域i的覆盖深度TDi,其中,m和i为自然数,i为区域编号,1≤i≤m,10<m。
在本发明的一个具体实施方式中,区域i的覆盖深度的计算公式为
Figure BDA0001203857430000061
或者
Figure BDA0001203857430000062
其中,i表示区域的编号。读段落到基因组上位置可以通过序列比对确定,比对可使用各种比对软件,例如SOAP(Short Oligonucleotide Analysis Package),bwa(Burrows-Wheeler Aligner),samtools,GATK(Genome Analysis Toolkit)等。
(三)基于目标样本基因组区域i的覆盖深度与k个参照样本的区域i的覆盖深度的差异程度,判断目标样本区域i的结构变异的发生,其中,k为自然数,k≥2。
在本发明的一个具体实施方式中,目标样本基因组区域i的覆盖深度与k个参照样本的区域i的覆盖深度的差异程度的比较,是通过比较目标样本和参照样本的基因组区域i的覆盖深度系数来实现的,目标样本基因组区域i的覆盖深度系数Ri的确定包括以下步骤,(a)对TDi进行第一校正以获得TDai,第一校正是通过对包含区域i在内的2n个连续区域的覆盖深度值进行线性回归实现的,其中,n为自然数,10<n≤m/2,在本发明的一个具体实施方式中,经第一校正线性回归获得的
Figure BDA0001203857430000063
其中,TDj表示n个连续区域中的第j个区域的覆盖深度,j为自然数,1≤j≤n;(b)在获得区域i的第一校正覆盖深度TDai后,进一步对TDai进行均一化获得
Figure BDA0001203857430000064
进而获得
Figure BDA0001203857430000065
在本发明一个具体实施方式中,对区域i的第一校正覆盖深度TDai进行均一化获得的
Figure BDA0001203857430000066
在本发明的一个具体实施方式中,在获得目标样本的Ri后进一步包括对Ri进行第二校正以获得ri
Figure BDA0001203857430000067
其中,Rai为k个参照样本基因组区域i的覆盖深度系数的平均值,
Figure BDA0001203857430000071
y为自然数表示参照样本编号,Ri,y表示参照样本y基因组区域i的覆盖深度系数。
在本发明的另一个具体实施方式中,在获得目标样本的Ri后进一步包括对Ri进行第二校正以获得ri
Figure BDA0001203857430000072
其中,Rai为k个参照样本和一个目标样本的基因组区域i的覆盖深度系数的平均值,
Figure BDA0001203857430000073
y为自然数表示参照样本编号,Ri,y表示参照样本y基因组区域i的覆盖深度系数。
上述计算处理目标样本基因组区域i的覆盖深度系数Ri的过程中,对中间数值的的校正、均一化等处理能减少因实验条件的波动、样品间本身的差异等带来的误差,使最后的ri能真实反映Ri且围绕1的波动幅度比Ri小,且多个样本的ri符合正态分布;上述实施方式中对TDi进行第一校正,接着对第一校正后的数值进行均一化,相当于两次求均值的过程,即在打算以包含区域i的n个连续区域的覆盖深度均值代表区域i的覆盖深度之前,n个区域中每个区域的覆盖深度值的计算都是利用以该区域为第一个区域的n个连续区域的覆盖深度均值表示的,这样相当于利用包含目标区域i的2n个连续区域的覆盖深度值来校正TDi,能使连续区域的覆盖深度保持稳定。需要说明的是,本领域人员可以利用其它校正或求平均值处理使相邻几个区域的覆盖深度值保持稳定,比如以与目标区域间隔多少个的几个区域的平均覆盖深度来校正目标区域覆盖深度,均属于本发明的构思。参照样本基因组区域i的覆盖深度系数的计算处理可以参考目标样本基因组区域i的覆盖深度系数的计算处理过程,参照样本数据可以预先计算处理好备用,也可以与目标样本的计算处理过程同步进行而获得。
在本发明的一个具体实施方式中,目标样本基因组区域i的覆盖深度与k个参照样本的区域i的覆盖深度的差异程度的判断,是通过t检验二者的覆盖深度系数的差异是否显著来实现的。在本发明的一个具体实施方式中,目标样本基因组区域i的t检验统计量的计算公式为
Figure BDA0001203857430000074
其中,
Figure BDA0001203857430000075
表示k个参照样本的ri,y的平均值,ri,y为参照样本y基因组区域i的经第二校正的覆盖深度系数,
Figure BDA0001203857430000076
为k个参照样本标准差,
Figure BDA0001203857430000077
基于目标样本基因组区域i的ti值,获得显著水平Pi,当Pi<0.05,判定所述区域i发生结构变异;反之,则判定所述区域i不发生结构变异。在本发明的另一个具体实施方式中,基于目标样本基因组区域i的ti值和预先确定的显著水平Pi0,获得ti理论值ti0,当ti≥ti0,判定所述区域i发生结构变异,反之,则判定所述区域i不发生结构变异,预先确定的Pi0≤0.05。根据t检验的t值表,预定Pi0后可查得对应的ti0
在本发明的一个实施方式中,为检测更大的CNV或***缺失,在进行步骤(三)之后,将同方向且连续的W个区域合并,获得一级合并区域,合并两个一级合并区域当两个一级合并区域是同方向的并且之间的跨度不超过L个区域,获得二级合并区域,检测二级合并区域的结构变异;其中,同方向区域指区域的覆盖深度的t统计量都大于0或者都小于0的区域,W和L均为自然数,W≥2,L-W≤1。要进一步检测更大的结构变异,可依次类推,如进一步合并符合条件的二级合并区域,合并条件可类似的为两个二级合并区域同方向且之间的在参考基因组上的距离不超过L个区域或L个二级合并区域。
在本发明的一个具体实施方式中,检测二级合并区域的结构变异,是基于目标样本基因组的所述二级合并区域的覆盖深度与多个参照样本基因组上对应的区域的覆盖深度的差异程度,来判断该二级合并区域是否发生结构变异,或者说来判断发生在区域i的结构变异是否横跨W个区域。参照样本基因组上对应的二级合并区域的覆盖深度的获得、目标样本基因组上的二级合并区域覆盖深度的t统计量的计算及结构变异判断过程可参见前面相对小的区域i的结构变异的计算判断过程。
根据本发明的再一个实施方式,提供一种适用于检测基因组结构变异中的杂合性丢失的的方法,包括以下步骤:
(1)对目标样本基因组核酸进行测序,以便获得基因组测序结果,所说的基因组测序结果由多个读段构成,基因组测序结果可以通过全基因测序获得,比如通过提取基因组DNA,依据现有高通量平台的指导手册,比如利用Illumina Hiseq2000/2500、Roche 454、Life technologies Ion Torrent、单分子或纳米孔测序平台等进行文库构建及上机测序获得读段(reads);或者通过探针捕获所述目标样本的基因组并进行测序获得,探针可以通过本发明一方面提供的探针的确定方法进行设计确定,接着按照现有的方法合成或制备而得的。
(2)将参考基因组分成m’个区域,基于测序结果中落在参考基因组区域i中的读段信息和群体区域i数据,获得目标样本基因组区域i和群体区域i共有的SNP集,分别计算目标样本和群体的共有SNP集中的各个SNP位点所在片段的杂合度,获得目标样本基因组区域i的杂合度集Ui,和群体区域i的杂合度集U0i,比较目标样本Ui和群体U0i以确定目标样本区域i杂合性丢失是否发生;其中,所述共有SNP集中的每个SNP的等位基因频率都大于0.1,所说的共有SNP集中的一个SNP位点所在片段是以与该SNP相邻的上下游两个SNP为边界点的,m’和i为自然数,m’≥i≥1,m’≥6。
在本发明的一个具体实施方式中,一个SNP位点所在片段的杂合度是以该SNP位点的次等位基因频率系数表示的,所述SNP位点的次等位基因频率系数Rhet=MAF/(1-MAF),MAF为该高频SNP的次等位基因频率。
在本发明的一个具体实施方式中,比较目标样本Ui和群体U0i以确定目标样本区域i杂合性丢失是否发生,包括利用F检验判断Ui的方差
Figure BDA0001203857430000091
和U0i的方差
Figure BDA0001203857430000092
是否有显著性差异,,若Ui和U0i的方差差异显著,则判定所述目标样本区域i存在杂合性丢失,反之,则判定所述目标样本区域i没有存在杂合性丢失。
在本发明的一个具体实施方式中,F检验包括分别计算Ui和Ui0的方差,利用所得目标样本Ui的方差
Figure BDA0001203857430000093
和群体Ui0的方差
Figure BDA0001203857430000094
计算获得两个互为倒数的统计量Fupper和Funder,利用互为所说的倒数的统计量获得显著水平pF,比较pF与预定显著水平pF0的大小,pF≤pF0说明两方差差异显著,反之则差异不显著,F检验包含计算公式,
Figure BDA0001203857430000095
Figure BDA0001203857430000096
Figure BDA0001203857430000097
pF=pupper+(1-punder),其中,v为目标样本基因组区域i和群体区域i共有SNP集中SNP的编号,q为目标样本基因组区域i和群体区域i共有SNP集中SNP的个数,Rhet,i,v为目标样本基因组区域i的共有SNP集中的第v个SNP的次等位基因频率系数,
Figure BDA0001203857430000098
为目标样本基因组区域i的共有SNP集中的q个SNP的次等位基因频率系数的平均值,Rhet,i0,v群体样本基因组区域i的共有SNP集中的第v个SNP的次等位基因频率系数,
Figure BDA0001203857430000099
为群体样本基因组区域i的共有SNP集中的q个SNP的次等位基因频率系数的平均值,pupper和punder分别根据Fupper和Funder获得,pF0≤0.05。pF0可以取通常设置的值、或者根据所掌握的已知信息、对检测准确性的要求等调整设置。
在本发明的一个实施方式中,为检测更大的LOH,在步骤(2)之后,将W’个发生杂合性丢失且连续的区域合并,获得三级合并区域,合并两个三级合并区域当所述两个三级合并区域之间的跨度不超过L’个区域时,获得四级合并区域,分别获得目标样本四级合并区域的杂合度集和群体同样区域的杂合度集,比较两个杂合度集,以确定目标样本四级合并区域是否发生杂合性丢失,其中,W’和L’均为自然数,W’≥2,W’/2≥L’。在本发明的一个具体实施方式中,W’≥4。要检测更大区域发生的LOH,可依次类推,比如进一步合并符合条件的四级合并区域,合并条件可类似的为两个四级合并区域之间的在参考基因组上的距离不超过L’个区域或L’个三级合并区域。
根据本发明的再一个实施方式,提供一种检测单亲二倍体的方法,当某目标样本基因组区域存在杂合性丢失时,计算这个区域的拷贝数,当这个区域拷贝数与同物种正常基因组上该区域的拷贝数一样时,判定所述目标样本的这个基因组区域存在UPD;基因组区域是否存在LOH可通过前面本发明披露的一方面的LOH检测方法进行。
本领域普通技术人员可以理解,上述实施方式中各种方法的全部或部分步骤可以通过程序指令相关硬件完成,该程序可以存储于一计算机可读存储介质中,存储介质可以包括:只读存储器、随机存储器、磁盘或光盘等。
根据本发明的最后一个实施方式,还提供一种检测基因组结构变异的装置,包括:数据输入单元,用于输入数据;数据输出单元,用于输出数据;存储单元,用于存储数据,其中包括可执行的程序;处理器,与上述数据输入单元、数据输出单元及存储单元数据连接,用于执行存储单元中存储的可执行的程序,程序的执行包括完成上述实施方式中各种方法的全部或部分步骤。
以下结合具体目标个体对依据本发明的具体探针设计方法及结构变异检测方法的运行结果进行详细的描述。下述过程涉及的名称定义或具体参数设置选择为:
1、将设计的探针称为选择目标区域探针(Seleted Target Region Primers,SeTR);
2、下文中的“覆盖深度”、“测序深度”和“深度”,可替换使用;下文中的“区域”和“目标区域”可替换使用;
3、文库构建、测序依据Hiseq 2000平台提供的小片段文库构建操作说明及上机测序说明来操作,文库的大小为300bp-350bp,双端测序(pair-end sequencing),读段长91bp(测序类型为PE91+8+91);
4、比对选择的参考基因组或参考序列为人类参考基因组(hg19,Build 37)。
实施例中未注明具体技术或条件的,按照本领域内的文献所描述的技术或条件(例如参考J.萨姆布鲁克等著,黄培堂等译的《分子克隆实验指南》,第三版,科学出版社)或者按照产品说明书进行。所用试剂或仪器未注明生产厂商者,均为可以通过市购获得的常规产品,例如可以采购自Illumina公司。
实施例1:芯片设计、制备、测试
通常,高(>60%)或者低(<35%)GC含量和高杂合度容易给其DNA片段在PCR或者探针捕获过程中是带来不利的影响,为了避免此种现象,我们设计了特殊的探针,我们将其称为SeTR.在设计SeTR探针的时候,遵循以下几个原则:a)探针序列的唯一性和稳定性较高,要求具有低杂合性和中等的GC(35%~60%)含量,b)含有离散型的高频SNP标记(SNPmarker),各SNP的等位基因频率(allele frequency,0.9>AF>0.1)以便更好检测全基因组的LOH,c)最终的目标区域呈现出相对均匀的分布。
SeTR探针设计或者说目标区域的选择流程如下:
1)基于千人基因数据库(ftp://ftp.ncbi.nih.gov/1000genomes/ftp/release),挑选出等位基因频率(Allele frequency,AF)为10%~90%的候选SNP集,然后再在SNP集中去掉两个SNP之间物理距离小于100pb的一个SNP,从而构成SNP maker1集。
2)以SNP maker1集的每一个SNP为中点,在其上下游各截取参考基因组序列50pb,构成100bp的理论探针序列集,然后将此探针序列集比回到参考序列上。如果某一探针序列的最佳比对没有错配,且其次佳比对也只有小于5%的错配,那么其对应的SNP则被保留,从而构成SNP maker2集。
3)基于SNP maker2集,我们挑选出在参考基因组物理位置上均匀的SNPmaker为最终的SNP maker集。在我们的研究中,我们选择了物理距离大约为10Kbp的SNP maker集。
4)如果在最终的SNP maker集中,有两个临近的SNP之间的距离大于了10Kbp,在选用他们之间的短重复序列(short tandem repeat,STR)来填补均匀。
设计完SeTR探针后,我们委托Roche来产出SeTR液相芯片。SeTR液相芯片含有278800个探针,总大小为41,795,106bp,其覆盖了有效的全基因组(2.89G)1.45%的区域。SeTR探针平均长度达到了149bp,相邻两两探针之间的平均物理距离为10.6kbp,如表1和图1所示。
表1 SeTR探针在每条染色体上的分布
Figure BDA0001203857430000111
Figure BDA0001203857430000121
用3个质检合格的DNA样品,YH(炎黄样本,中国人基因组DNA)、HG00537(千人基因组项目中的一个样本)和GM50275(获自柯瑞尔医学研究所Coriell Institute forMedical Research的人成纤维母细胞样本),来测试SeTR芯片的可用性,以保证此探针芯片能用于后续的检测研究。三个样本都利用SeTR捕获建库测序,获得测序序列(reads)。首先我们去掉被接头(adapter)污染和质量较低的比如平均质量值低于20的reads后,称剩下的reads为干净reads(clean reads),将干净reads比对到参考序列hg19上,得到了98.13%~99.29的reads比对到参考基因组上,其中比对到目标区域的达到了67.43%~67.87%,此外,有99.73%~99.95的目标区域至少被一条reads覆盖了,有超过99%的区域至少被覆盖到了10次,如表2所示,这些表现都要优于同类型的外显子组捕获(exome capture)芯片,比如Roche Nimblegen公司生产的外显子组液相芯片。此外,目标区域的深度分布,如图4A所示,类似于泊松分布(Poisson distribution),图4B显示目标区域内的绝大部分的高杂合位点的非参考序列碱基型(the non-reference allele)的reads支持数与参考序列碱基型(reference allele)的reads支持数几乎相当,即比对时某高杂合位点获得的正负reads支持数相当(正负reads分别来源两条同源染色体),这些都显示此探针无明显的单倍体型(常见为参考序列碱基型,即ref型)捕获的偏向性,以及对目标区域捕获均一性较优。
表2 三个样品的比对结果
Figure BDA0001203857430000122
Figure BDA0001203857430000131
实施例2:目标区域文库构建、测序
1、试验材料、试剂、仪器
样本:15例目标gDNA样本(人基因组DNA,样本编号见以下表3,“GM”“开头的都是人成纤维母细胞),24个参照DNA样本。
主要试剂仪器:PCR仪、移液器、离心机、舒适型恒温混匀仪、DNA打断仪、涡旋振荡器、磁力架、电泳仪、Hiseq2000测序仪、Nanodrop紫外分光光度计等,所用试剂或仪器未注明生产厂商者,均为可以通过市购获得的常规产品。
探针设计及合成:通过实施例一获得,在人的全基因组范围内,选取约41M的目标区域,从罗氏公司(Roche)定制NimbleGen SeqCap EZ液相探针,该探针集能够捕获对应的所设计的目标区域。
2、文库构建
1)基因组DNA提取
使用QIAGEN DNA提取试剂盒(DNA Mini Kit),并按照试剂盒说明书,从目标样本中提取基因组DNA约3-5μg,用于后续实验。将提取好的DNA(30-100ng)跑电泳检测,判断是否完整以及降解程度。
2)基因组DNA打断及纯化
使用covaris E210仪器对基因组DNA进行打断(参照仪器使用说明进行操作)。将DNA打断成200-250bp。使用QIAquick PCR Purification kit(250)试剂盒,按照试剂盒说明书操作,对打断后的DNA片段进行纯化,电泳检测主带大小是否符合要求,即主带大小是否为200-250bp。
3)末端修复、末端加A、加接头、预扩增
按建库要求,按双末端标签文库构建说明书步骤及其列明的试剂、反应条件等,对上述断裂纯化后的DNA片段进行末端修复,并进行纯化;加个碱基A于经末端修复纯化后的DNA片段的两端,纯化末端加A产物;在末端加A产品的两端连接测序接头,并利用能与测序接头互补结合的磁珠纯化带接头的DNA片段。配制PCR反应体系,扩增带接头的DNA片断,磁珠纯化PCR产物,电泳检测扩增产物主带大小是否在300-350bp;用Nanodrop紫外分光光度计检测DNA量,总量需大于1.0μg。
4)SeTR探针杂交及洗脱,扩增
依照市售的NimbleGen SeqCap EZ杂交洗脱试剂盒说明书进行,购买或配置试剂盒说明书中的杂交、洗脱相关试剂。准备1.5mL离心管,加入Cot-1DNA,通用封闭序列(Block),标签的封闭序列(index N Block)和经步骤3)后的DNA样品。然后离心1min,60℃真空浓缩干燥,然后加入杂交缓冲液等,震荡离心,放到95℃的金属干浴锅里变性10min,震荡后高速离心。在离心管中加入4.5ul探针,在PCR仪上杂交(47℃,64-72hours)。杂交完成后进行洗脱。然后按照文库构建说明书最后的扩增步骤进行PCR,按要求配制PCR反应体系,将杂交洗脱获得的DNA,聚合酶、底物、PCR反应缓冲液,Flowcell引物(依测序仪的测序芯片flowcell上带有的固定序列设计的引物)等反应物混合均匀。PCR程序为94℃预变性2min,94℃变性15s,58℃退火30s,72℃延伸30s,反应15个循环后,再72℃延伸5min。PCR完成后,取出PCR产物,离心,磁珠纯化,获得目标区域文库。用Nanodrop紫外分光光度计测文库的浓度,准备上机测序。
3、Hiseq2000高通量测序
质检合格的DNA文库,根据Hiseq2000操作说明进行上机测序。每个样本的数据量约4.5G,平均测序深度达到100X,但由于捕获芯片的效率很难达到100%,通过分析,最终的目标区域的有效测序深度为30X~45X。
实施例3:CNV、LOH和UPD的检测
总体流程参见图3。测序完成之后,下机数据为fastq文件格式。然后将过滤后得到高质量的reads与参考基因组(Hg19,Build 37)采用BWA软件进行比对得到SAM格式的比对文件,之后使用samtools软件将SAM比对文件格式转换成二进制的BAM文件,并对比对结果进行去重复和排序处理,接着,将再次使用samtools软件,将BAM格式转换成PILEUP格式具体详情请见生物信息分析策略部分。
一、测序数据过滤、比对
先将上述实施例illumina Hiseq2000下机的测序数据进行简单的数据过滤,将被adapter污染,含N比例高于5%,平均质量值低于Q20的reads进行去除。然后使用bwa比对软件将过滤后的数据比对到人类参考基因组上(hg19,Build 37),输出序列比对结果即SAM(sequence alignment/map)格式的比对文件(简称SAM文件),接着使用Samtools软件将SAM文件转换成二进制的BAM文件、去除掉PCR引起的重复(PCR duplicates)和进行排序处理,使用GATK软件对比对结果进行重比对和重校正。
二、计算目标区域的经第二校正的覆盖深度系数ri和片段的杂合度Rhet
根据上述过滤比对后获得的探针区域文件包含的信息计算出每个目标区域的ri和片段杂合度Rhet值。根据ri值,采用t检验预测CNV,根据RHet,采用F检验预测LOH和UPD。
三、检测CNV,LOH和UPD的分析
1、CNV检测
1.1计算每个目标区域的深度系数(Ri)
计算目标区域的深度,并用TDi表示(如公式1),为了保持连续的几个目标区域TD的稳定性,采用了公式2的方法来校正TDi,即利用第i区域后面的n‘个区域的深度来校正TDi,得到TDai,然后利用公式3和4对TDai进行均一化,此时得到每个目标区域的深度系数Ri
公式1:TDi=Tibase/Tilen,
公式2:
Figure BDA0001203857430000151
公式3:
Figure BDA0001203857430000152
公式4:
Figure BDA0001203857430000153
Tibase:比对到目标区域i的碱基数;Tilen:目标区域i的长度。
1.2利用多个参照样本(k=24)数据创建基准线,校正Ri获得ri
由于每次实验条件的波动和样品间本身的差异导致每次捕获的效率也存在一定的波动,进而引起Ri的波动,容易导致出现CNV假信号。因此,根据多个样的波动情况,创建统一的一个基准线是非常有益的。图4很好的体现出创建基准线利于这个检测,前体Ri(preRi)
Figure BDA0001203857430000154
的分布如图波动很大,而Ri波动相对小些,Ri通过基准线的校正后得到ri,其波动更小,更敏感更易检测CNV的发生。理论上,认为不同的样品中,不发生CNV的情况下,在同一个目标区域内,Ri值理论上是符合泊松分布的,并且都围绕各自特有的值相对稳定的上下波动,为了保持各自特有值的稳定性,通过调查多个样品的同一区域的Ri值,采用Ri平均值(mean Ri)来代替这个各自特有的值,为每个目标区域构建一个各自特有的基准线(robust baseline)。基于每个目标区域的Ri是围绕着mean Ri值上下波动的假设,我们将Ri除以mean Ri转化成ri,进而使得ri围绕1上下波动的正态分布。
1.3检测目标区域的CNV
理论上,来自多个样品的同一目标区域的ri值都应符合正态分布,因此调查某个样品的目标区域i时,可以通过比较多个样品此区域的ri值,利用T检验,t统计量的计算公式如下,来检测这个样品目标区域i的拷贝数变化情况。
Figure BDA0001203857430000155
公式中各参数下标中的1代表目标样本,2代表多个参照样本,
Figure BDA0001203857430000156
表示n1个待测样本的ri的平均数,
Figure BDA0001203857430000161
为n2个参照样本的ri的平均数,μ1为理论上所有待测样本的ri平均数,μ2理论上全部参照样本ri2平均数,S1和S2分别为待测样本和参照样本的标准差,df为自由度,df=n1+n2-2。
当待测样本为1即n1=1,待测样本和参照样本理论均值相同,上式化简为:
Figure BDA0001203857430000162
通过上面的简化公式,每个目标区域都对应一个可检测CNV的t值,进而得到P值(置信度),当某区域的P<0.05的时候,此区域则为一个发生CNV的区域。
1.4检测大CNV
基于单个区域t检验的p值,为每个区域附上一个伪信号值来表征是否被下一步CNV区域连接所考虑,再沿着染色体,将可能具有一致CNV的目标区域连接成块,从而确定CNV最终的大小及拷贝数。
伪信号值的标记规则为,当至少四个连续目标区域的测量值同方向(t值同时大于或者同时小于0)即偏离参照样品的相应区域时,若有3个区域的P值小于第一阈值(如0.05,常用的显著水平阈值),而且第四个不超过第二阈值(0.2,第一阈值的四倍),则四个区域均标记为偏离方向(比如偏大标为+,偏小标为-),合并成一个块;这里连续且同方向的区域个数和第一、第二阈值数值都是可调整的。如若一个块与另一个块之间的距离不超过5个区域的跨度,则这两个块进行合并为一个大块,依此类推,最后获得区块;参考前面1.3的方法公式,这个区块的r值以其所包含的所有区域的ri的平均值表示,对待测样本和参照样本的该区块域的r值进行t检验,计算该区块的P值。当该区块的P<0.05,此区块发生CNV,从而确定该区块的边界与大小,获得大CNV的边界和大小。
通过对目标15例样品的分析,我们得到的CNV结果与已知的验证结果(SNP-array结果)高度一致,并且不存在假阳性和假阴性,见表3。再者,我们模拟了8个30X的全基因组数据,其中包括5个正常样品,3个含有CNV的样品,通过对这8个模拟数据进行CNV检测分析,比较了当前已报到的exome区域CNV预测软件CONTRA(Li J,Lupat R,et al,CONTRA:copynumber analysis for targeted resequencing,Bioinformatics.2012 May 15;28(10):1307-13),结果显示,我们的方法敏感度和特异性均达到了100%,且各自的拷贝数也被精确的检测出来,对CNV的检测精度可达到500Kb且能精确定位,而CONTRA的敏感度为88.9%,特异性只为66.7%,拷贝数未给出,如表4所示。
表3
Figure BDA0001203857430000163
Figure BDA0001203857430000171
表4
Figure BDA0001203857430000172
Figure BDA0001203857430000181
2、LOH检测
2.1全基因组各区域的杂合状态检测
在待测样本基因组某区域内,找出在千人数据中基因频率(AF)为0.1~0.9的SNP位点,并按以下公式计算出千人中和待测样本中的这些SNP位点所在区域的RHet值。当待测样本中区域i为绝对杂合状态时,则Rhet=1,反之,为绝对纯合的时候,Rhet=0。
Rhet=MAF/(1-MAF),MAF(minor allele frequency)为次等位基因频率。
在待测样品中,以某区域内任意一个SNP位点m作为起始点,向后连续取n个SNP位点作为该区域内的杂合度集Sm,即Sm={Rhet,sm,Rhet,s(m+1),...,Rhet,s(m+n)},以同样的方式,在千人数据库中,取相同位置的SNP位点,构成杂合度集Pm,即Pm={Rhet,phetm,Rhet,p(m+1),...,Rhet,p(m+n)},F检验两个杂合度集的方差是否相等,具体的,分别计算待测样本该区域的杂合度集的方差
Figure BDA0001203857430000182
和千人样本的同样区域的杂合度集的方差
Figure BDA0001203857430000183
以及待测样本该区域杂合度集Sm的p值。
Figure BDA0001203857430000191
Figure BDA0001203857430000192
Figure BDA0001203857430000193
H0s=σp
HAs≠σp
Figure BDA0001203857430000194
Figure BDA0001203857430000195
p=pupper+(1-punder)
当p≤0.01的时候,我们接受HA,判断杂合度集Sm失去了群体中的杂合性,即集合Sm所在区域发生LOH。
2.2检测大的LOH
结合2.1的结果,采用检测大CNV步骤1.4的方式,记录连续的4个失去杂合状态的子集为一个最小单元。如若两个单元之间不超过2个子集跨度,则将两个单元合并成更大的单元,依此类推,最后连接成区块,此时,再根据待测样品和千人参考集之间的RHet值进行F检验,计算区块的p值,当p≦0.01的时候,我们则认为此区块发生LOH,否则为非LOH区块。
或者,为更准确地检测对合并条件可设置更严格,如为避免某些随机误差导致的假阳性,定义至少大于5M的区域才可能为一个真实的LOH。在此基础上,设置区块容错为1(即允许区块中1个子集的p值大于0.01)的条件下,将在2.1中p≤0.01的子集附近的p≤0.01的连续的子集与之合并。最后,对合并以后的区域内的RHet再进行了一次F检验,若其p值小于0.01,则认为该区块是一个真实的LOH。
3、UPD检测
结合上述全基因组的CNV和LOH检测结果,根据孟德尔遗传规律,进行UPD检测。
如果某一DNA区域在千人数据中显示为杂合状态,即RHet=1,而在实际检测中,其杂合状态消失,即RHet趋近于0,则判定此区域为发生了LOH,而如果在这个区域同时发生有CNV且有两个拷贝(CN=2),即拷贝数没有发生变化(本实施例的样本是二倍体样本,正常二倍体样本基因组各区域都是两个拷贝),则判定此区域发生了单亲二倍体(UPD)。
在15个样品的13个中,检测出了10个大于5M的LOH和4个大于5M的UPD,结果请见表5,LOH和UPD的检测在没有配对样本的情况下(一般是拿自身病变的组织和正常的组织进行比较的,这是配对样本,即有某种关联的样本,而本实施方式方检测LOH和UPD是把目标样本和多个参照样本集合做比较的,目标样本和参照样本集合没有相关性,所以不是配对样本),≥5M的LOH检测结果与CN=1的CNV结果一致(可利用CNV检测结果验证LOH检测结果的准确性),本发明方案检测LOH、UPD的准确性高,且可达到5M级别的精度。
Circos图(图5)综合展示了GM50275样本的CNV、LOH和UPD检测结果。
表5
Figure BDA0001203857430000201
工业实用性
本发明的基于参考序列确定探针序列的方法,能够有效用于确定探针序列,并且获得的探针,用于杂交捕获基因组获得多个基因组局部区域,捕获得的多个局部区域能够代表全基因组、能够反映全基因组变异信息,用于发现全基因范围的结构变异的发生。
尽管本发明的具体实施方式已经得到详细的描述,本领域技术人员将会理解。根据已经公开的所有教导,可以对那些细节进行各种修改和替换,这些改变均在本发明的保护范围之内。本发明的全部范围由所附权利要求及其任何等同物给出。
在本说明书的描述中,参考术语“一个实施例”、“一些实施例”、“示意性实施例”、“示例”、“具体示例”、或“一些示例”等的描述意指结合该实施例或示例描述的具体特征、结构、材料或者特点包含于本发明的至少一个实施例或示例中。在本说明书中,对上述术语的示意性表述不一定指的是相同的实施例或示例。而且,描述的具体特征、结构、材料或者特点可以在任何的一个或多个实施例或示例中以合适的方式结合。

Claims (32)

1.一种基于参考序列确定探针序列的方法,其特征在于,包括:
(1)基于多个离散高频SNP位点,构建第一候选探针集,其中,所述第一候选探针集由多个候选探针构成,其中,所述多个候选探针中的每一个均对应至少一个所述的离散高频SNP位点,所述多个离散高频SNP位点的每一个的等位基因频率分别为至少10%;
(2)将所述第一候选探针集中的所述多个候选探针与参考序列进行比对,以便获得比对结果;
(3)基于所述比对结果,对所述第一候选探针集进行第一筛选,以便获得由多个候选探针构成的第二候选探针集,
其中,所述第一筛选包括保留满足下列条件至少之一的候选探针:
与所述参考序列唯一比对的候选探针;
比对到所述参考序列的多个位置,并且所述多个位置中的至少两个位置的错配比例均小于10%的候选探针;
(4)将所述参考序列划分为多个分别具有预定长度的窗口,分别将所述第二候选探针集中的多个候选探针分配至各自匹配的窗口,以便确定所述多个候选探针各自的位置信息;
(5)基于所述位置信息以及所述离散高频SNP位点的等位基因频率,对所述第二候选探针集进行第二筛选,以便确定所述探针序列,
其中,按照下列步骤,确定所述探针:
(a)如果存在多个所述候选探针位于同一个窗口,则确定对应的离散高频SNP位点的等位基因频率最高的候选探针;
(b)如果同一个窗口仅存在一个对应的离散高频SNP位点的等位基因频率最高的候选探针,则选择该对应的离散高频SNP位点的等位基因频率最高的候选探针作为所述探针,如果同一个窗口存在多个对应的离散高频SNP位点的等位基因频率最高的候选探针,则选择所述多个对应的离散高频SNP位点的等位基因频率最高的候选探针中距离所述窗口中心最近的候选探针作为所述探针。
2.根据权利要求1所述的方法,其特征在于,所述多个离散高频SNP位点的每一个的等位基因频率分别为不超过90%。
3.根据权利要求1所述的方法,其特征在于,所述多个离散高频SNP位点中任意两个相邻离散高频SNP位点在所述参考序列上的物理距离不小于所述候选探针的长度。
4.根据权利要求1所述的方法,其特征在于,所述候选探针的长度为50~250mer。
5.根据权利要求4所述的方法,其特征在于,所述候选探针的长度为100mer。
6.根据权利要求1所述的方法,其特征在于,所述候选探针对应一个所述离散高频SNP位点,并且所述离散高频SNP位点对应于所述候选探针的中段。
7.根据权利要求6所述的方法,其特征在于,所述离散高频SNP位点对应于所述候选探针的中点。
8.根据权利要求1所述的方法,其特征在于,所述候选探针是从所述参考序列截取的。
9.根据权利要求1所述的方法,其特征在于,在进行所述比对之前,基于所述候选探针的GC含量以及单碱基重复数的至少之一,预先对所述第一候选探针集进行预筛选;
所述预筛选包括保留满足下列至少之一的候选探针:
GC含量为35%-65%;以及
单碱基重度小于7。
10.根据权利要求1所述的方法,其特征在于,在步骤(4)中,将所述参考序列划分为多个分别具有相同预定长度的窗口。
11.根据权利要求10所述的方法,其特征在于,将所述参考序列划分为多个长度为10Kb的窗口。
12.根据权利要求1所述的方法,其特征在于,对第二候选探针集进行第二筛选之后,当第二候选探针集中的分别落入相邻两个窗口的相邻两条候选探针在参考基因组上的距离大于相邻两窗口中任一窗口的长度时,进一步将参考基因组上的位于相邻两条候选探针之间的短串连重复序列或者短串联重复序列的一部分添加到经第二筛选后的第二候选探针集中,一起构成探针序列。
13.根据权利要求1所述的方法,其特征在于,所述参考序列为参考基因组或其一部分。
14.一种检测基因组结构变异的方法,所述基因组结构变异包括染色体非整倍性、拷贝数变异和***缺失的至少之一,所述方法用于非诊断目的,其特征在于,所述方法包括,
(1)对目标样本基因组核酸进行测序,以便获得基因组测序结果,所述基因组测序结果由多个读段构成,其中,所述测序包括采用探针进行筛选,其中,所述探针是通过权利要求1~13任一项所述的方法获得的;
(2)将参考基因组分为m个区域,利用落入区域i的所述读段的数目,计算区域i的覆盖深度TDi,m和i为自然数,i表示区域的编号,1≤i≤m,10<m;
(3)基于所述区域i的覆盖深度与k个参照样本的区域i的覆盖深度的差异程度,确定所述区域i是否存在结构变异,其中,k为自然数,k≥2。
15.根据权利要求14所述的方法,其特征在于,利用下列公式确定所述区域i的覆盖深度:
Figure FDA0003082386560000021
Figure FDA0003082386560000022
其中,i表示区域的编号。
16.根据权利要求14所述的方法,其特征在于,所述目标样本基因组区域i的覆盖深度与k个参照样本的区域i的覆盖深度的差异程度的检验,是通过t检验进行的。
17.根据权利要求14所述的方法,其特征在于,所述区域i的覆盖深度与k个参照样本的区域i的覆盖深度的差异程度的比较,是通过比较目标样本和参照样本的基因组区域i的覆盖深度系数进行的,其中,所述区域i的覆盖深度系数Ri的确定包括以下步骤,
(a)对TDi进行第一校正以获得第一校正覆盖深度TDai,所述第一校正是通过对包含区域i在内的2n个连续区域的覆盖深度值进行线性回归实现的,其中,n为自然数,10<n≤m/2;
(b)对TDai进行均一化获得
Figure FDA0003082386560000031
进而获得
Figure FDA0003082386560000032
18.根据权利要求17所述的方法,其特征在于,在步骤(a)中,基于下列公式,确定第一校正覆盖深度TDai
Figure FDA0003082386560000033
其中,TDj表示所述n个连续区域中的第j个区域的覆盖深度,j为自然数,1≤j≤n。
19.根据权利要求18所述的方法,其特征在于,在步骤(b)中,基于下列公式,对TDai进行均一化获得
Figure FDA0003082386560000034
Figure FDA0003082386560000035
20.根据权利要求16~19任一所述的方法,其特征在于,在获得目标样本的Ri后进一步包括对Ri进行第二校正以获得ri
Figure FDA0003082386560000036
其中,Rai为k个参照样本基因组区域i的覆盖深度系数的平均值,
Figure FDA0003082386560000037
y为自然数表示参照样本编号,Ri,y表示参照样本y基因组区域i的覆盖深度系数。
21.根据权利要求16~19任一所述的方法,其特征在于,在获得目标样本的Ri后进一步包括对Ri进行第二校正以获得ri
Figure FDA0003082386560000038
其中,Rai为k个参照样本和一个目标样本的基因组区域i的覆盖深度系数的平均值,
Figure FDA0003082386560000039
y为自然数表示参照样本编号,Ri,y表示参照样本y基因组区域i的覆盖深度系数。
22.根据权利要求20所述的方法,其特征在于,进行所述t检验,目标样本基因组区域i的t统计量的计算公式为
Figure FDA00030823865600000310
其中,
Figure FDA00030823865600000311
表示k个参照样本的ri,y的平均值,ri,y为参照样本y基因组区域i的经所述第二校正的覆盖深度系数,
Figure FDA00030823865600000312
S为k个参照样本标准差,
Figure FDA0003082386560000041
23.根据权利要求22所述的方法,其特征在于,基于目标样本基因组区域i的ti值,获得显著水平Pi,当Pi<0.05,判定所述区域i存在结构变异;反之,则判定所述区域i不存在结构变异。
24.根据权利要求22所述的方法,其特征在于,基于目标样本基因组区域i的ti值和预先确定的显著水平Pi0,获得ti理论值ti0,当ti≥ti0,判定所述区域i存在结构变异,反之,则判定所述区域i不存在结构变异;所述预先确定的Pi0≤0.05。
25.根据权利要求14~19任一所述的方法,其特征在于,在进行步骤(3)之后,将同方向且连续的W个区域合并,获得一级合并区域,当所述两个一级合并区域是同方向的并且之间的跨度不超过L个区域,合并两个一级合并区域,获得二级合并区域,基于目标样本基因组的所述二级合并区域的覆盖深度与多个参照样本基因组上对应的区域的覆盖深度的差异程度,来检测二级合并区域的结构变异;其中,同方向区域指区域的t统计量都大于0或者都小于0的区域,W和L均为自然数,W≥2,L-W≤1。
26.一种检测杂合性丢失的方法,所述方法用于非诊断目的,其特征在于,包括,
(1)对目标样本基因组核酸进行测序,以便获得基因组测序结果,所述基因组测序结果由多个读段构成,其中,所述测序包括采用探针进行筛选,其中,所述探针是通过权利要求1~13任一项所述的方法获得的;
(2)将参考基因组分成m’个区域,基于所述基因组测序结果中落在区域i中的读段信息和群体区域i的数据,获得目标样本基因组区域i和群体区域i共有的SNP位点构成共有SNP集,分别计算目标样本和群体的共有SNP集中的各个SNP位点所在片段的杂合度,获得目标样本基因组区域i的杂合度集Ui,和群体区域i的杂合度集U0i,比较目标样本Ui和群体U0i以确定目标样本区域i是否存在杂合性丢失;其中,所述SNP位点所在片段是以与该SNP相邻的上下游两个SNP为边界点的,m’和i为自然数,m’≥i≥1,m’≥6。
27.根据权利要求26所述的方法,其特征在于,所述共有SNP集中的每个SNP的等位基因频率都大于0.1。
28.根据权利要求26所述的方法,其特征在于,所述SNP位点所在片段的杂合度是以该SNP位点的次等位基因频率系数表示的,所述SNP位点的次等位基因频率系数Rhet=MAF/(1-MAF),MAF为该SNP的次等位基因频率。
29.根据权利要求28所述的方法,其特征在于,所述比较目标样本Ui和群体U0i以确定目标样本区域i杂合性丢失是否发生,包括利用F检验判断Ui的方差
Figure FDA0003082386560000043
和U0i的方差
Figure FDA0003082386560000042
是否有显著差异,若Ui和U0i的方差差异显著,则判定所述目标样本区域i存在杂合性丢失,反之,则判定所述目标样本区域i不存在杂合性丢失。
30.根据权利要求29所述的方法,其特征在于,所述F检验包括分别计算Ui和Ui0的方差,利用所得目标样本Ui的方差
Figure FDA0003082386560000051
和群体Ui0的方差
Figure FDA0003082386560000052
计算获得两个互为倒数的统计量Fupper和Funder,利用所述互为倒数的统计量获得显著水平pF,比较pF与预定显著水平pF0的大小,包含计算公式,
Figure FDA0003082386560000053
Figure FDA0003082386560000054
Figure FDA0003082386560000055
pF=pupper+(1-punder),其中,v为目标样本基因组区域i和群体区域i共有的高频SNP集中SNP的编号,q为目标样本基因组区域i和群体区域i共有的高频SNP集中SNP的个数,Rhet,i,v为目标样本基因组区域i的共有高频SNP集中的第v个SNP的次等位基因频率系数,
Figure FDA0003082386560000056
为目标样本基因组区域i的共有高频SNP集中的q个SNP的次等位基因频率系数的平均值,Rhet,i0,v为群体样本基因组区域i的共有高频SNP集中的第v个SNP的次等位基因频率系数,
Figure FDA0003082386560000057
为群体样本基因组区域i的共有高频SNP集中的q个SNP的次等位基因频率系数的平均值,pupper和punder分别根据Fupper和Funder获得,pF0≤0.05。
31.根据权利要求26~30任一所述的方法,其特征在于,在步骤(2)之后,将W’个发生杂合性丢失且连续的区域合并,获得三级合并区域,当所述两个三级合并区域之间的跨度不超过L’个区域时,合并两个三级合并区域,获得四级合并区域,分别获得目标样本四级合并区域的杂合度集和群体同样区域的杂合度集,比较两个杂合度集,以确定目标样本四级合并区域是否发生杂合性丢失,其中,W’和L’均为自然数,W’≥2,W’/2≥L’。
32.一种检测单亲二倍体的方法,所述方法用于非诊断目的,其特征在于,当检测目标样本基因组区域存在杂合性丢失时,计算这个基因组区域的拷贝数,当这个基因组区域的拷贝数与同物种正常基因组该区域的拷贝数一样时,判定所述目标样本基因组区域为单亲二倍体;目标样本基因组区域杂合性丢失的确定是通过权利要求26~30任一所述方法进行的。
CN201480080426.0A 2014-07-04 2014-07-04 确定探针序列的方法和基因组结构变异的检测方法 Active CN106715711B (zh)

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
PCT/CN2014/081686 WO2016000267A1 (zh) 2014-07-04 2014-07-04 确定探针序列的方法和基因组结构变异的检测方法

Publications (2)

Publication Number Publication Date
CN106715711A CN106715711A (zh) 2017-05-24
CN106715711B true CN106715711B (zh) 2021-09-17

Family

ID=55018343

Family Applications (1)

Application Number Title Priority Date Filing Date
CN201480080426.0A Active CN106715711B (zh) 2014-07-04 2014-07-04 确定探针序列的方法和基因组结构变异的检测方法

Country Status (2)

Country Link
CN (1) CN106715711B (zh)
WO (1) WO2016000267A1 (zh)

Families Citing this family (22)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
WO2018214010A1 (zh) * 2017-05-23 2018-11-29 深圳华大基因研究院 一种基于测序数据的变异检测方法、装置和存储介质
CN112739828B (zh) * 2018-06-11 2024-04-09 深圳华大生命科学研究院 确定待测样本类型的方法及***
CN110872618B (zh) * 2018-09-04 2022-04-19 北京果壳生物科技有限公司 一种基于Illumina人全基因组SNP芯片数据判断被检样本性别的方法及用途
CN109584963A (zh) * 2018-09-30 2019-04-05 南京派森诺基因科技有限公司 一种高通量测序数据的多样化抽取方法
CN111383714B (zh) * 2018-12-29 2023-07-28 安诺优达基因科技(北京)有限公司 模拟目标疾病仿真测序文库的方法及其应用
CN110079589A (zh) * 2019-05-21 2019-08-02 中国农业科学院农业基因组研究所 一种精准获得全基因组范围内结构变异的方法
CN110600078B (zh) * 2019-08-23 2022-03-18 北京百迈客生物科技有限公司 一种基于纳米孔测序检测基因组结构变异的方法
CN110592208B (zh) * 2019-10-08 2022-05-03 北京诺禾致源科技股份有限公司 地中海贫血症三类亚型的捕获探针组合物及其应用方法和应用装置
CN112662767B (zh) * 2020-11-25 2021-08-06 深圳华大基因股份有限公司 用于衡量基因组不稳定性的试剂盒、探针及其应用
CN112522382B (zh) * 2020-12-22 2024-03-22 广州深晓基因科技有限公司 一种基于液相探针捕获的y染色体测序方法
CN112885410B (zh) * 2021-01-28 2022-09-09 陈晓熠 用于cnv结构变异检测的基因分型芯片
CN113593644B (zh) * 2021-06-29 2024-03-26 广东博奥医学检验所有限公司 基于家系的低深度测序检测染色体单亲二体的方法
WO2023030233A1 (zh) * 2021-08-30 2023-03-09 广州燃石医学检验所有限公司 一种拷贝数变异的检测方法及其应用
CN113971986B (zh) * 2021-10-12 2023-03-21 江苏先声医疗器械有限公司 一种通过序列相似性排查测序样本交叉污染的方法
CN114220481B (zh) * 2021-11-25 2023-09-08 深圳思勤医疗科技有限公司 基于全基因组测序完成待测样本的核型分析的方法、***和计算机可读介质
CN114678067B (zh) * 2022-03-21 2023-03-14 纳昂达(南京)生物科技有限公司 构建多人群非外显子区snp探针集合的方法及装置
CN114582427B (zh) * 2022-03-22 2023-04-07 成都基因汇科技有限公司 一种渐渗区段鉴定方法及计算机可读存储介质
CN115101128B (zh) * 2022-06-29 2023-09-15 纳昂达(南京)生物科技有限公司 一种杂交捕获探针脱靶危险性评估的方法
CN115713971B (zh) * 2022-09-28 2024-01-23 上海睿璟生物科技有限公司 靶向序列捕获探针设计策略选择方法、***及终端
CN115631789B (zh) * 2022-10-25 2023-08-15 哈尔滨工业大学 一种基于泛基因组的群体联合变异检测方法
CN115713967B (zh) * 2022-11-17 2023-08-29 纳昂达(南京)生物科技有限公司 探针池的设计方法及相关装置
CN116144794B (zh) * 2023-03-09 2023-12-19 华中农业大学 牛12k sv液相芯片及其设计方法和应用

Citations (2)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN101213312A (zh) * 2005-06-30 2008-07-02 先正达参股股份有限公司 用于筛选基因特异杂交多态性(gshp)的方法及其在遗传作图和标记开发中的用途
WO2014099979A2 (en) * 2012-12-17 2014-06-26 Virginia Tech Intellectual Properties, Inc. Methods and compositions for identifying global microsatellite instability and for characterizing informative microsatellite loci

Family Cites Families (12)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
US20020086289A1 (en) * 1999-06-15 2002-07-04 Don Straus Genomic profiling: a rapid method for testing a complex biological sample for the presence of many types of organisms
US20050042654A1 (en) * 2003-06-27 2005-02-24 Affymetrix, Inc. Genotyping methods
JPWO2005001091A1 (ja) * 2003-06-27 2006-08-10 オリンパス株式会社 核酸の変異及び多型の検出用プローブセット、それを固相化したdnaアレイ、並びにそれらを用いた変異及び多型の検出方法
US20050136417A1 (en) * 2003-12-19 2005-06-23 Affymetrix, Inc. Amplification of nucleic acids
WO2007056825A1 (en) * 2005-11-21 2007-05-24 Simons Haplomics Limited Method and probes for identifying a nucleotide sequence
JP5237126B2 (ja) * 2006-03-01 2013-07-17 キージーン ナムローゼ フェンノートシャップ ライゲーションアッセイを用いてハイスループットシークエンスに基づき遺伝子関連配列を検出する方法
US7901882B2 (en) * 2006-03-31 2011-03-08 Affymetrix, Inc. Analysis of methylation using nucleic acid arrays
CN101790731B (zh) * 2007-03-16 2013-11-06 纳特拉公司 用于清除遗传数据干扰并确定染色体拷贝数的***和方法
CN101712959A (zh) * 2008-10-08 2010-05-26 中国人民解放军军事医学科学院放射与辐射医学研究所 一种新的人细胞生长抑制基因thap11及其应用
US20130157873A1 (en) * 2010-05-19 2013-06-20 Translational Genomics Research Institute Methods of assessing a risk of developing necrotizing meningoencephalitis
WO2012034251A2 (zh) * 2010-09-14 2012-03-22 深圳华大基因科技有限公司 一种基因组结构性变异检测方法和***
CN102127819B (zh) * 2010-11-22 2014-08-27 深圳华大基因科技有限公司 Mhc区域核酸文库的构建方法及用途

Patent Citations (2)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN101213312A (zh) * 2005-06-30 2008-07-02 先正达参股股份有限公司 用于筛选基因特异杂交多态性(gshp)的方法及其在遗传作图和标记开发中的用途
WO2014099979A2 (en) * 2012-12-17 2014-06-26 Virginia Tech Intellectual Properties, Inc. Methods and compositions for identifying global microsatellite instability and for characterizing informative microsatellite loci

Non-Patent Citations (1)

* Cited by examiner, † Cited by third party
Title
一种增强MLPA检测SNP位点的特异性方法;胡佳莉等;《贵阳医学院学报》;20120630;第37卷(第3期);231-234 *

Also Published As

Publication number Publication date
WO2016000267A1 (zh) 2016-01-07
CN106715711A (zh) 2017-05-24

Similar Documents

Publication Publication Date Title
CN106715711B (zh) 确定探针序列的方法和基因组结构变异的检测方法
CN107708556B (zh) 诊断方法
US20240170094A1 (en) Noninvasive prenatal molecular karyotyping from maternal plasma
De Roeck et al. NanoSatellite: accurate characterization of expanded tandem repeat length and sequence through whole genome long-read sequencing on PromethION
CN107077537B (zh) 用短读测序数据检测重复扩增
Tsai et al. Discovery of rare mutations in populations: TILLING by sequencing
JP5972448B2 (ja) コピー数変異を検出する方法及びシステム
WO2016176091A1 (en) Error suppression in sequenced dna fragments using redundant reads with unique molecular indices (umis)
US20200286586A1 (en) Sequence-graph based tool for determining variation in short tandem repeat regions
CN110770840A (zh) 用于对来自已知或未知基因型的多个贡献者的dna混合物分解和定量的方法和***
WO2021037016A1 (en) Methods for detecting absence of heterozygosity by low-pass genome sequencing
KR20230117036A (ko) 게놈의 반복 영역들에서의 짧은 판독물들을 시각화하기 위한 방법들 및 시스템들
Ahsan et al. A survey of algorithms for the detection of genomic structural variants from long-read sequencing data
US20180142300A1 (en) Universal haplotype-based noninvasive prenatal testing for single gene diseases
De Roeck et al. Accurate characterization of expanded tandem repeat length and sequence through whole genome long-read sequencing on PromethION
Deleye et al. Massively parallel sequencing of micro-manipulated cells targeting a comprehensive panel of disease-causing genes: A comparative evaluation of upstream whole-genome amplification methods
Huszar et al. Mitigating the effects of reference sequence bias in single-multiplex massively parallel sequencing of the mitochondrial DNA control region
WO2019132010A1 (ja) 塩基配列における塩基種を推定する方法、装置及びプログラム
CN110993024B (zh) 建立胎儿浓度校正模型的方法及装置与胎儿浓度定量的方法及装置

Legal Events

Date Code Title Description
PB01 Publication
PB01 Publication
SE01 Entry into force of request for substantive examination
GR01 Patent grant
GR01 Patent grant