CN113593639A - 一种用于病毒基因组变异分析、监测方法和系统 - Google Patents

一种用于病毒基因组变异分析、监测方法和系统 Download PDF

Info

Publication number
CN113593639A
CN113593639A CN202110896978.0A CN202110896978A CN113593639A CN 113593639 A CN113593639 A CN 113593639A CN 202110896978 A CN202110896978 A CN 202110896978A CN 113593639 A CN113593639 A CN 113593639A
Authority
CN
China
Prior art keywords
sequence
mutation
frequency
site
amino acid
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
CN202110896978.0A
Other languages
English (en)
Other versions
CN113593639B (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.)
Hunan University
Original Assignee
Hunan 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 Hunan University filed Critical Hunan University
Priority to CN202110896978.0A priority Critical patent/CN113593639B/zh
Publication of CN113593639A publication Critical patent/CN113593639A/zh
Application granted granted Critical
Publication of CN113593639B publication Critical patent/CN113593639B/zh
Active legal-status Critical Current
Anticipated expiration legal-status Critical

Links

Images

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
    • 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/50—Mutagenesis
    • 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
    • G16B40/00—ICT specially adapted for biostatistics; ICT specially adapted for bioinformatics-related machine learning or data mining, e.g. knowledge discovery or pattern finding
    • 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

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

Abstract

本发明属于生物信息学和生物医药领域，具体涉及一种用于病毒基因组变异分析、监测方法和系统。该方法通过序列比对后进行全基因注释获得编码基因信息，再进行翻译比对并根据遗传密码子表分析突变类型，获得基因组的突变位点、类型及频率，统计高频率突变位点，通过采集不同时间和区域的病毒株，分析高频率突变位点在时间和地区上的分布特征，并筛选位于免疫表位及其附近的高频率突变位点，以实现病毒基因组变异分析及监测。

Description

