CN108660200A - 一种检测短串联重复序列扩张的方法 - Google Patents

一种检测短串联重复序列扩张的方法 Download PDF

Info

Publication number
CN108660200A
CN108660200A CN201810499329.5A CN201810499329A CN108660200A CN 108660200 A CN108660200 A CN 108660200A CN 201810499329 A CN201810499329 A CN 201810499329A CN 108660200 A CN108660200 A CN 108660200A
Authority
CN
China
Prior art keywords
sequence
short tandem
ref
reference gene
gene group
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.)
Granted
Application number
CN201810499329.5A
Other languages
English (en)
Other versions
CN108660200B (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.)
Grandomics Biosciences Co ltd
Xiangya Hospital of Central South University
Original Assignee
Grandomics Biosciences Co ltd
Xiangya Hospital of Central South University
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 Grandomics Biosciences Co ltd, Xiangya Hospital of Central South University filed Critical Grandomics Biosciences Co ltd
Priority to CN201810499329.5A priority Critical patent/CN108660200B/zh
Publication of CN108660200A publication Critical patent/CN108660200A/zh
Application granted granted Critical
Publication of CN108660200B publication Critical patent/CN108660200B/zh
Active legal-status Critical Current
Anticipated expiration legal-status Critical

Links

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 (13)

1.一种检测短串联重复序列扩张的方法，其包括如下步骤：
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.根据权利要求1所述的检测短串联重复序列扩张的方法，其特征在于，所述步骤2)中序列比对软件为NGMLR。
3.根据权利要求1所述的检测短串联重复序列扩张的方法，其特征在于，所述三代测序数据为NA12878三代测序数据。
4.根据权利要求1所述的检测短串联重复序列扩张的方法，其特征在于，所述三代测序数据为用Pacbio Sequel平台测序得到的数据。
5.根据权利要求1所述的检测短串联重复序列扩张的方法，其特征在于，所述参考基因组为人hg19参考基因组。
6.根据权利要求1所述的检测短串联重复序列扩张的方法，其特征在于，在步骤3)中，当测序深度小于100X时，判断短串联重复区域是否存在重复单元扩张的方法为：
比较短串联重复区域的每一条reads上的重复单元数目ri与所述参考基因组上重复单元数目R，如果它们之间的碱基数目差di大于或等于阈值α，那么记为存在重复单元扩张的reads；如果重复单元扩张的reads的数目N与短串联重复区域平均深度的比值大于阈值β，则认为所述短串联重复区域存在重复单元扩张。
7.根据权利要求6所述的检测短串联重复序列扩张的方法，其特征在于，所述阈值α的值为10。
8.根据权利要求1所述的检测短串联重复序列扩张的方法，其特征在于，所述阈值β的值为0.1。
9.根据权利要求1所述的检测短串联重复序列扩张的方法，其特征在于，步骤4)中所述插入序列si的参考基因组位置包括染色体编号、开始位置以及结束位置。
10.根据权利要求1所述的检测短串联重复序列扩张的方法，其特征在于，所述步骤4)中阈值γ的值为10。
11.根据权利要求1至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。
12.根据权利要求11所述的检测短串联重复序列扩张的方法，其特征在于，计算插入序列在参考基因组上的位置和插入序列的长度的具体方法为：
如果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为所述片段在参考基因组的结束位置；
如果插入序列的长度大于阈值δ，则将其记录。
13.根据权利要求12所述的检测短串联重复序列扩张的方法，其特征在于，所述阈值δ的值为10。
CN201810499329.5A 2018-05-23 2018-05-23 一种检测短串联重复序列扩张的方法 Active CN108660200B (zh)

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 true CN108660200A (zh) 2018-10-16
CN108660200B 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)

Cited By (5)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN110066862A (zh) * 2019-05-22 2019-07-30 中南大学 一种基于高通量测序读数的重复dna序列识别方法
CN111378738A (zh) * 2018-12-29 2020-07-07 北京希望组生物科技有限公司 Htt和jph3基因的检测引物、试剂盒及方法
WO2020147657A1 (zh) * 2019-01-15 2020-07-23 广州柿宝生物科技有限公司 一种用于长链分子的数学序列重建方法
CN111863133A (zh) * 2019-12-30 2020-10-30 上海交通大学医学院附属瑞金医院 一种高通量测序数据的分析方法、试剂盒及分析系统
CN115273984A (zh) * 2022-09-30 2022-11-01 北京诺禾致源科技股份有限公司 鉴定基因组串联重复区域的方法及装置

