CN115331736B - 基于文本匹配延伸高通量测序基因的拼接方法 - Google Patents

基于文本匹配延伸高通量测序基因的拼接方法 Download PDF

Info

Publication number
CN115331736B
CN115331736B CN202210856831.3A CN202210856831A CN115331736B CN 115331736 B CN115331736 B CN 115331736B CN 202210856831 A CN202210856831 A CN 202210856831A CN 115331736 B CN115331736 B CN 115331736B
Authority
CN
China
Prior art keywords
sequence
sequences
splicing
sequencing
thousand
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
CN202210856831.3A
Other languages
English (en)
Other versions
CN115331736A (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.)
Foshan University
Original Assignee
Foshan 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 Foshan University filed Critical Foshan University
Priority to CN202210856831.3A priority Critical patent/CN115331736B/zh
Publication of CN115331736A publication Critical patent/CN115331736A/zh
Application granted granted Critical
Publication of CN115331736B publication Critical patent/CN115331736B/zh
Active legal-status Critical Current
Anticipated expiration legal-status Critical

Links

Classifications

    • G—PHYSICS
    • G16—INFORMATION AND COMMUNICATION TECHNOLOGY [ICT] SPECIALLY ADAPTED FOR SPECIFIC APPLICATION FIELDS
    • G16B—BIOINFORMATICS, i.e. INFORMATION AND COMMUNICATION TECHNOLOGY [ICT] SPECIALLY ADAPTED FOR GENETIC OR PROTEIN-RELATED DATA PROCESSING IN COMPUTATIONAL MOLECULAR BIOLOGY
    • G16B30/00—ICT specially adapted for sequence analysis involving nucleotides or amino acids
    • G16B30/20—Sequence assembly
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F16/00—Information retrieval; Database structures therefor; File system structures therefor
    • G06F16/30—Information retrieval; Database structures therefor; File system structures therefor of unstructured textual data
    • G06F16/33—Querying
    • G06F16/3331—Query processing
    • G06F16/334—Query execution
    • G—PHYSICS
    • G16—INFORMATION AND COMMUNICATION TECHNOLOGY [ICT] SPECIALLY ADAPTED FOR SPECIFIC APPLICATION FIELDS
    • G16B—BIOINFORMATICS, i.e. INFORMATION AND COMMUNICATION TECHNOLOGY [ICT] SPECIALLY ADAPTED FOR GENETIC OR PROTEIN-RELATED DATA PROCESSING IN COMPUTATIONAL MOLECULAR BIOLOGY
    • G16B30/00—ICT specially adapted for sequence analysis involving nucleotides or amino acids
    • G16B30/10—Sequence alignment; Homology search
    • G—PHYSICS
    • G16—INFORMATION AND COMMUNICATION TECHNOLOGY [ICT] SPECIALLY ADAPTED FOR SPECIFIC APPLICATION FIELDS
    • G16B—BIOINFORMATICS, i.e. INFORMATION AND COMMUNICATION TECHNOLOGY [ICT] SPECIALLY ADAPTED FOR GENETIC OR PROTEIN-RELATED DATA PROCESSING IN COMPUTATIONAL MOLECULAR BIOLOGY
    • G16B50/00—ICT programming tools or database systems specially adapted for bioinformatics
    • G16B50/10—Ontologies; Annotations
    • G—PHYSICS
    • G16—INFORMATION AND COMMUNICATION TECHNOLOGY [ICT] SPECIALLY ADAPTED FOR SPECIFIC APPLICATION FIELDS
    • G16B—BIOINFORMATICS, i.e. INFORMATION AND COMMUNICATION TECHNOLOGY [ICT] SPECIALLY ADAPTED FOR GENETIC OR PROTEIN-RELATED DATA PROCESSING IN COMPUTATIONAL MOLECULAR BIOLOGY
    • G16B25/00—ICT specially adapted for hybridisation; ICT specially adapted for gene or protein expression
    • G16B25/10—Gene or protein expression profiling; Expression-ratio estimation or normalisation
    • G—PHYSICS
    • G16—INFORMATION AND COMMUNICATION TECHNOLOGY [ICT] SPECIALLY ADAPTED FOR SPECIFIC APPLICATION FIELDS
    • G16B—BIOINFORMATICS, i.e. INFORMATION AND COMMUNICATION TECHNOLOGY [ICT] SPECIALLY ADAPTED FOR GENETIC OR PROTEIN-RELATED DATA PROCESSING IN COMPUTATIONAL MOLECULAR BIOLOGY
    • G16B25/00—ICT specially adapted for hybridisation; ICT specially adapted for gene or protein expression
    • G16B25/20—Polymerase chain reaction [PCR]; Primer or probe design; Probe optimisation
    • Y—GENERAL TAGGING OF NEW TECHNOLOGICAL DEVELOPMENTS; GENERAL TAGGING OF CROSS-SECTIONAL TECHNOLOGIES SPANNING OVER SEVERAL SECTIONS OF THE IPC; TECHNICAL SUBJECTS COVERED BY FORMER USPC CROSS-REFERENCE ART COLLECTIONS [XRACs] AND DIGESTS
    • Y02—TECHNOLOGIES OR APPLICATIONS FOR MITIGATION OR ADAPTATION AGAINST CLIMATE CHANGE
    • Y02A—TECHNOLOGIES FOR ADAPTATION TO CLIMATE CHANGE
    • Y02A90/00—Technologies having an indirect contribution to adaptation to climate change
    • Y02A90/10—Information and communication technologies [ICT] supporting adaptation to climate change, e.g. for weather forecasting or climate simulation

