CN115881221A - 一种冠状病毒转录组鉴定分析方法及分析系统 - Google Patents

Abstract

本发明公开了一种冠状病毒转录组鉴定分析方法，所述方法包括使用bwa软件将短reads档案库SRA中获得的冠状病毒二代测序原始fastq文件与NCBI参考基因组进行序列比对，生成BAM文件；通过bcftools调用和过滤单核苷酸多态性SNP，并用vcf‑annotator对SNP进行注释；对BAM文件进行CIGAR字符串解析及断点识别操作，进行sgRNA鉴定；将sgRNA鉴定结果生成表达矩阵并进行转录组分析。本发明还公开了一种冠状病毒转录组鉴定分析系统，所述系统包括储存器和处理器；所述存储器上存储有计算机程序，当所述计算机程序被所述处理器执行时，实现如上所述的方法。本发明所述的冠状病毒转录组鉴定分析方法及系统为冠状病毒生物学提供了新的认知手段，并为未来的治疗方法的发展提供了宝贵的资源。

Description

一种冠状病毒转录组鉴定分析方法及分析系统

技术领域

本发明属于生物信息学技术领域，涉及一种冠状病毒转录组鉴定分析方法及系统。

背景技术

新型冠状病毒肺炎(COVID-19)的病原体是严重急性呼吸综合征冠状病毒-2(SARS-COV-2)，一种类似于MERS-CoV(唯一的另一种仍在流行的β冠状病毒)的β冠状病毒。MERS-CoV是中东呼吸综合征(MERS)的病原体，与SARS-CoV-2相比，其毒力更强，但传染性较弱，与SARS-CoV-2的系统发生距离稍远(氨基酸序列同源性小于90％)。这两种病毒都有一个约30000碱基的单正链多聚腺苷化RNA基因组，编码4种结构蛋白(刺突蛋白 (S)、膜蛋白(M)、包膜蛋白(E)和核衣壳蛋白(N))，这些蛋白在此二类病毒中功能十分相近。但是，这两种病毒在侵入细胞的受体、附加蛋白(accessory proteins)和16种非结构蛋白(nsp1至nsp16)的特异功能方面存在分歧。Nsp是通过病毒蛋白酶切割ORF1a和ORF1b 编码的两个大的多聚蛋白而产生的。ORF1最接近5’端，是进入宿主细胞后由基因组RNA直接翻译而来，根据核糖体跳跃式机制，可将其分为ORF1a和ORF1b^[1]。MERS-CoV至少编码 5种附加蛋白(ORF3、ORF4a、ORF4b、ORF5和ORF8b)，而SARS-CoV-2至少编码6种附加蛋白(ORF3a、ORF6、ORF7a、ORF7b、ORF8和ORF10)^[2]。所有除ORF1a和ORF1b 蛋白质，必须由sgRNA翻译而来^[3-4]。SgRNA通过一种称为非连续延伸的机制，借由称为转录调控序列(TRS)的长度可变的存在于各个基因之间的短序列(通常是6到12个核苷酸(nt)) 与基因组首部70碱基位置进行配对，而后将负链延伸到正链的5端，产生短负链sgRNA中间体。然后复制RNA中间体，生成编码病毒蛋白的正链sgRNA^[5]。

病毒生物学的基础是病毒基因及其功能的鉴定和注释。由于正链RNA病毒序列本身的特性以及亚基因组RNA(sgRNA)的存在，获取冠状病毒转录组信息本身就是一项困难的工作。病毒转录组的注释是了解病毒生物学的基础，而病毒生物学是阻止病毒传播、复制和致病的关键。此前的冠状病毒疫情，如2003年爆发的严重急性呼吸系统综合症(SARS)和2012年开始并仍在持续的MERS疫情^[6-7]导致了对这些人畜共患冠状病毒的研究的增加，以便了解这些病毒的起源。通过比较不同冠状病毒的转录组变化，可以揭示其独特的致病性和传染性背后的机制，并可能解释种间传播背后的分子机制。系统地注释隐藏在宏转录组数据中的冠状病毒转录谱的差异，有助于进一步了解病毒的传播性和毒力。然而，目前对这些病毒的体外转录谱系统比较仍缺少相应研究。

对于新出现的SARS-CoV-2病毒，测序在菌株进化的诊断和监测中发挥着至关重要的作用^[2,8]。然而，目前来说，SARS-CoV-2和MERS-CoV的测序数据集限于体外细胞系感染或感染模式生物过程中产生的病毒和宿主转录本。由于缺乏合适的分析工具，来自人类不同毒株的病毒转录组的分析则被忽视了。

序列同源性在病毒基因功能注释中起着至关重要的作用。然而，单靠序列同源性并不能保证蛋白质的表达，因为快速突变的RNA病毒可能包含序列的改变，从而导致新的ORF产生或原有ORF无法转录。因此，直接分析病毒RNA是了解哪些病毒基因可以表达的重要一环。就SARS-CoV-2而言，近期有研究利用牛津纳米孔技术对培养细胞系中产生的病毒RNA进行了分析，确定了典型和非典型病毒转录本的存在。这些研究都使用分离的病毒株感染从非洲绿猴肾上皮细胞中分离出来的Vero细胞系，该细胞系在感染时不会启动干扰素(IFN) 反应。虽然这些研究确定了病毒转录组的基本特征，但这些独立的研究只描述一个毒株的转录组且无法确定病毒应对最基本免疫反应时表达谱的改变(例如IFNxx)^[9-11]。

发明内容

为了解决现有技术存在的不足，本发明的目的是提供一种冠状病毒转录组鉴定分析方法，所述方法采用冠状病毒注释器(coronavirus annotator)量化病毒基因表达，并在众多公开的宏转录组数据集中识别真实可靠的sgRNA。本发明除了概括sgRNA谱的变化及其相对表达外，还可以确定几种不同冠状病毒的新型sgRNA；还能够提出SARS-CoV和 SARS-CoV-2之间共享的核心sgRNA序列，以及MERS-CoV特有的sgRNA序列。此外，在蝙蝠和穿山甲体内发现的相关冠状病毒中，SARS-CoV-2和MERS-CoV的一个新sgRNA 亚群在进化上似乎是保守的。同时特定的sgRNA的转录在体内、体外以及不同冠状病毒之间存在显著差异。

本发明提供了一种冠状病毒转录组鉴定分析方法，所述方法具体包括如下步骤：

步骤一、使用bwa软件将短reads档案库SRA中获得的冠状病毒二代测序原始fastq文件与NCBI参考基因组进行序列比对，生成BAM文件；

步骤二、对步骤一获得的BAM文件通过bcftools调用和过滤单核苷酸多态性SNP，并用vcf-annotator对SNP进行注释；随后，根据某特定SNP的有无将收集的病毒株进行分组，此步的分组，将用于步骤四sgRNA表达量构建的表达矩阵进行转录组分析，该分析内容详见步骤四。此外，还将根据SNP所指示的碱基变化，在参考基因组的基础上对特定碱基进行更换，生成所谓同义基因组序列。

步骤三、对步骤一中获得的BAM文件进行CIGAR字符串解析及断点识别操作，进行sgRNA鉴定；

步骤四、将步骤三中的sgRNA鉴定结果构建表达矩阵并进行转录组分析，分析过程中，根据SNP有无，分离来源(来自体内实验抑或体外实验)，项目来源(根据项目识别码)进行分组，对不同分组的病毒株的各sgRNA表达量进行统计学比较。

步骤一中，所述冠状病毒包括冠状病毒科Coronaviridae下的所有病毒种。

步骤一中，所述序列比对是指将通过核酸测序技术获得的核酸序列，与参考序列进行比较，从而获得该核酸序列与参考序列碱基的有无以及排列差异的过程，用于筛选存在断点的核酸序列。

