CN110066862B - 一种基于高通量测序读数的重复dna序列识别方法 - Google Patents

一种基于高通量测序读数的重复dna序列识别方法 Download PDF

Info

Publication number
CN110066862B
CN110066862B CN201910428254.6A CN201910428254A CN110066862B CN 110066862 B CN110066862 B CN 110066862B CN 201910428254 A CN201910428254 A CN 201910428254A CN 110066862 B CN110066862 B CN 110066862B
Authority
CN
China
Prior art keywords
frequency
reading
mer
throughput sequencing
mers
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
CN201910428254.6A
Other languages
English (en)
Other versions
CN110066862A (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.)
Central South University
Original Assignee
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 Central South University filed Critical Central South University
Priority to CN201910428254.6A priority Critical patent/CN110066862B/zh
Publication of CN110066862A publication Critical patent/CN110066862A/zh
Application granted granted Critical
Publication of CN110066862B publication Critical patent/CN110066862B/zh
Active legal-status Critical Current
Anticipated expiration legal-status Critical

Links

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

本发明公开了一种基于高通量测序读数的重复DNA序列识别方法，包括：由高通量测序的读数得到高频k‑mer集合，根据高频k‑mer集合对读数进行筛选，使得包含高频k‑mer较多的读数保留下来，成为高频读数；使用序列组装工具组装高频读数，得到contigs序列；对contigs序列进行筛选，保留下的所有contigs序列即为重复DNA序列。本发明可以从高通量测序读数中识别重复DNA序列，而无需物种参考序列，可以适用于参考序列未知的物种的重复DNA序列识别，并且本发明是通过组装高频读数得到重复DNA序列，相对于组装高频k‑mer得到重复DNA序列，提高了识别重复DNA序列的准确率。

Description