一种用于病毒基因组变异分析、监测方法和系统
技术领域
本发明属于生物信息学和生物医药领域，具体涉及一种用于病毒基因组变异分析、监测方法和系统。
背景技术
2019年底，由SARS-CoV-2引起的新冠疫情爆发，给世界公共卫生、经济、社会可持续发展带来了巨大危害及挑战。随着SARS-CoV-2在人群中的广泛流行，其基因组不断进化和变异。一些突变会赋予病毒新的遗传特征，如影响感染性和疫苗疗效等。前期相关研究表明，一些突变位点会影响病毒的感染特性以及对中和抗体甚至疫苗的免疫保护作用。
然而，疫苗的开发周期较长，从SARS-CoV-2出现到第一批在疫苗投入使用，大约历经了一年。目前，即将投入使用的疫苗是基于疫情初期的毒株开发的。在SARS-CoV-2大流行的这一年内，SARS-CoV-2基因组上经历了较多的变异，其中一些突变位点如S蛋白上的可导致感染性增强的D614G成为了优势变异且在全世界大范围流行。
因此，在病毒性疫情大流行期间，基于大尺度基因组序列数据来对病毒进行变异研究和监测，筛选并评估具有潜在功能的突变，对疫情的防控是至关重要的。
发明内容
基于此，本发明通过序列比对后进行全基因注释获得编码基因信息，再进行翻译比对并根据遗传密码子表分析突变类型，获得基因组的突变位点、类型及频率，统计高频率突变位点，通过采集不同时间和区域的病毒株，分析高频率突变位点在时间和空间上的分布特征，并筛选位于免疫表位及其附近的高频率突变位点，以实现病毒基因组变异分析及监测。
本发明提供了一种用于病毒基因组变异分析方法，所述分析方法具体包括：
获取参考基因序列和待分析基因组序列，将所述待分析基因组序列进行质量控制后获得测序质量良好的全基因组序列集；
将所述测序质量良好的序列集进行多序列比对，获得alignment序列文件；
将所述alignment序列文件以字符串完全匹配方法匹配编码基因的起始字符串，并返回编码基因在全基因组序列中的位置下标，并以alignment序列文件中每个基因对应的位置下标，生成基因位置表；
根据所述基因位置表，循环遍历每条序列，截取每条序列对应的编码基因片段并存储；
根据所述编码基因序列进行翻译成氨基酸序列，并对氨基酸序列进行多序列比对，获得对齐的氨基酸序列文件，与对齐前的编码基因序列进行匹配，采用“密码子出现的顺序”作为映射关系进而进行“回译”，将对齐的氨基酸序列“回译”成对齐的核苷酸；
将所述对齐的核苷酸序列以每三个碱基扫描方式遍历序列的每个密码子位点，识别和记录序列中存在连续三个碱基插入位点的位置“---”，标记为“插入位点”，并删除“插入位点”，获取不含插入位点的突变分析序列；
将参考序列对应的密码子和/或其对应的氨基酸与所述突变分析序列的密码子和/或翻译氨基酸按照预设突变分析方法进行分析获得待分析基因组序列突变位点以及变异类型。
进一步的，所述质量控制具体包括：
循环扫描每条序列，统计位置碱基的数量，当序列中含有10个以上的连续未知碱基，则将对应的序列进行删除。
进一步的，所述预设突变分析方法具体包括：
将参考序列对应密码子变量名记作qury_seq_conding，突变分析序列的密码子变量名记作s；
如果s和qury_seq_conding相同，当二者均不为字符串“---”，不变位点数加1，当二者都为字符串“---”，则忽略；
如果s和qury_seq_conding不同，其中有一个为字符串“---”，当qury_seq_conding为字符串“---”，标记为碱基插入位点并计数，当s为第字符串“---”，标记为碱基删除位点并计数；
如果s和qury_seq_conding不同，且都不为字符串“---”，将二者翻译成氨基酸translate_s和translate_qury_seq，当氨基酸相同，则标记为同义突变并计数，当translate_s为字符“？”，标记为未知突变并计数，当translate_s为字符“*”，标记为提前终止并计数，其他属于非同义突变位点并计数；
所述非同义突变位点，比较并记录氨基酸性质的变化。
本发明实施例还提供了一种用于病毒基因组变异的监测方法，所述监测方法具体包括：
采集不同时间和区域的待分析病毒基因组，按照上述的突变分析方法进行分析获得所有的突变位点和突变类型；统计各突变位点的突变频率，将所述突变频率高于预设频率阈值对应的突变位点记为高频突变位点；
获取所述高频突变位点对应的基因组序列的采集时间计算高频突变在所有基因组中的占比，获得折线图，并进行拟合，从而获得具有增长趋势的突变位点和突变毒株；
获取所述高频突变位点对应的区域，构建聚类热图，用于不同地区疫苗设计提供参考；
针对所述高频突变位点进行免疫表位筛选，所述免疫表位筛选具体包括：B细胞表位预测和T细胞表位预测。
进一步的，所述预设频率阈值大于等于0.5。
进一步的，所述B细胞表位预测具体包括：
连续B细胞表位预测：预测氨基酸序列的线性B细胞表位，将存在6-25个氨基酸的线性表位肽进行抗原评估，获得评分超过0.5对应的表位肽序列为连续B细胞表位；
离散B细胞表位预测：获取氨基酸序列的三级结构，并预测所述三级结构中的不连续B细胞表位，获得倾向性评分超过-3.7的为离散B细胞表位。
进一步的，所述T细胞表位预测具体包括：
MHC-I(CD8 T细胞)结合表位预测：根据人类HLAI类等位基因预测MHC-I结合表位，获取评分超过0.85且VaxiJen得分超过0.5的表位；所述人类HLAI类等位基因为HLA-A01:01，HLA-A02:01，HLA-A03:01，HLA-A11:01，HLA-A24:02，HLA-B07:02,HLA-B08:01,HLA-B40:01；
MHC-Ⅱ(CD4 T细胞)结合表位预测：根据7个人类DRB类等位基因预测MHC-Ⅱ结合表位，获得adjustedrank值低于1且VaxiJen评分超过0.5的结合表位，所述人类DRB类等位基因为DRB103:01，DRB107:01，DRB115:01，DRB301:01，DRB302:02，DRB401:01和DRB501:01。
基于同一发明构思的，本发明实施例还提供了一种用于病毒基因组变异分析、检测系统，所述系统具体包括：
筛选模块，用于筛选待分析基因组序列，将所述待分析基因组序列进行质量控制获得测序质量良好的全基因组序列集；
序列比对模块，用于将所述测序质量良好的基因组序列集进行多序列比对，获得alignment序列文件；
基因组注释模块，用于将所述alignment序列文件以字符串完全匹配方法匹配编码基因的起始字符串，并返回编码基因在全基因组序列中的位置下标，并以alignment序列文件中每个基因对应的位置下标，生成基因位置表；根据所述基因位置表，循环遍历每条序列，截取每条序列对应的编码基因片段并存储；
翻译比对模块，用于根据所述编码基因序列进行翻译成氨基酸序列，并对氨基酸序列进行多序列比对，获得对齐的氨基酸序列文件，与对齐前的编码基因序列进行匹配，采用“密码子出现的顺序”作为映射关系进而进行“回译”，将对齐的氨基酸序列“回译”成对齐的核苷酸；
插入位点标记模块；将所述对齐的核苷酸序列以每三个碱基扫描方式遍历序列的每个密码子位点，识别和记录序列中存在连续三个碱基插入位点“---”的位置，标记为“插入位点”，获得并删除“插入位点”，获取不含插入位点的突变分析序列；
编码基因突变分析模块，用于将参考序列对应的密码子和/或其对应的氨基酸与所述突变分析序列的密码子和/或翻译氨基酸按照预设突变分析方法进行分析获得待分析基因组序列突变位点以及变异类型；
监测模块，用于获取不同采集时间和区域的突变分析结果，统计高频率突变位点，并分析所述高频率突变位点在时间和区域上的规律以及免疫表位情况。
进一步的，所述监测模块具体包括：
高频率突变位点分析子模块，用于统计各突变位点的突变频率，将所述突变频率高于预设频率阈值对应的突变位点记为高频突变位点；
增长趋势监测子模块，用于获取所述高频突变位点对应的基因组序列的采集时间计算高频突变在所有基因组中的占比，获得折线图，并进行拟合，从而获得具有增长趋势的突变位点和突变毒株；
地域监测子模块，获取所述高频突变位点对应的区域，构建聚类热图，用于不同地区疫苗设计提供参考；
免疫表位筛选子模块，针对所述高频突变位点进行免疫表位筛选，所述免疫表位筛选具体包括：B细胞表位预测和T细胞表位预测。
有益效果：
该方法通过序列比对后以字符串完全匹配方式匹配编码基因的启示字符串并返回编码基因在全基因组序列中的位置下标，获得基因位置表，可以实现大量全基因组序列的快速批量注释，并采用编码基因翻译比对方法保留对齐后编码基因对应的核苷酸序列中的遗传信息，进行多种突变类型的分析，实现了病毒基因组变异的全面分析，并经过数据统计和免疫表位筛选进行突变监测，为疫苗设置提供了数据参考。
应当理解的是，以上的一般描述和后文的细节描述仅是示例性和解释性的，并不能限制本公开。
附图说明
为了更清楚地说明本发明实施例或现有技术中的技术方案，下面将对实施例或现有技术描述中所需要使用的附图作简单地介绍，显而易见地，下面描述中的附图仅仅是本发明的一些实施例，对于本领域普通技术人员来讲，在不付出创造性劳动的前提下，还可以根据这些附图获得其他的附图。
图1为本发明实施例提供的一种用于病毒基因组突变分析、监测方法的流程图；
图2为本发明实施例提供的新冠病毒的S基因序列(S_codon.fasta)示例的碱基位点突变频率分布图；
图3为本发明实施例提供的D614G+A222V，D614G+L18F+A222V突变位点在种群中的突变频率趋势图；
图4为本发明实施例提供的40种高频突变位点的地区差异分布图。
具体实施方式
为了使本发明的目的、技术方案及优点更加清楚明白，以下结合附图及实施例，对本发明进行进一步详细说明。应当理解，此处所描述的具体实施例仅用以解释本发明，并不用于限定本发明。
如图1所示，在本发明实施例中，提出了一种用于病毒基因组突变分析、监测方法的流程图，针对于突变分析过程具体包括以下步骤：
步骤S101，获取参考基因序列和待分析基因组序列，将所述待分析基因组序列进行质量控制后获得测序质量良好的全基因组序列集。在本发明实施例中，对获取的基因组序列循环扫描，统计未知碱基(N字符串)的数量，如果序列中含有10个及以上连续的未知碱基(N)，则将该条序列剔除掉，不纳入后续的分析。过滤掉一些测序质量不好的序列。
步骤S102，将所述测序质量良好的基因组序列集进行多序列比对，获得alignment序列文件；在本发明实施例中，采用MAFFT程序进行多序列比对，获得alignment序列文件。
步骤S103，将所述alignment序列文件以字符串完全匹配方法匹配编码基因的起始字符串，并返回编码基因在全基因组序列中的位置下标，并以alignment序列文件中每个基因对应的位置下标，生成基因位置表，根据所述基因位置表，循环遍历每条序列，截取每条序列对应的编码基因片段并存储。
在本发明实施例中，以字符串完全匹配的方式匹配编码基因的起始字符串，返回编码基因在全基因组序列中的位置下标，进而对基因组进行序列注释。在进行注释前，需要先运行检查机制，判断用于识别的字符串是否在序列中存在。判断模式为：如果字符串存在，则继续运行下一步，否则输出提示信息，提示更换字符串。基于用户提供的参考序列的每个编码基因的信息文件，该文件包含基因名称，起始字符串，结束字符串，分别以逗号隔开。使用find函数，基于完全匹配(两个字符串相等)的方式，获得alignment序列文件中每个基因的对应的下标值(在序列中的坐标)，生成基因位置表。再循环遍历每条序列，基于基因位置表，截取出每条序列对应的编码基因片段，分别储存以每个编码基因名称命名的文件里。同时，对于每条序列，都输出一个*.gb格式的注释文件，里面包含了该全基因组序列的字符串信息，每个基因的基因名称和起始位置。该步骤可以实现大量全基因组序列数据的快速注释。
步骤S104，根据所述编码基因序列进行翻译比对，先将核苷酸序列翻译成氨基酸序列，并对氨基酸序列进行多序列比对，获得对齐的氨基酸序列文件，再“回译”成对齐的核苷酸序列。该比对过程可采用MAFFT程序。
步骤S105，根据所述的对齐的核苷酸序列。编码基因进行翻译比对后，对齐的序列文件中可能存在的“插入位点”(作为突变的一种)，会在参考序列中引入gap(-)，因为gap的位置不确定性，在不同次分析中会影响后续对突变位置的统一标记，因此，需要将“插入位点”单独进行统计和记录，并删除插入位点再用于后续分析。
将所述对齐的核苷酸序列以每三个碱基扫描方式遍历序列的每个密码子位点，识别和记录序列中存在连续三个碱基插入位点gap(---)的位置，标记为“插入位点”，并返回插入位点处其他非参考序列的密码子和翻译后的氨基酸，同时返回不含插入位点的突变分析序列。
步骤S106，将参考序列对应的密码子和/或其对应的氨基酸与所述突变分析序列密码子和/或翻译氨基酸按照预设突变分析方法进行分析获得待分析基因组序列突变位点以及变异类型。
在本发明实施例中，使用两个嵌套的循环，外循环的迭代数为密码子位点数量(总碱基数量除以3)，内循环迭代数为序列的条数。对于每个内循环，每条序列都与序列集中的第一条序列(参考序列)进行比较，从而识别5种不同的突变类型，包括同义替换、非同义替换、提前终止、碱基插入和删除。具体识别模式和比较方法为：
1)将参考序列对应密码子变量名记作qury_seq_conding，用于突变分析的序列的密码子变量名记作s；
2)如果s和qury_seq_conding相同。a.如果qury_seq_conding和s都不为“---”，不变位点数加1；b.如果s和qury_seq_conding都为“---”，则忽略；
3)如果s和qury_seq_conding不相同，其中有一个为“---”。a.如果qury_seq_conding为“---”，标记为碱基插入位点并计数；b.如果s为"---"，标记为碱基删除位点并计数；
4)如果s和qury_seq_conding不相同，且两个都不为“---”。将s和qury_seq_conding分别翻译成氨基酸，记作translate_s和translate_qury_seq。a.如果translate_s和translate_qury_seq相同，标记为同义突变并计数；b.如果translate_s为字符“？”，标记为未知突变并计数；c.如果translate_s为字符“*”，标记为提前终止并计数；d.如果都不符合上述条件，则标记为非同义突变并计数。
5)对于非同义突变位点，调用BioAider里内置的氨基酸性质表，判断并记录氨基酸的性质是否改变，包括极性和带电性。
6)对于上述突变分析步骤，写入生成详细的日志文件(log文件)以及汇总文件(summary文件)。log文件记录每条序列每个位点上的变异，汇总文件记录突变的密码子位点，突变的碱基，突变类型，突变频率，突变氨基酸是否改变等汇总信息。
7)根据用户指定的突变频数分布表，生成同义或非同义非同义突变位点对应的频数分布图，以便了解序列数据集整体的变异情况，并以位图*.jpg和矢量图.pdf的格式输出。
对于关联突变分析，可逐一扫描所有核苷酸位点与参考序列进行对比，汇总每条序列所有的突变位点信息，生成汇总报告(summary文件)，包括突变位点及突变频率、突变热点(高于指定的突变频率阈值)，并生成详细的日志记录文件(log文件)，详细记录每条序列每个位点上的变异，便于查询。
将上述编码基因进行关联突变分析，以每条序列作为一个整体标记所有的突变位点，如“L18F+D614G+A222V”，每条序列生成一个记录，并对它们进行统计得到这些关联突变位点的突变频率，最后根据指定的突变频率阈值输出高于阈值的突变热点。
所述监测方法具体包括：
采集不同时间和区域的待分析病毒基因组，按照上述突变分析方法进行分析获得所有的突变位点和突变类型；统计各突变位点的突变频率，将所述突变频率高于预设频率阈值对应的突变位点记为高频突变位点。如图2所示，上述基于密码子方式比对后的新冠病毒的S基因序列(S_codon.fasta)示例的碱基位点突变频率分布图。经关联突变分析获得40中突变频率大于0.1％的关联突变位点，如表1所示。
表1突变频率大于0.1％的关联突变位点
Figure BDA0003198346530000111
Figure BDA0003198346530000121
读取summary文件中突变位点汇总信息，默认选取突变频率(突变序列的数量占比)高于0.5％的突变位点作为用于时间分析的高频率突变位点。根据序列的采样时间，默认按月份为时间间隔，计算在每个时间间隔内这些高频率突变在群体中的占比，绘制折线图。此外，采用线性回归的方法进行拟合，计算回归系数，从而识别出具有增长趋势的突变位点和突变毒株。如图3示例中D614G+A222V，D614G+L18F+A222V突变位点在种群中的突变频率趋势图。
读取summary文件中突变位点汇总信息，默认选取突变频率(突变序列的数量占比)高于0.1％的突变位点作为用于地区分布分析的高频率突变位点。构建高频率突变位点及其在每个地区的数量矩阵，行为地区，列为突变位点，构建聚类热图。聚类热图采用欧式距离对行和列进行聚类，聚集在一起的行表示这些地区具有相似流行的突变毒株，聚集在一起的列表示这些突变位点具有相似地区分布。因为疫苗设计需要充分考虑现有毒株的流行情况，因此，对毒株基于高频率突变位点进行地理聚类，能为不同地区的疫苗设计提供参考，上述示例中的40种高频突变位点的地区差异分布图如图4所示。
对所述突变频率高于0.5％的高频突变位点进行免疫表位筛选，所述免疫表位筛选具体包括：B细胞表位预测和T细胞表位预测。
在本发明实施例中，连续B细胞表位预测：使用IEDB(http://www.iedb.org/)中的BepiPred 2.0程序来预测特定蛋白氨基酸序列的线性B细胞表位，使用默认阈值0.5。对于超过6个但小于25个氨基酸的线性表位肽，使用VaxiJen2.0在线服务器进行抗原评估，仅考虑使用VaxiJen预测得分超过0.5的表位肽。如果表位肽的长度超过25，则只考虑VaxiJen评分超过0.5的子片段序列。离散B细胞表位预测：对于特定蛋白的氨基酸序列，如果存在该蛋白质序列的三级结构，则直接使用。如果不存在，通过同源建模(SWISS-MODEL)或者从头预测(AlphaFold或RoseTTAFold)方式获得其三级结构。再通过IEDB中的DiscoTope 2.0程序预测特定蛋白三级结构中的不连续B细胞表位，并且只考虑倾向性评分(propensity)和discotope超过默认阈值-3.7的氨基酸残基。
在本发明实施例中，所述T细胞表位预测具体包括：MHC-I(CD8 T细胞)结合表位预测：使用IEDB中的NetMHCpan4.1方法基于8个最常见的人类HLAI类等位基因(HLA-A01:01，HLA-A02:01，HLA-A03:01，HLA-A11:01，HLA-A24:02，HLA-B07:02,HLA-B08:01,HLA-B40:01)预测了MHC-I结合表位,并且只考虑了评分超过0.85且VaxiJen得分超过0.5的表位。MHC-Ⅱ(CD4 T细胞)结合表位预测：基于常见的7个人类等位基因(包括DRB103:01，DRB107:01，DRB115:01，DRB301:01，DRB302:02，DRB401:01和DRB5*01:01)和IEDB推荐的2.22算法，表位肽的长度设置为15。类似地，MHC-Ⅱ的表位预测仅考虑adjusted rank值低于1且VaxiJen评分超过0.5的结合表位。以新冠病毒S蛋白氨基酸序列为例，采用上述方法进行B细胞表位，其结果如表2、3所示。
区域 开始 终止 免疫表位(BepiPred) 长度 VaxiJen得分
NTD 13 37 SQCVNLTTRTQLPPAYTNSFTRGVY 25 0.6860
NTD 59 81 FSNVTWFHAIHVSGTNGTKRFDN 23 0.6767
NTD 138 154 DPFLGVYYHKNNKSWME 17 0.5821
NTD 177 189 MDLEGKQGNFKNL 13 1.2592
NTD 206 221 KHTPINLVRDLPQGFS 16 0.6403
NTD 253 259 DSSSGWT 7 0.6067
S1 304 322 KSFTVEKGIYQTSNFRVQP 19 0.5729
RBD 351 360 YAWNRKRISN 10 0.5855
RBD 369 393 YNSASFSTFKCYGVSPTKLNDLCFT 25 1.4031
RBD 404 426 GDEVRQIAPGQTGKIADYNYKLP 23 1.1017
RBD 441 448 LDSKVGGN 8 0.8773
RBD 459 464 SNLKPF 6 0.5943
RBD 473 478 YQAGST 6 0.5812
RBD 487 492 NCYFPL 6 0.996
RBD 528 535 KKSTNLVK 8 0.658
S1 555 562 SNKKFLPF 8 1.3952
S1 627 632 DQLTPT 6 0.7329
S1 656 666 VNNSYECDIPI 11 0.6124
S1,S2 680 687 SPRRARSV 8 0.6844
S2 695 710 YTMSLGAENSVAYSNN 16 0.6434
S2 809 814 PSKPSK 6 1.1271
S2 1035 1043 GQSKRVDFC 9 1.779
S2 1110 1118 YEPQIITTD 9 0.8297
S2 1154 1169 KYFKNHTSPDVDLGDI 16 0.7333
CP 1255 1267 KFDEDDSEPVLKG 13 0.5066
表2连续B细胞表位
由表2可知，突变频率高于0.5％的高频率突变位点中，突变位点L18F、L18F和S477N突变热点位于连续B细胞免疫表位上，A222V突变热点位于连续B细胞免疫表位附近。L18F，A222V，S477N这3个高频率突变，值得后续进一步研究其突变对中和抗体及疫苗保护效果的影响。
表3离散的B细胞表位
Figure BDA0003198346530000141
Figure BDA0003198346530000151
Figure BDA0003198346530000161
由表3可知，突变频率高于0.5％的高频率突变位点中，突变位点N439K，S477N位于离散B细胞免疫表位上。N439K，S477N这2个高频率突变，值得后续进一步研究其突变对中和抗体及疫苗保护效果的影响。
本发明实施例还提供了一种用于病毒基因组变异分析、检测系统，所述系统具体包括：
筛选模块，用于筛选待分析基因组序列，将所述待分析基因组序列进行质量控制获得测序质量良好的全基因组序列集。
序列比对模块，用于将所述测序质量良好的基因组序列集进行多序列比对，获得alignment序列文件；
基因组注释模块，用于将所述alignment序列文件以字符串完全匹配方法匹配编码基因的起始字符串，并返回编码基因在全基因组序列中的位置下标，并以alignment序列文件中每个基因对应的位置下标，生成基因位置表；根据所述基因位置表，循环遍历每条序列，截取每条序列对应的编码基因片段并存储；
翻译比对模块，用于根据所述编码基因序列进行翻译成氨基酸序列，并对氨基酸序列进行多序列比对，获得对齐的氨基酸序列文件，与对齐前的编码基因序列进行匹配，采用“密码子出现的顺序”作为映射关系进而进行“回译”，将对齐的氨基酸序列“回译”成对齐的核苷酸；
插入位点标记模块；将所述对齐的核苷酸序列以每三个碱基扫描方式遍历序列的每个密码子位点，识别和记录序列中存在连续三个碱基插入位点“---”的位置，标记为“插入位点”，获得并删除“插入位点”，获取不含插入位点的突变分析序列；
编码基因突变分析模块，用于将参考序列对应的密码子和/或其对应的氨基酸与所述突变分析序列的密码子和/或翻译氨基酸按照预设突变分析方法进行分析获得待分析基因组序列突变位点以及变异类型；
监测模块，用于获取不同采集时间和区域的突变分析结果，统计高频率突变位点，并分析所述高频率突变位点在时间和区域上的规律以及免疫表位情况。
所述监测模块具体包括：高频率突变位点分析子模块，用于统计各突变位点的突变频率，将所述突变频率高于预设频率阈值对应的突变位点记为高频突变位点；增长趋势监测子模块，用于获取所述高频突变位点对应的基因组序列的采集时间计算高频突变在所有基因组中的占比，获得折线图，并进行拟合，从而获得具有增长趋势的突变位点和突变毒株；地域监测子模块，获取所述高频突变位点对应的区域，构建聚类热图，用于不同地区疫苗设计提供参考；免疫表位筛选子模块，对所述高频突变位点进行免疫表位筛选，所述免疫表位筛选具体包括：B细胞表位预测和T细胞表位预测。
以上所述实施例仅表达了本发明的几种实施方式，其描述较为具体和详细，但并不能因此而理解为对本发明专利范围的限制。应当指出的是，对于本领域的普通技术人员来说，在不脱离本发明构思的前提下，还可以做出若干变形和改进，这些都属于本发明的保护范围。因此，本发明专利的保护范围应以所附权利要求为准。
本领域技术人员在考虑说明书及实践这里公开的发明后，将容易想到本公开的其它实施方案。本申请旨在涵盖本公开的任何变型、用途或者适应性变化，这些变型、用途或者适应性变化遵循本公开的一般性原理并包括本公开未公开的本技术领域中的公知常识或惯用技术手段。说明书和实施例仅被视为示例性的，本公开的真正范围和精神由权利要求指出。
应该理解的是，虽然本发明各实施例的流程图中的各个步骤按照箭头的指示依次显示，但是这些步骤并不是必然按照箭头指示的顺序依次执行。除非本文中有明确的说明，这些步骤的执行并没有严格的顺序限制，这些步骤可以以其它的顺序执行。而且，各实施例中的至少一部分步骤可以包括多个子步骤或者多个阶段，这些子步骤或者阶段并不必然是在同一时刻执行完成，而是可以在不同的时刻执行，这些子步骤或者阶段的执行顺序也不必然是依次进行，而是可以与其它步骤或者其它步骤的子步骤或者阶段的至少一部分轮流或者交替地执行。
以上所述实施例的各技术特征可以进行任意的组合，为使描述简洁，未对上述实施例中的各个技术特征所有可能的组合都进行描述，然而，只要这些技术特征的组合不存在矛盾，都应当认为是本说明书记载的范围。

Claims (9)

1.一种用于病毒基因组变异分析方法，其特征在于，所述分析方法具体包括：
获取参考基因序列和待分析基因组序列，将所述待分析基因组序列进行质量控制后获得测序质量良好的全基因组序列集；
将所述测序质量良好的序列集进行多序列比对，获得alignment序列文件；
将所述alignment序列文件以字符串完全匹配方法匹配编码基因的起始字符串，并返回编码基因在全基因组序列中的位置下标，并以alignment序列文件中每个基因对应的位置下标，生成基因位置表；
根据所述基因位置表，循环遍历每条序列，截取每条序列对应的编码基因片段并存储；
根据所述编码基因序列进行翻译成氨基酸序列，并对氨基酸序列进行多序列比对，获得对齐的氨基酸序列文件，与对齐前的编码基因序列进行匹配，采用“密码子出现的顺序”作为映射关系进而进行“回译”，将对齐的氨基酸序列“回译”成对齐的核苷酸；
将所述对齐的核苷酸序列以每三个碱基扫描方式遍历序列的每个密码子位点，识别和记录序列中存在连续三个碱基插入位点的位置“---”，标记为“插入位点”，并删除“插入位点”，获取不含插入位点的突变分析序列；
将参考序列对应的密码子和/或其对应的氨基酸与所述突变分析序列的密码子和/或翻译氨基酸按照预设突变分析方法进行分析获得待分析基因组序列突变位点以及变异类型。
2.根据权利要求1所述的用于病毒基因组变异分析方法，其特征在于，所述质量控制具体包括：
循环扫描每条序列，统计位置碱基的数量，当序列中含有10个以上的连续未知碱基，则将对应的序列进行删除。
3.根据权利要求1所述的用于病毒基因组变异分析方法，其特征在于，所述预设突变分析方法具体包括：
将参考序列对应密码子变量名记作qury_seq_conding，突变分析序列的密码子变量名记作s；
如果s和qury_seq_conding相同，当二者均不为字符串“---”，不变位点数加1，当二者都为字符串“---”，则忽略；
如果s和qury_seq_conding不同，其中有一个为字符串“---”，当qury_seq_conding为字符串“---”，标记为碱基插入位点并计数，当s为第字符串“---”，标记为碱基删除位点并计数；
如果s和qury_seq_conding不同，且都不为字符串“---”，将二者翻译成氨基酸translate_s和translate_qury_seq，当氨基酸相同，则标记为同义突变并计数，当translate_s为字符“？”，标记为未知突变并计数，当translate_s为字符“*”，标记为提前终止并计数，其他属于非同义突变位点并计数；
所述非同义突变位点，比较并记录氨基酸性质的变化。
4.一种用于病毒基因组变异的监测方法，其特征在于，所述监测方法具体包括：
采集不同时间和区域的待分析病毒基因组，按照权利要求1-3任意所述的突变分析方法进行分析获得所有的突变位点和突变类型；统计各突变位点的突变频率，将所述突变频率高于预设频率阈值对应的突变位点记为高频突变位点；
获取所述高频突变位点对应的基因组序列的采集时间计算高频突变在所有基因组中的占比，获得折线图，并进行拟合，从而获得具有增长趋势的突变位点和突变毒株；
获取所述高频突变位点对应的区域，构建聚类热图，用于不同地区疫苗设计提供参考；
针对所述高频突变位点进行免疫表位筛选，所述免疫表位筛选具体包括：B细胞表位预测和T细胞表位预测。
5.根据权利要求4所述的用于病毒基因组变异的监测方法，其特征在于，所述预设频率阈值大于等于0.5。
6.根据权利要求4所述的用于病毒基因组变异的监测方法，其特征在于，所述B细胞表位预测具体包括：
连续B细胞表位预测：预测氨基酸序列的线性B细胞表位，将存在6-25个氨基酸的线性表位肽进行抗原评估，获得评分超过0.5对应的表位肽序列为连续B细胞表位；
离散B细胞表位预测：获取氨基酸序列的三级结构，并预测所述三级结构中的不连续B细胞表位，获得倾向性评分超过-3.7的为离散B细胞表位。
7.根据权利要求4所述的用于病毒基因组变异的监测方法，其特征在于，所述T细胞表位预测具体包括：
MHC-I(CD8 T细胞)结合表位预测：根据人类HLAI类等位基因预测MHC-I结合表位，获取评分超过0.85且VaxiJen得分超过0.5的表位；所述人类HLAI类等位基因为HLA-A01:01，HLA-A02:01，HLA-A03:01，HLA-A11:01，HLA-A24:02，HLA-B07:02,HLA-B08:01,HLA-B40:01；
MHC-Ⅱ(CD4 T细胞)结合表位预测：根据7个人类DRB类等位基因预测MHC-Ⅱ结合表位，获得adjusted rank值低于1且VaxiJen评分超过0.5的结合表位，所述人类DRB类等位基因为DRB103:01，DRB107:01，DRB115:01，DRB301:01，DRB302:02，DRB401:01和DRB501:01。
8.一种用于病毒基因组变异分析、检测系统，其特征在于，所述系统具体包括：
筛选模块，用于筛选待分析基因组序列，将所述待分析基因组序列进行质量控制获得测序质量良好的全基因组序列集；
序列比对模块，用于将所述测序质量良好的基因组序列集进行多序列比对，获得alignment序列文件；
基因组注释模块，用于将所述alignment序列文件以字符串完全匹配方法匹配编码基因的起始字符串，并返回编码基因在全基因组序列中的位置下标，并以alignment序列文件中每个基因对应的位置下标，生成基因位置表；根据所述基因位置表，循环遍历每条序列，截取每条序列对应的编码基因片段并存储；
翻译比对模块，用于根据所述编码基因序列进行翻译成氨基酸序列，并对氨基酸序列进行多序列比对，获得对齐的氨基酸序列文件，与对齐前的编码基因序列进行匹配，采用“密码子出现的顺序”作为映射关系进而进行“回译”，将对齐的氨基酸序列“回译”成对齐的核苷酸；
插入位点标记模块；将所述对齐的核苷酸序列以每三个碱基扫描方式遍历序列的每个密码子位点，识别和记录序列中存在连续三个碱基插入位点“---”的位置，标记为“插入位点”，获得并删除“插入位点”，获取不含插入位点的突变分析序列；
编码基因突变分析模块，用于将参考序列对应的密码子和/或其对应的氨基酸与所述突变分析序列的密码子和/或翻译氨基酸按照预设突变分析方法进行分析获得待分析基因组序列突变位点以及变异类型；
监测模块，用于获取不同采集时间和区域的突变分析结果，统计高频率突变位点，并分析所述高频率突变位点在时间和区域上的规律以及免疫表位情况。
9.根据权利要求8所述的用于病毒基因组变异分析、监测系统，其特征在于，所述监测模块具体包括：
高频率突变位点分析子模块，用于统计各突变位点的突变频率，将所述突变频率高于预设频率阈值对应的突变位点记为高频突变位点；
增长趋势监测子模块，用于获取所述高频突变位点对应的基因组序列的采集时间计算高频突变在所有基因组中的占比，获得折线图，并进行拟合，从而获得具有增长趋势的突变位点和突变毒株；
地域监测子模块，获取所述高频突变位点对应的区域，构建聚类热图，用于不同地区疫苗设计提供参考；
免疫表位筛选子模块，针对所述高频突变位点进行免疫表位筛选，所述免疫表位筛选具体包括：B细胞表位预测和T细胞表位预测。
CN202110896978.0A 2021-08-05 2021-08-05 一种用于病毒基因组变异分析、监测方法和系统 Active CN113593639B (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN202110896978.0A CN113593639B (zh) 2021-08-05 2021-08-05 一种用于病毒基因组变异分析、监测方法和系统

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN202110896978.0A CN113593639B (zh) 2021-08-05 2021-08-05 一种用于病毒基因组变异分析、监测方法和系统

Publications (2)

Publication Number Publication Date
CN113593639A true CN113593639A (zh) 2021-11-02
CN113593639B CN113593639B (zh) 2023-08-25

Family

ID=78255462

Family Applications (1)

Application Number Title Priority Date Filing Date
CN202110896978.0A Active CN113593639B (zh) 2021-08-05 2021-08-05 一种用于病毒基因组变异分析、监测方法和系统

Country Status (1)

Country Link
CN (1) CN113593639B (zh)

Cited By (8)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN114276422A (zh) * 2021-11-09 2022-04-05 中国人民解放军总医院 新型冠状病毒s蛋白多肽抗原及其应用
CN114550816A (zh) * 2022-03-01 2022-05-27 上海图灵智算量子科技有限公司 基于光子芯片病毒变异概率的预测方法
CN115312122A (zh) * 2022-10-12 2022-11-08 之江实验室 一种CRISPR-Cas酶可突变位点推荐方法和装置
CN115798578A (zh) * 2022-12-06 2023-03-14 中国人民解放军军事科学院军事医学研究院 一种分析与检测病毒新流行变异株的装置及方法
CN116343923A (zh) * 2023-03-21 2023-06-27 哈尔滨工业大学 一种基因组结构变异同源性识别方法
CN116741268A (zh) * 2023-04-04 2023-09-12 中国人民解放军军事科学院军事医学研究院 筛选病原体关键突变的方法、装置及计算机可读存储介质
WO2023180962A1 (en) * 2022-03-22 2023-09-28 Waters Technologies Ireland Limited Programmatic processing of protein or nucleic acid sequences to identify mutations at programmatically determined subsequences
CN117373527A (zh) * 2023-12-07 2024-01-09 中国科学院微生物研究所 Hiv序列质控方法、设备及存储介质

Citations (5)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN1813061A (zh) * 2003-04-25 2006-08-02 免疫医疗疫苗公司 重组副流感病毒表达系统以及包含源自间质肺病毒的异种抗原的疫苗
CN107122624A (zh) * 2017-05-01 2017-09-01 杨永臣 人类基因突变的hgvs名称生成及分析系统的实现方法
WO2019123398A1 (en) * 2017-12-21 2019-06-27 New Zealand Health Innovation Hub Management Limited Method of analysis of mutations in the hepatitis b virus and uses thereof
US20190237158A1 (en) * 2016-08-31 2019-08-01 Medgenome, Inc. Methods to analyze genetic alterations in cancer to identify therapeutic peptide vaccines and kits therefore
CN111445955A (zh) * 2020-04-10 2020-07-24 广州微远基因科技有限公司 新型冠状病毒变异分析方法及应用

Patent Citations (5)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN1813061A (zh) * 2003-04-25 2006-08-02 免疫医疗疫苗公司 重组副流感病毒表达系统以及包含源自间质肺病毒的异种抗原的疫苗
US20190237158A1 (en) * 2016-08-31 2019-08-01 Medgenome, Inc. Methods to analyze genetic alterations in cancer to identify therapeutic peptide vaccines and kits therefore
CN107122624A (zh) * 2017-05-01 2017-09-01 杨永臣 人类基因突变的hgvs名称生成及分析系统的实现方法
WO2019123398A1 (en) * 2017-12-21 2019-06-27 New Zealand Health Innovation Hub Management Limited Method of analysis of mutations in the hepatitis b virus and uses thereof
CN111445955A (zh) * 2020-04-10 2020-07-24 广州微远基因科技有限公司 新型冠状病毒变异分析方法及应用

Non-Patent Citations (1)

* Cited by examiner, † Cited by third party
Title
雷永良;王晓光;陶晓燕;李浩;孟胜利;陈秀英;柳付明;叶碧峰;唐青;: "浙江地区鼬獾和犬源狂犬病病毒分离株全基因组测序与分析", 病毒学报, no. 01 *

Cited By (12)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN114276422A (zh) * 2021-11-09 2022-04-05 中国人民解放军总医院 新型冠状病毒s蛋白多肽抗原及其应用
WO2023083092A1 (zh) * 2021-11-09 2023-05-19 中国人民解放军总医院 新型冠状病毒s蛋白多肽抗原及其应用
CN114550816A (zh) * 2022-03-01 2022-05-27 上海图灵智算量子科技有限公司 基于光子芯片病毒变异概率的预测方法
WO2023180962A1 (en) * 2022-03-22 2023-09-28 Waters Technologies Ireland Limited Programmatic processing of protein or nucleic acid sequences to identify mutations at programmatically determined subsequences
CN115312122A (zh) * 2022-10-12 2022-11-08 之江实验室 一种CRISPR-Cas酶可突变位点推荐方法和装置
CN115312122B (zh) * 2022-10-12 2022-12-16 之江实验室 一种CRISPR-Cas酶可突变位点推荐方法和装置
CN115798578A (zh) * 2022-12-06 2023-03-14 中国人民解放军军事科学院军事医学研究院 一种分析与检测病毒新流行变异株的装置及方法
CN116343923A (zh) * 2023-03-21 2023-06-27 哈尔滨工业大学 一种基因组结构变异同源性识别方法
CN116343923B (zh) * 2023-03-21 2023-12-08 哈尔滨工业大学 一种基因组结构变异同源性识别方法
CN116741268A (zh) * 2023-04-04 2023-09-12 中国人民解放军军事科学院军事医学研究院 筛选病原体关键突变的方法、装置及计算机可读存储介质
CN116741268B (zh) * 2023-04-04 2024-03-01 中国人民解放军军事科学院军事医学研究院 筛选病原体关键突变的方法、装置及计算机可读存储介质
CN117373527A (zh) * 2023-12-07 2024-01-09 中国科学院微生物研究所 Hiv序列质控方法、设备及存储介质

Also Published As

Publication number Publication date
CN113593639B (zh) 2023-08-25

Similar Documents

Publication Publication Date Title
CN113593639A (zh) 一种用于病毒基因组变异分析、监测方法和系统
US10777301B2 (en) Hierarchical genome assembly method using single long insert library
EP2834762B1 (en) Sequence assembly
Yuan et al. CNV_IFTV: an isolation forest and total variation-based detection of CNVs from short-read sequencing data
CA2424031C (en) System and process for validating, aligning and reordering genetic sequence maps using ordered restriction map
CN108830044B (zh) 用于检测癌症样本基因融合的检测方法和装置
CN110379458A (zh) 致病性变异位点判定方法、装置、计算机设备及存储介质
CN110621785B (zh) 基于三代捕获测序对二倍体基因组单倍体分型的方法和装置
WO2018218788A1 (zh) 一种基于全局种子打分优选的三代测序序列比对方法
CN110289047B (zh) 基于测序数据的肿瘤纯度及绝对拷贝数预测方法及系统
CN115083521B (zh) 一种单细胞转录组测序数据中肿瘤细胞类群的鉴定方法及系统
Peng et al. A novel codon-based de Bruijn graph algorithm for gene construction from unassembled transcriptomes
CN107292129A (zh) 易感基因型检测方法
Kearse et al. The Geneious 6.0. 3 read mapper
Vasquez-Gross et al. A haplotype-phased genome of wheat stripe rust pathogen Puccinia striiformis f. sp. tritici, race PST-130 from the Western USA
He et al. Phylogenomics reveal extensive phylogenetic discordance due to incomplete lineage sorting following the rapid radiation of alpine butterflies (Papilionidae: Parnassius)
CN112489727B (zh) 一种快速获取罕见病致病位点的方法和系统
McLay et al. Phylogenomics reveals extreme gene tree discordance in a lineage of dominant trees: Hybridization, introgression, and incomplete lineage sorting blur deep evolutionary relationships despite clear species groupings in Eucalyptus subgenus Eudesmia
CN109308935A (zh) 一种基于支持向量机预测非编码dna的方法及应用平台
Gu et al. SVLR: genome structural variant detection using Long-read sequencing data
CN110476215A (zh) 用于多序列文件的签名-散列
Krause et al. Sensitive and error-tolerant annotation of protein-coding DNA with BATH
Spang et al. Sequence database search using jumping alignments.
JP5433894B2 (ja) 立体構造データ帰属方法、立体構造データ帰属プログラム及び立体構造データ帰属装置
Bianchetti et al. vALId: validation of protein sequence quality based on multiple alignment data

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