所述筛选，是指在比对结果的BAM文件中，查找CIGAR字符串存在字母“H”或“S”的read，并将其信息记录下来的过程，此处“H”指代“hard clip”，即该read存在一端与参考序列对应位置完全不匹配的现象，“S”指代“soft clip”，即该read虽有一端存在与参考序列不一致的现象，但并非完全不匹配，仍有部分碱基与参考序列一致。

步骤二中，所述调用是指使用bcftools读取步骤一生成的BAM文件，对被比对序列与参考序列的碱基差异以及差异碱基的位置进行描述的过程。

步骤二中，所述过滤是指将同一位置占比超过1％的差异碱基的信息保留下来的过程。

步骤二中，所述注释，是指明确筛选出的差异碱基，是否改变该碱基所在位点的密码子以及氨基酸，并记录该改变是由何碱基改变至何碱基，何密码子改变至何密码子，何氨基酸改变至何氨基酸的过程。

步骤二中，所述更换特定碱基指将参考序列上存在差异碱基的位置所对应的碱基替换为差异碱基的过程。

在步骤二生成的同义基因组序列，用于与已完成全基因组测序的其它冠状病毒进行多序列比对，从而构建冠状病毒毒株的进化树，该过程免去了常规的基因组拼接的过程。

步骤三中，所述字符串解析是指，获取BAM文件中被比对序列与参考序列的比对情况的CIGAR字符串，对其进行信息提取，明确被比对序列分别与参考序列完全一致，不完全一致或完全不一致的位点位置，并获取所述三种比对情况在被比对序列上延伸长度的过程。

步骤三中，所述断点识别，是指上述三种比对情况中任意两种同时存在于被比对序列时，提取两种比对情况所在区域的接合处在被比对序列以及参考序列上的位置并记录的过程。

步骤三中，所述sgRNA鉴定，是指根据被比对序列的断点两侧序列在参考序列上的碱基位置，推测该被比对序列原属于该病毒何种基因表达产物的过程。

步骤四中，所述表达矩阵构建是指，计算样品中每一种步骤三所鉴定的sgRNA总数，并计算每一种sgRNA占该样品所有sgRNA的比例，并以样品编号为行，sgRNA类型为列，或sgRNA为行，样品编号为列，列表描述sgRNA在各样品中所占比例的过程。

本发明中使用的冠状病毒转录组数据集主要来自于公开可用的数据库，进行搜索和整理后，在NCBI的SRA数据库中筛选了19个包含原始reads的数据集，所述数据集包含了588个SARS-CoV-2以及相关冠状病毒样本；同时还使用了最近发表样本的数据集^[10]。

本发明所述的冠状病毒转录组鉴定分析方法利用了高精度第二代测序技术产生的序列，可以从单个read中识别转录调控TRS序列。

具体来说，在本发明方法中，原始read首先与各自的参考病毒基因组序列进行比对，所述参考病毒基因组序列包括SARS-CoV-2(GeneBank ID NC_045512.2)、SARS(GeneBank ID NC_004718.3)、MERS(GeneBank ID NC_019843.3)或其他冠状病毒的参考序列，如表1。

表1本发明加入的生物项目信息表

通过比对断点分析，可以推测出特异的sgRNA，识别出跨越5个前导序列(TRS)和更多远端基因组序列之间连接的reads(图2A)。特定的sgRNA的相对丰度等同于基因的相对表达。为确定病毒基因型和病毒起源(例如体内和体外)如何影响病毒基因表达，本发明还构建了热图(图3)。

本发明中的冠状病毒注释器被设计用来描述所有可能的断点。然而，为了辨别真正可靠的sgRNA，本发明去除了罕见的断点和样本间不一致的断点。一个完整的断点包含两个独立的基因组位置(图2A)。此外，本发明还分析了非sgRNA断点，这些序列的5’端没有包含 leader TRS(LTRS)。数据表明，非sgRNA断点非常罕见(通常低于总sgRNA断点的0.05％)，而且存在样本间不一致的现象，由于这些断点仅在一项研究中被发现，因此，本发明专注于研究那些由传统5’端引导转录调节序列(LTRS)以及3’端转录调节序列(BTRS)形成的 sgRNA。

基于上述所述的方法，本发明还提出了一种冠状病毒转录组鉴定分析系统，所述系统包括储存器和处理器；所述存储器上存储有计算机程序，当所述计算机程序被所述处理器执行时，实现本发明冠状病毒转录组鉴定分析方法。

本发明还提供了所述的冠状病毒转录组鉴定分析方法、新发冠状病毒TRS序列的鉴定中的应用。

本发明还提供了sgRNA，其基因序列如SEQ ID NO.1，SEQ ID NO.2所示，该序列为参考基因组中，对应sgRNA序列，其余冠状病毒中对应序列的一致性应在70％以上。

本发明还提供了所述sgRNA在制备鉴定冠状病毒的产品中的应用。

所述冠状病毒包括冠状病毒科Coronaviridae下的所有病毒种。

本发明还提供了一种诊断试剂，所述诊断试剂包括如上所述的sgRNA。

本发明的有益效果包括：本发明所述的冠状病毒转录组鉴定分析方法为冠状病毒生物学提供了新的认知手段，并为未来的治疗方法的发展提供了宝贵的资源。通过本发明的方法进行鉴定分析，结果显示，SARS-CoV、SARS-CoV-2和MERS-CoV都具备通过经典途径生成的核心sgRNA表达谱。此外，在其中几种冠状病毒中发现了编码进化保守结构多肽的新型sgRNA，且这些sgRNA在体外和体内样品中都会表达，这极大地增加了已预测的冠状病毒蛋白数量。此外，两种新发现的多肽，通过序列可预测分别存在IFN反应性和阻断IL17E(IL25)信号的能力，暗示可能与疾病有直接的功能相关性。此外，S蛋白的sgRNAs体内表达水平明显高于体外水平，而核衣壳蛋白的水平则相反，这可能与冠状病毒的传染性和传播性相关。

附图说明

图1为本发明冠状病毒转录组鉴定分析方法的流程图。

图2为研究总览以及sgRNA图谱。图2A为研究流程总览。图2B为SARS-CoV-2经典断点图。图2C为三种冠状病毒开放阅读框注释以及比较。

图3为SARS-CoV-2的sgRNA表达热图以及对应菌株的SNP注释。

图4为SARS-CoV-2和MERS新发现sgRNA及其对应表达产物。图4A为SARS-CoV-2 新发现sgRNA的断点图。图4B、C、D为新发现sgRNA对应的蛋白产物在不同新冠病毒中的保守性分析结果。

图5为新发现sgRNA断点以及对应的转录调控序列。图5A为SARS-CoV-2的pORF2b的两个转录调控序列。图5B为在SARS-CoV-2新发现的M蛋白转录调控序列。图5C为 SARS-CoV-2中新发现tORF7b的转录调控序列。图5D为MERS-CoV中新发现pORF8c的转录调控序列。

图6为体内和体外条件下，各sgRNA表达量的比较。图6A为SARS-CoV-2的sgRNA 在体内与体外sgRNA表达差异的比较，图6B为MERS-CoV的sgRNA在体内与体外sgRNA 表达差异的比较，图6C为冠状病毒S以及N基因的表达以及这些冠状病毒的进化关系。

图7为体内以及体外数据集reads覆盖度图，其中图7A为体外数据，图7B为体内数据。

图8为除SARS-CoV-2外的四种冠状病毒断点图谱，其中图8A为SARS-CoV，图8B为MERS-CoV，图8C为一种来自穿山甲的冠状病毒，图8D为HKU1。

图9为来自穿山甲的冠状病毒的pORF2b与IL-17受体的同源部分的示意图。

图10为SARS-CoV-2和MERS-CoV详细的基因表达图谱。

图11为格列卫和IFN-β治疗下病毒载量的变化。

具体实施方式