一种基于高通量测序读数的重复DNA序列识别方法
技术领域
本发明涉及分子生物学领域，具体而言，涉及一种基于高通量测序读数的重复DNA序列识别方法。
背景技术
脱氧核糖核酸(DNA)是一种由双链结构组成的分子，两条链相互缠绕形成双螺旋结构。DNA携带着生物用于引导生长发育和生命机能运作的遗传指令并控制着所有生物的生长和繁殖。两条DNA链被称为多核苷酸，每条链由核苷酸这种更小的单元线性链接而成，其中每个核苷酸单元由胞嘧啶(C)，鸟嘌呤(G)，腺嘌呤(A)，胸腺嘧啶(T)四种含氮碱基之一、脱氧核糖以及磷酸基团组成。核苷酸通过一个核苷酸的核糖与下一个核苷酸的磷酸之间的共价键在链中彼此连接，产生交替的核糖——磷酸骨架。根据碱基配对规则(A与T和C与G)，氢键将两条分开的多核苷酸链的含氮碱基结合在一起，形成双链DNA。基于DNA的结构信息，现代分子生物学通常用由A、T、C、G组成的字符串来代表DNA分子的内容信息，并将这些字符串存储在存储介质中以便进行科学研究。DNA测序技术是指确定DNA中核苷酸顺序的一种生物技术。近年来，对低成本测序的高需求推动了高通量测序技术(High-throughputsequencing)，也被称为下一代测序技术(Next-generation sequencing)的发展，这种技术使得测序过程可以并行化进行，可以一次对几十万到几百万条DNA分子进行序列测定。高通量测序技术降低了DNA测序难度，使得DNA测序数据的规模成指数型增长。据统计，美国国立生物技术信息中心(NCBI)建立的公共DNA序列数据库GenBank在1983年时仅有7万条DNA序列，数据规模为680338个碱基对(base pair)；而截至到2018年底，GenBank有超过两亿条DNA序列，数据规模超过两千亿个碱基对。如何有效地分析这些数据并发现其中的规律以指导生物学研究和实验，是当今分子生物学和生物信息学的研究重点和难点。
重复DNA序列是指在DNA序列中以相同或相似的方式出现多次的序列片段。在许多物种的基因组DNA序列中，重复DNA序列占有很高的比例，例如人类基因组中超过60％、玉米基因组中超过80％和果蝇基因组中超过20％都是由重复DNA序列构成。重复DNA序列对目前生物信息学的许多基础的分析任务，例如序列从头组装、序列比对、序列纠错等的进行带来了一定的挑战，重复DNA序列也已经被证实对多重真核生物的进化规律和疾病的产生起到重要影响。
目前已经有了诸如基于同源性分析、基于结构相似性、基于高通量测序读数三类重复DNA序列的识别方法或工具，其中基于高通量测序读数的识别方法由于其不需要重复DNA序列结构和物种的参考序列等先验知识，具有更加灵活的特点。目前技术中，RepARK和REPdenovo是两种被普遍采用的基于高通量测序读数的重复DNA序列识别方法。2014年，Koch等人提出了从测序得到的序列中识别高频k-mer并将k-mer进行序列组装得到重复DNA序列的RepARK方法。2016年，Chu等人在RepARK的基础上，将组装得到的contigs序列进行合并得到重复DNA序列。上述两种方法都是将高频k-mer直接进行序列组装得到重复DNA序列，但是直接使用高频k-mer进行序列组装过程存在两点不足：第一，常见的序列组装工具在序列组装过程中需要达到一定测序深度的序列作为输入，才会得到可靠的组装结果，互不相同的高频k-mer显然不能满足测序深度的要求；第二，k-mer的长度较短，常见的序列组装工具并不适用于组装如此短的序列，这会导致序列组装工具在这种情况下不能得出最准确的组装结果，进而影响最终识别的重复DNA序列的质量。因此，采用序列组装方法或工具直接对k-mer进行组装，得到的组装结果不够可靠，识别重复DNA序列的准确率不高，还有进一步提升的空间。
发明内容
本发明所提供的技术问题是，针对现有技术的不足，提供一种基于高通量测序读数的重复DNA序列识别方法，准确率高。
本发明所提供的技术方案为：
本发明提供了一种基于高通量测序读数的重复DNA序列识别方法，包括以下步骤：
步骤1、基于高通量测序读数(reads)，得到高频k-mer集合；
步骤2、根据高频k-mer集合，对高通量测序读数进行筛选，保留其中包含高频k-mer多于设定阈值的读数，作为高频读数；
步骤3、使用序列组装工具组装高频读数，得到contigs序列，即重复DNA序列。
进一步地，所述步骤1包括以下步骤：
步骤1.1、将高通量测序读数(reads)中的每一条DNA序列分别分成长度为k的子序列，即k-mer，得到k-mer列表L，由L中互不相同的所有k-mer构成k-mer集合；对k-mer集合中的每一条k-mer，分别统计其在高通量测序读数中出现的频次，得到每条k-mer的频次信息；
步骤1.2、确定高频阈值，去除掉k-mer集合中频次低于高频阈值的k-mer，将保留下的k-mer组合成高频k-mer集合。
进一步地，所述步骤1.2中，将高通量测序读数的测序深度D与覆盖系数c相乘得到高频阈值τ1，c的大小可以设置为大于1的正数，通常设置为1.5～3之间，c值越大表示对高频k-mer的选取越严格。在一种可选的方案中，高通量测序读数的测序深度D是已知的，那么可以跳过之前的测序深度估算部分，使用D与覆盖系数相乘c也可得出高频阈值τ1。遍历k-mer集合，将其中出现频次低于高频阈值τ1的k-mer去除，余下的k-mer组成为高频k-mer集合。
进一步地，有的高通量测序读数数据集在下载时会给出测序深度D，若没有给出，则可以通过以下方法估算测序深度D：
1)将高通量测序读数中的每一条DNA序列分别分成长度为k1的子序列，即k1-mer，得到k1-mer列表L1，由L1中互不相同的所有k1-mer构成集合S1；k1为用于估算测序深度使用的子序列长度；
2)对于S1中的每一条k1-mer，分别统计其在高通量测序读数中出现的频次，将S1中的第i条k1-mer mi在高通量测序读数中出现的频次记为Fi，i∈1，2，...，n，n为S1中k1-mer的个数n；并由所有(mi，Fi)构成集合Sk1；
3)基于集合Sk1，统计其中出现特定频次的k1-mer的个数，将集合Sk1中频次为t的k1-mer的个数记为f(t)；
4)将t和f(t)的关系绘制为“频次——k1-mer数目”曲线图，其中横坐标表示频次t，纵坐标频次为t的k1-mer的个数f(t)；在高通量测序读数中测序错误较少时，曲线通常会呈现近似的泊松分布或高斯分布，随着频次从0开始增加，曲线在开始会有一个斜度较大的下降，随后会有明显的上升，且会上升到一个峰值点，这个峰值点处的横坐标为p；
5)根据p的值估计出高通量测序读数的测序深度D，计算公式如下：
Figure BDA0002068168080000031
式中，p为曲线图中峰值点处横坐标，l为高通量测序得到的所有读数的平均长度。
进一步地，k1的数值设置在15至20之间。
进一步地，如果将k长度设置的过长会对长度较短的高通量测序读数的k-mer频次统计造成影响，如果设置得过短会使测序错误对k-mer的频次统计结果的影响加大，所以本发明将k的数值大小设置为31。
进一步地，所述步骤2中，对每一条读数，根据以下方法确定其是否为高频读数：
将该读数分成长度为k的子序列，得到的所有k-mer依次记为s1，s2，...，sq，q为得到的k-mer总数；依次检查每一条k-mer是否包含在高频k-mer集合中，记录该读数包含在高频k-mer集合中的k-mer数目m，如果该读数首尾两条k-mer(s1和sq)均包含在高频k-mer集合中，并且m/q超过一定的阈值τ2，那么将该读数归为高频读数；对于高通量测序得到的每一条读数，若其为单端读数，则根据以上方法确定其是否为高频读数；若其为双端读数，则根据以上方法分别确定其两条读数是否都为高频读数，若是，则该双端读数为高频读数；将所有高频读数及其对应的质量分数(在组装过程中组装工具会用到质量分数信息)以fastq文件格式存储于文件中，并将此文件称为高频读数文件。
进一步地，所述步骤3中，对于所有contigs序列，滤除掉在组装过程中构图覆盖度(该参数是序列组装构图过程产生的一个可信度数据)较低的序列，保留下的所有contigs序列即为重复DNA序列。
有益效果：
本发明提出的基于高通量测序读数的重复DNA序列识别方法，包括：基于高通量测序读数进行k-mer频次统计，得到在所有在序列中出现的k-mer的集合及每条k-mer的频次信息；根据k-mer的频次结果估算测序序列的测序深度；依据已知的测序深度(有的高通量测序读数数据集在下载时会给出测序深度)或者估算出的测序深度确定高频阈值，去除掉k-mer集合中频次低于高频阈值的k-mer，将保留下的k-mer组合成高频k-mer集合；根据高频k-mer集合，采用一种筛选规则对高通量测序读数进行筛选，使得包含高频k-mer较多的读数保留下来，成为高频读数；使用序列组装工具组装高频读数，得到contigs序列；对contigs序列根据一定的方法进行筛选，保留下的所有contigs序列即为重复DNA序列。容易注意到的是，本发明实施例可以从高通量测序读数中识别重复DNA序列，而无需物种参考序列，可以适用于参考序列未知物种的高通量测序读数。另一方面，大部分序列组装方法或工具都是基于读数组装而不是k-mer组装而设计的，因此本发明先通过高频k-mer选择出高频读数，再通过序列组装方法或工具组装高频读数，从而得到重复DNA序列，而不是直接对高频k-mer进行序列组装，从而提高了识别重复DNA序列的准确率；在采用序列组装方法或工具组装高频读数的过程之后，滤除了组装过程中构图覆盖度较低的序列，进一步提升了识别重复DNA序列的准确率。
附图说明
图1为本发明实施例的重复DNA序列识别方法的流程示意图。
图2为本发明实施例“频次——k-mer数目”曲线示意图。
图3为本发明实施例筛选高频读数流程示意图。
具体实施方式
为了使本技术领域的人员更好地理解本发明方案，下面将结合本发明实施例中的附图，对本发明实施例中的技术方案进行清楚、完整地描述。
首先对本发明实施例中出现的技术名词进行解释如下：
高通量测序(high-throughput sequencing)技术：又称“下一代”测序技术("Next-generation"sequencing，简称NGS)，是对应于传统的桑格测序技术(SangerSequencing)的新一代测序技术，以能一次并行对几十万到几百万条DNA分子进行序列测定和一般读长较短等为标志。目前具有代表性的高通量测序平台有454测序仪、Illumina测序仪、Polonator测序仪和HeliScope测序仪等等。高通量测序技术目前已经在物种的从头测序(de novo sequencing)，小分子RNA测序(small RNA sequencing)，遗传病诊断、基因表达水平的测量、构建进化树、环境基因组学等多个分子生物学领域得到广泛应用。高通量测序分为单端测序(single-end sequencing)和双端测序(paired-end sequencing)两种。由于双端测序的准确率更高，相对于单端测序，双端测序被更广泛的使用。双端测序的流程是，先在序列的一端进行测序，然后反过来测互补序列的另一端，得到的两个序列就成为双端读数。本发明实施例所需的高通量测序读数可以为单端测序或双端测序中的任一种读数。
测序深度：DNA测序的深度(或覆盖度)是描述与已知参考序列成功比对或“覆盖”上的平均读数数目。测序深度可以用测序得到的碱基总量(bp)与基因组大小(Genome)的比值计算得出，是评价测序量的指标之一。
序列组装：序列组装是指将短DNA片段通过比对和拼接等方法重建原始长序列的过程。高通量测序中经常需要序列组装过程，因为高通量测序技术无法一次性读取全基因组，而是根据所使用的技术读取几十到几百个碱基的片段。contigs序列是指通过序列拼接而获得的结果序列。
k-mer：通过DNA测序获得的读数包含的所有长度为k的子序列。给定一条长度为L的测序读数，一共可以得到(L-k+1)条k-mer。k-mer经常用于生物信息学的各种基础处理方法中，例如序列从头组装(de novo assembly)或者序列比对中。本发明实施例中，将k-mer作为最小的序列重复单元来进行频次统计和根据“频次——k-mer数目”曲线图规律对高通量测序读数的测序深度进行估计。
fastq：fastq格式是一种保存生物序列(通常为核酸序列)及其测序质量得分信息的文本格式。序列与质量得分皆由单个ASCII字符表示。目前，fastq格式已经成为了保存高通量测序结果的标准文件格式。
本发明实施例提供一种基于高通量测序读数的重复DNA序列识别方法。
本发明提供了一种基于高通量测序读数的重复DNA序列识别方法，具体步骤如图1所示，包括：
本发明实施例中输入的高通量测序读数为从NCBI SRA数据库中下载得到的果蝇(Drosophila)的全基因组测序读数，在NCBI SRA数据中的文库编号为SRX040484和SRX040486，下文称这两个文库分别为文库一和文库二。
步骤1：估算测序深度，确定高频阈值；
估算输入高通量测序读数的测序深度，计算出区别高频k-mer和非高频k-mer的阈值；
步骤1.1：从测序读数中获取15-mer并统计15-mer频次；
将每一条读数的所有长度为15的子序列收集起来成为15-mer列表L1，将L1中互不相同的所有15-mer组合为集合S1。
计算S1中15-mer的个数n，对于S1中的每一条的15-mer mi(i∈1，2，...，n1)统计其在高通量测序读数中出现的频次记为Fi，并将每一条15-mer和它频次的映射(mi，Fi)组合起来为集合Sk1。
这里将文库一fastq文件中的五对共十条序列取出对步骤1.1的过程进行举例说明。
取出的十条序列如表1所示：
表1
Figure BDA0002068168080000061
十条序列的长度均为74bp，可产生共(L-k+1)*10条k-mer，列表L1共包含600条k-mer。统计L1中每条k-mer的出现次数，集合S1中共含有584条互不相同的k-mer，Sk1集合在高通量测序读数中出现频次较高的部分结果如表2所示。
表2
k-mer 频次
GTGTGCGAAATACGT 11
TGTGCGAAATACGTG 7
TGTGCGAAATACGTC 6
GTGCGAAATACGTCG 5
TGGTGTGCGAAATAC 5
GGTGTGCGAAATACG 4
TGCGAAATACGTGCG 4
GTGCGAAATACGTGC 4
CTGTGCGAAATACGT 3
GTGCGAAATACGTGA 3
TGTGTGCGAAATACG 3
步骤1.2：绘制“频次——15-mer数目”曲线图，确定峰值处横坐标；
基于集合Sk1的结果，统计数据的15-mcr的频次频数直方图数据。具体来讲，对于特定的频次t，用f(t)表示出现频次为t的15-mer的个数，计算公式如下：
Figure BDA0002068168080000071
如果Fi＝t，则Cit＝1，否则Cit＝0
将t和f(t)的关系绘制为“频次——15-mcr数目”曲线图，其中横坐标表示频次t，纵坐标频次为t的15-mer的个数f(t)。在高通量测序读数中测序错误较少时，曲线通常会呈现近似的泊松分布或高斯分布，随着频次从0开始增加，曲线在开始会有一个斜度较大的下降，随后会有明显的上升，且会上升到一个峰值点，这个峰值点处的横坐标为p。
本实施例在果蝇测序数据的基础上绘制的“频次——15-mer数目”曲线图，如图2所示，可见在最高峰值处的横坐标p约为33。
步骤1.3：估算测序深度；
根据p的值估计出高通量测序的测序深度D，计算公式如下：
Figure BDA0002068168080000072
式中，p为曲线图中峰值点处横坐标，l为高通量测序得到的所有序列的平均长度。需要指出的是，k1的数值设置在15至20之间可以得出良好的估算效果。
本发明实施例中，p值为33，l值为82，k1值为15，计算出的测序深度D取整后约为40。
步骤1.4：根据测序深度计算高频阈值；
将估算深度的D与覆盖系数c相乘得出高频阈值τ1，c的大小可以设置为大于1的数，通常设置为1.5～3之间是较为合理的，c值越大表示对高频k-mer的选取越严格。在一种可选的方案中，高通量测序的测序深度D是已知的，那么可以无需执行步骤1.1至步骤1.3，使用D与覆盖系数相乘就可得出高频阈值τ1。
本发明实施例中，将c值设置为2.2，计算出高频阈值τ1为88。
步骤2：从测序文件中获取31-mer并统计31-mer频次；
与步骤1.1相似，将每一条读数序列中所有长度为31的子序列收集起来成为31-mer列表L2，将L2中互不相同的所有31-mer组合为集合S2。
计算S2中31-mer的个数n，对于S2中的每一条的31-mer mi(i∈1，2，...，n2)统计其在高通量测序读数中出现的频次记为Fi，并将每一条k-mer及其频次的映射(mi，Fi)组合起来为集合Sk2。
在本发明实施例中，Sk2集合中共包括757881473条31-mer及其频次的映射。
步骤3：根据高频阈值获得高频31-mer集合；
遍历集合S2将其中出现频次低于高频阈值τ1的31-mer从集合S2中去除，余下的31-mer组成为高频31-mer集合Sh。
在本发明实施例中，Sh中共2251900条高频31-mer。
步骤4：从测序序列中根据高频31-mer集合筛选出高频读数；
对于任一条读数，设包含的所有31-mer序列分别为s1，s2，...，sq，检查每一条31-mer是否包含在高频31-mer集合Sh中，记录包含在高频31-mer集合Sh中的数目m，如果该读数首尾的两条31-mer(s1和sq)都包含在高频31-mer集合Sh中，并且m/q超过一定的阈值τ2，那么将序列归为高频读数。对于单端读数，将所有高频读数及其对应的质量分数以fastq文件格式存储于文件中；对于双端读数，当一对读数的两条序列都是高频读数时，将两条读数及其质量分数以fastq文件格式存储于文件中，并将此文件称为高频读数文件。此步骤的流程示意图如图3所示。
在本发明实施例中，将t2设置为80％，共得到1475722条高频读数，高频读数文件大小为730Mb。
步骤5：组装高频读数得到contigs序列；
使用序列组装工具SPAdes对高频读数文件进行序列组装，得到组装结果contigs序列。需要说明的是，在本发明实施例中，可以采用现有技术中提供的序列组装工具或方法进行序列组装，本发明实施例对此不做具体限定。
在本实施例中，使用SPAdes工具对高频读数进行序列组装，得到52791条contigs序列。
步骤6：筛选contigs序列得到重复DNA序列
对于contigs序列中的每一条序列，滤除掉在组装过程中构图覆盖度较低的序列，得到识别的重复DNA序列。
本发明实施例中将组装过程中构图覆盖度低于高频阈值两倍(τ1*2)的contigs序列去除，得到50680条重复DNA序列。
实验结果评估
本发明实施例在使用本发明所述方法的情况下，将所得重复DNA序列与RepARK方法和REPdenovo方法进行比较。具体而言，使用本发明所述方法与RepARK方法和REPdenovo方法在使用相同高通量测序读数作为输入的前提下，对识别的重复DNA序列的序列属性、参考序列比对率和多比对率、覆盖Repbase数据库中序列的比例等方面进行比较。Repbase(https：//www.girinst.org/repbase/)是一个重复DNA序列数据库，它提供了部分常见物种的重复DNA序列参考，覆盖Repbase数据库中序列的比例越大，说明重复DNA序列识别方法得到的重复DNA序列越精确。通过以上几种方法，对本发明实施例识别的重复DNA序列质量进行评估。
以下分别有五个本发明实施例产生的重复DNA序列与RepARK方法和REPdenovo方法识别的重复DNA序列进行比较。五个本发明实施例所使用的高通量测序读数分别来自NCBI SRA数据库(https：//www.ncbi.nlm.nih.gov/sra/)和GAGE网站(http：//gage.cbcb.umd.edu/)，分别为果蝇(Drosophila)、酵母菌(Saccharomyces cerevisiae)、切叶蚁(Acromyrmex echinatior)、小鼠(Mus musculus)和人类第十四条染色体(HumanChromosome 14)。五个高通量测序读数数据集的具体信息如表3所示。
表3高通量测序读数数据集信息
Figure BDA0002068168080000091
表4五个数据集上的重复DNA序列比较结果
Figure BDA0002068168080000092
Figure BDA0002068168080000101
表5重复DNA序列在参考序列上的覆盖率比较结果
Figure BDA0002068168080000102
表6重复DNA序列在Repbase上的覆盖率比较结果
Figure BDA0002068168080000103
表4中包含五个物种的本发明实施例从重复DNA序列的重复DNA序列数目、文库大小、最大/最小长度、N50、N90、比对率和多比对率几个方面对重复DNA序列的比较结果。其中N50和N90是指将重复DNA序列按照长度从长至短排序，并将序列长度从长至短相加，当得到的序列长度和达到所有序列总长度的50％和90％时，最后一个相加的序列长度，即为重复DNA序列的N50和N90。比对率和多比对率分别是指将所有重复DNA序列通过序列比对工具Bowtie2比对到其所属物种的参考序列上，可以比对到参考序列上一个以上位置和两个以上位置的重复DNA序列数目占重复DNA序列总数的比例。使用RepeatMasker工具对本发明实施例生成的重复DNA序列和RepARK以及REPdenovo的重复DNA序列结果在部分物种的RepBase和参考序列上的覆盖情况进行比较，结果包含在表5和表6中。从结果中可以看出，本发明实施例中得到的重复DNA序列相比于RepARK工具和REPdenovo工具来讲，主要有几个优势：长度更长，体现在序列的N50、N90和最长序列长度上；准确性更高，体现在比对率和多比对率上，尤其是多比对率，多次比对到物种的参考序列上是体现重复DNA序列准确性的重要指标之一；更全面，体现在文库大小、覆盖参考序列的比例、覆盖RepBase数据库中序列的比例更大。需要说明的是，在本发明实施例中，可以采用现有技术中提供的序列比对算法或工具对序列的比对率和多比对率进行比较，本发明实施例对此不做具体限定。
需要强调的是，本发明所述的实例是说明性的，而不是限定性的，因此本发明不限于具体实施方式中所述的实例，凡是由本领域技术人员根据本发明的技术方案得出的其他实施方式，不脱离本发明宗旨和范围的，不论是修改还是替换，同样属于本发明的保护范围。

Claims (7)

1.一种基于高通量测序读数的重复DNA序列识别方法，其特征在于，包括以下步骤：
步骤1、基于高通量测序读数，得到高频k-mer集合；
步骤2、根据高频k-mer集合，对高通量测序读数进行筛选，保留其中包含高频k-mer多于设定阈值的读数，作为高频读数，其中，所述步骤2中，对每一条读数，根据以下方法确定其是否为高频读数：
将该读数分成长度为k的子序列，得到的所有k-mer依次记为s1，s2，...，sq，q为得到的k-mer总数；依次检查每一条k-mer是否包含在高频k-mer集合中，记录该读数包含在高频k-mer集合中的k-mer数目m，如果s1和sq均包含在高频k-mer集合中，并且m/q超过阈值τ2，那么将该读数归为高频读数；
对于高通量测序得到的每一条读数，若其为单端读数，则根据以上方法确定其是否为高频读数；若其为双端读数，则根据以上方法分别确定其两条读数是否都为高频读数，若是，则该双端读数为高频读数；
步骤3、使用序列组装工具组装高频读数，得到contigs序列，从而得到重复DNA序列。
2.根据权利要求1所述的基于高通量测序读数的重复DNA序列识别方法，其特征在于，所述步骤1包括以下步骤：
步骤1.1、将高通量测序读数中的每一条DNA序列分别分成长度为k的子序列，即k-mer，得到的所有k-mer构成k-mer集合；对k-mer集合中的每一条k-mer，分别统计其在高通量测序读数中出现的频次，得到每条k-mer的频次信息；
步骤1.2、确定高频阈值，去除掉k-mer集合中频次低于高频阈值的k-mer，保留下的所有k-mer构成高频k-mer集合。
3.根据权利要求2所述的基于高通量测序读数的重复DNA序列识别方法，其特征在于，所述步骤1.2中，将高通量测序读数的测序深度D与覆盖系数c相乘得到高频阈值。
4.根据权利要求3所述的基于高通量测序读数的重复DNA序列识别方法，其特征在于，通过以下方法估算高通量测序读数的测序深度D：
1)将高通量测序读数中的每一条DNA序列分别分成长度为k1的子序列，即k1-mer，得到的所有k1-mer构成集合S1；
2)对于S1中的每一条k1-mer，分别统计其在高通量测序读数中出现的频次，将S1中的第i条k1-mer mi在高通量测序读数中出现的频次记为Fi，i∈1，2，...，n，n为S1中k1-mer的个数n；并由所有(mi，Fi)构成集合Sk1；
3)基于集合Sk1，统计其中出现特定频次的k1-mer的个数，将集合Sk1中频次为t的k1-mer的个数记为f(t)；
4)将t和f(t)的关系绘制为“频次——k1-mer数目”曲线图，其中横坐标表示频次t，纵坐标频次为t的k1-mer的个数f(t)；将曲线图中峰值点处的横坐标记为p；
5)根据p的值估计出高通量测序读数的测序深度D，计算公式如下：
Figure FDA0002637560200000021
式中，l为高通量测序得到的所有读数的平均长度。
5.根据权利要求4所述的基于高通量测序读数的重复DNA序列识别方法，其特征在于，k1的数值设置在15至20之间。
6.根据权利要求1～5中任一项所述的基于高通量测序读数的重复DNA序列识别方法，其特征在于，k的数值设置为31。
7.根据权利要求1所述的基于高通量测序读数的重复DNA序列识别方法，其特征在于，所述步骤3中，对于所有contigs序列，滤除掉在组装过程中构图覆盖度低于设定阈值的序列，保留下的所有contigs序列即为重复DNA序列。
CN201910428254.6A 2019-05-22 2019-05-22 一种基于高通量测序读数的重复dna序列识别方法 Active CN110066862B (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN201910428254.6A CN110066862B (zh) 2019-05-22 2019-05-22 一种基于高通量测序读数的重复dna序列识别方法

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN201910428254.6A CN110066862B (zh) 2019-05-22 2019-05-22 一种基于高通量测序读数的重复dna序列识别方法

Publications (2)

Publication Number Publication Date
CN110066862A CN110066862A (zh) 2019-07-30
CN110066862B true CN110066862B (zh) 2021-02-12

Family

ID=67371132

Family Applications (1)

Application Number Title Priority Date Filing Date
CN201910428254.6A Active CN110066862B (zh) 2019-05-22 2019-05-22 一种基于高通量测序读数的重复dna序列识别方法

Country Status (1)

Country Link
CN (1) CN110066862B (zh)

Citations (4)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN105637099A (zh) * 2013-08-23 2016-06-01 考利达基因组股份有限公司 使用短读段的长片段从头组装
CN106021985A (zh) * 2016-05-17 2016-10-12 杭州和壹基因科技有限公司 一种基因组数据压缩方法
CN108491687A (zh) * 2018-03-22 2018-09-04 中南大学 基于contig质量评估分类及图优化的scaffolding方法
CN108660200A (zh) * 2018-05-23 2018-10-16 北京希望组生物科技有限公司 一种检测短串联重复序列扩张的方法

Family Cites Families (1)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
US10198454B2 (en) * 2014-04-26 2019-02-05 Bonnie Berger Leighton Quality score compression for improving downstream genotyping accuracy

Patent Citations (4)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN105637099A (zh) * 2013-08-23 2016-06-01 考利达基因组股份有限公司 使用短读段的长片段从头组装
CN106021985A (zh) * 2016-05-17 2016-10-12 杭州和壹基因科技有限公司 一种基因组数据压缩方法
CN108491687A (zh) * 2018-03-22 2018-09-04 中南大学 基于contig质量评估分类及图优化的scaffolding方法
CN108660200A (zh) * 2018-05-23 2018-10-16 北京希望组生物科技有限公司 一种检测短串联重复序列扩张的方法

Non-Patent Citations (3)

* Cited by examiner, † Cited by third party
Title
Comparison of the two major classes of assembly algorithms: overlap^layout^consensus and de-bruijn-graph;Zhenyu Li等;《BRIEFINGS IN FUNCTIONAL GENOMICS》;20111219;第II卷(第I期);25-37 *
DNA 片段拼接中重复序列算法研究;王磊等;《计算机科学》;20061231;第33卷(第7期);164-170 *
Improving de novo Assembly Based on Read Classification;Xingyu Liao等;《IEEE/ACM TRANSACTIONS ON COMPUTATIONAL BIOLOGYAND BIOINFORMATICS》;20180731;第17卷(第1期);177-188 *

Also Published As

Publication number Publication date
CN110066862A (zh) 2019-07-30

Similar Documents

Publication Publication Date Title
US11817180B2 (en) Systems and methods for analyzing nucleic acid sequences
CN114999573B (zh) 一种基因组变异检测方法及检测系统
CN110997937A (zh) 具有可变长度非随机独特分子标识符的通用短衔接子
CN112466404B (zh) 一种宏基因组重叠群无监督聚类方法及系统
EP2923293B1 (en) Efficient comparison of polynucleotide sequences
JP6141310B2 (ja) 強固な変異体特定および検証
US20190177719A1 (en) Method and System for Generating and Comparing Reduced Genome Data Sets
CN110088840B (zh) 校正核酸序列读数的重复区域中的碱基调用的方法、系统和计算机可读媒体
CN115719616A (zh) 一种病原物种特异性序列的筛选方法及系统
CN110066862B (zh) 一种基于高通量测序读数的重复dna序列识别方法
CN110232951B (zh) 判断测序数据饱和的方法、计算机可读介质和应用
CN116246703A (zh) 一种核酸测序数据的质量评估方法
WO2019132010A1 (ja) 塩基配列における塩基種を推定する方法、装置及びプログラム
CN111798922B (zh) 基于重测序数据中多态性位点密度鉴定小麦育种的基因组选择利用区间的方法
CN111094591A (zh) 用于对生物分子进行测序的方法
US20180121600A1 (en) Methods, Systems and Computer Readable Storage Media for Generating Accurate Nucleotide Sequences
Aldawiri et al. A Novel Approach for Mapping Ambiguous Sequences of Transcriptomes
CN114420214A (zh) 核酸测序数据的质量评估方法和筛选方法
RU2799778C9 (ru) Способ определения показателя, коррелированного с вероятностью того, что два мутированных прочтения последовательности происходят от одной и той же содержащей мутации последовательности
KR102111731B1 (ko) 핵산 시퀀스를 분석하는 방법 및 장치
CN114300055B (zh) 优化的宏基因组纳米孔测序数据定量方法
CN112634983B (zh) 病原物种特异pcr引物优化设计方法
CN110600083B (zh) 基于无拼接组装wgs数据的醋酸钙—鲍曼不动杆菌复合群鉴定方法
RU2799778C1 (ru) Способ определения показателя, коррелированного с вероятностью того, что два мутированных прочтения последовательности происходят от одной и той же содержащей мутации последовательности
CN116230086B (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