CN108660200B - 一种检测短串联重复序列扩张的方法 - Google Patents
一种检测短串联重复序列扩张的方法 Download PDFInfo
- Publication number
- CN108660200B CN108660200B CN201810499329.5A CN201810499329A CN108660200B CN 108660200 B CN108660200 B CN 108660200B CN 201810499329 A CN201810499329 A CN 201810499329A CN 108660200 B CN108660200 B CN 108660200B
- Authority
- CN
- China
- Prior art keywords
- sequence
- ref
- short tandem
- tandem repeat
- reads
- 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
Links
- 108091092878 Microsatellite Proteins 0.000 title claims abstract description 129
- 238000000034 method Methods 0.000 title claims abstract description 24
- 238000003780 insertion Methods 0.000 claims abstract description 89
- 230000037431 insertion Effects 0.000 claims abstract description 89
- 238000001514 detection method Methods 0.000 claims abstract description 41
- 238000007671 third-generation sequencing Methods 0.000 claims abstract description 36
- 238000002864 sequence alignment Methods 0.000 claims abstract description 4
- 210000000349 chromosome Anatomy 0.000 claims description 35
- 239000012634 fragment Substances 0.000 claims description 26
- 238000012163 sequencing technique Methods 0.000 claims description 23
- 230000010261 cell growth Effects 0.000 claims description 3
- 230000035945 sensitivity Effects 0.000 description 14
- 206010039101 Rhinorrhoea Diseases 0.000 description 7
- 238000005516 engineering process Methods 0.000 description 6
- 238000012216 screening Methods 0.000 description 5
- 238000012360 testing method Methods 0.000 description 4
- 238000011144 upstream manufacturing Methods 0.000 description 4
- 230000000875 corresponding effect Effects 0.000 description 3
- 238000012937 correction Methods 0.000 description 2
- 230000003247 decreasing effect Effects 0.000 description 2
- 238000011156 evaluation Methods 0.000 description 2
- 238000007672 fourth generation sequencing Methods 0.000 description 2
- 230000035772 mutation Effects 0.000 description 2
- 108090000623 proteins and genes Proteins 0.000 description 2
- 230000005945 translocation Effects 0.000 description 2
- 108020004414 DNA Proteins 0.000 description 1
- 108091028043 Nucleic acid sequence Proteins 0.000 description 1
- 238000013459 approach Methods 0.000 description 1
- 238000003556 assay Methods 0.000 description 1
- 230000037429 base substitution Effects 0.000 description 1
- 230000009286 beneficial effect Effects 0.000 description 1
- 230000008859 change Effects 0.000 description 1
- 238000012217 deletion Methods 0.000 description 1
- 230000037430 deletion Effects 0.000 description 1
- 238000010586 diagram Methods 0.000 description 1
- 229940079593 drug Drugs 0.000 description 1
- 239000003814 drug Substances 0.000 description 1
- 230000000694 effects Effects 0.000 description 1
- 238000001914 filtration Methods 0.000 description 1
- 230000014509 gene expression Effects 0.000 description 1
- 239000000203 mixture Substances 0.000 description 1
- 238000012986 modification Methods 0.000 description 1
- 230000004048 modification Effects 0.000 description 1
- 238000007481 next generation sequencing Methods 0.000 description 1
- 239000002773 nucleotide Substances 0.000 description 1
- 125000003729 nucleotide group Chemical group 0.000 description 1
- 230000035790 physiological processes and functions Effects 0.000 description 1
- 230000008569 process Effects 0.000 description 1
- 230000003252 repetitive effect Effects 0.000 description 1
- 230000002441 reversible effect Effects 0.000 description 1
- 238000010200 validation analysis Methods 0.000 description 1
Images
Classifications
-
- C—CHEMISTRY; METALLURGY
- C12—BIOCHEMISTRY; BEER; SPIRITS; WINE; VINEGAR; MICROBIOLOGY; ENZYMOLOGY; MUTATION OR GENETIC ENGINEERING
- C12Q—MEASURING 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/00—Measuring or testing processes involving enzymes, nucleic acids or microorganisms; Compositions therefor; Processes of preparing such compositions
- C12Q1/68—Measuring or testing processes involving enzymes, nucleic acids or microorganisms; Compositions therefor; Processes of preparing such compositions involving nucleic acids
- C12Q1/6869—Methods for sequencing
Landscapes
- Life Sciences & Earth Sciences (AREA)
- Chemical & Material Sciences (AREA)
- Proteomics, Peptides & Aminoacids (AREA)
- Organic Chemistry (AREA)
- Zoology (AREA)
- Wood Science & Technology (AREA)
- Health & Medical Sciences (AREA)
- Engineering & Computer Science (AREA)
- Microbiology (AREA)
- Immunology (AREA)
- Biotechnology (AREA)
- Molecular Biology (AREA)
- Biophysics (AREA)
- Analytical Chemistry (AREA)
- Physics & Mathematics (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
本发明提供了一种检测短串联重复序列扩张的方法，其包括如下步骤：1)序列比对；2)RepeatHMM检测三代测序数据短串联重复；3)inScan检测短串联重复区域的序列插入；4)计算RepeatHMM检测结果与短串联重复区域的序列插入检测结果的交集。本发明结合序列插入检测和RepeatHMM短串联重复检测结果，提高了检测短串联重复序列扩张的特异性。
Description
技术领域
本发明属于基因测序技术领域，具体涉及短串联重复序列(STR，short tandemrepeat)扩张检测方法。
背景技术
短串联重复指DNA序列中的多个核苷酸(重复单元数目大于或等于2小于或等于6)前后首尾相连而构成的重复序列，重复单元数目的变化会对基因组结构造成重要影响，进而可能会影响基因的表达、修饰和相应的生理功能；短串联重复单元数目增多时称为短串联重复扩张。
三代测序指的是单条DNA/RNA分子测序技术，目前商用的三代测序技术有Pacbio公司的单分子实时测序技术和Nanopore公司的纳米孔测序技术。Pacbio公司的单分子实时测序技术测得的reads平均长度为10Kb，部分可以达到100kbp，Nanopore公司的纳米孔测序技术测得的reads平均长度也为10Kb，部分可以达到2.3Mbp。三代测序与二代测序相比的优势是读长更长、无GC偏好性，缺点是序列的错误率较高(约15％的错误率)。三代测序产生的长reads可以跨越短串联重复，从而准确检测重复单元的数目，同时可以检测二代测序无法检测的大尺度的短串联重复扩张(扩张长度大于二代测序的读长(100-300bp))。
现有用于三代测序检测短串联重复的方法是RepeatHMM(Liu,Q.,Zhang,P.,Wang,D.,Gu,W.&Wang,K.Interrogating the“unsequenceable”genomic trinucleotide repeatdisorders by long-read sequencing.Genome Medicine,65(2017))，该方法主要应用隐马尔可夫模型进行短串联重复单元的识别。但是RepeatHMM在检测短串联重复扩张时存在较高的假阳性。该技术方案的主要步骤如下：
1)选择感兴趣的短串联重复：从参考基因组序列中选择感兴趣的短串联重复，记录该短串联重复在参考基因组上的位置(染色体编号、起始位置、终止位置)、重复单元(如：CGG)和重复单元数目；
2)将三代测序的长reads比对到参考基因组：首先使用TRF软件(tandem repeatfinder)检测长reads上是否存在步骤1)中预设的串联重复，如果存在就将长reads切割为多段侧翼区序列和重复区序列；然后使用bwa mem和特定参数将侧翼区的片段比对到参考基因组，如果这些有序的侧翼区序列都成功比对，那么利用这些比对信息确定长reads上重复区域的起始位置和终止位置。这种方式被定义为“切分-重比对”策略；对于“切分-重比对”策略不能识别的长reads，使用bwa mem直接将其比对到参考基因组，如果长reads不能比对到参考基因组则丢弃；
3)确定reads上的重复区域：利用包含重复区的长reads的上游和下游序列(默认为18bp，用户可自定义)信息；分别使用bwa mem将上游和下游序列比对到参考基因组，如果它们都有高的比对一致性，那么在重复区和上下游序列之间加上一定数目的N，用来确保下面的过程不会将上下游序列识别为重复区序列；
4)长reads错误校正：首先构造一段比参考基因组上短串联重复长50％的完美短串联重复(重复单元完全一致的短串联重复)序列，例如，参考基因组上短串联重复的重复单元为CTG，长度为30个单元，新构造的完美短串联重复的重复单元也为CTG，长度为45个单元。然后使用一种非对称的比对算法UnsymSeqAlg，将长reads比对到新构造的完美短串联重复序列，根据比对结果进行序列校正；
5)检测重复单元数目：区域内的每一条长reads都通过隐马尔可夫模型识别重复单元数目；
6)重复单元数目峰值识别：应用混合高斯模型识别所有reads上重复单元数目的峰值，并计算样本短串联重复的基因型(genotype)；
7)检测短串联重复扩张：如果样本短串联重复单元数目大于参考基因组短串联重复单元数目，那么样本中存在短串联重复扩张。
然而，现有技术有如下缺陷：
1)对于短间隔相邻短串联重复，RepeatHMM容易产生假阳性检测结果。如：重复单元为AT，重复数目为n的短串联重复，(AT)n，和重复单元为AT，重复数目为m的短串联重复，(AT)m之间有大小为5bp的间隔碱基序列GCAGT，由于三代测序错误率较高，间隔碱基序列可能发生碱基替换、插入、缺失等突变，RepeatHMM可能错误地将(AT)n和(AT)m识别为一个连续的短串联重复，产生假阳性短串联重复扩张检测结果；
2)短串联重复扩张是一种特殊类型的序列插入，现有用于三代测序检测基因组结构变异的工具并不是为检测短串联重复区域的序列插入而设定的，检测这类序列插入的敏感度低；
3)RepeatHMM与现有用于三代测序检测基因组结构变异的工具结合使用可以提高短串联重复检测的特异性，但是会造成较大的敏感性损失。
发明内容
本发明的目的在于提高利用三代测序检测基因组序列插入的敏感性，并降低RepeatHMM检测短串联重复扩张的假阳性，为此本发明提供了一种检测短串联重复序列扩张的方法，其包括如下步骤：
1)获得三代测序数据；
2)序列比对
使用序列比对软件将三代测序数据比对到参考基因组；
3)RepeatHMM检测三代测序数据短串联重复
使用RepeatHMM检测短串联重复单元数目，判断短串联重复区域是否存在重复单元扩张；
4)inScan检测短串联重复区域的序列插入
对三代测序数据比对结果，提取目标区域内的reads；
计算reads片段内的插入序列si的参考基因组位置和长度，如果si的长度大于或等于阈值γ，那么记录si；
检测reads片段间插入序列，设一条reads在比对时切分为n条片段Fr1至Frn，所述片段按照其在reads上的开始位置read_start，从小到大进行排序得到片段组成的数组Fr，数组的长度为n，组合其中两个reads片段，计算所述两个reads片段的相对位置，判断两个reads片段之间是否存在插入序列，计算插入序列在参考基因组上的位置和插入序列的长度；
5)计算RepeatHMM检测结果与短串联重复区域的序列插入检测结果的交集
对于一个短串联重复区域，如果RepeatHMM检测到该短串联重复区域存在重复单元扩张，同时检测到该短串联重复区域存在序列插入，则所述短串联重复区域称为RepeatHMM检测结果与短串联重复区域的序列插入检测结果的交集。
根据本发明的实施方式，所述步骤2)中序列比对软件为NGMLR。
根据本发明的实施方式，所述三代测序数据为NA12878三代测序数据。
根据本发明的实施方式，所述三代测序数据为用Pacbio Sequel平台测序得到的数据。
根据本发明的实施方式，所述参考基因组为人hg19参考基因组。
根据本发明的实施方式，在步骤3)中，当测序深度小于100X时，判断短串联重复区域是否存在重复单元扩张的方法为：
比较短串联重复区域的每一条reads上的重复单元数目ri与所述参考基因组上重复单元数目R，如果它们之间的碱基数目差di大于或等于阈值α，那么记为存在重复单元扩张的reads；如果重复单元扩张的reads的数目N与短串联重复区域平均深度的比值大于阈值β，则认为所述短串联重复区域存在重复单元扩张。
根据本发明的实施方式，所述阈值α的值为10。
根据本发明的实施方式，所述阈值β的值为0.1。
根据本发明的实施方式，步骤4)中所述插入序列si的参考基因组位置包括染色体编号、开始位置以及结束位置。
根据本发明的实施方式，所述步骤4)中阈值γ的值为10。
根据本发明的实施方式，所述步骤4)中判断片段之间是否存在插入序列的具体方法为：
对于片段Fr[i]与Fr[j]，其中i>＝1且i<＝n-1，j>i且j<＝n，如果Fr[i]与Fr[j]比对到同一条染色体、比对方向相同且它们在reads上的距离drij大于它们在参考基因组上的距离dfij，那么Fr[i]与Fr[j]之间存在序列插入；如果Fr[i]与Fr[j]之间存在序列插入，那么i＝i+1；如果Fr[i]与Fr[j]之间不存在序列插入且Fr[i]与Fr[j]在同一条染色体上，那么i＝i+1；
计算插入序列在参考基因组上的位置和插入序列的长度，如果Fr[i]与Fr[j]之间存在插入序列，那么分3种情况计算插入序列在参考基因组上的位置和插入序列的长度：
a.INS/INDEL类型的序列插入，如果Fr[j].ref_start>＝Fr[i].ref_start，则插入序列的长度insert_lenght＝Fr[j].read_start–Fr[i].read_end，插入序列在参考基因组染色体上的起始位置insert_ref_start＝Fr[i].ref_end，插入序列在参考基因组染色体上的终止位置insert_ref_end＝Fr[j].ref_start；
b.TANDEM_DUP类型的序列插入，如果Fr[j].ref_start<Fr[i].ref_start且Fr[j].ref_end>＝Fr[i].ref_end，则插入序列的长度insert_lenght＝(Fr[j].read_start–Fr[i].read_end)–(Fr[j].ref_start–Fr[i].ref_end)，插入序列在参考基因组染色体上的起始位置insert_ref_start＝Fr[i].ref_end，插入序列在参考基因组染色体上的终止位置insert_ref_end＝Fr[i].ref_end；
c.DUP类型的序列插入，如果Fr[j].ref_start<Fr[i].ref_start且Fr[j].ref_end<Fr[i].ref_end，则插入序列的长度insert_lenght＝Fr[j].read_end–Fr[i].read_end，插入序列在参考基因组染色体上的起始位置insert_ref_start＝Fr[i].ref_end，插入序列在参考基因组染色体上的终止位置insert_ref_end＝Fr[i].ref_end；
其中，read_start为所述片段在reads上的起始位置，read_end为所述片段在reads上的结束位置，ref为所述片段比对到的参考基因组染色体，ref_start为所述片段在参考基因组的开始位置，ref_end为所述片段在参考基因组的结束位置；
如果插入序列的长度大于阈值δ，则将其记录。
根据本发明的实施方式，所述阈值δ的值为10。
本发明技术方案带来的有益效果
本技术方案可以有效地降低利用三代测序数据检测短串联重复序列扩张的假阳性，同时不会对检测的敏感性造成大的损失：
1.我们对RepeatHMM检测低测序深度(小于100X)样本的短串联重复序列扩张方法进行了改进，直接比较每一条reads的重复单元数目与参考基因组重复单元数目的差异并利用测序深度信息进行过滤(设置阈值)来检测短串联重复序列扩张；
2.结合序列插入检测和RepeatHMM短串联重复检测结果，可以提高检测短串联重复序列扩张的特异性；
3.开发专用于三代测序检测基因组序列插入的方法，提高检测序列插入的敏感性，使得本技术方案不会对短串联重复扩张检测的敏感性造成大的损失。
附图说明
图1为本发明主要方法的流程图；
图2为本发明检测复杂序列插入的示意图；
图3为本发明可以检测的4种类型的序列插入示意图；
图4为短串联重复区域序列扩张长度与序列插入长度的相关性图。
具体实施方式
下面将结合本发明实施例中的附图，对本发明实施例中的技术方案进行清楚、完整地描述，显然，所描述的实施例仅仅是本发明一部分实施例，而不是全部的实施例。基于本发明中的实施例，本领域普通技术人员在没有做出创造性劳动前提下所获得的所有其他实施例，都属于本发明保护的范围。
除非另有定义，本文所使用的所有的技术和科学术语与属于本发明的技术领域的技术人员通常理解的含义相同。本文中在本发明的说明书中所使用的术语只是为了描述具体的实施例的目的，不是旨在于限制本发明。
实施例
二代测序数据读长短(100-300bp)，但是测序准确度高(>99％)，适用于检测短重复扩张长度小于其读长的重复扩张。本实施例首先使用HipSTR(Willems,T.et al.Genome-wide profiling of heritable and de novo STR variations.Nature Methods14,590(2017))检测NA12878二代测序数据(illuminaHiSeq2500，150bp paired-end，30X)中的短串联重复数目，随机抽取1000个存在短串联重复扩张的区域和4000个不存在短串联重复扩张的区域作为测试数据集，然后利用该测试数据集测试不同检测方法检测NA12878三代测序数据(Pacbio Sequel测序平台，约44X)短串联重复扩张的效果。本技术方案的主要方法流程图如图1所示：
1、构建验证数据集
随机抽取人hg19参考基因组上的50000个短串联重复区域，对NA12878的二代测序数据使用HipSTR检测短串联重复单元数目和基因型；
二代测序短reads完全跨过变异区域的检测结果可信度高，按照该条件筛选得到高可信度的短串联重复区域集合；
如果样本短串联重复序列的长度与参考基因组上相应的短串联重复序列长度的差值(d)大于或等于10bp，那么定义该短串联重复序列存在重复扩张。从上述高可信度的短串联重复区域集合中，按照d大于或等于10进行筛选，然后随机抽取1000个存在重复扩张的短串联重复区域；
如果样本短串联重复序列的长度与参考基因组上相应的短串联重复序列长度的差值(d)等于0，那么该短串联重复序列与参考基因组上短串联重复序列相同。从高可信度的短串联重复区域集合中，按照d等于0进行筛选，然后随机抽取4000个与参考基因组上短串联重复序列相同的短串联重复区域；
合并随机抽取的1000个存在重复扩张的短串联重复区域和随机抽取的4000个与参考基因组上短串联重复序列相同的短串联重复区域，得到包含5000个短串联重复区域的测试数据集。
2、序列比对
使用NGMLR(Sedlazeck,F.J.et al.Accurate detection of complexstructural variations using single molecule sequencing.bioRxiv(2017).)序列比对软件(-X Pacbio)将约44X的NA12878三代测序数据(Pacbio Sequel平台)比对到人参考基因组hg19。
3、RepeatHMM检测三代测序数据短串联重复
1)使用RepeatHMM检测测试数据集的5000个短串联重复单元数目；
2)RepeatHMM检测低深度的测序样本(小于100X)时存在无法判断基因型的情况，因此对于低深度样本，比较短串联重复区域的每一条reads上的重复单元数目(ri)与参考基因组上重复单元数目(R)，如果它们之间的碱基数目差di大于或等于阈值α(默认10)，那么记为存在重复单元扩张的reads；如果重复单元扩张的reads的数目(N)与短串联重复区域平均深度的比值大于阈值β(默认0.1)，那么认为该短串联重复区域存在重复单元扩张。
4、inScan检测短串联重复区域的序列插入
1)对NA12878三代测序数据(Pacbio Sequel)NGMLR比对结果，提取目标区域内的reads；
2)计算reads片段内的插入序列si的参考基因组位置(染色体编号、开始位置、结束位置)和长度，如果si的长度大于或等于阈值γ(默认为10)，那么记录si；
3)检测reads片段间插入序列：reads的每一个比对片段包含5个主要属性：在reads上的开始位置read_start,在read上的结束位置read_end，比对到的参考基因组染色体ref，在参考基因组的开始位置ref_start，在参考基因组的结束位置ref_end。设一条reads在比对时切分为n条片段(Fr1,Fr2,Fr3…Frn)，这些片段按照它们在reads上的开始位置read_start，从小到大进行排序可以得到片段组成的数组Fr(数组的长度为n)。组合两个reads片段，计算它们的相对位置，可以推测片段之间是否存在插入序列、插入序列在参考基因组上的位置和插入序列的长度。判断插入序列是否存在：对于片段Fr[i](i>＝1且i<＝n-1)与Fr[j](j>i且j<＝n)，如果Fr[i]与Fr[j]比对到同一条染色体、比对方向相同且它们在reads上的距离drij大于它们在参考基因组上的距离dfij，那么Fr[i]与Fr[j]之间存在序列插入。如果Fr[i]与Fr[j]之间存在序列插入，那么i＝i+1。如果Fr[i]与Fr[j]之间不存在序列插入且Fr[i]与Fr[j]在同一条染色体上，那么i＝i+1。例如：图2所示，一条三代测序reads分成5个片段，以一种复杂的方式比对到参考基因组。其中，Fr[1]比对到chr1染色体，Fr[2]反向互补后比对到chr1染色体，Fr[3]比对到chr1染色体，Fr[4]比对到chr20染色体，Fr[5]比对到chr1染色体。inScan在第四步时发现Fr[3]和Fr[5]之间存在插入片段Fr[4]。而sniffles和nanosv等结构变异检测软件将检测出chr1于chr20之间存在两个易位(Translocation)断点(Breakpoint)。计算在参考基因组上的位置和插入序列的长度：如果Fr[i]与Fr[j]之间存在插入序列，那么分3种情况计算在参考基因组上的位置和插入序列的长度：a.INS/INDEL类型的序列插入，如果Fr[j].ref_start>＝Fr[i].ref_start，那么插入序列的长度insert_lenght＝Fr[j].read_start–Fr[i].read_end，插入序列在参考基因组染色体上的起始位置insert_ref_start＝Fr[i].ref_end，插入序列在参考基因组染色体上的终止位置insert_ref_end＝Fr[j].ref_start；b.TANDEM_DUP类型的序列插入，如果Fr[j].ref_start<Fr[i].ref_start且Fr[j].ref_end>＝Fr[i].ref_end，那么插入序列的长度insert_lenght＝(Fr[j].read_start–Fr[i].read_end)–(Fr[j].ref_start–Fr[i].ref_end)，插入序列在参考基因组染色体上的起始位置insert_ref_start＝Fr[i].ref_end，插入序列在参考基因组染色体上的终止位置insert_ref_end＝Fr[i].ref_end；c.DUP类型的序列插入，如果Fr[j].ref_start<Fr[i].ref_start且Fr[j].ref_end<Fr[i].ref_end，那么插入序列的长度insert_lenght＝Fr[j].read_end–Fr[i].read_end(插入序列长度的下限)，插入序列在参考基因组染色体上的起始位置insert_ref_start＝Fr[i].ref_end，插入序列在参考基因组染色体上的终止位置insert_ref_end＝Fr[i].ref_end。如果插入序列的长度大于阈值δ(默认为10)，那么将其记录。4种序列插入的示意图如图3所示。
4)对检测结果按照插入序列长度小于150bp、支持reads数目大于或等于2进行筛选。
作为对照，用Sniffles检测短串联重复区域的序列插入，其步骤为：
1)对NA12878三代测序数据(Pacbio Sequel)NGMLR比对结果，使用Sniffles(检测的最小变异长度设置为10bp，最小支持reads数目设置为1)软件检测全基因组结构变异；
2)按照变异类型为插入(insertion)或重复(duplication)、插入序列长度小于150bp(HipSTR利用二代测序数据检测短串联重复时无法检出大于reads长度的重复单元扩张)、支持的reads数目大于或等于2进行筛选，得到Sniffles检测到的序列插入集合；
3)使用bedtools intersect计算短串联重复区域与Sniffles检测到的序列插入的overlap，计算短串联重复区域内序列插入长度的中值。
5、计算RepeatHMM检测结果与短串联重复区域的序列插入检测结果的交集
对于一个短串联重复区域，如果RepeatHMM检测到该短串联重复区域存在重复单元扩张，同时检测到该短串联重复区域存在序列插入，那么这种短串联重复区域称为RepeatHMM检测结果与短串联重复区域的序列插入检测结果的交集。
6、结果评估
对于利用二代测序数据检出存在重复单元扩张的1000个短串联重复区域，使用inScan检测三代测序数据(Pacbio Sequel平台)发现短串联重复区域内序列插入长度的中值与短串联重复扩张长度具有很高的相关性(P值小于2.2X10-16，R2等于0.695，图4)。在三代测序数据中，存在重复单元扩张的1000个短串联重复区域中的993个包含插入序列。单独使用RepeatHMM检测随机选择的5000个短串联重复区域(1000个存在重复扩张，4000个重复单元数目不变)，共发现1136个短串联重复区域存在重复扩张，其中934个真阳性结果，202个假阳性结果，敏感度为93.4％，阳性率为82.2％，假阳性率为17.8％；结合使用RepeatHMM和Sniffles检测随机选择的5000个短串联重复区域，发现498个短串联重复区域存在重复扩张，其中458个真阳性结果，40个假阳性结果，敏感度为45.8％，阳性率为92％，假阳性率为8％；结合使用RepeatHMM和inScan检测随机选择的5000个短串联重复区域，发现980个短串联重复区域存在重复扩张，其中920个真阳性结果，60个假阳性结果，敏感度为92％，阳性率为93.9％，假阳性率为6.1％；说明结合使用RepeatHMM和inScan可以有效提高利用三代测序数据检测短串联重复扩张的特异性(假阳性率下降11.7％，从17.8％下降到6.1％)，同时不会对检测的敏感度造成大的损失(敏感度仅下降1.4％，从93.4％下降到92％)。结合使用RepeatHMM和Sniffles也可以提高利用三代测序数据检测短串联重复扩张的特异性(假阳性率下降9.8％，从17.8％下降到8％)，但是有非常大的敏感度损失(敏感度下降47.6％，从93.4％下降到45.8％)(表1)。综上所述，结合使用RepeatHMM和inScan的短串联重复序列扩张检测结果表现最佳。
表1不同检测方法的结果评估
| 检测方法 | 敏感度(％) | 阳性率(％) | 假阳性率(％) |
| RepeatHMM | 93.4 | 82.2 | 17.8 |
| RepeatHMM+Sniffles | 45.8 | 92.0 | 8.0 |
| RepeatHMM+inScan | 92.0 | 93.9 | 6.1 |
Claims (8)
1.一种检测短串联重复序列扩张的方法，其包括如下步骤：
1)获得三代测序数据；
2)序列比对
使用序列比对软件将所述三代测序数据比对到参考基因组；
3)RepeatHMM检测所述三代测序数据短串联重复
使用RepeatHMM检测短串联重复单元数目，判断短串联重复区域是否存在重复单元扩张；
在步骤3)中，当测序深度小于100X时，判断短串联重复区域是否存在重复单元扩张的方法为：
比较短串联重复区域的每一条reads上的重复单元数目ri与所述参考基因组上重复单元数目R，如果它们之间的碱基数目差di大于或等于阈值α，那么记为存在重复单元扩张的reads；如果重复单元扩张的reads的数目N与短串联重复区域平均深度的比值大于阈值β，则认为所述短串联重复区域存在重复单元扩张；
4)inScan检测短串联重复区域的序列插入
对三代测序数据比对结果，提取目标区域内的reads；
计算reads片段内的插入序列si的参考基因组位置和长度，如果si的长度大于或等于阈值γ，那么记录si；
检测reads片段间插入序列，设一条reads在比对时切分为n条片段Fr1至Frn，所述片段按照其在reads上的开始位置read_start，从小到大进行排序得到片段组成的数组Fr，数组的长度为n，组合其中两个reads片段，计算所述两个reads片段的相对位置，判断所述两个reads片段之间是否存在插入序列，计算插入序列在参考基因组上的位置和插入序列的长度；
所述步骤4)中判断片段之间是否存在插入序列的具体方法为：
对于片段Fr[i]与Fr[j]，其中i>＝1且i<＝n-1，j>i且j<＝n，如果Fr[i]与Fr[j]比对到同一条染色体、比对方向相同且它们在reads上的距离drij大于它们在参考基因组上的距离dfij，那么Fr[i]与Fr[j]之间存在序列插入；如果Fr[i]与Fr[j]之间存在序列插入，那么i＝i+1；如果Fr[i]与Fr[j]之间不存在序列插入且Fr[i]与Fr[j]在同一条染色体上，那么i＝i+1；
计算插入序列在参考基因组上的位置和插入序列的长度的具体方法为：
如果Fr[i]与Fr[j]之间存在插入序列，那么分3种情况计算插入序列在参考基因组上的位置和插入序列的长度：
a.INS/INDEL类型的序列插入，如果Fr[j].ref_start>＝Fr[i].ref_start，则插入序列的长度insert_lenght＝Fr[j].read_start–Fr[i].read_end，插入序列在参考基因组染色体上的起始位置insert_ref_start＝Fr[i].ref_end，插入序列在参考基因组染色体上的终止位置insert_ref_end＝Fr[j].ref_start；
b.TANDEM_DUP类型的序列插入，如果Fr[j].ref_start<Fr[i].ref_start且Fr[j].ref_end>＝Fr[i].ref_end，则插入序列的长度insert_lenght＝(Fr[j].read_start–Fr[i].read_end)–(Fr[j].ref_start–Fr[i].ref_end)，插入序列在参考基因组染色体上的起始位置insert_ref_start＝Fr[i].ref_end，插入序列在参考基因组染色体上的终止位置insert_ref_end＝Fr[i].ref_end；
c.DUP类型的序列插入，如果Fr[j].ref_start<Fr[i].ref_start且Fr[j].ref_end<Fr[i].ref_end，则插入序列的长度insert_lenght＝Fr[j].read_end–Fr[i].read_end，插入序列在参考基因组染色体上的起始位置insert_ref_start＝Fr[i].ref_end，插入序列在参考基因组染色体上的终止位置insert_ref_end＝Fr[i].ref_end；
其中，read_start为所述片段在reads上的起始位置，read_end为所述片段在reads上的结束位置，ref为所述片段比对到的参考基因组染色体，ref_start为所述片段在参考基因组的开始位置，ref_end为所述片段在参考基因组的结束位置；
如果插入序列的长度大于阈值δ，则将其记录；
步骤4)中所述插入序列si的参考基因组位置包括染色体编号、开始位置以及结束位置；
所述步骤4)中阈值γ的值为10；
5)计算RepeatHMM检测结果与短串联重复区域的序列插入检测结果的交集
对于一个短串联重复区域，如果RepeatHMM检测到该短串联重复区域存在重复单元扩张，同时检测到该短串联重复区域存在序列插入，则所述短串联重复区域称为RepeatHMM检测结果与短串联重复区域的序列插入检测结果的交集。
2.根据权利要求1所述的检测短串联重复序列扩张的方法，其特征在于，所述步骤2)中序列比对软件为NGMLR。
3.根据权利要求1所述的检测短串联重复序列扩张的方法，其特征在于，所述三代测序数据为NA12878三代测序数据。
4.根据权利要求1所述的检测短串联重复序列扩张的方法，其特征在于，所述三代测序数据为用Pacbio Sequel平台测序得到的数据。
5.根据权利要求1所述的检测短串联重复序列扩张的方法，其特征在于，所述参考基因组为人hg19参考基因组。
6.根据权利要求1所述的检测短串联重复序列扩张的方法，其特征在于，所述阈值α的值为10。
7.根据权利要求1所述的检测短串联重复序列扩张的方法，其特征在于，所述阈值β的值为0.1。
8.根据权利要求1所述的检测短串联重复序列扩张的方法，其特征在于，所述阈值δ的值为10。
Priority Applications (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN201810499329.5A CN108660200B (zh) | 2018-05-23 | 2018-05-23 | 一种检测短串联重复序列扩张的方法 |
Applications Claiming Priority (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN201810499329.5A CN108660200B (zh) | 2018-05-23 | 2018-05-23 | 一种检测短串联重复序列扩张的方法 |
Publications (2)
| Publication Number | Publication Date |
|---|---|
| CN108660200A CN108660200A (zh) | 2018-10-16 |
| CN108660200B true CN108660200B (zh) | 2022-10-18 |
Family
ID=63777486
Family Applications (1)
| Application Number | Title | Priority Date | Filing Date |
|---|---|---|---|
| CN201810499329.5A Active CN108660200B (zh) | 2018-05-23 | 2018-05-23 | 一种检测短串联重复序列扩张的方法 |
Country Status (1)
| Country | Link |
|---|---|
| CN (1) | CN108660200B (zh) |
Families Citing this family (5)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN111378738B (zh) * | 2018-12-29 | 2022-07-01 | 北京希望组生物科技有限公司 | Htt和jph3基因的检测引物、试剂盒及方法 |
| CN109801679B (zh) * | 2019-01-15 | 2021-02-02 | 广州柿宝生物科技有限公司 | 一种用于长链分子的数学序列重建方法 |
| CN110066862B (zh) * | 2019-05-22 | 2021-02-12 | 中南大学 | 一种基于高通量测序读数的重复dna序列识别方法 |
| CN111863133B (zh) * | 2019-12-30 | 2023-07-18 | 上海交通大学医学院附属瑞金医院 | 一种高通量测序数据的分析方法、试剂盒及分析系统 |
| CN115273984B (zh) * | 2022-09-30 | 2022-11-29 | 北京诺禾致源科技股份有限公司 | 鉴定基因组串联重复区域的方法及装置 |
Citations (2)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN106295250A (zh) * | 2016-07-28 | 2017-01-04 | 北京百迈客医学检验所有限公司 | 二代测序短序列快速比对分析方法及装置 |
| CN106951731A (zh) * | 2017-03-28 | 2017-07-14 | 上海至本生物科技有限公司 | 一种大片段插入或缺失的预测方法及系统 |
Family Cites Families (1)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| US20090023190A1 (en) * | 2007-06-20 | 2009-01-22 | Kai Qin Lao | Sequence amplification with loopable primers |
-
2018
- 2018-05-23 CN CN201810499329.5A patent/CN108660200B/zh active Active
Patent Citations (2)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN106295250A (zh) * | 2016-07-28 | 2017-01-04 | 北京百迈客医学检验所有限公司 | 二代测序短序列快速比对分析方法及装置 |
| CN106951731A (zh) * | 2017-03-28 | 2017-07-14 | 上海至本生物科技有限公司 | 一种大片段插入或缺失的预测方法及系统 |
Non-Patent Citations (1)
| Title |
|---|
| Interrogating the "unsequenceable" genomic trinucleotide repeat disorders by long-read sequencing;Qian Liu et al.;《Genome Med》;20170718;第9卷(第1期);第1-16页 * |
Also Published As
| Publication number | Publication date |
|---|---|
| CN108660200A (zh) | 2018-10-16 |
Similar Documents
| Publication | Publication Date | Title |
|---|---|---|
| CN108660200B (zh) | 一种检测短串联重复序列扩张的方法 | |
| CN110299185B (zh) | 一种基于新一代测序数据的插入变异检测方法及系统 | |
| US11694764B2 (en) | Method for large scale scaffolding of genome assemblies | |
| CN106951731B (zh) | 一种大片段插入或缺失的预测方法及系统 | |
| CN108350495B (zh) | 对分隔长片段序列进行组装的方法和装置 | |
| CN110010193A (zh) | 一种基于混合策略的复杂结构变异检测方法 | |
| CN110600078A (zh) | 一种基于纳米孔测序检测基因组结构变异的方法 | |
| CN106715711A (zh) | 确定探针序列的方法和基因组结构变异的检测方法 | |
| CN107480470B (zh) | 基于贝叶斯与泊松分布检验的已知变异检出方法和装置 | |
| KR20140140122A (ko) | 복제 수 변이를 검측하기 위한 방법 및 시스템 | |
| CN110033829A (zh) | 基于差异snp标记物的同源基因的融合检测方法 | |
| CN108304694B (zh) | 基于二代测序数据分析基因突变的方法 | |
| CN106355045B (zh) | 一种基于扩增子二代测序小片段插入缺失检测的方法及装置 | |
| CN111326212A (zh) | 一种结构变异的检测方法 | |
| CN103902852A (zh) | 基因表达的定量方法及装置 | |
| CN111755072A (zh) | 一种同时检测甲基化水平、基因组变异和插入片段的方法及装置 | |
| CN112365922A (zh) | 用于检测msi的微卫星位点、其筛选方法及应用 | |
| CN105483210A (zh) | 一种rna编辑位点的检测方法 | |
| CN115896256A (zh) | 基于二代测序技术的rna插入缺失突变的检测方法、装置、设备和存储介质 | |
| CN109920480B (zh) | 一种校正高通量测序数据的方法和装置 | |
| KR20160039386A (ko) | Itd 검출 장치 및 방법 | |
| McShane et al. | Characterizing nascent transcription patterns of PROMPTs, eRNAs, and readthrough transcripts in the ENCODE4 deeply profiled cell lines | |
| WO2019132010A1 (ja) | 塩基配列における塩基種を推定する方法、装置及びプログラム | |
| CN113416770B (zh) | 一种染色体结构变异断点的定位方法及装置 | |
| CN115831222A (zh) | 一种基于三代测序的全基因组结构变异鉴定方法 |
Legal Events
| Date | Code | Title | Description |
|---|---|---|---|
| PB01 | Publication | ||
| PB01 | Publication | ||
| SE01 | Entry into force of request for substantive examination | ||
| SE01 | Entry into force of request for substantive examination | ||
| GR01 | Patent grant | ||
| GR01 | Patent grant |