结合以下具体实施例和附图，对本发明作进一步的详细说明。实施本发明的过程、条件、实验方法等，除以下专门提及的内容之外，均为本领域的普遍知识和公知常识，本发明没有特别限制内容。

本发明实施例中所有测序数据来源于NCBI Short ReadsArchive(SRA)。一些纳米孔数据集从各自研究的文章中描述的在线存储位置下载^[10]。通过关键词“冠状病毒(coronavirus)”搜索和人工筛选确定生物项目，只保存宏转录组数据，使用包含wget的自定义脚本从SRA 下载原始测序文件，使用SRAtoolkit从下载的SRA文件生成压缩格式的fastq文件。用bwa 对SARS-CoV、SARS-CoV-2或MERS-CoV的参考基因组序列进行初始序列比对后，滤除病毒reads过少的样本。CORONATATOR仅使用第二代测序技术(Illumina)生成的reads，纳米孔数据用于比较。

本发明所提供的冠状病毒注释器包括一系列用于分析冠状病毒RNA-Seq数据的perl和 bash脚本，所述冠状病毒注释器包括3个主要步骤，包括预处理、断点识别、sgRNA读取和分析，详细信息如下：

1、预处理

以SARS-CoV、SARS-CoV-2和MERS-CoV的基因组序列作为参考序列进行比对，生成预处理BAM文件，蝙蝠病毒和穿山甲病毒对应的基因组则以NCBI上的拼接序列作为参考。通过bcftools调用和过滤SNP^[12]，使用vcf-annotator进行SNP注释^[13]。此外，通过筛选出的SNP生成同义基因组序列以供进一步分析。

2、断点标识

从soft-clip和hard-clip比对中识别出断点，这些比对都是部分比对，主要是由重组接合位置的reads引起的，这是由冠状病毒产生sgRNA的机制产生。在这一步中，生成了reads 信息、断点位置、CIGAR字符串以及可能的TRS序列的矩阵。

3、sgRNA读取和分析

经典sgRNA是通过对参考基因组序列的两个断点坐标进行识别和定义的，这些断点是通过从部分比对中提取(一个是初级比对，另一个是次级比对，这样结果会在被比对的一对read 同时比对到参考序列的两个区域时出现，此时比对匹配度高的比对结果称之为初级比对，次之的为次级比对)得到的。为了识别可能的TRS模式，从之前生成的同义基因组序列中提取断点对之间的序列。然后，通过人工比较已知病毒基因起始密码子与其断点之间的距离，鉴定出相应的sgRNA的基因。保留sgRNA计数超过20的生物样本进行进一步分析，在这些样本中，sgRNA按照对应的基因计数并使用总sgRNA计数归一化得到转录组矩阵。

本发明在实施过程中还进行了新开放阅读框的识别：使用Prodigal^[14]预测潜在的ORF，并使用-s参数输出所有潜在的基因。还使用了一个python脚本来标识非常短的ORF。然后，对于多个生物项目支持的sgRNA，计算并排序它们的断点与所有已识别的起始密码子位点之间的距离。将起始距离上游断点最近的ORF被标记，并人工检查以进行验证，通过标识新开放阅读框，可以为新发现的sgRNA提供其可靠性生物学依据，当新发现的sgRNA断点下游存在新开放阅读框，尤其是距离非常接近时，新发现的sgRNA将很有可能表达出来并有生物学功能。

本发明在实施过程中还进行了序列比对和系统发育分析：采用SARS-CoV-2、SARS-CoV、 MERS-CoV的同义基因组序列，以及蝙蝠、穿山甲或其他人类冠状病毒生物样本中sgRNA 大于20个的序列进行系统发育分析。多序列比对采用MAFFT^[15]，最大似然共识树采用 IQ-TREE^[16]构建，进行1000次bootstrap。

本发明在实施过程中还将纳米孔sgRNA比例转化为短reads形式：Kim等人^[10]的数据集包括纳米孔数据和reads数据。两者之间的比值用于将其他纳米孔数据集转换为与本研究中其他数据集可相互比较的比例。

本发明中涉及到的热图(图3)用于显示基因表达谱。使用heatmap.plus包，用ggplot2 包制作sgRNA表达点图和箱形图，比较不同样本来源间基因表达差异，采用t检验和wilcoxon 检验进行统计学分析。

