【7.4.2】 HUMAnN2--宏基因组和宏转录组的物种水平功能分析
HUMANN2 是一个 Python2/3 兼容包。最新版本可以通过 pip 或 HomeBrew 安装(或通过 http://huttenhower.sph.harvard.edu/humann2 从源代码安装)。HUMAnN2 还捆绑为 bioBakery 虚拟机的一部分,可用作 Vagrant Box、Google Cloud 映像和 Amazon Web Services AMI(通过 http://huttenhower.sph.harvard.edu/biobakery )
微生物群落的功能概况通常是使用全面的宏基因组或宏转录组序列读取搜索生成的,这很耗时,容易出现虚假映射,并且通常仅限于群落水平的量化。我们开发了 HUMAnN2,这是一种分层搜索策略,可以对宿主相关和环境群落进行快速、准确和物种解析的功能分析。HUMAnN2 识别社区的已知物种,将读数与其泛基因组对齐,对未分类的读数执行翻译搜索,最后量化基因家族和通路。相对于纯翻译搜索,HUMAnN2 速度更快,产生更准确的基因家族谱。我们应用 HUMAnN2 研究海洋代谢的临床变异、人类微生物组途径中的生态贡献模式、物种基因组与转录贡献的差异,以及菌株分析。此外,我们引入了“贡献多样性”来解释不同微生物群落类型的生态组装模式(ecological assembly)。
前言
从宏基因组学和元转录组学 (‘meta’omic’) 测序数据中分析微生物群落功能是微生物生态学中一个至关重要的挑战。它有可能表征在许多社区1 中观察到的广泛的生化“暗物质” ,并将特定的分子活动与环境2和健康相关的3表型联系起来。与分类学分析相比,功能分析旨在量化由已知和未表征的社区成员4贡献的基因和代谢途径的内容。虽然分类学分析可以在元组学测序读数的最大信息子集上进行5 , 6,全面的功能分析必须考虑所有读取和它们可能源自的巨大基因空间,从而增加了相当大的分析复杂性。
对于宏基因组的功能谱存在几种方法7,8,9,它的一个子集已应用于metatranscriptomes 10,11,12,13。其中包括 HUMAnN 14,我们在人类微生物组计划 (HMP) 15期间开发了它,用于与宿主相关和环境相关的元组功能分析。与后来的方法一样,HUMAnN 解释元组学测序读数的翻译搜索以重建代谢功能。尽管现有方法受益于翻译搜索的最新进展16 , 17 , 18,它们仍然比核苷酸水平分析慢得多。此外,虽然一些功能分析方法结合了用于数据库细化7或目标量化9 的分类学概念,但大多数仅限于报告社区水平的丰度,而不是每个生物的贡献。同样,功能分析落后于微生物群落19、20、21菌株水平分析的努力,尽管对物种内菌株特定功能的认识越来越多。
我们开发了 HUMAnN2 以将分类信息与功能配置文件相结合,并通过将分层方法与核苷酸级搜索、加速翻译搜索和通路重建组件相结合来限制翻译搜索瓶颈。HUMAnN2 超过了纯翻译搜索策略的准确性和性能。此外,由 HUMAnN2 量化的基因和通路丰度自动分为已知和未表征物种的贡献。这为解释宿主相关和环境群落元组提供了以前无法获得的细节。
一、结果
算法概述
HUMAnN2 实施了“分层搜索”策略,以在物种级分辨率下快速准确地分析元组的功能内容(图1a,补充图1,方法),其结果也可用于应变分析. 在第一层中,HUMAnN2 通过使用 MetaPhlAn2 筛选 DNA 或 RNA 读数来快速识别样本中的已知微生物物种(参考文献22)。HUMAnN2 然后通过合并已识别物种的预先构建的、功能注释的泛基因组来构建特定于样本的数据库23. 在第二层中,HUMAnN2 将所有样本读数与样本的泛基因组数据库进行核苷酸级映射。相对于全面的翻译搜索,针对相关泛基因组的核苷酸级映射可以快速解释大部分读数,而虚假比对的机会较少。在第三层也是最后一层,对与已识别物种的泛基因组不一致的读取进行针对综合蛋白质数据库(默认情况下,UniRef90 或 UniRef50(参考24))的加速翻译搜索。
a,HUMAnN2 用于元组学功能分析的分层搜索算法概述(在补充图1 中进行了扩展)。b,在合成肠道宏基因组上评估的 HUMAnN2 的分层搜索与纯翻译搜索。c、d、敏感性、精确度和整体准确度(1 - Bray-Curtis 不相似性)计算了 ( c ) 基因家族和 ( d ) 通路丰度谱相对于全社区水平(“整体”)的金标准和对于每个分层。e,HUMAnN2 在量化社区总 COG 丰度的任务中与其他方法进行了比较。运行时反映了 8 个 CPU 内核上的多线程执行。
分层搜索生成元组学读数到具有已知或模糊分类法的基因序列的映射。这些映射按质量和序列长度加权,以估计每个生物体和社区总基因家族的丰度,这些丰度可以重新组合到其他功能系统(例如,COG 25、KOs 26、Pfam 域27和 GO 术语28)。最后,进一步分析注释到代谢酶的基因家族,以重建和量化社区和每个生物体中完整的代谢途径(默认情况下,MetaCyc 29)。
分层搜索优于纯翻译搜索 Tiered search outperforms pure translated search
我们通过分析合成宏基因组(方法)来评估 HUMAnN2 的准确性。我们首先模拟了一个包含 1000 万个 100 核苷酸 (nt) DNA 读数 (1 Gnt) 的人类肠道宏基因组,这些读数来自 HMP 粪便样本15 中20 种最丰富的细菌物种。物种的丰度在几何上从 ~0.1x 到 ~70x 基因组覆盖范围交错(图1b),包括拟杆菌属的 9 个成员——这两个挑战都对每个物种的准确分析提出了挑战。我们使用 HUMAnN2 的分层搜索和纯翻译搜索策略分析了这个合成宏基因组(补充说明1和补充图2;对 100 名成员、非人类相关社区的平行分析)。
对于群落基因家族(UniRef90)丰度,HUMAnN2 分层搜索的灵敏度、精确度和总体准确度(1 – Bray-Curtis 不相似性)分别为 86%、90% 和 89%(图1c)。因此,HUMAnN2 (i) 检测到社区中大多数预期的基因家族,(ii) 仅报告一小部分虚假检测的家族,以及 (iii) 将绝大多数读数正确分配给其源家族。通过纯翻译搜索分析的基因家族不太准确(总体准确度为 67%),部分原因是在将所有样本读数与综合蛋白质数据库进行比对时,虚假比对的可能性更大。
HUMAnN2 分层搜索的每个物种准确性对于以 1 倍或更高基因组覆盖率存在的 14 个物种仍然很高,包括九个拟杆菌属物种。低于 1 倍的覆盖率,灵敏度和整体准确度会随着覆盖率的下降而下降,因为该域中更多的基因家族样本不足。然而,低覆盖率物种的精确度始终保持在高水平,表明它们的泛基因组没有招募大量不相关的读数。传递到翻译搜索层并映射到蛋白质的读数的小子集 (1.4%) 产生了“未分类”分层,对总体错误的贡献很小(补充说明2)。
HUMAnN2 分层搜索的准确率趋势在通路水平上相似,通路精度普遍超过基因家族精度(图1d)。这是因为虚假匹配完整通路的难度更大,这需要多个不同的反应(基因家族)来虚假招募读取。同时,HUMAnN2 对检测完整(或接近完整)路径的要求会导致灵敏度和整体准确度随着覆盖率的降低而衰减得更快。对于特征不太好的样本,纯翻译搜索策略固有的基因级错误在通路量化过程中趋于“平滑”,尽管来自纯翻译搜索的通路配置文件仍然不如来自 HUMAnN2 分层搜索的通路配置文件准确(87% 98%)。
在综合评估中,HUMAnN2 的分层搜索也比纯翻译搜索快 3 倍(运行时间 <1 小时;图1e和补充说明 1和3)。我们进一步对跨越六个身体部位的 397 个 HMP 宏基因组的分层搜索的性能进行了基准测试(补充说明4)。在典型的样本中,在泛基因组搜索期间映射了约 60% 的读数,在翻译搜索期间映射了额外的约 20%(补充图3)。因此,对于表征良好的真实世界宏基因组,HUMAnN2 解释了快速泛基因组搜索期间的大部分样本读取,使其比纯翻译搜索策略更有效。
与现有方法的比较
我们将 HUMAnN2 与基于纯翻译搜索的现有功能分析方法进行了比较:HUMAnN1(参考文献14)、COGNIZER 10、MEGAN 12和 ShotMAP 13(图1e))。这种比较基于对直系同源群 (COG) 丰度的社区级集群的估计,这是所有方法 (方法) 共有的输出格式。我们基于 UniRef90 为 ShotMAP 构建了一个自定义搜索数据库,并使用了其他三种方法的推荐数据库。我们注意到这三种方法的 COG 定义系统可能与我们基于 UniProt 的黄金标准不同,这可能会影响它们相对于 HUMAnN2 和 ShotMAP 的准确性。然而,在这些评估中采样的分离基因组早于除 HUMAnN1 之外的所有方法,这限制了由于数据库覆盖而导致的潜在偏差。
HUMAnN2 的分层搜索的总体准确率最高 (97%),其次是 HUMAnN2 的纯翻译搜索 (83%)、ShotMAP (72%)、HUMAnN1 (59%)、MEGAN (56%) 和 COGNIZER (43%)。HUMAnN2 纯翻译搜索准确性的提高可能归因于我们旨在最大化特异性的事后对齐过滤和加权(方法;补充图4和5)。HUMAnN2 的分层搜索在 45 分钟内分析了 1000 万次读取的合成宏基因组。这与使用加速翻译搜索16(42 分钟)的HUMAnN1 类似,但 HUMAnN2 提供了更详细的输出,并考虑了大约 20 倍的序列空间。HUMAnN2 比所有其他方法快 3 倍以上。
在元转录组和非参考物种上的表现
我们对 HUMAnN2 进行了广泛的额外评估。在分析广泛定义的基因家族(UniRef50;补充说明3)或合成肠道宏转录组(补充说明5和补充图6)时,HUMAnN2 保持准确和有效。至关重要的是,HUMAnN2 在包含已知物种和新物种的新分离株的宏基因组上表现出色(后者由翻译的搜索层进行分析)。这是通过分析一个复杂的(100 个成员)合成群落同时保留 HUMAnN2 泛基因组数据库的一部分来模拟新物种来实现的(补充说明1和补充图2)并且通过施加HUMAnN2等方法来分离物基因组的社区该后日期方法数据库(补充说明6和补充图7 - 10)。
我们另外将 HUMAnN2 与合成宏基因组的宏基因组组装进行了比较(补充说明7)。该评估扩展了之前对真实世界人类宏基因组30的组装和基于参考的方法的比较,产生了非常相似的域级功能多样性排名。虽然组装有利于揭示深度测序的人类宏基因组中的新序列多样性,但 HUMAnN2 在具有适度测序深度的宏基因组中鉴定了更多已知域。这一优势源于通过组装31检测低覆盖率宏基因组序列的已知挑战,这也在我们的综合评估中观察到
。。。。
二、方法
这些方法详细介绍了 HUMAnN2 算法、其数据库的构建、我们对合成宏基因组的评估以及贡献多样性计算。与我们的 HUMAnN2 应用相关的方法(HMP 宏基因组、红海宏基因组和配对 IBDMDB 元组的分析,the analyses of HMP metagenomes, Red Sea metagenomes, and paired IBDMDB meta’omes)在补充说明4中提供。与我们的评估上合成metatranscriptomes,新颖分离物的基因组,并组装宏基因组的方法是在提供补充说明 5,6,和7分别。方法详细信息也可以在自然研究报告摘要(Nature Research Reporting)中找到。
2.1 算法概述
HUMAnN2 是一个系统,用于加速对宿主相关和环境相关微生物群落的鸟枪法宏基因组和宏转录组学 (meta’omic) 测序进行功能分析。HUMAnN2 实施分层搜索策略,包括三个搜索阶段(层)。
- 在第一层搜索中,元组学样本被快速筛选以识别潜在群落中的已知物种。然后,通过连接检测到的物种的预先计算的、功能注释的泛基因组,该信息用于构建样本的自定义基因序列数据库。
- 在第二个搜索层中,整个样本与该数据库对齐,产生 (i) 每个物种、每个基因的比对统计数据和 (ii) 未映射读取的集合。
- 在最后的搜索层中,未映射的读取通过翻译搜索与用户指定的(通常是全面的和非冗余的)蛋白质数据库进行比对,产生 (i) 分类未分类的每个基因比对统计数据和 (ii) 新读取的集合。每个基因的比对统计数据根据比对质量、覆盖率和序列长度进行加权,以产生 (i) 社区的基因丰度值和 (ii) 根据每个物种和“未分类”贡献分层。基因丰度值最终应用于代谢网络重建,以识别和量化社区中的途径(也根据每个物种和“未分类”贡献进行分层)。这些过程,包括基础数据库和搜索参数,将在下面详细展开。
2.2 基因和通路参考数据作为 HUMAnN2 的固定输入
综合蛋白质数据库:
HUMAnN2 使用 UniRef90 和 UniRef50(参考文献24)作为综合的、非冗余的蛋白质序列数据库。
- 简而言之,UniRef90 代表 UniProt 45中所有非冗余蛋白质序列的聚类,使得簇中的每个序列与簇中最长序列(簇种子)的 90% 同一性和 80% 覆盖率对齐。每个结果簇由单个序列表示(通常是簇中注释最好的成员,不一定是种子)。
- UniRef50 是通过对所有 UniRef90 代表性序列进行聚类来构建的,以使聚类与 50% 的氨基酸序列同一性和 80% 的簇种子覆盖率对齐。
我们使用 UniRef90 和 UniRef50 簇 (i) 作为描述微生物基因组中基因家族结构的基础和 (ii) 作为翻译元组学搜索的综合数据库(见下文)。HUMAnN2 使用的蛋白质注释(例如,酶委员会 (EC) 编号、COG 25、KO 26、Pfam 域27, 和 GO term 28 assignments) 是从代表性 UniProt 序列的注释中推断出来的。
ChocoPhlAn 泛基因组:
HUMAnN2 中的核苷酸水平搜索是使用物种泛基因组的集合进行的。我们在 HUMAnN2 中将这个集合称为“ChocoPhlAn”。(早期版本的 ChocoPhlAn 发布为 MetaRef 46;纳入 HUMAnN2 的 ChocoPhlAn 版本与底层 MetaPhlAn2 及其标记数据库22 相同。)物种的泛基因组是物种蛋白质编码潜力的非冗余表示。为了构建给定物种的泛基因组,我们从 NCBI GenBank 和/或 RefSeq 下载该物种的所有可用分离基因组,以及相关的编码序列 (CDS) 注释。使用 PhyloPhlAn 47分析每个分离株基因组以确认正确的分类位置。使用 UCLUST 48,然后我们将来自给定物种的高质量分离基因组的所有 CDS 以 97% 的核苷酸同一性进行聚类。保存来自每个集群的一个代表性(质心)序列。这些质心序列构成了物种的泛基因组。这些步骤是在 MetaPhlAn2 开发过程中进行的。
为了使用 ChocoPhlAn 进行功能分析,我们将每个泛基因组质心序列注释为 UniRef90 和 UniRef50,方法是 (i) 翻译质心以产生氨基酸序列,然后 (ii) 对 UniRef90 执行蛋白质级搜索。如果 UniRef90 中质心的最佳命中符合纳入相应 UniRef90 集群的标准(>90% 氨基酸同一性和 >80% 覆盖率),则该质心被注释到 UniRef90 集群并继承其相应的 UniRef50 注释。如果不是,则质心被标记为“UniRef90_unknown”,并对 UniRef50 进行类似的搜索(需要与 UniRef50 序列具有 >50% 的同一性)。如果此搜索也失败,则质心被标记为“UniRef50_unknown”。ChocoPhlAn 包括 >4,000 种细胞微生物(细菌、古细菌和真菌)的泛基因组,其中包括超过 1800 万个基因簇。HUMAnN2 v0.9.6 增加了对 >3,000 个病毒泛基因组的支持,其中包括 >100,000 个基因簇。
将 UniRef90/50 基因家族与 MetaCyc 反应相关联:
HUMAnN2 生成的所有比对都折叠为 UniRef90 或 UniRef50 基因家族,它们构成了该方法解析度最高的主要输出。在代谢途径重建之前,基因家族必须进一步分解为酶/反应丰度。这需要生成一个映射,将 UniRef90 和 UniRef50 标识符与 MetaCyc 反应联系起来。这些链接是通过两种方式建立的。
- 首先,MetaCyc 反应与 UniProt 中的蛋白质子集相关,这些蛋白质由 UniProt 登录号 (AC) 标识。由于 UniProt 中的每个蛋白质都与 UniRef90 簇(并且通过扩展为 UniRef50 簇)相关联,因此反应-AC 关联被转换为反应-UniRef90 和反应-UniRef50 关联以用于 HUMAnN2。
- 其次,MetaCyc 反应与酶委员会 (EC) 目录中的条目相关联,这是酶活性的四级分层描述。 UniProt 条目(以及,通过扩展,UniRef 条目)也与 EC 编号相关联。 这种关系使 MetaCyc 反应和 UniRef90/50 标识符的附加传递关联使用 EC 注释作为桥梁。 为了保持特异性,在此过程中仅使用了最高特异性的 EC 注释(例如,与 EC 1.1.1 关联的 UniRef90 条目不会链接到与 EC 1.1.1.1 关联的 MetaCyc RXN,反之亦然 允许映射)。 具有至少一个 UniRef90(或 UniRef50)关联的 MetaCyc RXN 在 HUMAnN2 中被认为是“可量化的”。
MetaCyc 对通路作图的反应
HUMAnN1(参考文献14)结合了 KEGG 的结构化通路句法26以提高通路重建和量化的准确性。此句法指定 (i) 必须满足以完成途径的反应,以及 (ii) 可能通过途径的替代途径(可通过不同的反应组合满足)。我们通过解析 MetaCyc 的通路定义文件为 MetaCyc 通路生成了相应的结构。更具体地说,每个途径都被解析为连接初始反应物和最终产物的有向无环图。(MetaCyc 的“超级路径”被解析为它们各自的子路径,并删除了递归路径。)
路径中的每个反应节点都被注释以描述它是否通过 AND 或 OR 关系与其他节点连接(例如,表明反应 1 和 2 都需要将 A 转换为 B,或者反应 1 或 2 可以执行转换)。如果存在一条从初始反应物到最终产物的路径,该路径仅通过在给定元组学样本(见下文)中检测到(非零丰度)的反应节点,则称该路径是满足的。路径被排除 (i) 如果它们包含少于四个可量化的反应(与 4 级 EC 编号相关的反应,而后者又与 UniRef90 和 UniRef50 家族相关)或 (ii) 如果它们包含 >10% 的不可量化反应(在否则可接受的路径在结构化路径语法中被标记为“可选”)。
2.3 通过分层搜索量化基因家族
分类预筛选
HUMAnN2 将质量控制(including host-read-depleted)元组作为 FASTA 或 FASTQ 文件(带有可选的 GZIP 压缩)作为输入。DNA/RNA 读数最初使用带有默认参数的 MetaPhlAn2 进行筛选(生成的 MetaPhlAn2 输出保存为 HUMAnN2 中的临时输出)。MetaPhlAn2 检测到的高于目标相对丰度阈值的微生物种类将传递到下一个搜索层(泛基因组搜索)。使用 0.0001 (0.01%) 相对丰度的宽松检测阈值作为默认值,这相当于 10-Gnt 宏基因组中 5-Mbp 微生物基因组的 0.1 倍覆盖率,其中 50% 的读数映射到已测序的分离基因组.
泛基因组搜索
HUMAnN2 接下来将预筛选中检测到的物种的泛基因组连接为单个 FASTA 文件,然后将其作为输入用于构建 Bowtie 2 索引49。然后,在“非常敏感”模式下使用 Bowtie 2 对所有样本读数(如上所述)进行分析。由于 HUMAnN2 与孤立的编码序列比对,因此在评估 Bowtie 2 比对质量时不考虑读取末端配对关系。
翻译搜索
未能与泛基因组数据库对齐的读取通过针对用户指定的蛋白质数据库的翻译搜索进行映射。有四种选择:UniRef90 和 UniRef50 的完整版本,以及 UniRef90 和 UniRef50 的简化版本,仅包含与 MetaCyc 反应相关的蛋白质(在补充说明3 中进一步讨论)。HUMAnN2 可以调用三个已翻译的搜索二进制文件来完成此任务:DIAMOND 16、RAPSearch2(参考文献17)和 USEARCH 48. DIAMOND 是推荐的默认值。HUMAnN2 会根据用户是针对 UniRef90 集群还是更广泛(更具包容性)的 UniRef50 集群进行映射来调整翻译搜索的参数。例如,当使用 DIAMOND 对 UniRef50 进行翻译搜索时,会调用“–sensitive”搜索标志。翻译搜索的最终输出是读取与蛋白质比对统计的表格报告(表格 BLAST 格式)。
对齐后处理
HUMAnN2 中的比对经过后处理以考虑作图质量和数据库序列长度。如果一个 read 与不同的数据库序列有两个或多个高质量的比对,则该 read 的单个计数将按照平方比对同一性按比例划分到相应的序列中。这是 HUMAnN1 中实现的默认对齐加权程序的更通用版本,它基于对齐E值(在某些对齐软件中缺乏严格等效的统计数据,例如 Bowtie 2)。值得注意的是,在 HUMAnN1 评估期间发现各种类似的加权方案同样好,并且都明显优于朴素的最佳命中映射14。
序列的加权计数通过数据库序列的可比对长度(以千碱基为单位)进一步标准化,以产生以每千碱基 (RPK) 为单位的读数计数。(可比对长度是数据库序列的总长度减去读取的比对长度加 1:与数据库序列等效比对可能已经开始的位置数。)这些程序适用于针对 ChocoPhlAn 泛基因组的核苷酸水平比对并根据 UniRef90/UniRef50 翻译对齐。根据 UniRef90/UniRef50 注释(或 UniRef90_unknown/UniRef50_unknown,如果不存在注释),对 ChocoPhlAn 泛基因组中序列的加权命中在物种内求和。在翻译搜索期间对 UniRef90/UniRef50 家族的加权直接命中被求和并分配到“未分类”物种箱。
HUMAnN2 的翻译搜索使用全面的(而不是特定于样本的)序列数据库,这会导致出现更多虚假比对的机会。为了弥补这一点,HUMAnN2 在应用上面概述的一般加权程序之前,以两种额外的方式过滤翻译对齐结果。首先,如果在比对中使用了大部分读段(可调默认值:90% 查询覆盖率),我们说读段与蛋白质“对齐”。这迫使读取的翻译对齐更接近于 Bowtie 2 的非局部对齐模式(如在泛基因组搜索中使用)。接下来,读取的权重仅分布在序列被对齐读取“良好覆盖”的蛋白质上(可调默认值:覆盖位置的 50%)。如果没有这样的过滤器,对于小、5)。从未“良好对齐”或仅与覆盖不良的蛋白质对齐的读取与未对齐的读取一起导出,用于 HUMAnN2 外部的下游分析(例如,新基因序列的组装)。
2.4 量化通路丰度和覆盖率
使用上述 UniRef50/UniRef90 到 MetaCyc 反应映射,反应的丰度计算为映射到反应的所有基因家族的丰度总和。这些总和是针对每个物种、“未分类”层和整个群落计算的,与 HUMAnN2 的基因水平丰度报告一致。
HUMAnN2 计算通路丰度(拷贝数)和覆盖率(检测置信度)的程序主要按照 HUMAnN1 中的描述和基准计算(参考文献14),并添加了修改以解释(i)从 KEGG 到基于 MetaCyc 的通路定义和 (ii) 需要以分层(按物种)的方式以及社区范围内的方式计算值。从一组反应丰度开始,HUMAnN2 首先执行(可选)间隙填充步骤,以解决明显耗尽的反应或注释不足的问题。HUMAnN2 中的默认间隙填充用次最少丰度反应的丰度替换了通路中最少丰度的反应。在间隙填充计算中不考虑可选反应。接下来,MinPath 50用于确定一组简洁的途径来解释观察到的反应。然后按照 HUMAnN1 的结构化(默认)或非结构化路径定义方法计算每个路径的丰度和覆盖率。对于结构化途径,丰度计算为反应丰度的调和平均值(在优化替代子途径和可选反应之后);对于非结构化途径,丰度计算为途径中前 50% 丰度最高的反应的平均值。将反应丰度转换为反应检测置信度的度量后,覆盖率的计算方法类似。这些程序针对在每个物种中检测到的反应、“未分类”反应丰度值和群落总丰度值进行。
参考资料
- 2018,Species-level functional profiling of metagenomes and metatranscriptomes。 https://www.nature.com/articles/s41592-018-0176-y
个人公众号,比较懒,很少更新,可以在上面提问题,如果回复不及时,可发邮件给我: tiehan@sina.cn