CN114582419B - 一种基于滑动窗口的基因序列多聚腺苷酸尾巴提取方法 - Google Patents

一种基于滑动窗口的基因序列多聚腺苷酸尾巴提取方法 Download PDF

Info

Publication number
CN114582419B
CN114582419B CN202210110546.7A CN202210110546A CN114582419B CN 114582419 B CN114582419 B CN 114582419B CN 202210110546 A CN202210110546 A CN 202210110546A CN 114582419 B CN114582419 B CN 114582419B
Authority
CN
China
Prior art keywords
tail
sequence
poly
base
value
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
CN202210110546.7A
Other languages
English (en)
Other versions
CN114582419A (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.)
Suzhou University
Original Assignee
Suzhou 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 Suzhou University filed Critical Suzhou University
Priority to CN202210110546.7A priority Critical patent/CN114582419B/zh
Publication of CN114582419A publication Critical patent/CN114582419A/zh
Application granted granted Critical
Publication of CN114582419B publication Critical patent/CN114582419B/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
    • G16B20/00—ICT specially adapted for functional genomics or proteomics, e.g. genotype-phenotype associations
    • G16B20/30—Detection of binding sites or motifs

Landscapes

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

Abstract

本发明公开了一种基于滑动窗口的基因序列多聚腺苷酸尾巴提取方法，特点是包括输入基因序列文件查找连续N个A碱基片段和连续N个T碱基片段作为潜在的多聚腺苷酸尾巴的起始位置的步骤；初始化设置滑动窗口参数并移动滑动窗口直到滑动窗口到达序列最右或左端，或者错配惩罚值达到阈值5则停止滑动，获得多聚腺苷酸尾巴序列和尾巴长度值的步骤；若提供该物种的参考基因组序列，将待比对序列通过序列比对工具与参考基因组序列进行比对或者使用过滤条件过滤多聚腺苷酸尾巴，获得精确度更高的多聚腺苷酸尾巴，并根据多聚腺苷酸尾巴数量确定序列类型的步骤；最后确定多聚腺苷酸尾巴为类型的步骤，优点是精度高、运算速度快且用户友好。

Description