Landscapes

  • Engineering & Computer Science (AREA)
  • Physics & Mathematics (AREA)
  • Life Sciences & Earth Sciences (AREA)
  • Theoretical Computer Science (AREA)
  • Health & Medical Sciences (AREA)
  • Biotechnology (AREA)
  • Medical Informatics (AREA)
  • Biophysics (AREA)
  • Bioinformatics & Cheminformatics (AREA)
  • Bioinformatics & Computational Biology (AREA)
  • Spectroscopy & Molecular Physics (AREA)
  • Evolutionary Biology (AREA)
  • General Health & Medical Sciences (AREA)
  • Analytical Chemistry (AREA)
  • Proteomics, Peptides & Aminoacids (AREA)
  • Chemical & Material Sciences (AREA)
  • Databases & Information Systems (AREA)
  • Bioethics (AREA)
  • Computational Linguistics (AREA)
  • Data Mining & Analysis (AREA)
  • General Engineering & Computer Science (AREA)
  • General Physics & Mathematics (AREA)
  • Measuring Or Testing Involving Enzymes Or Micro-Organisms (AREA)

Abstract

本发明公开了一种基于文本匹配延伸高通量测序基因的拼接方法，涉及生物信息学领域。该拼接方法包括：获取测序序列，将其编号、组合，并选取种子序列；将种子序列第一侧的预设长度的序列隔开，得到查询序列，进行文本匹配，找到的相同序列至少两条或以上，然后合并。再用查询序列进行首尾拼接，选取最长的拼接，比较拼接序列与原序列头部，确定不是相同序列则进行新一轮首尾拼接。将序列反转互补后进行相同的拼接。本发明提供了一种精确、快速简易的高通量测序序列拼接的有效方法，可广泛用于转录组序列拼接、可变剪接、全转录谱等。

Description