Citations (3)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
US20140235845A1 (en) * 2007-06-20 2014-08-21 Life Technologies Corporation Sequence amplification with loopable primers
CN106295250A (zh) * 2016-07-28 2017-01-04 北京百迈客医学检验所有限公司 二代测序短序列快速比对分析方法及装置
CN106951731A (zh) * 2017-03-28 2017-07-14 上海至本生物科技有限公司 一种大片段插入或缺失的预测方法及系统

Patent Citations (3)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
US20140235845A1 (en) * 2007-06-20 2014-08-21 Life Technologies Corporation Sequence amplification with loopable primers
CN106295250A (zh) * 2016-07-28 2017-01-04 北京百迈客医学检验所有限公司 二代测序短序列快速比对分析方法及装置
CN106951731A (zh) * 2017-03-28 2017-07-14 上海至本生物科技有限公司 一种大片段插入或缺失的预测方法及系统

Non-Patent Citations (1)

* Cited by examiner, † Cited by third party
Title
QIAN LIU ET AL.: "Interrogating the "unsequenceable" genomic trinucleotide repeat disorders by long-read sequencing", 《GENOME MED》 *

Cited By (9)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN111378738A (zh) * 2018-12-29 2020-07-07 北京希望组生物科技有限公司 Htt和jph3基因的检测引物、试剂盒及方法
CN111378738B (zh) * 2018-12-29 2022-07-01 北京希望组生物科技有限公司 Htt和jph3基因的检测引物、试剂盒及方法
WO2020147657A1 (zh) * 2019-01-15 2020-07-23 广州柿宝生物科技有限公司 一种用于长链分子的数学序列重建方法
CN110066862A (zh) * 2019-05-22 2019-07-30 中南大学 一种基于高通量测序读数的重复dna序列识别方法
CN110066862B (zh) * 2019-05-22 2021-02-12 中南大学 一种基于高通量测序读数的重复dna序列识别方法
CN111863133A (zh) * 2019-12-30 2020-10-30 上海交通大学医学院附属瑞金医院 一种高通量测序数据的分析方法、试剂盒及分析系统
CN111863133B (zh) * 2019-12-30 2023-07-18 上海交通大学医学院附属瑞金医院 一种高通量测序数据的分析方法、试剂盒及分析系统
CN115273984A (zh) * 2022-09-30 2022-11-01 北京诺禾致源科技股份有限公司 鉴定基因组串联重复区域的方法及装置
CN115273984B (zh) * 2022-09-30 2022-11-29 北京诺禾致源科技股份有限公司 鉴定基因组串联重复区域的方法及装置

Also Published As

Publication number Publication date
CN108660200B (zh) 2022-10-18

Similar Documents

Publication Publication Date Title
CN108660200A (zh) 一种检测短串联重复序列扩张的方法
CN110010193B (zh) 一种基于混合策略的复杂结构变异检测方法
US20180225416A1 (en) Systems and methods for visualizing a pattern in a dataset
CN103403182B (zh) 与癌症相关的遗传或分子畸变的检测
CN104298892B (zh) 基因融合的检测装置和方法
CN113160882B (zh) 一种基于三代测序的病原微生物宏基因组检测方法
CN104794371B (zh) 检测逆转座子插入多态性的方法和装置
Larsson et al. Comparative microarray analysis
CN108830044A (zh) 用于检测癌症样本基因融合的检测方法和装置
CN110299185A (zh) 一种基于新一代测序数据的插入变异检测方法及系统
WO2019222757A1 (en) Inferring selection in white blood cell matched cell-free dna variants and/or in rna variants
CN113674803A (zh) 一种拷贝数变异的检测方法及其应用
Williams et al. Plant microRNA prediction by supervised machine learning using C5. 0 decision trees
Couvillion et al. Transcription elongation is finely tuned by dozens of regulatory factors
CN110219054A (zh) 一种核酸测序文库及其构建方法
CN108460248B (zh) 一种基于Bionano平台检测长串联重复序列的方法
CN104293892A (zh) 检测核基因组中与表型形状相关基因的方法
CN116434843A (zh) 一种碱基测序质量评估方法
US20070203653A1 (en) Method and system for computational detection of common aberrations from multi-sample comparative genomic hybridization data sets
Bell-Glenn et al. A novel framework for the identification of reference dna methylation libraries for reference-based deconvolution of cellular mixtures
CN116348957A (zh) 检测测序数据中的交叉污染
CN108304693B (zh) 利用高通量测序数据分析基因融合的方法
Dou et al. Differential expression analysis in RNA-Seq by a naive Bayes classifier with local normalization
Palejev Comparison of RNA-seq differential expression methods
Fay et al. Detecting hitchhiking from patterns of DNA polymorphism

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