一种基于滑动窗口的基因序列多聚腺苷酸尾巴提取方法
技术领域
本发明涉及一种多聚腺苷酸尾巴提取方法，尤其是涉及一种基于滑动窗口的基因序列多聚腺苷酸尾巴提取方法。
背景技术
随着测序技术的不断发展，第三代测序技术即单分子实时测序技术逐渐成熟，该技术同时可以测量上亿序列模板，1秒可以测量十个碱基且测序准确率高达99.9％，这种高通量、高速度、高精度的全长DNA测序技术迅速产生了大量基因序列数据。这些基因序列数据，包含着丰富的信息，可用于研究选择性多聚腺苷酸化(Alternative Polyadenylation,APA)，特别是是多聚腺苷酸尾巴(poly(A)tail)长度。poly(A)尾巴是位于mRNA末端的一连串主要由腺苷酸A组成的序列，但是高通量测序技术测得的基因序列中poly(A)尾巴可能处于序列中部，且可能因为测序错误含有非A碱基，需要开发算法进行提取。目前大量学者已证明poly(A)尾巴在许多生物过程中起着重要作用，且其长度和翻译效率、稳定性关系密切。但目前仍缺少高效准确且灵活易用的算法工具能从序列中识别并提取poly(A)尾巴。
现有的能从序列中识别并提取poly(A)尾巴提取方法，包括
(1)PAIso-seq(Liu,Y.,et al.Poly(A)inclusive RNA isoform sequencing(PAIso-seq)reveals wide-spread non-adenosine residues within RNA poly(A)tails.Nature Communications 2019；10(1):5292)提供的流程安装依赖较多，仅适用于PAIso-seq测序产生的三代全长测序数据，不适用于基于短测序策略产生的数据；PAIso-seq采用经验条件过滤比对后被剪切掉的序列来获取尾巴的具体碱基构成和尾巴长度。具体来说，就是将原始测序数据比对到参考基因组，未能成功比对到参考基因组上的序列被认为是潜在的尾巴，再进一步根据以下条件过滤：1)序列长度不小于15nt；2)序列中至少包含5个连续A碱基；3)序列中的非A碱基个数小于20个；4)序列中的非A碱基比例少于50％。满足以上四个条件便被PAIso-seq认定为poly(A)尾巴。该方法仅能提取满足固定过滤条件的尾巴，仅适用于全长测序数据，且需要通过序列比对才能获得尾巴，不适用于无参考基因组的物种，且一条序列只能找到一条尾巴。
(2)FLAMseq(Legnini,I.,et al.FLAM-seq:full-length mRNA sequencingreveals principles of poly(A)tail length control.Nat.Methods 2019；16(9):879-886.)开发的寻找尾巴的工具FLAMAnalysis采用投票的方式确定尾巴长度：首先使用人工标记序列(FLAMseq测序结果中会包含的一段已知的序列)确定尾巴的起始点，再从该起始点往后取固定长度的序列，判断序列中的尾巴非A碱基比例是否满足要求，满足则继续向后搜索直到不满足条件。该方法在确定尾巴时，采用多组参数，比如将每次扫描的序列长度设为L，A碱基的比例设为N，则通过多组L-N组合来进行结果投票，票数高的尾巴结果定为最终结果。该方法使用多轮投票方式，投票算法涉及到的参数多且不直观，导致其计算速度慢，难以确定参数值。此外，该方法仅适用于FLAMseq产生的带有固定结构的测序序列，不适用于其它测序方法产生的数据。
(3)专门针对Nanopore开发的尾巴定量工具tailfindr(Krause,M.,etal.tailfindr:alignment-free poly(A)length measurement for Oxford Nanopore RNAand DNA sequencing.2019；25(10):1229-1241)仅能从Nanopore测序数据(即电平信号)中近似定量尾巴长度而不能准确定量尾巴的碱基组成。由于Nanopore测序方法采用不同碱基导电率的差异进行基因测序，ACTG四种碱基作为电阻在相同电压下产生的电流大小不同，即当不同碱基穿过电压时会产生不同强度的电信号，通过有监督的学***信号产生的速度很高，因此当大量相同碱基出现之后(比如poly(A)尾巴)，很难准确预估序列真实长度，且完全无法辨识其中少量混杂的其它碱基产生的电信号。
(4)FLEP-seq(Long,Y.,et al.FLEP-seq:simultaneous detection of RNApolymerase II position,splicing status,polyadenylation site and poly(A)taillength at genome-wide scale by single-molecule nascent RNA sequencing.NatureProtocols 2021；16(9):4355-4381.)针对PacBio(全长测序)和Nanopore两种不同的全长测序方法开发了两个不同的脚本来定量poly(A)尾巴。FLEP-seq中针对PacBio数据采用积分方式寻找尾巴，即从序列起点开始，如果是A碱基则积分+1，否则积分-1.5，到某个碱基时积分为正则认为这里是尾巴起点继续向后扫描，如果为负则等积分重新变为正开始认为是尾巴起点，直到扫完为止，记录下积分最大的点为尾巴结束点，并截取尾巴。该方法只能得到序列的尾巴长度信息，并不能得到尾巴的具体序列或位置等信息。
(5)Poly(A)-seq(Yu,F.,et al.Poly(A)-seq:A method for direct sequencingand analysis of the transcriptomic poly(A)-tails.PLOS ONE 2020；15(6):e0234696)中提到的名为PA-finder的流程或PAIso-seq2均未提供任何可用的代码。Poly(A)-seq采用复杂的方式进行尾巴定位：首先将测序过程中添加的人工接头去掉；然后通过A*9模式匹配尾巴的起点，允许最高一个错配，并从起点处切割序列；找到满足条件的序列之后通过两个经验条件：1)长度至少10nt，2)非A碱基含量小于等于4nt，过滤出尾巴。该方法只能找到至少包含9个连续A碱基的尾巴，且仅能提取满足固定过滤条件的尾巴，且一条序列只能找到一条尾巴。此外，该方法的流程仅在论文中简单描述，作者并未提供代码或工具包。
发明内容
本发明所要解决的技术问题是提供一种精度高、运算速度快且用户友好的基于滑动窗口的基因序列多聚腺苷酸尾巴提取方法。
本发明解决上述技术问题所采用的技术方案为：一种基于滑动窗口的基因序列多聚腺苷酸尾巴提取方法，包括以下步骤：
(1)在输入FastQ格式的基因序列文件中查找连续N个A碱基片段和连续N个T碱基片段作为潜在的多聚腺苷酸尾巴的起始位置；
(2)初始化设置滑动窗口参数；
(3)若潜在的多聚腺苷酸尾巴为A碱基片段，则以片段末端为起始位置，移动滑动窗口直到滑动窗口到达序列最右端，或者错配惩罚值达到阈值则停止滑动，获得多聚腺苷酸尾巴序列和尾巴长度值；
(4)若潜在的多聚腺苷酸尾巴为T碱基片段，则以片段首位碱基为起始位置，移动滑动窗口直到滑动窗口到达序列最左端，或者错配惩罚值达到阈值则停止滑动，获得多聚腺苷酸尾巴序列和尾巴长度值；
(5)将步骤(3)及步骤(4)获得的多聚腺苷酸尾巴过滤后，获得精确度更高的多聚腺苷酸尾巴，并根据序列中提取的多聚腺苷酸尾巴数量确定步骤(1)输入的基因序列类型；
(6)若提供接头序列，则在步骤(5)获得的多聚腺苷酸尾巴末尾位置的3’端查找接头序列，若找到的接头序列的起始位置位于多聚腺苷酸尾巴末尾位置的3’端的10个碱基范围内，则将该多聚腺苷酸尾巴类型标记为“structural”类型，否则将其标记为“nonstructural”类型；
(7)若未提供接头序列，则计算步骤(5)多聚腺苷酸尾巴末尾位置到序列3’端终点的碱基数，如果超过25个碱基，则标记多聚腺苷酸尾巴类型为“nonstructural”类型，否则标记为“structural”类型。
进一步，步骤(1)具体为：输入FastQ格式的基因序列文件，在输入的基因序列中查找连续N个A碱基片段和连续N个T碱基片段，将查找到的A碱基片段和T碱基片段作为潜在的多聚腺苷酸尾巴的起始位置，若都没找到，则认为该基因序列中没有多聚腺苷酸尾巴，其中N为初始尾巴长度，默认值为8。
进一步，步骤(2)具体为：将步骤(1)获得的每条潜在的多聚腺苷酸尾巴设置滑动窗口，滑动窗口的初始大小为该条潜在多聚腺苷酸尾巴的长度，初始位置为该条潜在多聚腺苷酸尾巴起点位置，滑动距离固定为1，滑动窗口滑动方向为序列的5’端到3’端，即若是A碱基片段尾巴，则滑动方向为从左到右，若是T碱基片段尾巴，则滑动方向为从右到左。
进一步，步骤(3)具体为：若潜在的多聚腺苷酸尾巴为A碱基片段，则以片段末位碱基为起始位置，将初始化后的滑动窗口向右移动1个碱基，如果该碱基为A且错配惩罚值为0，则尾巴长度计数值加1，错配惩罚值不变；如果该碱基为A且错配惩罚值不为0，则尾巴长度计数值加1，错配惩罚值减1，如果此时错配惩罚值为负，则将其值重新置0；如果该碱基为非A，则长度计数值和错配惩罚值都加1，重复该过程，直到滑动窗口到达序列最右端，或者错配惩罚值达到阈值则停止滑动，获得多聚腺苷酸尾巴序列和尾巴长度值，多聚腺苷酸尾巴序列为尾巴最左端到错配惩罚值最后一次增加1的左边一个碱基之间的这段序列。
所述的多聚腺苷酸尾巴长度计数值初始值为初始尾巴长度N，N值为8，错配惩罚值表示A碱基片段尾巴中非A碱基个数，初始值为0；错配惩罚阈值表示A碱基片段尾巴中出现的最长连续非A碱基个数，默认值为5。
进一步，步骤(4)具体为：若潜在的多聚腺苷酸尾巴为T碱基片段，则以该片段首位碱基为起始位置，将初始化后的滑动窗口向左移动1个碱基，如果该碱基为T且错配惩罚值为0，则尾巴长度计数值加1，错配惩罚值不变；如果该碱基为T且错配惩罚值不为0，则尾巴长度计数值加1，错配惩罚值减1，如果此时错配惩罚值为负，则将其值重新置0；如果该碱基为非T，则长度计数值和错配惩罚值都加1，重复该过程，直到滑动窗口到达序列最左端，或者错配惩罚值达到阈值则停止滑动，获得多聚腺苷酸尾巴序列和尾巴长度值，多聚腺苷酸尾巴序列为尾巴最右端到错配惩罚值最后一次增加1的右边一个碱基之间的这段序列。
所述的多聚腺苷酸尾巴长度计数值初始值为初始尾巴长度N，N值为8，错配惩罚值表示T碱基片段尾巴中非T碱基个数，初始值为0；错配惩罚阈值表示T碱基片段尾巴中出现的最长连续非T碱基个数，默认值为5。
进一步，步骤(5)中若提供该物种的参考基因组序列，则将步骤(3)和步骤(4)获得的每条多聚腺苷酸尾巴序列及其5’端外侧的200个碱基序列作为待比对序列，若该条多聚腺苷酸尾巴序列5’端外侧的序列长度不足200个碱基，则将该条多聚腺苷酸尾巴序列及其5’端外侧的全部序列作为待比对序列，将待比对序列通过序列比对工具与参考基因组序列进行比对，将待比对序列中的多聚腺苷酸尾巴序列中不能比对上参考基因组序列的片段提取过滤后，获得精确度更高的多聚腺苷酸尾巴，并根据序列中提取的多聚腺苷酸尾巴数量确定步骤(1)输入的基因序列类型；若多聚腺苷酸尾巴序列与参考基因组序列碱基无差异，则表明该尾巴序列实际是来自参考基因组的序列，并非真正的多聚腺苷酸尾巴。
进一步，步骤(5)中所述的过滤具体如下：
A.若多聚腺苷酸尾巴长度少于12个碱基长度，则舍弃；
B.若一条物种全长基因序列识别到多个多聚腺苷酸尾巴，则将长度低于所有多聚腺苷酸尾巴平均长度值的三分之二长度的尾巴舍弃；
C.将A碱基尾巴中非A碱基含量高于90％的多聚腺苷酸尾巴和T碱基尾巴中非T碱基含量高于90％的多聚腺苷酸尾巴舍弃。
所述的序列比对工具为minimap2、Bowtie或STAR软件。
与现有技术相比，本发明的优点在于：本发明提出一种基于滑动窗口的基因序列多聚腺苷酸尾巴提取方法，首先扫描序列中的多聚A或T碱基片段得到潜在的多聚腺苷酸尾巴位置，再对每个位置通过滑动窗口算法扫描完整尾巴，并通过严格的条件进行过滤尾巴，最终获得精确的多聚腺苷酸尾巴位置及长度。具有以下优点：1)因其使用多聚A或T的识别及滑动窗口法，运行速度快、精度高；2)同时适用于短读段或全长测序数据，且对无参考基因组的物种也适用。
综上所述，本发明通过多聚A或T的识别及滑动窗口进行多聚腺苷酸尾巴定位和识别，获得精确的尾巴位置及长度。这种方法可以直接在原始序列准确提取尾巴，无需进行序列比对，也可在提供参考基因组序列的情况下进一步修正尾巴，适用于高通量短读段测序数据或全长测序数据。
具体实施方式
以下结合实施例对本发明作进一步详细描述。
在下文中，仅简单地描述了某些示例性实施例。正如本领域技术人员可认识到的那样，在不脱离本发明申请实施例的精神或范围的情况下，可通过各种不同方式修改所描述的实施例。因此，描述被认为本质上是示例性的而非限制性的。
一、具体实施例
由于基因序列是双链的，在实际数据分析过程中，由于未知序列的具体的方向，在计算时，需要同时考虑+向和-向的情况。即，针对尾巴提取问题，将同时考虑A碱基的尾巴以及T碱基的尾巴。
参数初始化：
接头序列(adapterSeq)：接头(adapter)是测序仪在测序过程中添加的一段已知的短核苷酸序列，有些测序数据的接头序列可能已被去除(例如，测序公司会对测序仪产出的数据进行预处理，去除接头后再返回给客户)，有些则仍被保留。例如，FLAM-seq测序会在尾巴结束位置***一段特定序列(9个G碱基)以标记尾巴的结束位置。通过参数adapterSeq传入这里可以根据用户数据的实际情况，指定接头序列。若不存在则可以不初始化。adapterSeq默认值为“”。
初始尾巴长度(anchorLen)：用于确定尾巴在序列中的起始位置，在序列中首先找到anchorLen个连续的A碱基或T碱基作为初始尾巴。默认值为8。
非A碱基个数上限(dropoff)：指定尾巴中可能出现的最长连续非A碱基数(若是A碱基的尾巴)或非T碱基数(若是T碱基的尾巴)，默认值为5。
物种的参考基因组序列：指定要扫描尾巴的序列数据对应的物种的参考基因组序列，序列文件格式为FASTA。默认为“”，即表示不提供。
一种基于滑动窗口的基因序列多聚腺苷酸尾巴提取方法，包括以下步骤：
1、在输入FastQ格式的基因序列文件中查找连续N个A碱基片段和连续N个T碱基片段作为潜在的多聚腺苷酸尾巴的起始位置，具体如下：输入FastQ格式的基因序列文件，在输入的基因序列中查找连续N个A碱基片段和连续N个T碱基片段，将查找到的A碱基片段和T碱基片段作为潜在的多聚腺苷酸尾巴的起始位置，若都没找到，则认为该基因序列中没有多聚腺苷酸尾巴，其中N为初始尾巴长度，默认值为8。
2、初始化设置滑动窗口参数，具体如下：将步骤1获得的每条潜在的多聚腺苷酸尾巴设置滑动窗口，滑动窗口的初始大小为该条潜在多聚腺苷酸尾巴的长度，初始位置为该条潜在多聚腺苷酸尾巴起点位置，滑动距离固定为1，滑动窗口滑动方向为序列的5’端到3’端，即若是A碱基片段尾巴(序列方向5’端到3’端)，则滑动方向为从左到右，若是T碱基片段尾巴(序列方向3’端到5’端)，则滑动方向为从右到左。
3、若潜在的多聚腺苷酸尾巴为A碱基片段，则以该片段末位碱基为起始位置，移动滑动窗口直到滑动窗口到达序列最右端，或者错配惩罚值达到阈值则停止滑动，获得多聚腺苷酸尾巴序列和尾巴长度值，具体如下：若潜在的多聚腺苷酸尾巴为A碱基片段，则以该片段末位碱基为起始位置，将初始化后的滑动窗口向右移动1个碱基，如果该碱基为A且错配惩罚值为0，则尾巴长度计数值加1，错配惩罚值不变；如果该碱基为A且错配惩罚值不为0，则尾巴长度计数值加1，错配惩罚值减1，如果此时错配惩罚值为负，则将其值重新置0；如果该碱基为非A，则长度计数值和错配惩罚值都加1，重复该过程，直到滑动窗口到达序列最右端，或者错配惩罚值达到阈值则停止滑动，获得多聚腺苷酸尾巴序列和尾巴长度值，多聚腺苷酸尾巴序列为尾巴最左端到错配惩罚值最后一次增加1的左边一个碱基之间的这段序列。其中多聚腺苷酸尾巴长度计数值初始值为初始尾巴长度N，N值为8，错配惩罚值表示A碱基片段尾巴中非A碱基个数，初始值为0；错配惩罚阈值表示A碱基片段尾巴中出现的最长连续非A碱基个数，默认值为5。
4、若潜在的多聚腺苷酸尾巴为T碱基片段，则以该片段首位碱基为起始位置，移动滑动窗口直到滑动窗口到达序列最左端，或者错配惩罚值达到阈值则停止滑动，获得多聚腺苷酸尾巴序列和尾巴长度值，具体如下：若潜在的多聚腺苷酸尾巴为T碱基片段，则以该片段首位碱基为起始位置，将初始化后的滑动窗口向左移动1个碱基，如果该碱基为T且错配惩罚值为0，则尾巴长度计数值加1，错配惩罚值不变；如果该碱基为T且错配惩罚值不为0，则尾巴长度计数值加1，错配惩罚值减1，如果此时错配惩罚值为负，则将其值重新置0；如果该碱基为非T，则长度计数值和错配惩罚值都加1，重复该过程，直到滑动窗口到达序列最左端，或者错配惩罚值达到阈值则停止滑动，获得多聚腺苷酸尾巴序列和尾巴长度值，多聚腺苷酸尾巴序列为尾巴最右端到错配惩罚值最后一次增加1的右边一个碱基之间的这段序列，其中多聚腺苷酸尾巴长度计数值初始值为初始尾巴长度N，N值为8，错配惩罚值表示T碱基片段尾巴中非T碱基个数，初始值为0；错配惩罚阈值表示T碱基片段尾巴中出现的最长连续非T碱基个数，默认值为5。
5、将步骤3及步骤4获得的多聚腺苷酸尾巴过滤后，获得精确度更高的多聚腺苷酸尾巴，并根据序列中提取的多聚腺苷酸尾巴数量确定步骤1输入的基因序列类型，其中过滤具体如下：
A.若多聚腺苷酸尾巴长度少于12个碱基长度，则舍弃；
B.若一条物种全长基因序列识别到多个多聚腺苷酸尾巴，则将长度低于所有多聚腺苷酸尾巴平均长度值的三分之二长度的尾巴舍弃；
C.将A碱基尾巴中非A碱基含量高于90％的多聚腺苷酸尾巴和T碱基尾巴中非T碱基含量高于90％的多聚腺苷酸尾巴舍弃。
若提供该物种的参考基因组序列，则将步骤3和步骤4获得的每条多聚腺苷酸尾巴序列及其5’端外侧的200个碱基序列作为待比对序列，若该条多聚腺苷酸尾巴序列5’端外侧的序列长度不足200个碱基，则将该条多聚腺苷酸尾巴序列及其5’端外侧的全部序列作为待比对序列，将待比对序列通过序列比对工具与参考基因组序列进行比对，将待比对序列中的多聚腺苷酸尾巴序列中不能比对上参考基因组序列的片段提取过滤后，获得精确度更高的多聚腺苷酸尾巴，并根据序列中提取的多聚腺苷酸尾巴数量确定步骤(1)输入的基因序列类型；若多聚腺苷酸尾巴序列与参考基因组序列碱基无差异，则表明该尾巴序列实际是来自参考基因组的序列，并非真正的多聚腺苷酸尾巴。序列比对是绝大多数生物数据分析的例行步骤，有许多被广泛接受的比对工具，如minimap2、Bowtie、STAR等，用户可以自行从网络上下载并安装一种比对工具，可使用默认参数比对。
6、若提供接头序列，则在步骤5获得的多聚腺苷酸尾巴末尾位置的3’端(如果是A碱基尾巴，则是在尾巴的右端；若是T碱基尾巴，则是在尾巴的左端)查找接头序列，若找到的接头序列的起始位置位于多聚腺苷酸尾巴末尾位置的3’端的10个碱基范围内，则将该多聚腺苷酸尾巴类型标记为“structural”类型，否则将其标记为“nonstructural”类型。
7、若未提供接头序列，则计算步骤(5)多聚腺苷酸尾巴末尾位置到序列3’端终点的碱基数，如果超过25个碱基，则标记多聚腺苷酸尾巴类型为“nonstructural”类型，否则标记为“structural”类型。尾巴一般都出现在序列的3’末端，“structural”类型表明该尾巴和测序结构相对应，更有可能是真物学上真正的尾巴。
最终输出尾巴列表，内容包括序列ID、序列方向(如果有提供参考基因组)、序列比对上的位置(如果有提供参考基因组)、尾巴完整序列、尾巴长度、非A/T碱基数量、非A/T碱基比例(非A/T碱基数量总数/尾巴总长度)、C碱基数量、C碱基比例、G碱基数量、G碱基比例、T/A碱基数量、T/A碱基比例、尾巴类型(structural或nonstructural)、序列类型“one-tail”(序列只有一条尾巴)、“two-tail”(序列有两条尾巴)、“multi-tail”(序列有超过2条尾巴)。
二、结果分析
表1依据本发明提取方法设计的软件包PolyAtailor与其它工具的运算速度比较
Figure BDA0003494956460000091
注：“无需比对”指未提供物种的参考基因组时，不通过序列比对，直接从测序序列中提取poly(A)尾巴。“基于比对”指提供物种的参考基因组时，通过序列比对，从测序序列中提取poly(A)尾巴。表格中的“——”表示相应的工具不提供该功能。
上述说明并非对本发明的限制，本发明也并不限于上述举例。本技术领域的普通技术人员在本发明的实质范围内，做出的变化、改型、添加或替换，也应属于本发明的保护范围。