本发明的功能注释中，新肽序列用EMBL在线工具FASTA (https://www.ebi.ac.uk/Tools/sss/fasta/)对UniProtKB/Swiss-Prot数据库进行默认参数比对。使用NCBI CD Blast在线服务鉴定蛋白结构域。

本发明还检验了病毒多肽序列的保守性：为了检测相关病毒种中推测的多肽序列的保守性，本发明建立了一个包含所有相关病毒基因组预测的开放码的参考数据库。使用DC MegaBlast(DisContinuous MegaBlast)搜索种间同源序列。参数设置如下：window_size 0, gapopen 0,gapextend 2,penalty-1,reward 1,num_alignments 1。然后利用MAFFT对一组同源 ORF进行多序列比对(MSA)。然后用CLUSTAO(Clustal Omega)计算MSA结果的矩阵。核苷酸和氨基酸序列也用相同的方法进行了检测。

实施例1证明sgRNA的丰度为冠状病毒ORF提供了不同水平

SARS-CoV-2的很多ORF是基于序列一致性进行注释的，其中一些注释结果在蛋白质组学和测序研究中存在争议^[2,9]。本发明实施例将常见的ORF划分为三类(核心、低支持度和无支持度)(图2B)。识别真正的sgRNA需要多研究和多样本分析，这是由于在文库制备或测序过程中往往会产生独特的人工序列^[17-18]。此外，许多低丰度的非经典sgRNA可能是随机产生的没有特定功能的异常转录本^[10]。因此，只有在多个研究和数据集中存在的sgRNA才能作为真正的sgRNA候选。为了对每个病毒基因进行分类，本发明实施例考虑了一些因素，如 sgRNA相对丰度、TRS保守和可能受起始密码子劫持(start-codon hijacking)影响的核糖体跳读(leaky ribosome scanning)^[19]。

本发明实施例首先通过查看特异性sgRNA测序证据的支持程度验证了SARS-CoV-2、SARS-CoV和MERS-CoV中最常被注释的ORF。对于SARS-CoV-2，去除sgRNA中read数小于20的样本后，剩余34份样本。为了在这些样本中识别核心的sgRNA，并将其归入本发明实施例第一个sgRNA类别，使用加权平均方法按照病毒种类识别了所有可能的sgRNA，并记录了它们的相对丰度。在相对丰度为0.5％时，发现其中8个典型断点(breakpoint)分别对应了8种sgRNA，这些sgRNA包含SARS-CoV-2的8个已被深入研究的ORF：S、E、M、 N、ORF3a、ORF6、ORF7a和ORF8(图2B-C，表2，表3，图8A)。

表2三种冠状病毒注释信息表

表2-1SARS-CoV-2

#a:此处呈现ORF3b在SARS-CoV-2的一个常用注释，另一个注释在括号里显示。

#b:ORF7b在SARS-CoV-2的表达相对丰度相比其他经典sgRNA明显偏低且在SARS-CoV中无同源蛋白。

表2-2SARS-CoV

表2-3MERS-CoV

表3样品sgRNA表达矩阵

这些ORF的sgRNA断点位于起始密码子上游9～162nt之间。N是最丰富的核心sgRNA，占所有样品中被鉴定核心sgRNA的54％。E sgRNA含量最低，为1.5％，是最近蛋白质组学研究中唯一未被鉴定出的核心蛋白^[9,20]。此外，ORF7a、M、ORF3a、S、ORF8和ORF6 的检出率分别为10.6％、8.4％、6.9％、6.1％、5.9和2.7％。这8种核心sgRNA合计占总sgRNA 的70％-100％，具体数值取决于样本类型(体内和体外)、毒株和read覆盖度(图3)。

这8个核心sgRNA除了具有较高的相对丰度外，还有一个相同的TRS，其核心序列为ACGAAC，为这组sgRNA所独有。该核心序列可能是sgRNA形成的必要和充分条件。同样的8个核心sgRNAs以及核心TRS序列也被SARS-CoV所共享(图8A)。此外，除了N (其TRS核心序列为ACGAA)，MERS-CoV的7个核心sgRNA(S、E、M、N、ORF3、 ORF4a和ORF5)(图2C)也使用了这个核心序列。

第二类ORF的sgRNA通常相对丰度较低，且不使用这个核心序列。这类sgRNA包括SARS-CoV-2和SARS-CoV中的ORF7b，SARS中的ORF3b，MERS-CoV中的ORF4b和 ORF8b。对于SARS-CoV-2，E的平均相对丰度为1.5％，是核心sgRNA丰度中最低的，而 ORF7b仅为0.02％。这种sgRNA形成的低丰度或低效率可能是使用非经典TRS造成的。这组sgRNA不像核心sgRNA那样使用保守的核心TRS序列，即它们重组所依赖的序列总是以核心序列为中心偏移几个碱基，并且经常出现前导TRS和3’TRS之间的不匹配的现象。

其余的ORF属于第三类，这些ORF不存在支持的sgRNA。这类ORF可进一步分为两个子类。第一类是虽没有sgRNA支持，但可能被翻译。已有研究表明，一些冠状病毒的 ORF可以通过一种被称为核糖体跳读的机制表达^[19]。SARS-CoV-2的ORF9b就属于这一类。事实上，最近多项蛋白质组学研究表明，ORF9b蛋白产物在SARS-CoV-2中是存在的^[20]。它在SARS-CoV中的同源基因也被命名为ORF9b，也在这一类别。SARS-CoV-2和SARS- CoV的ORF7b在之前的研究中也被提及应属于这一类ORF^[19]，而且ORF7b的起始密码子与其前部的ORF7a之间存在的很长一段序列(SARS-CoV-2中为362bp；SARS-CoV为365bp) 是没有额外的起始密码子的。然而，这些ORF在本发明中仍能形成自己的sgRNA，只是含量很低。

第二类是最可疑的ORF，除无sgRNA支持外，在它们和最接近的sgRNA断点之间由于存在起始密码子，导致其不太可能表达。这一类别包括几个经常被注释的ORF，包括 SARS-CoV-2的ORF3b、ORF9c和ORF10、SARS-CoV的ORF8b。这些ORF和之前的ORF 之间存在多个起始密码子，此外，其相应的sgRNA以及蛋白质组学证据也均不存在^[10,20-21]。事实上，ORF10的存在在最近研究中^[9,11]也确实受到了质疑。上述证据表明，在假定的开放框中进行不含sgRNA或蛋白质组学支持的病毒产物实验是可能存在缺陷的。例如，最近的一项研究生成了SARS-CoV-2中预测的ORF3b截短版本的合成产物，推测SARS-CoV-2中假定的截短版本比SARS版本具有更强的抗IFN活性^[22]。

实施例2SARS-CoV-2、MERS-CoV和SARS-CoV中包含非经典TRS的新型sgRNA的鉴定

在核心sgRNA谱的形成过程中，包含最小核心序列的3’TRS将与先导TRS配对。对于每个特定的核心sgRNA，使用的两个TRS核心序列必须具有相同的长度和序列，尽管不同的sgRNA的TRS区域的长度可能不同(图2C)。本实施例中发现这些典型的SARS-CoV-2 TRS的平均长度为～9.6碱基。SARS-CoV使用了与SARS-CoV-2相同的核心序列，而MERS 则使用另一类具有不同核心序列的6个核苷酸的TRS(图8A，图8B)。

当寻找丰度超过0.2％的sgRNA转录本时，发现了至少在两个独立样本和研究中存在的三种新型sgRNA(图4A)。所有三种新的sgRNAs没有使用核心sgRNA使用的经典TRS 序列。但这些sgRNA的断点支持不连续延伸模型(discontinuous extension model)，因为来自3’的序列在最终转录本的TRS序列中被发现了(图5)。另外，对链特异性RNA文库序列分析确定了这些sgRNA负链形式的存在。这些新型sgRNA在之前的Nanopore测序研究中没有被描述^[9-11]。如前所述，对这些非经典断点序列的分析发现当前导TRS与3’TRS核心序列不互补时，则可能利用核心序列周边序列进行互补，而导致sgRNA的产生(例Firgure 2b,tORF7b和3pORF2b)。本实施例的分析也证实了TRS序列在亲缘关系较远病毒之间有很大的差别，同时也发现经典TRS序列长度在一些冠状病毒中可以超过30碱基(图8D)。

由三个新型TRS产生的三种新型sgRNA，称之为pORF2b，aM和tORF7b。最长的新型sgRNApORF2b位于S基因内，在22501附近有两个备选的TRS，编码一个短肽，其结构域在近缘的冠状病毒中是保守的(图4A，图4B)。第二个sgRNA断点位于M的经典断点下游31碱基即26494号碱基处。该sgRNA也将表达M蛋白，但其5’UTR不同于M的经典sgRNA(图5B)。三种新型sgRNA中最短的sgRNA其断点位于27761，编码截断的 ORF7b(tORF7b)。该截断去除了肽链胞外结构域和组成跨膜结构域的24个氨基酸中的14 个(图4A，图4C)。该sgRNA在体内和体外均有较高水平的表达，可能具有新的功能。

pORF2b翻译产物为36个氨基酸的肽链。PSIPRED预测^[23]其有一胞内蛋白结合结构域和两个短的位于跨膜结构域的alpha螺旋，其中一个alpha螺旋部分位于细胞外(图4A)。pORF2b在两个独立的研究的4个样本中被发现。在来自美国华盛顿州的患者样本(SRX7884411)中，pORF2b的表达量最高，占总sgRNA的11.1％。在同一项目的另一个患者样本(SRX7884409)中，占已鉴定sgRNA的1.2％(图3)。SNP分析显示感染这两名患者的病毒株只差一个核苷酸。来自同一研究的其他5个不同病毒株的患者样本(图3)没有产生pORF2b的sgRNA。这表明pORF2b转录本的水平可能与病毒株之间存在一定的相关性，同时说明pORF2b并非是测序随机产生的假阳性结果。SgRNApORF2b也在另一个独立研究(PRJNA615032)中被鉴定，该研究在两个体外样本中使用了不同于体内研究中鉴定的病毒株(图3)。

在其他Sarbecoviruses中确定了pORF2b的序列保守程度，包括SARS-CoV、HKU3(蝙蝠冠状病毒)、RaTG13(蝙蝠冠状病毒，被认为与SARS-CoV-2直接相关)和感染穿山甲的冠状病毒(SRX7732088)^[24]。四种病毒均有相应的ORF，其中RaTG13同源性最高，核苷酸同源性为91.85％(表4，图4B)。pORF2b以及C端延伸的穿山甲版本与人类IL17RB的配体结合域具有很高的相似性(图9)。

表4sgRNA保守程度信息表

表4-1SARS-CoV-2

表4-2MERS-CoV

第三个新的断点位于ORF7b内的27761位置，编码ORF7b的一个截断型(tORF7b)。在两个包括一个以上病毒株的独立研究中鉴定出了该转录本及其在体内和体外的相对丰度。该转录本最近也在一个被单一病毒株感染的VERO细胞系中被鉴定出来^[10]。这种新型sgRNA在体内和体外都有相对高水平的表达(图3)，而且在两个研究的几个SARS-CoV样本中也存在这种sgRNA的同源基因产物。ORF7b的这一截短版本缺失了胞内结构域和超过一半的跨膜结构域，但保留了亲水的胞外结构域(图4A)。ORF7b作为SARS-CoV-2ORF7b 的同源基因也存在于SARS-CoV病毒颗粒中^[19]。在ORF7b中由tORF7b编码的部分在SARS 中高度保守(表4，图4C)。然而，一项研究中的SARS ORF7b存在45nt缺失(该部分也存在于tORF7b缺失序列中)且导致病毒对β干扰素诱导力减弱，进而提供了复制的优势^[25]。未来的研究将揭示这种新的sgRNA是否编码了一种具有拮抗IFN功能的毒性肽段，同时破坏干扰素反应的启动。

本发明还获得了大量的MERS-CoV体内和体外序列数据集，能够识别大量的非典型sgRNA(图8)。sgRNA(推测为ORF8c或pORF8c)编码一个可翻译成51个氨基酸的肽段。该sgRNA在5个独立的体内和体外研究中被鉴定出来，其丰度从0.03％到1.0％不等。 PSIPRED提示这种新肽具有一个连接了细胞质螺旋结构的跨膜结构域。在其他Merbecovirus 包括HKU4、HKU5和一种Erinacus冠状病毒(Erinacus coronavirus)中确认了该蛋白的保守性。pORF8c在3个样本中均存在，保守性不尽相同，在Merbecovrius中，该蛋白的胞质 N端是最保守的，而在HKU5和Erinaceus中观察到C端延长版本(图4D)。

为了穷尽对新型的sgRNA的寻找，将筛选阈值降低到相对丰度0.01％，同时保持其他标准。该分析在多个研究发现了更多SARS-CoV-2、SARS-CoV和MERS-CoV新型sgRNA (表5)。进一步的测序和未来的实验将确定pORF2b、tORF7B和aM以及其他大量丰度极低的新sgRNA的重要性。

表5本发明新发现的sgRNA信息表

表5-1SARS-CoV-2新型sgRNA

#Kim_2020为一个单样本数据集，Kim等2020.

#*指无重叠序列

表5-2SARS-CoV新型sgRNA

#a:该sgRNA加权比例低于阈值10^-4，但该sgRNA可翻译为tORF7b蛋白产物，该蛋白肽链在 SARS-CoV-2保守且其存在被多项研究支持，因此包含在内。

表5-3MERS-CoV新型sgRNA

实施例3CORONATATOR检测实验诱导的新型pORF8c相对丰度变化

本实施例将验证本发明冠状病毒转录组鉴定分析方法的实验效用，并验证新的sgRNA 在实验刺激时与其它已研究基因类似的反应方式。为此使用了一个实验数据集，本实施例在体外治疗MERS-CoV感染过程中测试了格列卫和IFN-β对宿主基因表达的影响(PRJNA233943&amp；PRJNA233944)(图11)。本实施例分析了这些药物对病毒载量、病毒基因表达和新发现pORF8c表达的影响。最初的分析表明，病毒载量的降低导致单个病毒基因的表达的增加(图11)。重要的是，即使在病毒载量较低的情况下，N与S的比值仍然很高，这表明这个比值不受病毒丰度的影响，而是受体内外环境的影响。此外，格列卫和 IFN-β对病毒基因表达的影响不一致，对不同病毒基因的影响不同。在IFN-β和格列卫的应答中，pORF8c的表达与核心sgRNA，特别是N和E，随着病毒载量呈现相同的趋势。这说明在这一背景下，pORF8c在基因表达方面与一些核心sgRNA具有相同的生物学反应。

实施例4在体内，SARS-CoV-2的S sgRNA的相对丰度升高

在处理数据集时，本发明注意到SARS-CoV-2的测序数据在参考基因组上有两种不同的reads覆盖模式，这表明病毒的reads有两种来源。经过进一步检查，发现这两个模式的来源是体内和体外条件(图7)。前者由细胞外病毒颗粒和存在于BALF(人)、鼻洗液(雪貂)或肺匀浆(MERS)中的受感染宿主细胞组成，而后者由受感染细胞系组成，不受机体或可能存在的先天(如VERO细胞不产生IFN)的抗病毒反应的影响。来自SARS-CoV-2 的样品(主要是BALF)内的病毒序列一般覆盖整个病毒参考基因组，3’端相对于其他区域的差别较小。相比之下，在体外样本中，由于病毒转录过程中会形成sgRNA的嵌套，病毒基因组3’端覆盖度高。

SARS-CoV-2和MERS-CoV是唯二尚在流行的有毒力冠状病毒，且都有体内和体外来源的宏转录组数据集。本发明中分析了SARS-CoV-2和MERS-CoV体内和体外生成的 sgRNA的相对丰度。在比较病毒sgRNA在体内与体外生成的差异时，发现S与N的比值在体内明显高于体外，尤其是SARS-CoV-2(0.04体外vs 0.69体内，p＝0.0012，Wilcoxon ranksum test)(图6A-B，图10)。这种sgRNA水平的显著差异的一个一般解释是，环境压力的差异影响了病毒复制对这些sgRNA需求。选择压力可能促进有利于病毒增殖的基因转录。例如，S蛋白的主要功能围绕宿主细胞的识别和入侵，而N蛋白的功能以调节病毒RNA 促进病毒复制为中心，该过程通过结合病毒基因组3’端、病毒包装信号和TRS介导的^[26-28]。另一种可能同时存在的解释是，体内S/N比值的增加由病毒复制转录复合物(RTC)的内部改变导致，该过程使转录过程更有利于TRS通读，优先生成更长的sgRNA。这一全局的病毒转录调整可能涉及宿主因素的改变，为传染性支气管炎病毒(IBV)，在传染性支气管病毒(infectious bronchitisvirus)，一种γ冠状病毒，其N蛋白被胞内GSK3磷酸化后，招募解旋酶DDX1促进TRS通读从而促进长sgRNA的形成^[29]。在这一前提下，由于N蛋白是由最短的sgRNA产生的，而S蛋白是由最长的sgRNA产生的，故TRS通读将提高S的表达量。未来对病毒颗粒体内和体外的电子显微镜研究将确定SARS-CoV-2中S sgRNA的丰度是否与病毒颗粒表面的S蛋白水平相关。sgRNAs在体内和体外表达显著差异的其他例子包括SARS-CoV-2和MERS-CoV通过多种途径抑制免疫应答时的附属sgRNAs水平总体增加(图6A、B)^[30]。

为了更清楚地了解SARS-CoV-2的sgRNA在体内和体外与其他冠状病毒相比的相对丰度差异，以及确定是否遗漏了其他新的sgRNA，本实施例利用CARONTATOR分析了其他冠状病毒。该分析包括OC43、NL63、HKU1以及与SARS-COV-2高度同源的蝙蝠病毒和穿山甲病毒^[8,24](图6C，表1和图8)。一些数据集没有足够的断点来提供信息。例如， SARS-COV-2同源性最高的蝙蝠病毒RaTG13的分析只得到1个断点，因此在图6C中省略。

在不同冠状病毒中，SARS-COV-2的S sgRNA含量最高，特别是在体内(图6C)。分析表明，该现象不是由于菌株特异性导致的，因为不同菌株的体内样本中都能观测到高水平SsgRNA表达(图3)。高水平的S蛋白表达可能在病毒中发挥作用从而有利于跨越物种屏障从而表现出高传染性。同样，注意到S sgRNA的相对水平与冠状病毒的传染性正相关。 SsgRNAs在体内的病毒感染性和水平如下：SARS-CoV-2>HKU1>MERS^[31]。然而，单是S 蛋白水平不足以导致高水平的传播能力，因为S蛋白的稳定性、受体的亲和力^[32]和病毒颗粒的稳定性^[33]等因素也将影响病毒的传播。

实施例5RTC突变逆转体内外N和S sgRNA的表达

本发明中还实现了对病毒RTC组分的突变的观测，会改变S到N的表达谱，特别是，Kim等人^[10]使用的病毒株在RTC组分nsp3(编码papin蛋白酶，可结合N和M蛋白)中有一个独特的非同义突变(图3)。这个病毒毒株的体外转录组显示S与N的比例大幅增加，类似于体内条件下的病毒表达谱(图3)。一个体内实验的病毒毒株(SRX7852918)也有两个nsp3的非同义突变，除此之外还有nsp6和nsp12的突变，该病毒的转录谱与体外条件的表达谱，其S/N比值降低。

在两株基因表达改变的病毒中存在nsp3突变也引人深思。据报道，nsp3可结合病毒基因组3’末端TRS，该部位是病毒全局RNA包装信号同时编码了N和M蛋白^[27,34-35]。此外，N蛋白的磷酸化会改变其构象，从而更倾向结合病毒RNA，就像前述在IBV中，促进 TRS通读那样^[2,36]会促进长链sgRNA生成。这一现象初步表明，nsp3内的突变通过影响病毒全局构造来影响各sgRNA的相对丰度，并可能以上述IBV类似的机制作用。此外，在上述机制影响下，体内和体外N sgRNA相对丰度的改变，可能反馈其与nsp3的相互作用，并影响nsp3突变在体内和体外的功能(图6A-B)。

迄今为止，大量SARS-CoV-2序列数据主要用于新发毒株的分型和跟踪。虽然这很重要，但是会导致对宝贵的资料利用不足。通过开发CORONATATOR生物信息分析流程，在病毒株的鉴定，描述冠状病毒sgRNA表达和发现通过非经典TRS配对机制生成的保守 sgRNA新未知功能(图5)方面迈出了重要一步。这些新发现的蛋白质功能预测仍在进行中。初步发现穿山甲病毒SARS-CoV-2pORF2b同源物与人IL17RB s配体结合域具有广泛的相似性(图9)。冠状病毒通过该sgRNA可能产生一种短肽，理论上可以破坏IL17B和 IL17E(IL25)信号，后两者通常与促进或抑制特定情况下的炎症反应有关。未来需要进行蛋白质组学研究和/或核糖体测序研究，以验证本发明所鉴定的新型sgRNA编码的蛋白质产物的存在。

本发明的分析还表明，不同的SARS-CoV-2病毒株表达的各sgRNA的水平不同(表S3，图3)，尤其是新发现的sgRNA。本发明结果强调，只有对病毒处于选择性压力下的患者样本进行深度测序，才能真正了解sgRNA表达方面的病毒致病机制。这需要对样本进行深入了解，根据不同病毒株在新冠肺炎不同阶段进行彻底的测序和分析，这也将推动真正个性化治疗的产生。

尽管其他人畜共患的病毒可能在基因或基因组水平上与SARS-CoV-2具有广泛的序列相似性，但仅具有序列相似性不足以产生人类易感的具备致病性的病毒。在讨论人畜共患病毒起源时，人们一般不考虑病毒基因的特异性表达水平，如S蛋白，但这可能对跨越物种障碍很重要。比方说，由于未取样的人畜共患病毒的数量十分巨大，能够跨越物种屏障的S蛋白很可能早已存在，但尚未以足够高的水平表达以实现可持续的人际传播。低水平的S蛋白表达可以让该病毒由蝙蝠向人类零星传播，但在人类无法将传播过程持续下去，加之人类生活环境卫生状况较好，人与人之间的传播效率会很低。与此照应的是，生活在蝙蝠洞穴附近的人会携带病毒特异性抗体，但却从未经历过严重疾病^[37]。

本发明对宏转录组数据集的分析确定了许多来源的RNA，如宿主RNA和微生物RNA(尽管不是最好的分子捕获方式)。在目前还不清楚为什么有些人会死于SARS-CoV-2感染，而有些人则不会，这些有价值的序列不应被浪费，如果这些数据集能够携带更多的临床信息，就可以变得更有使用价值。理论上大多数关于SARS-CoV-2的GISAID条目都应当有支持它的宏转录组数据集。然而，目前GISAID中描述病毒基因组序列和毒株的条目远远超过在SRA中发现的原始reads数据，共享原始reads信息将极大地帮助科研人员研究这种病毒，并最终遏制它。

参考文献

1.Knipe,D.M.,Howley,P.M.,2013.Fields virology,sixth ed.Wolters KluwerHealth/Lippincott Williams &Wilkins,Philadelphia.

2.Wu,F.,Zhao,S.,Yu,B.,Chen,Y.M.,Wang,W.,Song,Z.G.,Hu,Y.,Tao,Z.W.,Tian,J.H.,Pei,Y.Y., Yuan,M.L.,Zhang,Y.L.,Dai,F.H.,Liu,Y.,Wang,Q.M.,Zheng,J.J.,Xu,L.,Holmes,E.C.,Zhang,Y.Z.,2020. A new coronavirus associatedwithhuman respiratory disease in China.Nature 579,265-+.

3.Brian,D.A.,Baric,R.S.,2005.Coronavirus genome structure andreplication.Curr Top Microbiol 287, 1-30.

4.Yount,B.,Curtis,K.M.,Fritz,E.A.,Hensley,L.E.,Jahrling,P.B.,Prentice,E.,Denison,M.R.,Geisbert, T.W.,Baric,R.S.,2003.Reverse genetics witha full-length infectious cDNA of severe acute respiratory syndromecoronavirus.P Natl Acad Sci USA 100,12995-13000.

5.Sola,I.,Almazan,F.,Zuniga,S.,Enjuanes,L.,2015.Continuous andDiscontinuous RNA Synthesis in Coronaviruses.Annu RevVirol 2,265-288.

6.Assiri,A.,McGeer,A.,Perl,T.M.,Price,C.S.,Al Rabeeah,A.A.,Cummings,D.A.T.,Alabdullatif,Z.N., Assad,M.,Almulhim,A.,Makhdoom,H.,Madani,H.,Alhakeem,R.,Al-Tawfiq,J.A.,Cotten,M.,Watson,S.J., Kellam,P.,Zumla,A.I.,Memish,Z.A.,Team,K.M.-C.I.,2013.Hospital Outbreak ofMiddle East RespiratorySyndrome Coronavirus.New Engl J Med 369,407-416.

7.Peiris,J.S.M.,Yuen,K.Y.,Osterhaus,A.D.M.E.,Stohr,K.,2003.Currentconcepts:The severe acute respiratory syndrome.New Engl JMed 349,2431-2441.

8.Zhou,P.,Yang,X.L.,Wang,X.G.,Hu,B.,Zhang,L.,Zhang,W.,Si,H.R.,Zhu,Y.,Li,B.,Huang,C.L., Chen,H.D.,Chen,J.,Luo,Y.,Guo,H.,Jiang,R.D.,Liu,M.Q.,Chen,Y.,Shen,X.R.,Wang,X.,Zheng,X.S., Zhao,K.,Chen,Q.J.,Deng,F.,Liu,L.L.,Yan,B.,Zhan,F.X.,Wang,Y.Y.,Xiao,G.F.,Shi,Z.L.,2020.A pneumonia outbreak associatedwith a new coronavirus of probable bat origin.Nature 579,270-273.

9.Davidson,A.D.,Williamson,M.K.,Lewis,S.,Shoemark,D.,Carroll,M.W.,Heesom,K.,Zambon,M., Ellis,J.,Lewis,P.A.,Hiscox,J.A.,Matthews,D.A.,2020.Characterisation of the transcriptome and proteome of SARS-CoV-2 usingdirect RNA sequencing and tandem mass spectrometry reveals evidence for acell passage induced in-frame deletion in the spike glycoprotein that removesthe furin-like cleavage site.bioRxiv.

10.Kim,D.,Lee,J.-Y.,Yang,J.-S.,Kim,J.W.,Kim,V.N.,Chang,H.,2020.TheArchitecture of SARS-CoV-2 Transcriptome.Cell 181,914-921.e910.

11.Taiaroa,G.,Rawlinson,D.,Featherstone,L.,Pitt,M.,Caly,L.,Druce,J.,Purcell,D.,Harty,L.,Tran,T., Roberts,J.,Catton,M.,Williamson,D.,Coin,L.,Duchene,S.,2020.Direct RNA sequencing and early evolution ofSARS-CoV-2.bioRxiv.

12.Li,H.,2011.A statistical framework for SNP calling,mutationdiscovery,association mapping and population genetical parameter estimationfrom sequencing data.Bioinformatics 27,2987-2993.

13.Petit,R.,2014.vcf-annotator,https://github.com/rpetit3/vcf-annotator.

14.Hyatt,D.,Chen,G.L.,Locascio,P.F.,Land,M.L.,Larimer,F.W.,Hauser,L.J.,2010.Prodigal: prokaryotic gene recognition and translation initiationsite identification.BMC Bioinformatics 11,119.

15.Katoh,K.,Misawa,K.,Kuma,K.,Miyata,T.,2002.MAFFT:a novel method forrapid multiple sequence alignment based on fast Fourier transform.NucleicAcids Res 30,3059-3066.

16.Minh,B.Q.,Schmidt,H.A.,Chernomor,O.,Schrempf,D.,Woodhams,M.D.,vonHaeseler,A.,Lanfear, R.,2020.IQ-TREE 2:New Models and Efficient Methods forPhylogenetic Inference in the Genomic Era.Mol Biol Evol 37,1530-1534.

17.Lebrigand,K.,Magnone,V.,Barbry,P.,Waldmann,R.,2020.High throughputerror corrected Nanopore single cell transcriptome sequencing.Nat Commun 11,4025.

18.Peng,Q.,Vijaya Satya,R.,Lewis,M.,Randad,P.,Wang,Y.,2015.Reducingamplification artifacts in high multiplex amplicon sequencing by usingmolecular barcodes.BMC Genomics 16,589.

19.Schaecher,S.R.,Mackenzie,J.M.,Pekosz,A.,2007.The ORF7b protein ofsevere acute respiratory syndrome coronavirus(SARS-CoV)is expressed in virus-infected cells and incorporated into SARS-CoV particles.J Virol 81,718-731.

20.Bojkova,D.,Klann,K.,Koch,B.,Widera,M.,Krause,D.,Ciesek,S.,Cinatl,J.,Münch,C.,2020. Proteomics ofSARS-CoV-2-infected host cells reveals therapytargets.Nature.

21.Gordon,D.E.,Jang,G.M.,Bouhaddou,M.,Xu,J.,Obernier,K.,White,K.M.,O’Meara,M.J.,Rezelj, V.V.,Guo,J.Z.,Swaney,D.L.,Tummino,T.A.,Huettenhain,R.,Kaake,R.M.,Richards,A.L.,Tutuncuoglu,B., Foussard,H.,Batra,J.,Haas,K.,Modak,M.,Kim,M.,Haas,P.,Polacco,B.J.,Braberg,H.,Fabius,J.M., Eckhardt,M.,Soucheray,M.,Bennett,M.J.,Cakir,M.,McGregor,M.J.,Li,Q.,Meyer,B.,Roesch,F.,Vallet,T.,Mac Kain,A.,Miorin,L.,Moreno,E.,Naing,Z.Z.C.,Zhou,Y.,Peng,S.,Shi,Y.,Zhang,Z.,Shen,W.,Kirby, I.T.,Melnyk,J.E.,Chorba,J.S.,Lou,K.,Dai,S.A.,Barrio-Hernandez,I.,Memon,D.,Hernandez-Armenta,C., Lyu,J.,Mathy,C.J.P.,Perica,T.,Pilla,K.B.,Ganesan,S.J.,Saltzberg,D.J.,Rakesh,R.,Liu,X.,Rosenthal,S.B., Calviello,L.,Venkataramanan,S.,Liboy-Lugo,J.,Lin,Y.,Huang,X.-P.,Liu,Y.,Wankowicz,S.A.,Bohn,M., Safari,M.,Ugur,F.S.,Koh,C.,Savar,N.S.,Tran,Q.D.,Shengjuler,D.,Fletcher,S.J.,O’Neal,M.C.,Cai,Y., Chang,J.C.J.,Broadhurst,D.J.,Klippsten,S.,Sharp,P.P.,Wenzell,N.A.,Kuzuoglu,D.,Wang,H.-Y.,Trenker, R.,Young,J.M.,Cavero,D.A.,Hiatt,J.,Roth,T.L.,Rathore,U.,Subramanian,A.,Noack,J.,Hubert,M.,Stroud,R.M.,Frankel,A.D.,Rosenberg,O.S.,Verba,K.A.,Agard,D.A.,Ott,M.,Emerman,M.,Jura,N.,von Zastrow, M.,Verdin,E.,Ashworth,A.,Schwartz,O.,d’Enfert,C.,Mukherjee,S.,Jacobson,M.,Malik,H.S.,Fujimori, D.G.,Ideker,T.,Craik,C.S.,Floor,S.N.,Fraser,J.S.,Gross,J.D.,Sali,A.,Roth,B.L.,Ruggero,D.,Taunton,J.,Kortemme,T.,Beltrao,P.,Vignuzzi,M.,García-Sastre,A.,Shokat,K.M.,Shoichet,B.K.,Krogan,N.J.,2020.A SARS-CoV-2 protein interaction map reveals targetsfor drug repurposing.Nature.

22.Konno,Y.,Kimura,I.,Uriu,K.,Fukushi,M.,Irie,T.,Koyanagi,Y.,Nakagawa,S.,Sato,K.,2020. SARS-CoV-2 ORF3b is a potent interferon antagonistwhose activity is further increased by a naturally occurring elongationvariant.bioRxiv.

23.Buchan,D.W.A.,Jones,D.T.,2019.The PSIPRED Protein AnalysisWorkbench:20 years on.Nucleic Acids Res 47,W402-W407.

24.Lam,T.T.-Y.,Shum,M.H.-H.,Zhu,H.-C.,Tong,Y.-G.,Ni,X.-B.,Liao,Y.-S.,Wei,W.,Cheung, W.Y.-M.,Li,W.-J.,Li,L.-F.,Leung,G.M.,Holmes,E.C.,Hu,Y.-L.,Guan,Y.,2020.Identifying SARS-CoV-2 related coronaviruses in Malayanpangolins.Nature.

25.Pfefferle,S.,

V.,Ditt,V.,Grywna,K.,Mühlberger,E.,Drosten,C.,2009.Reverse genetic characterization ofthe natural genomic deletion in SARS-Coronavirus strain Frankfurt-1 open reading frame 7b reveals an attenuatingfunction ofthe 7b protein in-vitro and in-vivo.Virol J 6,131-131.

26.Fan,H.,Ooi,A.,Tan,Y.W.,Wang,S.,Fang,S.,Liu,D.X.,Lescar,J.,2005.Thenucleocapsid protein of coronavirus infectious bronchitis virus:crystalstructure of its N-terminal domain and multimerization properties.Structure13,1859-1868.

27.Liang,Y.,Wang,M.-L.,Chien,C.-S.,Yarmishyn,A.A.,Yang,Y.-P.,Lai,W.-Y.,Luo,Y.-H.,Lin,Y.-T., Chen,Y.-J.,Chang,P.-C.,Chiou,S.-H.,2020.Highlight ofImmune Pathogenic Response and Hematopathologic Effect in SARS-CoV,MERS-CoV,and SARS-Cov-2 Infection.Frontiers in Immunology 11.

28.Molenkamp,R.,Spaan,W.J.,1997.Identification of a specificinteraction between the coronavirus mouse hepatitis virus A59nucleocapsidprotein andpackaging signal.Virology 239,78-86.

29.Wu,C.-H.,Chen,P.-J.,Yeh,S.-H.,2014.Nucleocapsid Phosphorylationand RNA Helicase DDX1 Recruitment Enables Coronavirus Transition fromDiscontinuous to Continuous Transcription.Cell Host& Microbe 16,462-472.

30.Canton,J.,Fehr,A.R.,Fernandez-Delgado,R.,Gutierrez-Alvarez,F.J.,Sanchez-Aparicio,M.T., Garcia-Sastre,A.,Perlman,S.,Enjuanes,L.,Sola,I.,2018.MERS-CoV 4b protein interferes with the NF-kappaB-dependent innateimmune response during infection.PLoS Pathog 14,e1006838.

31.Kissler,S.M.,Tedijanto,C.,Goldstein,E.,Grad,Y.H.,Lipsitch,M.,2020.Projecting the transmission dynamics of SARS-CoV-2through thepostpandemic period.Science 368,860-868.

32.Wrobel,A.G.,Benton,D.J.,Xu,P.,Roustan,C.,Martin,S.R.,Rosenthal,P.B.,Skehel,J.J.,Gamblin, S.J.,2020.SARS-CoV-2and bat RaTG13 spikeglycoprotein structures inform on virus evolution and furin-cleavageeffects.Nat Struct Mol Biol 27,763-767.

33.Aboubakr,H.A.,Sharafeldin,T.A.,Goyal,S.M.,2020.Stability of SARS-CoV-2and other coronaviruses in the environment and on common touch surfacesand the influence of climatic conditions:A review.Transbound Emerg Dis.

34.Hurst,K.R.,Koetzner,C.A.,Masters,P.S.,2013.Characterization of acritical interaction between the coronavirus nucleocapsid protein andnonstructural protein 3of the viral replicase-transcriptase complex.J Virol87,9159-9172.

35.Lei,J.,Kusov,Y.,Hilgenfeld,R.,2018.Nsp3 of coronaviruses:Structures and functions of a large multi-domain protein.Antiviral Research149,58-74.

36.Chang,C.-k.,Hou,M.-H.,Chang,C.-F.,Hsiao,C.-D.,Huang,T.-h.,2014.TheSARS coronavirus nucleocapsid protein–Forms and functions.Antiviral Research103,39-50.

37.Wang,N.,Li,S.-Y.,Yang,X.-L.,Huang,H.-M.,Zhang,Y.-J.,Guo,H.,Luo,C.-M.,Miller,M.,Zhu,G., Chmura,A.A.,Hagan,E.,Zhou,J.-H.,Zhang,Y.-Z.,Wang,L.-F.,Daszak,P.,Shi,Z.-L.,2018.Serological Evidence of Bat SARS-Related CoronavirusInfection in Humans,China.Virol Sin 33,104-107.

本发明的保护内容不局限于以上实施例。在不背离本发明构思的精神和范围下，本领域技术人员能够想到的变化和优点都被包括在本发明中，并且以所附的权利要求书为保护范围。

SEQUENCE LISTING

<110> 上海市免疫学研究所

<120> 一种冠状病毒转录组鉴定分析方法及分析系统

<160> 2

<170> PatentIn version 3.5

<210> 1

<211> 111

<212> RNA

<213> pORF2b

<400> 1

augcuuggaa caggaagaga aucagcaacu guguugcuga uuauucuguc cuauauaauu 60

ccgcaucauu uuccacuuuu aaguguuaug gagugucucc uacuaaauua a 111

<210> 2

<211> 63

<212> RNA

<213> tORF7b

<400> 2

augcuuauua ucuuuugguu cucacuugaa cugcaagauc auaaugaaac uugucacgcc 60

uaa 63

Claims (10)

1.一种冠状病毒转录组鉴定分析方法，其特征在于，所述方法包括如下步骤：

步骤一、使用bwa软件将短reads档案库SRA中获得的冠状病毒二代测序原始fastq文件与NCBI参考基因组进行序列比对，生成BAM文件；

步骤二、对步骤一获得的BAM文件通过bcftools调用和过滤单核苷酸多态性SNP，并用vcf-annotator对SNP进行注释；随后，根据某特定SNP的有无将收集的病毒株进行分组，用于后续转录组分析；同时，根据SNP所指示的碱基变化，在参考基因组的基础上对特定碱基进行更换，生成同义基因组序列；

步骤三、对所述步骤一中获得的BAM文件进行CIGAR字符串解析及断点识别操作，进行sgRNA鉴定；

步骤四、将所述步骤三中的sgRNA鉴定结果构建表达矩阵并进行转录组分析，分析过程中，根据SNP有无，分离来源，项目来源进行分组，对不同分组的病毒株的各sgRNA表达量进行统计学比较。

2.如权利要求1所述的方法，其特征在于，步骤一中，所述冠状病毒包括冠状病毒科Coronaviridae下的所有病毒种。

3.如权利要求1所述的方法，其特征在于，步骤一中，所述序列比对是指将通过核酸测序技术获得的核酸序列，与参考序列进行比较，从而获得该核酸序列与所述参考序列碱基的有无以及排列差异的过程，用于筛选存在断点的核酸序列。

4.如权利要求3所述的方法，其特征在于，所述筛选是指在比对结果的BAM文件中，查找CIGAR字符串存在字母“H”或“S”的read，并将其信息记录下来的过程；所述字母“H”指代“hard clip”，即该read存在一端与参考序列对应位置完全不匹配的现象；所述字母“S”指代“soft clip”，即该read虽有一端存在与参考序列不一致的现象，但并非完全不匹配，仍有部分碱基与参考序列一致。

5.如权利要求1所述的方法，其特征在于，步骤二中，所述调用是指使用bcftools读取所述步骤一生成的BAM文件，对被比对序列与参考序列的碱基差异以及差异碱基的位置进行描述的过程；所述过滤是指将同一位置占比超过1％的差异碱基的信息保留下来的过程；所述注释，是指明确筛选出的差异碱基，记录该碱基所在位点的密码子以及氨基酸是否改变，和/或如何改变；所述更换特定碱基指将参考序列上存在差异碱基的位置所对应的碱基替换为差异碱基的过程。

6.如权利要求1所述的方法，其特征在于，在步骤二生成的同义基因组序列，用于与已完成全基因组测序的其它冠状病毒进行多序列比对，从而构建冠状病毒毒株的进化树。

7.如权利要求1所述的方法，其特征在于，步骤三中，所述字符串解析是指获取BAM文件中被比对序列与参考序列的比对情况的CIGAR字符串，对其进行信息提取，明确被比对序列分别与参考序列完全一致，不完全一致或完全不一致的位点位置，并获取所述三种比对情况在被比对序列上延伸长度的过程；所述断点识别，是指上述三种比对情况中任意两种同时存在于被比对序列时，提取两种比对情况所在区域的接合处在被比对序列以及参考序列上的位置并记录的过程；所述sgRNA鉴定，是指根据被比对序列的断点两侧序列在参考序列上的碱基位置，推测该被比对序列原属于该病毒何种基因表达产物的过程。

8.如权利要求1所述的方法，其特征在于，步骤四中，所述构建表达矩阵是指计算样品中每一种所述步骤三所鉴定的sgRNA总数，并计算每一种sgRNA占该样品所有sgRNA的比例，并以样品编号为行，sgRNA类型为列，或sgRNA为行，样品编号为列，列表描述sgRNA在各样品中所占比例的过程。

9.一种冠状病毒转录组鉴定分析系统，其特征在于，所述系统包括储存器和处理器；所述存储器上存储有计算机程序，当所述计算机程序被所述处理器执行时，实现如权利要求1-8之任一项所述的方法。

10.如权利要求1-8之任一项所述的方法、或如权利要求9所述的系统在冠状病毒转录组鉴定分析中的应用。

Images