基于文本匹配延伸高通量测序基因的拼接方法
技术领域
本发明涉及生物信息学领域，尤其涉及一种基于文本匹配延伸高通量测序基因的拼接方法。
背景技术
深入了解组织细胞内的基因活动首先要对RNA表达水平进行分析，高通量RNA测序(RNA-Seq)开展十多年来为了解基因调控提供了十分高效的途径。测序提供了大量RNA表达的短序列，长度通常为100-300个核苷酸(不同平台有差异)，而从这些短序列获得基因表达水平并将短序列拼接得到长RNA序列成为大多数研究人员开发的重要目标。序列组装通常基于两种策略，即基于参考基因组(Mapping-first，ab initio)的拼接和从头(Assembling-first,de novo)拼接，随着转录组测序的广泛开展，从头组装转录谱担当了重要角色。目前主流拼接方法有：Bridger、Oases、Trinity等，其基本上都是采用kmer(k＝20)比对进行匹配拼接。申请人在前序申请中也采用kmer序列进行了文本匹配分析了RNA-Seq序列的基因表达量(CN108388772B)。
RNA-Seq测序提供的短序列不管是单端或者双端测序，RNA测序建库进行反转录获得互补DNA，任意片段(包括其互补链)被切割和测序的几率是相同的，表达水平高的RNA其被测序几率就较高，可以获得拼接延伸机会就大。由于高通量的测序错误率不低，只是利用kmer节点进行重叠拼接延伸是不够的，利用较长的拼接重叠区是必要的。此外，序列拼接之后应可以追溯原始序列来源，另外可以根据匹配的原始序列高度相似性设计荧光定量PCR扩增的引物，不同组织的拼接序列进行比对可以发现可变剪接。多数从头拼接程序仍然受限于kmer算法，可变剪接体，测序错误率，不同表达水平转录本，计算内存等影响，转录组拼接需要进一步提高序列拼接准确率。
发明内容
本发明所要解决的技术问题在于，提供一种基于文本匹配延伸高通量测序基因的拼接方法，其拼接效率高，简单快捷，拼接过程可追溯，可用于进一步PCR扩增、可变剪接分析等。
为了解决上述技术问题，本发明提供了一种基于文本匹配延伸高通量测序基因的拼接方法，其包括：
S1：通过高通量测序平台获取待分析样品的测序序列；
S2：将测序序列加编号，打散，随机组合；
S3：选取10万条测序序列作为种子序列；
S4：将种子序列第一侧预设长度的序列隔开，得到查询序列；
S5：将查询序列与一组或多组100万条测序序列进行文本匹配，并保留至少匹配两次的多条序列作为匹配序列；
S6：采用查询序列作为拼接点，将匹配序列与种子序列第一侧拼接，并选取拼接后长度最长的序列作为拼接序列；
S7：比较拼接序列第一侧预设长度的序列与种子序列第一侧预设长度的序列；
S8：若两者不同，则以拼接序列作为新的种子序列，并进入步骤S4；直至所得拼接序列第一预设长度的序列与种子序列第一侧预设长度的序列相同，或进行至少5次拼接后进入S9；
S9：若两者相同，则将所得拼接序列反转后作为新的种子序列，并进入步骤S4，直至完成种子序列两侧的拼接。
作为上述技术方案的改进，步骤S4中，将每条种子序列左侧40mer每隔5mer切下5个20mer序列，得到查询序列。
作为上述技术方案的改进，步骤S5中，将查询序列与一组或多组100万条测序序列进行文本匹配，得到返回序列；每个查询序列得到的返回序列进行排序并计数重复数，保留超过2条相同序列以上的序列组作为匹配序列。
作为上述技术方案的改进，步骤S2中，采用逐级按照每10万→100万→5万条进行随机排序并按随机方式合并序列，其中，每个100万分成若干个目录单独进行每1万条的随机排序，再随机组合合并所有序列。
作为上述技术方案的改进，步骤S1中，所述高通量测序平台选用Roche的454测序平台、Illumina的HySeq4000测序平台和ABI的SOLiD测序平台中的一种或多种。
作为上述技术方案的改进，所述待分析样品为动物、植物、微生物的DNA或RNA提取物；或者，所述待分析样本为水、大气、土壤中的微生物的DNA或RNA提取物。
作为上述技术方案的改进，步骤S5中，采用开源代码系统内嵌shell命令进行文本过滤；
所述shell命令选用awk、sed、sort、grep、tr、split、comm、paste、cat中的一种或多种。
相应的，本发明还公开了一种基于文本匹配延伸高通量测序基因的拼接系统，其包括：
高通量测序平台，用于获取待分析样品的测序序列；
数据处理系统：所述的数据处理系统用于实现上述的拼接方法。
作为上述技术方案的改进，所述数据处理系统为开源代码系统。
相应的，本发明还公开了一种计算机可读介质，其上存储有计算机指令，其特征在于，该指令被处理执行时实现上述的拼接方法的步骤。
实施本发明，具有如下有益效果：
本发明采取开源系统文本过滤命令，对测序序列两侧各从头隔开预设长度得到的查询序列，采用查询序列与100万条测序序列进行文本匹配，找到的相同序列至少两条或以上然后合并。然后用查询序列分别进行首尾拼接，选取最长的拼接，比较拼接序列与原序列头部，确定不是相同序列则进行新一轮首尾拼接。将序列反转互补后进行相同的拼接。实施本发明，可利用150bp的RNA-Seq序列可获得最长约1.3kb拼接序列，约十分之一序列没有获得匹配延伸，若增加目标查找序列则部分可以获得延伸，10万条分开多个如10个1万条在不同窗口运行6天即可完成序列拼接。拼接序列与用Trinity拼接unigene序列进行比较能够很好匹配，利用匹配序列进行远程blast比对，返回了更多的高打分匹配注释。因而，本发明提供了一种精确、快速简易的高通量测序序列拼接的有效方法，可广泛用于转录组序列拼接、可变剪接、全转录谱等分析。
附图说明
图1是本发明实施例1中拼接方法的流程图；
图2是本发明实施例1中采用首尾40mer内5个20mer与仅仅采用首尾20mer获得拼接序列长度比较图；其中，三角形标志为采用首尾40mer内5个20mer；圆形标志为仅采用首尾20mer；
图3是本发明实施例1中10万条序列两端40mer匹配100万条序列首尾拼接序列长度分布图；
图4是本发明实施例1中匹配100万条没有获得延伸的50mer(首1万条中1147条)用300万条匹配延伸结果图；
图5是采用Trinity拼接序列长度分布图；
图6是本发明实施例1中拼接序列及Trinity unigene序列进行本地blast获得序列比对匹配图；其中，A为Trinity unigene序列，B为实施例1所得序列，C为实施例1序列比对Trinity unigene序列；
图7是随机选取一条拼接序列(15000018)与所有拼接序列和Trinity unigene序列进行本地blast匹配序列正负链Clustal-Omiga多序列比对结果图；
图8是本发明实施例1、Trinity所得拼接序列长度与注释数量的关系图，其中，三角形标志为实施例1，×标志为Trinity。
具体实施方式
为使本发明的目的、技术方案和优点更加清楚，下面将结合附图对本发明作进一步地详细描述。仅此声明，本发明在文中出现或即将出现的上、下、左、右、前、后、内、外等方位用词，仅以本发明的附图为基准，其并不是对本发明的具体限定。
实施例1拼接方法
1材料与方法
1.1材料
操作系统：开放源代码系统FreeBSD 12(https://www.freebsd.org/)，安装在ThinkSystem SR650服务器(内存64G，双CPU为Intel Xeon Silver 4110 8C 85W2.1GHz处理器)。
文本过滤命令：FreeBSD系统内嵌shell命令，如awk,sed,sort,grep,tr,split,comm,paste,cat等等。
RNA-Seq测序数据：由北京奥维森基因科技有限公司完成采用Illumina二代高通量测序平台，PE150测序策略获得的香叶天竺葵叶片精准转录组测序数据，长度150mer，fastq.gz格式的双端“clean data”序列数据大小约4G。公司提供了由Trinity程序获得的5个叶片样品测序的拼接序列(unigene)及其注释。
1.2拼接方法(具体流程参图1)
(1)序列预处理：采用逐级按照每10万→100万→5万条进行随机排序并按随机方式合并序列，其中，每个100万分成若干个目录单独进行每1万条的随机排序，再随机组合合并所有序列。测序公司提供的“clean data”程序，先去除超过10个的多聚A/C/G/T序列，每个测序样品序列按照每万条进行随机打散，分别切割为每100万条序列组，随机选取3个100万条进行进一步分析。
(2)算法分析：
选用两侧40mer序列作为延伸种子序列，将两侧的40mer隔开5mer切割得到了5个20mer，然后将其与目标100万条序列进行匹配查找。从返回匹配序列中可以进行首尾连接，若其中找到的匹配序列只是找到一条则认为其仍然只是从前面得到的5个20mer的一条找到的，以其进行延伸则显得匹配长度太短，故须至少使用2条或以上更多的匹配序列进行延伸。因而，两侧的40mer的5个20mer匹配长度范围包括25～40mer，即找到的两个相同匹配序列可以是从开始第1和第2个20mer到首尾即第1个和第5个20mer返回的匹配序列。如果查找返回的序列与原始序列头部20mer相同，则程序终止查找。将序列反转互补进行查找就可以获得右侧延伸。获得的拼接序列两侧可以继续进行多次延伸，获得更长的延伸拼接。
(3)程序运行：全部利用shell脚本进行序列拼接过程。种子序列为3个100万条序列第一组序列的开始10万条，目标匹配查找序列为第二组100万条序列。如果种子序列没有得到延伸，则可以考虑扩大利用3组或更多的100万条序列合并进行匹配。拼接延伸从左侧开始，10万条种子序列(含id)按顺序进行开始40mer进行切割获得5个20mer短序列，将其分别用于查找目标100万条序列获得匹配序列，每个20mer查找匹配序列得到返回序列进行排序并计数重复数，保留超过2条相同序列以上的序列组。5个20mer均用于首尾拼接，选择最长的拼接序列用于下一轮首尾拼接，5次首尾拼接后将拼接序列反转互补，也同样进行首尾拼接5次，最后反转互补恢复种子序列原向保存拼接序列于同一个文件(fasta格式)。
对实施例1的拼接方法进行比对与分析，具体如下：
(一)拼接效率
基于上述实施例1的拼接方法，利用每万条准随机打散的10万条序列查找100万条RNA-Seq测序序列，完成匹配查找和拼接过程耗时可以根据服务器的内存和CPU运算效率估算，如果分开每万条一个窗口运行(10个进程)，一周内可以完成上述拼接工作。服务器可以同时运行若干个窗口，50个窗口同时拼接对服务器计算资源调配不会受到太大影响(服务器内存64G)。
从左侧40mer内隔开5mer选取5个20mer序列，同时查找100万条序列内匹配5个20mer的重叠群序列，因而出现同一个序列被匹配至少一次(5个20mer只有1个得到匹配)，用匹配两个或多个20mer(返回匹配同一序列至少2次)与只是利用仅仅匹配一个20mer获得拼接序列相比则明显拓宽了拼接重叠区长度，拼接的精确率也大大提高，相当于是从25～40mer(开头两个20mer～全部5个20mer)用于拼接延伸。从拼接长度看，最长可以达到1.3kb(首尾均延伸5次)，约10％序列没有获得延伸。与单纯直接利用首尾20mer延伸进行了比较，40mer的拼接效率显著提升(图2)。
(二)拼接所得序列长度分析
进一步的，用了1个叶片样品RNA测序10万条150mer序列两端40mer内5个20mer与100万条目标序列匹配进行首尾拼接获得了拼接序列(分开10个终端窗口运行了6天)，长度分布情况见图3。没有获得拼接(150mer)的占了11.6％，151-300占了13.5％，400至最长的1351占了67.4％，即三分之二以上的拼接长度超过400mer，而其中600-1000之间的占了37.2％。
(三)文本匹配采用测序序列数量的影响
拼接长度较短的序列其基因表达量也较低，增加目标匹配序列数量将会获得更多和更长的拼接，部分序列可能由于表达量过低仍然无法得到较长延伸，图4为增加目标程序序列为300万条后上图没有延伸的部分序列(150mer，10万条的首1万条里没有延伸的1147条)拼接长度分布，仍有近一半150mer序列没有获得延伸，151-500长度也约占了半数，因而10万条序列没有拼接的150mer在增加目标序列为300万条后再次拼接仍无法获得拼接的数量将降至5％左右。增加目标匹配序列数量用于获得10万条种子序列延伸会延缓拼接进程，因而只是用短拼接序列及150mer没有延伸序列进行再次拼接就会较好节省计算资源(同样多终端运行)，同样最好也应对300mer及以下的拼接序列进行再次延伸。
(四)与Trinity拼接unigene序列的比较
所有Trinity程序拼接获得的unigene序列95805条，序列长度均为200mer以上，最长达到17kb，长度分布如图5。从图5可以看出，200-400mer短片段占了45％，随着拼接长度增加序列数量急剧下降，而上面图3本方法拼接出现了较集中的分布即600-1000mer的占了37％，而图5的只有19％。图5中最长拼接unigene达到了17kb，而从图中随着长度增加出现陡降的情况看，拼接的整体效果并不很令人满意。
为了进一步了解获得拼接序列相互之间比对关系，对Trinity拼接获得的unigene序列(95805条)进行本地blast比对(打分没有超过100分的，比对长度少于100个核苷酸)，57条没有比对结果，38625条(40％)没有找到匹配比对序列，匹配1条以上序列占60％(57124条)，最多匹配序列数100条，结果见图6。同时也对10万条序列进行本地blast比对，87条没有比对结果，14378条(14％)序列没有找到匹配(其中7440为没有延伸的150mer种子序列)，其他可以找到超过1条或以上匹配序列(86％，最高匹配序列数250条)(图6)，blast比对匹配序列数量多需要进一步进行匹配拼接序列合并。利用本方法10万条序列(拼接序列及没有拼接的种子序列)与Trinity拼接unigene序列进行本地blast比对，2114条序列没有获得比对(1289条为没有延伸的150mer序列，长度超过300的283条)，其他97886条匹配1条以上(图6)，匹配Trinity的unigene序列去除重复后共3.6万条。Trinity的unigene长度仅选取了200mer以上的(图5)且拼接源序列来自所有5个样品，本方法分析来源于1个样品测序序列(每个样品数据均应反映出各自的基因活动情况)，且拼接长度低于200mer的占了15％(图3)。
由于两者各自及相互之间仍然存在比对匹配关系，完全区分两者之间的基因数量仍需要进一步分析。随机对其中一条拼接序列与所有拼接序列和Trinity unigene序列进行本地blast获得的匹配序列进行正负链多序列Clustal-Omiga比对，结果见图7，高匹配部位(相似性超过90％)用红色实线表示。结果表明，两种拼接方法得到的拼接长序列仍然可以从拼接序列中找到匹配较长序列，但仍然存在多个不同长度即匹配不同部位的拼接序列，个别匹配序列多序列比对出现了匹配不佳的区段，需要进一步利用多序列比对工具Clustal-Omiga结合文本匹配打分进行合并匹配序列获得较长的唯一序列。
实施例2拼接序列比对及注释
将实施例1得到的拼接序列进行远程NCBIblast比对(https://blast.ncbi.nlm.nih.gov/Blast.cgi？PROGRAM＝blastn&PAGE_TYPE＝BlastSearch&LINK_LOC＝blasthome)，每次250条序列手动录入，多窗口同时比对，10万条1天即可完成比对。blast选择nt/nr核酸库，比对参数仅选择高度相似性即“megablast”比对，选取打分最高的第一个获取其注释信息。为了比较用Trinity得到的拼接序列与上面得到的拼接序列的拼接效果，两者均进行远程及本地blast比对分析，本地blast程序为NCBI的blast-2.2.26-FreeBSD。多序列比对工具为Clustal-Omiga(远程比对登录https://www.ebi.ac.uk/Tools/msa/clustalo/)，本地安装Clustal-Omega-1.2.2到FreeBSD服务器。
基于上述比对及注释发现，10万条拼接序列(包括没有获得延伸的种子序列)进行NCBI远程blast比对，获得注释的有74898条，没有注释的25497条，对返回注释结果选择第一个匹配基因的注释进行提取，去除重复后得到22733条不同基因注释(一些序列出现相同注释多次，与它们之间存在blast比对匹配关系有关)。
Trinity注释结果得到接近5万条注释基因(测序公司分析注释)，但是通过unigene序列的blast比对发现不同unigene之间有匹配关系(60％以上匹配了两条或以上)。为了更好比较两者的注释，利用公司提供的unigene序列(94361条)重新进行了NCBI远程blast比对注释，40996条获得注释(去除重复后得到30821条单一注释)，53365条没有获得注释(高度相似性megablast)，而很多序列在公司提供的注释是被注释了的，个别序列进一步采取“Somewhat similar sequences(blastn)”进行比对可以返回注释，从比对结果看只是部分区段获得了比对，即比对序列注释连续性有待考证。
如果拼接序列连续性较好，拼接序列越长进行NCBI远程BLAST比对获得注释机会会越高。比较了本发明获得的一个样品所有拼接序列与Trinity利用所有5个样品unigene序列注释与序列长度，见图8。由于Trinity的unigene基因长度呈现陡降态势即随着长度增加数量明显下降(图5)，注释数量也是呈现随长度快速下降(最后翘起部分为合并长序列的注释数)。而本发明利用一个样品获得的拼接序列其注释呈现居中分布，与序列拼接长度分布(图3)趋势一致，亦即更长序列获得注释是越多的。
基于上述实施例1、2可以看出：
(一)RNA-Seq测序序列两侧40mer首尾拼接延伸效率
本实施例利用测序序列两侧40mer区段尝试将其延伸拼接，更长的两端匹配(覆盖40mer)可增加拼接准确率，提高拼接长度可塑性，为了避免测序错误率带来的首尾不衔接，从两侧40mer进行隔开5mer切下5个20mer用于从测序序列中找到匹配序列进行延伸拼接，找到的匹配序列至少是由两个或以上20mer(可以覆盖25～40mer)查找返回即同一条序列被匹配两次以上，这样就可以减少由于测序错误产生错配而导致延伸中断(图2，仅仅首尾20mer延伸拼接效果不理想，40mer则效果很佳)。算法实现仅由开放源代码系统的文本匹配过滤等脚本命令即可完成(图1流程图)。本方法算法思路朴素、程序运行高效，序列拼接长度达到较满意，在增加测序目标序列数量达到300万条，用10万条种子序列进行匹配延伸，仅5％序列没有获得延伸(图3、图4)，拼接长度呈现居中分布，600-1000mer占了近40％。进行远程NCBIblast进行高相似性“megablast”比对，3/4拼接序列获得了注释(图8)，不同注释基因数量达到了2.2万条。如果利用多个10万条种子序列进行多次拼接，无疑能够更好的获得转录组测序的完整全长转录谱序列，也可以利用不同组织的测序数据进行可变剪接等分析。
本方法拼接可以追溯每一步拼接步骤的实现过程，如可以随时回找某一拼接序列的拼接来源序列，将匹配的源序列进行多序列比对就能够判断测序错误出现位置，这样便于开展定量PCR的引物设计以进一步验证基因表达情况。
(二)与Trinity拼接unigene序列的比较
Trinity从头拼接是目前公司测序使用最多的程序之一，本发明利用其获得的序Trinity本方法列拼接数据进行了一些比较分析，拼接长度方面本方法获得的长度分布居中且较长序列数量多(图3、图5)。用相同blast参数进行远程NCBI blast比对，本方法返回的注释数量也较多且分布与拼接长度居中呈现一致性(图8、图3)。与Trinity拼接相比，出现了更多的拼接序列之间相互匹配重叠关系(图6)，本方法拼接序列从Trinity拼接序列中找到的匹配unigene数量达到3.6万条单一序列，由于Trinity序列本身仍然存在相互之间的匹配重叠关系，还不能够区别两者的拼接序列的更多异质性，从图7多序列比对(Clustal-Omiga图)可以看出，不同的unigene序列可匹配同一条拼接序列的不同位置，反之亦然
(三)进一步将重叠匹配拼接序列合并延伸的必要性
从上面的分析可以看出，Trinity拼接的unigene与本方法获得的拼接序列可以较好相互匹配，Clustal-Omiga多序列比对可以明显看出匹配序列之间的重叠关系(图7)，而获得唯一序列对于了解组织基因转录水平等是十分必要的。虽然可以通过人工肉眼辨认多序列比对结果并确认唯一序列，但工作量将不堪重负。利用Clustal-Omiga比对结果，下一步将尝试结合文本匹配合并等方法合并重叠匹配拼接序列。
小结而言：
1、本发明利用开放源代码系统shell文本过滤命令，建立了高通量RNA-Seq测序序列拼接新方法并获得了较好从头拼接效果；
2、本发明的算法思路及实现简单有效：利用测序序列首尾40mer每隔开5mer切下的5个20mer查找100万条或更多的目标测序序列，将获得的每条多余两条及以上匹配序列用5个20mer进行首尾连接，保留最长的拼接序列；
3、本发明得到的拼接序列长度呈现居中分布，长度600-1000mer占了近40％，75％序列可以进行远程NCBI核酸数据库blast获得注释(megablast高相似性比对)，86％拼接序列本地blast比对可以匹配1条以上序列；
4、拼接过程完全可追溯原始拼接序列，可用于进行PCR引物设计、可变剪接分析、唯一序列多序列合并分析等。
以上所述是发明的优选实施方式，应当指出，对于本技术领域的普通技术人员来说，在不脱离本发明原理的前提下，还可以做出若干改进和润饰，这些改进和润饰也视为本发明的保护范围。

Claims (10)

1.一种基于文本匹配延伸高通量测序基因的拼接方法，其特征在于，包括：
S1：通过高通量测序平台获取待分析样品的测序序列；
S2：将测序序列加编号，打散，随机组合；
S3：选取10万条测序序列作为种子序列；
S4：将种子序列第一侧预设长度的序列隔开，得到查询序列；
S5：将查询序列与一组或多组100万条测序序列进行文本匹配，并保留至少匹配两次的多条序列作为匹配序列；
S6：采用查询序列作为拼接点，将匹配序列与种子序列第一侧拼接，并选取拼接后长度最长的序列作为拼接序列；
S7：比较拼接序列第一侧预设长度的序列与种子序列第一侧预设长度的序列；
S8：若两者不同，则以拼接序列作为新的种子序列，并进入步骤S4；直至所得拼接序列第一预设长度的序列与种子序列第一侧预设长度的序列相同，或进行至少5次拼接后进入S9；
S9：若两者相同，则将所得拼接序列反转后作为新的种子序列，并进入步骤S4，直至完成种子序列两侧的拼接。
2.如权利要求1所述的基于文本匹配延伸高通量测序基因的拼接方法，其特征在于，步骤S4中，将每条种子序列左侧40mer每隔5mer切下5个20mer序列，得到查询序列。
3.如权利要求1所述的基于文本匹配延伸高通量测序基因的拼接方法，其特征在于，步骤S5中，将查询序列与一组或多组100万条测序序列进行文本匹配，得到返回序列；每个查询序列得到的返回序列进行排序并计数重复数，保留超过2条相同序列以上的序列组作为匹配序列。
4.如权利要求1所述的基于文本匹配延伸高通量测序基因的拼接方法，其特征在于，步骤S2中，采用逐级按照每10万→100万→5万条进行随机排序并按随机方式合并序列，其中，每个100万分成若干个目录单独进行每1万条的随机排序，再随机组合合并所有序列。
5.如权利要求1所述的基于文本匹配延伸高通量测序基因的拼接方法，其特征在于，步骤S1中，所述高通量测序平台选用Roche的454测序平台、Illumina的HySeq4000测序平台和ABI的SOLiD测序平台中的一种或多种。
6.如权利要求1所述的基于文本匹配延伸高通量测序基因的拼接方法，其特征在于，所述待分析样品为动物、植物、微生物的DNA或RNA提取物；或者，所述待分析样本为水、大气、土壤中的微生物的DNA或RNA提取物。
7.如权利要求1所述的基于文本匹配延伸高通量测序基因的拼接方法，其特征在于，步骤S5中，采用开源代码系统内嵌shell命令进行文本过滤；
所述shell命令选用awk、sed、sort、grep、tr、split、comm、paste、cat中的一种或多种。
8.一种基于文本匹配延伸高通量测序基因的拼接系统，其特征在于，包括：
高通量测序平台，用于获取待分析样品的测序序列；
数据处理系统：所述的数据处理系统用于实现权利要求1-7任一项所述的拼接方法。
9.如权利要求8所述的基于文本匹配延伸高通量测序基因的拼接系统，其特征在于，所述数据处理系统为开源代码系统。
10.一种计算机可读介质，其上存储有计算机指令，其特征在于，该指令被处理执行时实现权利要求1-7任一项所述的拼接方法的步骤。
CN202210856831.3A 2022-07-20 2022-07-20 基于文本匹配延伸高通量测序基因的拼接方法 Active CN115331736B (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN202210856831.3A CN115331736B (zh) 2022-07-20 2022-07-20 基于文本匹配延伸高通量测序基因的拼接方法

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN202210856831.3A CN115331736B (zh) 2022-07-20 2022-07-20 基于文本匹配延伸高通量测序基因的拼接方法

Publications (2)

Publication Number Publication Date
CN115331736A CN115331736A (zh) 2022-11-11
CN115331736B true CN115331736B (zh) 2023-07-25

Family

ID=83917711

Family Applications (1)

Application Number Title Priority Date Filing Date
CN202210856831.3A Active CN115331736B (zh) 2022-07-20 2022-07-20 基于文本匹配延伸高通量测序基因的拼接方法

Country Status (1)

Country Link
CN (1) CN115331736B (zh)

Citations (5)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN105389481A (zh) * 2015-12-22 2016-03-09 武汉菲沙基因信息有限公司 一种三代全长转录组中可变剪切体的检测方法
CN108388772A (zh) * 2018-01-26 2018-08-10 佛山科学技术学院 一种利用文本比对分析高通量测序基因表达水平的方法
CN110257481A (zh) * 2019-05-09 2019-09-20 扬州大学 一种基于比较基因组学的转座子插入多态tip分子标记的挖掘方法
WO2021226558A1 (en) * 2020-05-08 2021-11-11 The Broad Institute, Inc. Methods and compositions for simultaneous editing of both strands of a target double-stranded nucleotide sequence
CN113891936A (zh) * 2019-03-19 2022-01-04 布罗德研究所股份有限公司 编辑核苷酸序列的方法和组合物

Patent Citations (7)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN105389481A (zh) * 2015-12-22 2016-03-09 武汉菲沙基因信息有限公司 一种三代全长转录组中可变剪切体的检测方法
CN108388772A (zh) * 2018-01-26 2018-08-10 佛山科学技术学院 一种利用文本比对分析高通量测序基因表达水平的方法
CN113891936A (zh) * 2019-03-19 2022-01-04 布罗德研究所股份有限公司 编辑核苷酸序列的方法和组合物
CN113891937A (zh) * 2019-03-19 2022-01-04 布罗德研究所股份有限公司 编辑核苷酸序列的方法和组合物
CN114127285A (zh) * 2019-03-19 2022-03-01 布罗德研究所股份有限公司 编辑核苷酸序列的方法和组合物
CN110257481A (zh) * 2019-05-09 2019-09-20 扬州大学 一种基于比较基因组学的转座子插入多态tip分子标记的挖掘方法
WO2021226558A1 (en) * 2020-05-08 2021-11-11 The Broad Institute, Inc. Methods and compositions for simultaneous editing of both strands of a target double-stranded nucleotide sequence

Non-Patent Citations (5)

* Cited by examiner, † Cited by third party
Title
Comparing DNA sequence collections by direct comparison of compressed text indexes;Cox A J等;《Algorithms in Bioinformatics: 12th International Workshop》;214-224 *
Perl语言环境下生物信息学的数据库技术;郭文久等;《安康学院学报》(第5期);80-84 *
Unix文本比对分析高通量RNA-Seq测序基因表达;宋东光等;《生物信息学》;第16卷(第2期);119-129 *
基于短序列分组和拼接策略的子序列快速查询算法;范纯龙等;《计算机应用研究》;第37卷(第6期);1702-1706 *
真实感汉语可视语音合成关键技术研究;赵晖;《中国博士学位论文全文数据库 (信息科技辑)》(第4期);I136-21 *

Also Published As

Publication number Publication date
CN115331736A (zh) 2022-11-11

Similar Documents

Publication Publication Date Title
US9334532B2 (en) Complexity reduction method
Kiethega et al. RNA-level unscrambling of fragmented genes in Diplonema mitochondria
US20210375397A1 (en) Methods and systems for determining fusion events
CN110692101A (zh) 用于比对靶向的核酸测序数据的方法
CN109337997B (zh) 一种山茶属多态性叶绿体基因组微卫星分子标记引物及筛选和甄别近缘种的方法
CN102899335A (zh) 一种高通量Small RNA测序获得番木瓜环斑病毒基因组序列的方法
Monger et al. Towards next generation CHO cell biology: Bioinformatics methods for RNA‐Seq‐based expression profiling
CN114708910A (zh) 一种利用单细胞测序数据计算池测序中细胞亚群富集分数的方法
AU2010329825B2 (en) RNA analytics method
CN111192636A (zh) 一种适用于oligodT富集的mRNA二代测序结果分析方法
CN108388772B (zh) 一种利用文本比对分析高通量测序基因表达水平的方法
CN107862177B (zh) 一种区分鲤群体的单核苷酸多态性分子标记集的构建方法
Pereira et al. RNA‐seq: applications and best practices
CN115331736B (zh) 基于文本匹配延伸高通量测序基因的拼接方法
US20240141425A1 (en) Correcting for deamination-induced sequence errors
US20200232010A1 (en) Methods, compositions, and systems for improving recovery of nucleic acid molecules
US20140364321A1 (en) Method for analyzing DNA methylation based on MspJI cleavage
US20200131566A1 (en) Methods, compositions and systems for calibrating epigenetic partitioning assays
WO2021108708A1 (en) Methods, compositions and systems for improving the binding of methylated polynucleotides
CN104951673A (zh) 一种基因组酶切图谱拼接方法及系统
CN106520961B (zh) 玉米微卫星标记位点开发方法与微卫星标记位点内的微卫星标记的长度检测方法
CN110684830A (zh) 一种石蜡切片组织rna分析方法
WO2023192568A1 (en) Methods and systems for detecting ribonucleic acids
JP2023534124A (ja) 遺伝子シーケンシング解析方法、装置、記憶媒体及びコンピュータ機器
II et al. Targeted capture of phylogenetically-informative Ves SINE insertions in genus Myotis

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