Claims (7)

1.一种基于滑动窗口的基因序列多聚腺苷酸尾巴提取方法，其特征在于包括以下步骤：
（1）在输入FastQ格式的基因序列文件中查找连续N个A碱基片段和连续N个T碱基片段作为潜在的多聚腺苷酸尾巴的起始位置；
（2）初始化设置滑动窗口参数；
（3）若潜在的多聚腺苷酸尾巴为A碱基片段，则以该片段末位碱基为起始位置，移动滑动窗口直到滑动窗口到达序列最右端，或者错配惩罚值达到阈值则停止滑动，获得多聚腺苷酸尾巴序列和尾巴长度值，具体为：若潜在的多聚腺苷酸尾巴为A碱基片段，则以该片段末位碱基为起始位置，将初始化后的滑动窗口向右移动1个碱基，如果该碱基为A且错配惩罚值为0，则尾巴长度计数值加1，错配惩罚值不变；如果该碱基为A且错配惩罚值不为0，则尾巴长度计数值加1，错配惩罚值减1，如果此时错配惩罚值为负，则将其值重新置0；如果该碱基为非A，则长度计数值和错配惩罚值都加1，重复该过程，直到滑动窗口到达序列最右端，或者错配惩罚值达到阈值则停止滑动，获得多聚腺苷酸尾巴序列和尾巴长度值，多聚腺苷酸尾巴序列为尾巴最左端到错配惩罚值最后一次增加1的左边一个碱基之间的这段序列；
（4）若潜在的多聚腺苷酸尾巴为T碱基片段，则以该片段首位碱基为起始位置，移动滑动窗口直到滑动窗口到达序列最左端或者错配惩罚值达到阈值则停止滑动，获得多聚腺苷酸尾巴序列和尾巴长度值，具体为：若潜在的多聚腺苷酸尾巴为T碱基片段，则以该片段首位碱基为起始位置，将初始化后的滑动窗口向左移动1个碱基，如果该碱基为T且错配惩罚值为0，则尾巴长度计数值加1，错配惩罚值不变；如果该碱基为T且错配惩罚值不为0，则尾巴长度计数值加1，错配惩罚值减1，如果此时错配惩罚值为负，则将其值重新置0；如果该碱基为非T，则长度计数值和错配惩罚值都加1，重复该过程，直到滑动窗口到达序列最左端，或者错配惩罚值达到阈值则停止滑动，获得多聚腺苷酸尾巴序列和尾巴长度值，多聚腺苷酸尾巴序列为尾巴最右端到错配惩罚值最后一次增加1的右边一个碱基之间的这段序列；
（5）将步骤（3）及步骤（4）获得的多聚腺苷酸尾巴过滤后，获得精确度更高的多聚腺苷酸尾巴，并根据序列中提取的多聚腺苷酸尾巴数量确定步骤（1）输入的基因序列类型，所述的过滤具体如下：
A．若多聚腺苷酸尾巴长度少于12个碱基长度，则舍弃；
B．若一条物种全长基因序列识别到多个多聚腺苷酸尾巴，则将长度低于所有多聚腺苷酸尾巴平均长度值的三分之二长度的尾巴舍弃；
C．将A碱基尾巴中非A碱基含量高于90%的多聚腺苷酸尾巴和T碱基尾巴中非T碱基含量高于90%的多聚腺苷酸尾巴舍弃；
（6）若提供接头序列，则在步骤（5）获得的多聚腺苷酸尾巴末尾位置的3’端查找接头序列，若找到的接头序列的起始位置位于多聚腺苷酸尾巴末尾位置的3’端的10个碱基范围内，则将该多聚腺苷酸尾巴类型标记为“structural”类型，否则将其标记为“nonstructural”类型；
（7）若未提供接头序列，则计算步骤（5）多聚腺苷酸尾巴末尾位置到序列3’端终点的碱基数，如果超过25个碱基，则标记多聚腺苷酸尾巴类型为“nonstructural”类型，否则标记为“structural”类型。
2.根据权利要求1所述的一种基于滑动窗口的基因序列多聚腺苷酸尾巴提取方法，其特征在于步骤（1）具体为：输入FastQ格式的基因序列文件，在输入的基因序列中查找连续N个A碱基片段和连续N个T碱基片段，将查找到的A碱基片段和T碱基片段作为潜在的多聚腺苷酸尾巴的起始位置，若都没找到，则认为该基因序列中没有多聚腺苷酸尾巴，其中N为初始尾巴长度，默认值为8。
3.根据权利要求2所述的一种基于滑动窗口的基因序列多聚腺苷酸尾巴提取方法，其特征在于步骤（2）具体为：将步骤（1）获得的每条潜在的多聚腺苷酸尾巴设置滑动窗口，滑动窗口的初始大小为该条潜在多聚腺苷酸尾巴的长度，初始位置为该条潜在多聚腺苷酸尾巴起点位置，滑动距离固定为1，滑动窗口滑动方向为序列的5’端到3’端，即若是A碱基片段尾巴，则滑动方向为从左到右，若是T碱基片段尾巴，则滑动方向为从右到左。
4.根据权利要求1所述的一种基于滑动窗口的基因序列多聚腺苷酸尾巴提取方法，其特征在于：所述的多聚腺苷酸尾巴长度计数值初始值为初始尾巴长度N，N值为8，错配惩罚值表示A碱基片段尾巴中非A碱基个数，初始值为0；错配惩罚阈值表示A碱基片段尾巴中出现的最长连续非A碱基个数，默认值为5。
5.根据权利要求1所述的一种基于滑动窗口的基因序列多聚腺苷酸尾巴提取方法，其特征在于：所述的多聚腺苷酸尾巴长度计数值初始值为初始尾巴长度N，N值为8，错配惩罚值表示T碱基片段尾巴中非T碱基个数，初始值为0；错配惩罚阈值表示T碱基片段尾巴中出现的最长连续非T碱基个数，默认值为5。
6.根据权利要求1所述的一种基于滑动窗口的基因序列多聚腺苷酸尾巴提取方法，其特征在于：步骤（5）中若提供该物种的参考基因组序列，则将步骤（3）和步骤（4）获得的每条多聚腺苷酸尾巴序列及其5’端外侧的200个碱基序列作为待比对序列，若该条多聚腺苷酸尾巴序列5’端外侧的序列长度不足200个碱基，则将该条多聚腺苷酸尾巴序列及其5’端外侧的全部序列作为待比对序列，将待比对序列通过序列比对工具与参考基因组序列进行比对，将待比对序列中的多聚腺苷酸尾巴序列中不能比对上参考基因组序列的片段提取过滤后，获得精确度更高的多聚腺苷酸尾巴，并根据序列中提取的多聚腺苷酸尾巴数量确定步骤（1）输入的基因序列类型；若多聚腺苷酸尾巴序列与参考基因组序列碱基无差异，则表明该尾巴序列实际是来自参考基因组的序列，并非真正的多聚腺苷酸尾巴。
7.根据权利要求6所述的一种基于滑动窗口的基因序列多聚腺苷酸尾巴提取方法，其特征在于：所述的序列比对工具为minimap2、Bowtie或STAR软件。
CN202210110546.7A 2022-01-29 2022-01-29 一种基于滑动窗口的基因序列多聚腺苷酸尾巴提取方法 Active CN114582419B (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN202210110546.7A CN114582419B (zh) 2022-01-29 2022-01-29 一种基于滑动窗口的基因序列多聚腺苷酸尾巴提取方法

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN202210110546.7A CN114582419B (zh) 2022-01-29 2022-01-29 一种基于滑动窗口的基因序列多聚腺苷酸尾巴提取方法

Publications (2)

Publication Number Publication Date
CN114582419A CN114582419A (zh) 2022-06-03
CN114582419B true CN114582419B (zh) 2023-02-10

Family

ID=81769566

Family Applications (1)

Application Number Title Priority Date Filing Date
CN202210110546.7A Active CN114582419B (zh) 2022-01-29 2022-01-29 一种基于滑动窗口的基因序列多聚腺苷酸尾巴提取方法

Country Status (1)

Country Link
CN (1) CN114582419B (zh)

Citations (2)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN102732629A (zh) * 2012-08-01 2012-10-17 复旦大学 利用高通量测序同时测定基因表达量和多聚腺苷酸加尾的方法
CN104711340A (zh) * 2013-12-17 2015-06-17 北京大学 一种转录组测序方法

Family Cites Families (11)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
US20110059453A1 (en) * 2009-08-23 2011-03-10 Affymetrix, Inc. Poly(A) Tail Length Measurement by PCR
EP3461914A1 (en) * 2010-10-22 2019-04-03 Cold Spring Harbor Laboratory Varietal counting of nucleic acids for obtaining genomic copy number information
US20180265912A1 (en) * 2011-08-23 2018-09-20 Rutgers, The State University Of New Jersey Modified 3' region extraction and deep sequencing of polydenylation sites and poly(a) tail length analysis
US10829804B2 (en) * 2015-03-23 2020-11-10 The University Of North Carolina At Chapel Hill Method for identification and enumeration of nucleic acid sequences, expression, splice variant, translocation, copy, or DNA methylation changes using combined nuclease, ligase, polymerase, terminal transferase, and sequencing reactions
CN105734053B (zh) * 2016-04-20 2018-12-14 武汉生命之美科技有限公司 一种高通量测序分析Poly(A)尾长度文库的构建方法
WO2019116346A1 (en) * 2017-12-15 2019-06-20 Novartis Ag Polya tail length analysis of rna by mass spectrometry
US20210277379A1 (en) * 2018-08-03 2021-09-09 Beam Therapeutics Inc. Multi-effector nucleobase editors and methods of using same to modify a nucleic acid target sequence
CA3112262A1 (en) * 2018-10-02 2020-04-09 Max-Delbruck-Centrum Fur Molekulare Medizin In Der Helmholtz-Gemeinschaft Full-length rna sequencing
WO2020177012A1 (zh) * 2019-03-01 2020-09-10 武汉华大医学检验所有限公司 用于rna直接建库的核酸序列、基于rna样本直接构建测序文库的方法及应用
CN110499356B (zh) * 2019-09-05 2021-06-08 中国科学院遗传与发育生物学研究所 一种待测样品中具有poly(A)尾巴的RNA的测序文库的构建方法
CN112481363B (zh) * 2020-03-09 2024-05-28 南京大学 突变Aerolysin单体在检测RNA碱基序列及RNA修饰方面的应用

Patent Citations (2)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN102732629A (zh) * 2012-08-01 2012-10-17 复旦大学 利用高通量测序同时测定基因表达量和多聚腺苷酸加尾的方法
CN104711340A (zh) * 2013-12-17 2015-06-17 北京大学 一种转录组测序方法

Also Published As

Publication number Publication date
CN114582419A (zh) 2022-06-03

Similar Documents

Publication Publication Date Title
US20170335390A1 (en) Method for genotyping clonotype profiles using sequence tags
KR102356323B1 (ko) 서열 변이체 콜에 대한 검증방법 및 시스템
KR101795124B1 (ko) 복제 수 변이를 검측하기 위한 방법 및 시스템
CN110211633B (zh) Mgmt基因启动子甲基化的检测方法、测序数据的处理方法及处理装置
WO2000000637A2 (en) Method for sequencing nucleic acids with reduced errors
CN104404160A (zh) 一种mit引物设计方法和利用高通量测序构建浮游动物条形码数据库的方法
CN115312121B (zh) 靶基因位点检测方法、装置、设备及计算机存储介质
CN115691672B (zh) 针对测序平台特征的碱基质量值矫正方法、装置、电子设备和存储介质
CN112011615A (zh) 用于人甲状腺癌基因融合试剂盒及检测方法
CN114582419B (zh) 一种基于滑动窗口的基因序列多聚腺苷酸尾巴提取方法
CN106021980B (zh) 一种dna及蛋白质水平突变分析***
CN111292806B (zh) 一种利用纳米孔测序的转录组分析方法
CN104789675A (zh) 一种检测荷斯坦奶牛瘤胃微生物的方法
Thanaraj et al. Prediction of exact boundaries of exons
CN106011313B (zh) 一种快速区分iltv、ibv、mg和ms的多重荧光免疫分析方法及试剂
CN115691679A (zh) 一种基于二代和三代测序技术的宏病毒组分析方法
WO2019132010A1 (ja) 塩基配列における塩基種を推定する方法、装置及びプログラム
JP7160349B2 (ja) 核酸をシークエンシングする方法および解析する方法
CN104951673B (zh) 一种基因组酶切图谱拼接方法及***
CN104769129B (zh) 一种主要组织相容性复合体mhc分型方法及其应用
CN110684830A (zh) 一种石蜡切片组织rna分析方法
CN112750501A (zh) 一种宏病毒组流程的优化分析方法
CN114171121B (zh) 一种mRNA 5’3’末端差异的快速检测方法
CN111826428A (zh) 一种基于二代测序的微卫星不稳定性检测的方法及***
CN113403371B (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