拿到金葡菌全基因组测序数据后该怎么分析这个问题几乎每个刚接触生物信息学的微生物研究者都会卡住。测序公司返回来一堆raw data看起来就是几GB的fastq文件而你要从中得到结论——这是什么序列型有没有耐药基因有没有毒力因子它在暴发溯源中和其它菌株是什么关系如果脑子里面没有一张清晰的分析路径图很容易陷在“装了软件但不知道下一步干嘛”的状态里。这篇文章就围绕金黄色葡萄球菌以下简称金葡菌全基因组测序数据的完整分析流程展开从拿到fastq开始到最终拿到分型、耐药、毒力、溯源结论为止把每一步怎么做、为什么这么做、有哪些坑全部串起来。无论你是刚开始接触WGS数据分析的研究生还是已经跑过几条流程但想理顺原理的从业者这篇文章都值得收藏。1. 数据质控与组装分析链路的第一道关卡1.1 原始数据质控的几个关键指标测序公司交付的数据一般是双端fastq文件命名常见为Sample_R1.fastq.gz和Sample_R2.fastq.gz。拿到数据的第一件事不是急着组装而是做质量评估和过滤。这一步直接决定后面所有分析的可信度省不掉。我习惯用fastp做质控一条命令能同时完成质量过滤、接头去除、长度过滤和统计报告生成。相比之下传统的FastQC Trimmomatic流程也能做但需要手动串联两个软件速度也慢不少。fastp的好处是自动识别并去除接头序列而且会对双端reads做paired-end overlap分析能顺带修正一部分测序错误。fastp -i Sample_R1.fastq.gz -I Sample_R2.fastq.gz \ -o clean_R1.fastq.gz -O clean_R2.fastq.gz \ --detect_adapter_for_pe \ --thread 16 \ --html Sample_fastp.html跑完看几个核心指标Q20/Q30比例分别代表碱基错误率低于1%和0.1%的碱基占比。金葡菌基因组GC含量约33%属于AT-rich基因组如果测序平台是Illumina NovaSeq一般Q30不低于90%。clean reads数过滤后至少保留原始reads数的80%以上低于这个比例要检查是否测序质量本身有问题。插入片段长度分布正常应该在300-500bp左右。如果双端reads的overlap比例异常高说明插入片段太短对后续组装连续性会有影响。还有一个容易被忽视的指标是duplication rate。如果是细菌基因组测序文库构建时DNA投入量不足容易导致PCR重复率高。fastp默认不自动去重对细菌基因组这种低复杂度样本我一般会在参数里加上--dedup否则重复read会显著影响后续的组装完整性和变异检测准确性。注意质控报告一定要保存好。写文章或报告时审稿人大概率会要求提供测序数据质量参数回头再补跑非常麻烦。1.2 组装工具选择与参数调整质控之后进入组装环节。金葡菌基因组大小约2.8-3.0Mb属于比较小的基因组Illumina双端150bp reads的组装难度不高主流工具都能拿到不错的结果。我用得最多的是SPAdes它基于de Bruijn graph算法对细菌基因组的效果非常稳定。spades.py -1 clean_R1.fastq.gz -2 clean_R2.fastq.gz \ --careful \ -k 21,33,55,77 \ --cov-cutoff auto \ -o spades_output这里几个参数值得解释一下--carefulSPAdes官方建议在组装完成后进行mismatch correction虽然会显著增加运行时间但对提高碱基准确率有帮助。金葡菌基因组不大增加的时间完全可以接受所以我每次都开。-kk-mer长度列表。细菌基因组一般用21到77这几个值就能覆盖如果测序读长是151bp最大k-mer不要超过读长减10否则k-mer无法有效覆盖reads末端区域。--cov-cutoff auto自动过滤低覆盖度的contig。金葡菌如果混入了少量人类宿主DNA或环境DNA污染auto模式能够帮助过滤掉这些低覆盖度的contig。组装完成后用QUAST评估结果quast.py contigs.fasta -r reference.fasta -o quast_output金葡菌组装的核心指标我一般这么判断contig数在50-150之间属于正常NovaSeq数据通常更少N50在10万以上说明组装结果不错完整度检查用CheckM或BUSCO确认基因组完整性在99%以上。如果contig数超过300要怀疑是不是有质粒、前噬菌体等重复区域把组装打碎了——这个后面在可移动元件分析里会详细展开。整个分析链路中组装质量是最容易被事后追溯验证的一环。之前有人送测的样本混了两种不同菌株组装出来contig数异常多N50很低但质控报告又完全正常。后来用Mash检查污染才发现样本本身就是混合的。所以组装指标异常时不要急着调参数先怀疑数据本身。2. 菌株鉴定与分子分型确定你的菌株身份2.1 物种确认别拿到数据就默认是金葡菌很多人拿到数据后跳过物种鉴定直接开始分析这是一个隐患。实验室分离的“金葡菌”不一定真是金黄色葡萄球菌——有可能是表皮葡萄球菌、溶血葡萄球菌等CoNS凝固酶阴性葡萄球菌或者干脆是混合污染。物种鉴定最快的方法是Kraken2配合标准数据库kraken2 --db standard --threads 16 \ --paired clean_R1.fastq.gz clean_R2.fastq.gz \ --report Sample_kraken_report.txt正常金葡菌样本的reads分类结果中Staphylococcus aureus的比例应该超过95%。如果这个比例偏低需要重新审视菌株纯度。更精确的做法是用组装后的contigs跑ANI平均核苷酸一致性和参考基因组比对。金葡菌种内ANI值一般在99%以上而它与表皮葡萄球菌的ANI通常在75%-80%区间差异非常明显。2.2 MLST分型和spa分型的原理与实操物种确认无误后第一个拿到的分子分型结果通常是MLST多位点序列分型。金葡菌的MLST方案基于7个管家基因arcC、aroE、glpF、gmk、pta、tpi、yqiL的等位基因编号每个菌株得到一个7位数的数字组合即ST型。用mlst工具可以直接从contigs文件跑mlst contigs.fasta输出结果示例contigs.fasta staphylococcus_aureus 398 1 3 1 1 620 1 1其中第二个字段表示物种第三个字段是ST型后面7个数字是7个管家基因的等位基因编号。spa分型针对的是金黄色葡萄球菌特有的spa基因protein A基因的X区该区域由一系列24bp的重复单元构成不同菌株的重复单元排列组合不同形成不同的spa型别。spaTyper可以从fastq或组装结果中直接分型spaTyper --contigs contigs.fastaMLST分型和spa分型在暴发调查中扮演的角色不太一样。MLST分辨率中等适合描述菌株的大背景——比如某个地区流行的是ST239还是ST5。spa分型分辨率略高一些在金葡菌院内感染暴发调查中spa型别相同的分离株往往提示可能存在克隆传播。但要注意spa分型只分析单个位点如果spa基因发生重组或缺失突变可能出现无法分型的情况。之前遇到一个菌株spa重复区发生了部分缺失spaTyper报了一个new allele后来用Sanger测序确认是spa基因的自然变异。3. 耐药基因与毒力因子鉴定挖掘表型背后的遗传基础3.1 耐药基因鉴定从数据库到临床意义金葡菌的耐药问题集中在MRSA耐甲氧西林金黄色葡萄球菌其分子基础是mecA基因或mecC等变异体编码的PBP2a蛋白导致β-内酰胺类抗生素失效。另外还有vanA基因介导的万古霉素耐药、erm基因介导的大环内酯类耐药、fosB基因介导的磷霉素耐药等。用ABRicate可以批量扫描耐药基因abricate --db resfinder contigs.fasta sample_resfinder.txt abricate --db card contigs.fasta sample_card.txtABRicate内置了多个数据库ResFinder、CARD、ARG-ANNOT等。实测下来CARD数据库对点突变型耐药基因的覆盖更全比如gyrA的喹诺酮耐药突变ResFinder对获得性耐药基因的检索更准。两个数据库的结果可以互相印证。有一个容易翻车的地方是mecA基因检出只代表存在耐药决定簇不代表体外药敏试验一定耐药。mecA的表达受mecI-mecR1调控系统影响有些菌株携带mecA但表达量很低表现为隐匿型MRSAoxacillin敏感。反过来也一样mecA阴性但苯唑西林耐药的情况也可能出现——这种表型通常由其他机制比如blaZ过表达或PBP突变介导。所以基因型报告和药敏结果不一致时不要急着质疑测序数据先把两边的证据再核对一遍。3.2 毒力因子鉴定金葡菌的“武器库”扫描金葡菌的毒力因子构成非常复杂主要包括黏附素如fnbA、fnbB、clfA、毒素如PVL杀白细胞素、TSST-1中毒性休克毒素、肠毒素sea-see等、免疫逃逸因子如spa、cIfB、以及生物膜形成相关基因如icaADBC操纵子。同样用ABRicate扫描VFDB数据库就行abricate --db vfdb contigs.fasta sample_vfdb.txt对于金葡菌有几个毒力因子的临床意义特别值得关注PVLlukF-PV/lukS-PV杀白细胞素社区获得性MRSACA-MRSA的重要标志。如果检测到PVL阳性基本可以判断菌株属于CA-MRSA谱系。TSST-1tst基因导致中毒性休克综合征的超抗原毒素。tst基因通常由前噬菌体携带散在分布于金葡菌基因组中。肠毒素基因簇sea-see食物中毒的重要致病因子金葡菌食物中毒暴发调查中必查。ABRicate默认的覆盖度阈值是80%有时候可能会漏掉部分截短的毒力基因。建议扫描完看一眼原始比对结果特别是那些“partial positive”的基因不要直接按二元变量处理要结合比对覆盖度判断这个基因是否完整存在。3.3 表型与基因型关联分析时的注意事项耐药基因检测结果要跟表型药敏数据做关联分析时有几个逻辑陷阱有些耐药基因不是组成型表达比如诱导型克林霉素耐药iMLSB基因检测只能告诉你“有这个基因”无法判断表达模式。染色体介导的耐药位点比如gyrA、grlA的QRDR区域点突变在ABRicate默认参数下可能检测不到因为这类突变不是基因的“存在/缺失”问题而是点突变问题。需要用PointFinder或其他工具做专门的点突变筛查。基因组测序检出的耐药基因如果与药敏结果明显矛盾考虑是否为杂合菌群或者样本污染必要时回到原始平板重新挑取单克隆做验证。4. 系统发育分析与暴发溯源SNP层级的分辨率4.1 为什么MLST不够用SNP才是金标准在院内感染暴发调查中同一个ST型别甚至同一个spa型别的菌株并不代表它们一定来自同一个传播链。MLST和spa分型的分辨率不足以区分短期暴发中的菌株差异。这时候需要做核心基因组SNP分析cgSNP——比较菌株间在核心基因组所有菌株共有基因中的单核苷酸多态性差异。不同菌株之间SNP差异数量的流行病学解释大致有个经验范围具体取决于菌株谱系和暴发时间跨度SNP差异范围流行病学推断参考0-12个SNP高度提示近期克隆传播12-30个SNP不排除传播关系结合流调判断30-100个SNP可能是较长时间或间接传播链100个SNP倾向于流行背景多样性非直接传播这个阈值不是写在教科书里的绝对标准不同研究机构的cutoff取值不完全一致但0-12的区间被广泛接受为“同一个暴发克隆”的参考阈值。4.2 cgSNP分析流程与工具实操做cgSNP分析我的标准流程如下选一个高质量的参考基因组。金葡菌常用的参考基因组有NCTC 8325ATCC 25923衍生株、Newman、USA300_FPR3757等。参考基因组的选择要与研究菌株的谱系尽量接近否则read mapping时比对效率会显著下降。用Snippy对每个样本独立mapping到参考基因组并鉴定变异snippy --outdir sample_snippy --ref reference.gbk --R1 clean_R1.fastq.gz --R2 clean_R2.fastq.gz用Snippy-core将多个样本的变异位点合并成core SNP matrixsnippy-core --ref reference.gbk sample1_snippy sample2_snippy sample3_snippy做重组区域过滤。金葡菌经常通过水平基因转移获取外源DNA片段这些重组区域的SNP不符合“垂直遗传”的假设如果不剔除会影响进化树结构。Gubbins可以识别并去除重组区域run_gubbins.py -p gubbins_out core.full.aln用IQ-TREE构建最大似然树iqtree -s gubbins_out.filtered_polymorphic_sites.fasta -m GTR -bb 1000 -nt 16-bb 1000表示做1000次ultrafast bootstrap这是目前构建细菌系统发育树的主流做法。最终得到的tree文件可以用iTOL在线工具进行可视化和注释把ST型、耐药表型、分离时间等meta信息做成热图或色条标在树旁边。4.3 核心基因组SNP数据的常见坑核心基因组定义的影响Snippy-core默认只分析在所有样本中都存在的位点。如果样本量太少例如只有3个菌株核心基因组会比较大但如果混入了亲缘很远的谱系core区域被大量过滤掉分辨率反而下降。样本量特别大时可以考虑用Roary或Panaroo先构建泛基因组再提取core基因做比对。重组过滤不是可选项金葡菌的mecA毒力岛和SCCmec元件经常发生基因组重组事件导致SNP密度出现局部偏高。不剔除重组区直接建树可能会得到错误的分支关系。这个步骤在发高分文章时几乎是必跑项。mapping质量与低覆盖区域Snippy默认的SNP筛选条件对低覆盖区域可能漏检。如果样本测序深度偏低低于30×建议单独检查候选SNP位点的比对情况防止因为reads错配产生假阳性SNP。之前处理一批来自甲醛固定石蜡包埋样本的测序数据时由于DNA降解严重SNP检测噪声非常大后来只能降低阈值并手动检查位点——这种类型的样本如果测序深度不够后面的分析质量很难保证。5. 可移动元件与MRSA分子流行病学进阶分析方向5.1 SCCmec分型把MRSA分得再细一点mecA基因位于一个称为SCCmec葡萄球菌染色体盒式元件的大型可移动遗传元件上。SCCmec本身还分为多个型别I型到XIII型不同型别与菌株来源高度相关HA-MRSA医院获得性MRSA以SCCmec II、III型为主CA-MRSA社区获得性MRSA以SCCmec IV、V型为主SCCmec分型可以用SCCmecFinder在线分析也可以用KMerFinder配合SCCmec数据库做本地扫描。如果只想快速判断是哪个大类可以先跑mlst并检查mecA基因是否检测到——ST239几乎都携带SCCmec III型USA300克隆则与SCCmec IVa型相关联这些谱系特征都是综合判断的依据。SCCmec元件的边界区域常含有ccr基因复合体位点特异性重组酶造成SCCmec有多种结构变异。测序组装如果只得到短contigSCCmec区域往往是断裂的不能仅凭mecA一个基因就报告“SCCmec IV型”需要结合ccr基因型别和周边结构才能正确分型。5.2 前噬菌体、质粒与毒力岛的快速筛查金葡菌基因组中包含大量可移动元件很多关键的毒力因子和耐药基因都偶联在这些元件上。tst、PVL、部分肠毒素都位于前噬菌体区域mecA位于SCCmec四环素耐药基因tetK可能存在于质粒上这些基因的移动性对流行病学传播有显著影响。常用工具有PHASTER在线预测前噬菌体区域。输出结果包含每个前噬菌体的位置、长度和完整性。plasmidSPAdes从全基因组测序数据中单独组装质粒。MobTyper或mlplasmids对质粒contig进行复制起始蛋白分型。拿到前噬菌体预测结果后可以关联检查tst基因是否刚好落在某个噬菌体区域内——如果是说明该毒力基因具有水平转移潜能。在某个金葡菌ST398菌株里检测到tst基因同时PHASTER显示该菌株存在一个完整的φSa3-like噬菌体tst正好位于噬菌体整合酶下游这就从结构上解释了为什么该菌株能获得TSST-1的编码能力。5.3 全基因组数据驱动的传播链研究设计做金葡菌暴发溯源研究时如果条件允许尽量把以下数据全部纳入综合判断核心基因组SNP系统发育树提供菌株间遗传距离耐药基因谱与药敏表型提供临床干预信息毒力因子谱提供致病性差异的解释流行病学信息病房分布、时间线、接触史基因测序数据提供的是“可能性”证据不能单独下定论。SNP差异很少但流行病学上完全没有联系时要考虑到可能存在一个未被发现的中间传播者或环境储菌库。SNP差异较大但流行病学指向明显的也不能直接排除传播关系——毕竟感染者在体内定植数周至数月期间菌株可能积累额外突变。6. 分析流程总览与工具速查完整走完一遍金葡菌WGS分析流程后整理一个速查表供日常参考分析环节核心工具关键参数/数据库判定参考原始数据质控fastpQ3090%dedup开启保留reads80%基因组组装SPAdes--careful, k-mer自动调整contig200, N5050kb组装质量评估QUAST/CheckM参考基因组或marker基因完整度99%物种确认Kraken2/Mashstandard数据库金葡菌reads95%MLST分型mlstpubmlst数据库获得ST型spa分型spaTyperspa数据库获得spa型耐药基因扫描ABRicateCARD/ResFinder覆盖度阈值默认80%关注与表型一致性毒力因子扫描ABRicateVFDB注意部分阳性基因关注PVL/TSST等系统发育分析SnippyGubbinsIQ-TREEcgSNP矩阵过滤重组0-12 SNP考虑传播SCCmec分型SCCmecFinder/KMerFinderSCCmec数据库区分HA/CA-MRSA前噬菌体预测PHASTER在线平台完整性三档判断质粒分析plasmidSPAdes/MobTyper质粒数据库耐药基因水平转移证据以上是金葡菌全基因组测序数据分析的标准路径。就拿我自己的经验来说刚上手的那段时间总觉得软件越多越好跑完这个跑那个结果一堆输出文件堆在服务器里真正能写进文章里的却没几个。后来逐渐明白分析流程的核心不是工具的数量而是问题导向——你手上有什么样本想回答什么临床或流行病学问题再倒推需要做哪些分析。如果只让我留三条建议给刚入门的同行第一质控和组装报告务必保存好这是所有下游分析可信度的基础第二拿到任何基因型结果后都要和表型数据做交叉验证数据库预测只是参考第三系统发育分析中重组过滤不可跳过否则结论很容易被金葡菌高频的同源重组误导。最后再分享一个小技巧整个分析过程中养成写分析日志的习惯把每一个软件版本、参数、数据库版本都记录下来这样不仅方便自己回溯写methods部分的时候也会省力得多。
金黄色葡萄球菌全基因组测序数据分析：从fastq到分型、耐药与溯源
拿到金葡菌全基因组测序数据后该怎么分析这个问题几乎每个刚接触生物信息学的微生物研究者都会卡住。测序公司返回来一堆raw data看起来就是几GB的fastq文件而你要从中得到结论——这是什么序列型有没有耐药基因有没有毒力因子它在暴发溯源中和其它菌株是什么关系如果脑子里面没有一张清晰的分析路径图很容易陷在“装了软件但不知道下一步干嘛”的状态里。这篇文章就围绕金黄色葡萄球菌以下简称金葡菌全基因组测序数据的完整分析流程展开从拿到fastq开始到最终拿到分型、耐药、毒力、溯源结论为止把每一步怎么做、为什么这么做、有哪些坑全部串起来。无论你是刚开始接触WGS数据分析的研究生还是已经跑过几条流程但想理顺原理的从业者这篇文章都值得收藏。1. 数据质控与组装分析链路的第一道关卡1.1 原始数据质控的几个关键指标测序公司交付的数据一般是双端fastq文件命名常见为Sample_R1.fastq.gz和Sample_R2.fastq.gz。拿到数据的第一件事不是急着组装而是做质量评估和过滤。这一步直接决定后面所有分析的可信度省不掉。我习惯用fastp做质控一条命令能同时完成质量过滤、接头去除、长度过滤和统计报告生成。相比之下传统的FastQC Trimmomatic流程也能做但需要手动串联两个软件速度也慢不少。fastp的好处是自动识别并去除接头序列而且会对双端reads做paired-end overlap分析能顺带修正一部分测序错误。fastp -i Sample_R1.fastq.gz -I Sample_R2.fastq.gz \ -o clean_R1.fastq.gz -O clean_R2.fastq.gz \ --detect_adapter_for_pe \ --thread 16 \ --html Sample_fastp.html跑完看几个核心指标Q20/Q30比例分别代表碱基错误率低于1%和0.1%的碱基占比。金葡菌基因组GC含量约33%属于AT-rich基因组如果测序平台是Illumina NovaSeq一般Q30不低于90%。clean reads数过滤后至少保留原始reads数的80%以上低于这个比例要检查是否测序质量本身有问题。插入片段长度分布正常应该在300-500bp左右。如果双端reads的overlap比例异常高说明插入片段太短对后续组装连续性会有影响。还有一个容易被忽视的指标是duplication rate。如果是细菌基因组测序文库构建时DNA投入量不足容易导致PCR重复率高。fastp默认不自动去重对细菌基因组这种低复杂度样本我一般会在参数里加上--dedup否则重复read会显著影响后续的组装完整性和变异检测准确性。注意质控报告一定要保存好。写文章或报告时审稿人大概率会要求提供测序数据质量参数回头再补跑非常麻烦。1.2 组装工具选择与参数调整质控之后进入组装环节。金葡菌基因组大小约2.8-3.0Mb属于比较小的基因组Illumina双端150bp reads的组装难度不高主流工具都能拿到不错的结果。我用得最多的是SPAdes它基于de Bruijn graph算法对细菌基因组的效果非常稳定。spades.py -1 clean_R1.fastq.gz -2 clean_R2.fastq.gz \ --careful \ -k 21,33,55,77 \ --cov-cutoff auto \ -o spades_output这里几个参数值得解释一下--carefulSPAdes官方建议在组装完成后进行mismatch correction虽然会显著增加运行时间但对提高碱基准确率有帮助。金葡菌基因组不大增加的时间完全可以接受所以我每次都开。-kk-mer长度列表。细菌基因组一般用21到77这几个值就能覆盖如果测序读长是151bp最大k-mer不要超过读长减10否则k-mer无法有效覆盖reads末端区域。--cov-cutoff auto自动过滤低覆盖度的contig。金葡菌如果混入了少量人类宿主DNA或环境DNA污染auto模式能够帮助过滤掉这些低覆盖度的contig。组装完成后用QUAST评估结果quast.py contigs.fasta -r reference.fasta -o quast_output金葡菌组装的核心指标我一般这么判断contig数在50-150之间属于正常NovaSeq数据通常更少N50在10万以上说明组装结果不错完整度检查用CheckM或BUSCO确认基因组完整性在99%以上。如果contig数超过300要怀疑是不是有质粒、前噬菌体等重复区域把组装打碎了——这个后面在可移动元件分析里会详细展开。整个分析链路中组装质量是最容易被事后追溯验证的一环。之前有人送测的样本混了两种不同菌株组装出来contig数异常多N50很低但质控报告又完全正常。后来用Mash检查污染才发现样本本身就是混合的。所以组装指标异常时不要急着调参数先怀疑数据本身。2. 菌株鉴定与分子分型确定你的菌株身份2.1 物种确认别拿到数据就默认是金葡菌很多人拿到数据后跳过物种鉴定直接开始分析这是一个隐患。实验室分离的“金葡菌”不一定真是金黄色葡萄球菌——有可能是表皮葡萄球菌、溶血葡萄球菌等CoNS凝固酶阴性葡萄球菌或者干脆是混合污染。物种鉴定最快的方法是Kraken2配合标准数据库kraken2 --db standard --threads 16 \ --paired clean_R1.fastq.gz clean_R2.fastq.gz \ --report Sample_kraken_report.txt正常金葡菌样本的reads分类结果中Staphylococcus aureus的比例应该超过95%。如果这个比例偏低需要重新审视菌株纯度。更精确的做法是用组装后的contigs跑ANI平均核苷酸一致性和参考基因组比对。金葡菌种内ANI值一般在99%以上而它与表皮葡萄球菌的ANI通常在75%-80%区间差异非常明显。2.2 MLST分型和spa分型的原理与实操物种确认无误后第一个拿到的分子分型结果通常是MLST多位点序列分型。金葡菌的MLST方案基于7个管家基因arcC、aroE、glpF、gmk、pta、tpi、yqiL的等位基因编号每个菌株得到一个7位数的数字组合即ST型。用mlst工具可以直接从contigs文件跑mlst contigs.fasta输出结果示例contigs.fasta staphylococcus_aureus 398 1 3 1 1 620 1 1其中第二个字段表示物种第三个字段是ST型后面7个数字是7个管家基因的等位基因编号。spa分型针对的是金黄色葡萄球菌特有的spa基因protein A基因的X区该区域由一系列24bp的重复单元构成不同菌株的重复单元排列组合不同形成不同的spa型别。spaTyper可以从fastq或组装结果中直接分型spaTyper --contigs contigs.fastaMLST分型和spa分型在暴发调查中扮演的角色不太一样。MLST分辨率中等适合描述菌株的大背景——比如某个地区流行的是ST239还是ST5。spa分型分辨率略高一些在金葡菌院内感染暴发调查中spa型别相同的分离株往往提示可能存在克隆传播。但要注意spa分型只分析单个位点如果spa基因发生重组或缺失突变可能出现无法分型的情况。之前遇到一个菌株spa重复区发生了部分缺失spaTyper报了一个new allele后来用Sanger测序确认是spa基因的自然变异。3. 耐药基因与毒力因子鉴定挖掘表型背后的遗传基础3.1 耐药基因鉴定从数据库到临床意义金葡菌的耐药问题集中在MRSA耐甲氧西林金黄色葡萄球菌其分子基础是mecA基因或mecC等变异体编码的PBP2a蛋白导致β-内酰胺类抗生素失效。另外还有vanA基因介导的万古霉素耐药、erm基因介导的大环内酯类耐药、fosB基因介导的磷霉素耐药等。用ABRicate可以批量扫描耐药基因abricate --db resfinder contigs.fasta sample_resfinder.txt abricate --db card contigs.fasta sample_card.txtABRicate内置了多个数据库ResFinder、CARD、ARG-ANNOT等。实测下来CARD数据库对点突变型耐药基因的覆盖更全比如gyrA的喹诺酮耐药突变ResFinder对获得性耐药基因的检索更准。两个数据库的结果可以互相印证。有一个容易翻车的地方是mecA基因检出只代表存在耐药决定簇不代表体外药敏试验一定耐药。mecA的表达受mecI-mecR1调控系统影响有些菌株携带mecA但表达量很低表现为隐匿型MRSAoxacillin敏感。反过来也一样mecA阴性但苯唑西林耐药的情况也可能出现——这种表型通常由其他机制比如blaZ过表达或PBP突变介导。所以基因型报告和药敏结果不一致时不要急着质疑测序数据先把两边的证据再核对一遍。3.2 毒力因子鉴定金葡菌的“武器库”扫描金葡菌的毒力因子构成非常复杂主要包括黏附素如fnbA、fnbB、clfA、毒素如PVL杀白细胞素、TSST-1中毒性休克毒素、肠毒素sea-see等、免疫逃逸因子如spa、cIfB、以及生物膜形成相关基因如icaADBC操纵子。同样用ABRicate扫描VFDB数据库就行abricate --db vfdb contigs.fasta sample_vfdb.txt对于金葡菌有几个毒力因子的临床意义特别值得关注PVLlukF-PV/lukS-PV杀白细胞素社区获得性MRSACA-MRSA的重要标志。如果检测到PVL阳性基本可以判断菌株属于CA-MRSA谱系。TSST-1tst基因导致中毒性休克综合征的超抗原毒素。tst基因通常由前噬菌体携带散在分布于金葡菌基因组中。肠毒素基因簇sea-see食物中毒的重要致病因子金葡菌食物中毒暴发调查中必查。ABRicate默认的覆盖度阈值是80%有时候可能会漏掉部分截短的毒力基因。建议扫描完看一眼原始比对结果特别是那些“partial positive”的基因不要直接按二元变量处理要结合比对覆盖度判断这个基因是否完整存在。3.3 表型与基因型关联分析时的注意事项耐药基因检测结果要跟表型药敏数据做关联分析时有几个逻辑陷阱有些耐药基因不是组成型表达比如诱导型克林霉素耐药iMLSB基因检测只能告诉你“有这个基因”无法判断表达模式。染色体介导的耐药位点比如gyrA、grlA的QRDR区域点突变在ABRicate默认参数下可能检测不到因为这类突变不是基因的“存在/缺失”问题而是点突变问题。需要用PointFinder或其他工具做专门的点突变筛查。基因组测序检出的耐药基因如果与药敏结果明显矛盾考虑是否为杂合菌群或者样本污染必要时回到原始平板重新挑取单克隆做验证。4. 系统发育分析与暴发溯源SNP层级的分辨率4.1 为什么MLST不够用SNP才是金标准在院内感染暴发调查中同一个ST型别甚至同一个spa型别的菌株并不代表它们一定来自同一个传播链。MLST和spa分型的分辨率不足以区分短期暴发中的菌株差异。这时候需要做核心基因组SNP分析cgSNP——比较菌株间在核心基因组所有菌株共有基因中的单核苷酸多态性差异。不同菌株之间SNP差异数量的流行病学解释大致有个经验范围具体取决于菌株谱系和暴发时间跨度SNP差异范围流行病学推断参考0-12个SNP高度提示近期克隆传播12-30个SNP不排除传播关系结合流调判断30-100个SNP可能是较长时间或间接传播链100个SNP倾向于流行背景多样性非直接传播这个阈值不是写在教科书里的绝对标准不同研究机构的cutoff取值不完全一致但0-12的区间被广泛接受为“同一个暴发克隆”的参考阈值。4.2 cgSNP分析流程与工具实操做cgSNP分析我的标准流程如下选一个高质量的参考基因组。金葡菌常用的参考基因组有NCTC 8325ATCC 25923衍生株、Newman、USA300_FPR3757等。参考基因组的选择要与研究菌株的谱系尽量接近否则read mapping时比对效率会显著下降。用Snippy对每个样本独立mapping到参考基因组并鉴定变异snippy --outdir sample_snippy --ref reference.gbk --R1 clean_R1.fastq.gz --R2 clean_R2.fastq.gz用Snippy-core将多个样本的变异位点合并成core SNP matrixsnippy-core --ref reference.gbk sample1_snippy sample2_snippy sample3_snippy做重组区域过滤。金葡菌经常通过水平基因转移获取外源DNA片段这些重组区域的SNP不符合“垂直遗传”的假设如果不剔除会影响进化树结构。Gubbins可以识别并去除重组区域run_gubbins.py -p gubbins_out core.full.aln用IQ-TREE构建最大似然树iqtree -s gubbins_out.filtered_polymorphic_sites.fasta -m GTR -bb 1000 -nt 16-bb 1000表示做1000次ultrafast bootstrap这是目前构建细菌系统发育树的主流做法。最终得到的tree文件可以用iTOL在线工具进行可视化和注释把ST型、耐药表型、分离时间等meta信息做成热图或色条标在树旁边。4.3 核心基因组SNP数据的常见坑核心基因组定义的影响Snippy-core默认只分析在所有样本中都存在的位点。如果样本量太少例如只有3个菌株核心基因组会比较大但如果混入了亲缘很远的谱系core区域被大量过滤掉分辨率反而下降。样本量特别大时可以考虑用Roary或Panaroo先构建泛基因组再提取core基因做比对。重组过滤不是可选项金葡菌的mecA毒力岛和SCCmec元件经常发生基因组重组事件导致SNP密度出现局部偏高。不剔除重组区直接建树可能会得到错误的分支关系。这个步骤在发高分文章时几乎是必跑项。mapping质量与低覆盖区域Snippy默认的SNP筛选条件对低覆盖区域可能漏检。如果样本测序深度偏低低于30×建议单独检查候选SNP位点的比对情况防止因为reads错配产生假阳性SNP。之前处理一批来自甲醛固定石蜡包埋样本的测序数据时由于DNA降解严重SNP检测噪声非常大后来只能降低阈值并手动检查位点——这种类型的样本如果测序深度不够后面的分析质量很难保证。5. 可移动元件与MRSA分子流行病学进阶分析方向5.1 SCCmec分型把MRSA分得再细一点mecA基因位于一个称为SCCmec葡萄球菌染色体盒式元件的大型可移动遗传元件上。SCCmec本身还分为多个型别I型到XIII型不同型别与菌株来源高度相关HA-MRSA医院获得性MRSA以SCCmec II、III型为主CA-MRSA社区获得性MRSA以SCCmec IV、V型为主SCCmec分型可以用SCCmecFinder在线分析也可以用KMerFinder配合SCCmec数据库做本地扫描。如果只想快速判断是哪个大类可以先跑mlst并检查mecA基因是否检测到——ST239几乎都携带SCCmec III型USA300克隆则与SCCmec IVa型相关联这些谱系特征都是综合判断的依据。SCCmec元件的边界区域常含有ccr基因复合体位点特异性重组酶造成SCCmec有多种结构变异。测序组装如果只得到短contigSCCmec区域往往是断裂的不能仅凭mecA一个基因就报告“SCCmec IV型”需要结合ccr基因型别和周边结构才能正确分型。5.2 前噬菌体、质粒与毒力岛的快速筛查金葡菌基因组中包含大量可移动元件很多关键的毒力因子和耐药基因都偶联在这些元件上。tst、PVL、部分肠毒素都位于前噬菌体区域mecA位于SCCmec四环素耐药基因tetK可能存在于质粒上这些基因的移动性对流行病学传播有显著影响。常用工具有PHASTER在线预测前噬菌体区域。输出结果包含每个前噬菌体的位置、长度和完整性。plasmidSPAdes从全基因组测序数据中单独组装质粒。MobTyper或mlplasmids对质粒contig进行复制起始蛋白分型。拿到前噬菌体预测结果后可以关联检查tst基因是否刚好落在某个噬菌体区域内——如果是说明该毒力基因具有水平转移潜能。在某个金葡菌ST398菌株里检测到tst基因同时PHASTER显示该菌株存在一个完整的φSa3-like噬菌体tst正好位于噬菌体整合酶下游这就从结构上解释了为什么该菌株能获得TSST-1的编码能力。5.3 全基因组数据驱动的传播链研究设计做金葡菌暴发溯源研究时如果条件允许尽量把以下数据全部纳入综合判断核心基因组SNP系统发育树提供菌株间遗传距离耐药基因谱与药敏表型提供临床干预信息毒力因子谱提供致病性差异的解释流行病学信息病房分布、时间线、接触史基因测序数据提供的是“可能性”证据不能单独下定论。SNP差异很少但流行病学上完全没有联系时要考虑到可能存在一个未被发现的中间传播者或环境储菌库。SNP差异较大但流行病学指向明显的也不能直接排除传播关系——毕竟感染者在体内定植数周至数月期间菌株可能积累额外突变。6. 分析流程总览与工具速查完整走完一遍金葡菌WGS分析流程后整理一个速查表供日常参考分析环节核心工具关键参数/数据库判定参考原始数据质控fastpQ3090%dedup开启保留reads80%基因组组装SPAdes--careful, k-mer自动调整contig200, N5050kb组装质量评估QUAST/CheckM参考基因组或marker基因完整度99%物种确认Kraken2/Mashstandard数据库金葡菌reads95%MLST分型mlstpubmlst数据库获得ST型spa分型spaTyperspa数据库获得spa型耐药基因扫描ABRicateCARD/ResFinder覆盖度阈值默认80%关注与表型一致性毒力因子扫描ABRicateVFDB注意部分阳性基因关注PVL/TSST等系统发育分析SnippyGubbinsIQ-TREEcgSNP矩阵过滤重组0-12 SNP考虑传播SCCmec分型SCCmecFinder/KMerFinderSCCmec数据库区分HA/CA-MRSA前噬菌体预测PHASTER在线平台完整性三档判断质粒分析plasmidSPAdes/MobTyper质粒数据库耐药基因水平转移证据以上是金葡菌全基因组测序数据分析的标准路径。就拿我自己的经验来说刚上手的那段时间总觉得软件越多越好跑完这个跑那个结果一堆输出文件堆在服务器里真正能写进文章里的却没几个。后来逐渐明白分析流程的核心不是工具的数量而是问题导向——你手上有什么样本想回答什么临床或流行病学问题再倒推需要做哪些分析。如果只让我留三条建议给刚入门的同行第一质控和组装报告务必保存好这是所有下游分析可信度的基础第二拿到任何基因型结果后都要和表型数据做交叉验证数据库预测只是参考第三系统发育分析中重组过滤不可跳过否则结论很容易被金葡菌高频的同源重组误导。最后再分享一个小技巧整个分析过程中养成写分析日志的习惯把每一个软件版本、参数、数据库版本都记录下来这样不仅方便自己回溯写methods部分的时候也会省力得多。