古代和现代DNA中的瘟疫
博文名称：Plague in Ancient and Modern DNA
博文链接：https://towardsdatascience.com/plague-in-ancient-and-modern-dna-f8e52e3d8a11
发表时间：Apr 3, 2021
老彼得·勃鲁盖尔《死亡的胜利》，马德里普拉多博物馆
这是《生命科学的数学统计和机器学习》专栏的第23篇文章，在这里，我用通俗易懂的语言讨论了计算生物学中常见的一些神秘的分析技术。DNA测序技术应用于考古材料，极大地丰富了我们对人类历史的认识。例如，分析历史墓葬中人类遗骸的DNA可以提供有关古代流行病（例如鼠疫耶尔森氏菌引起的瘟疫）的重要信息。然而，检测古代病原体的常用方法缺乏特异性，可能会导致错误的发现。
在本文中，我将展示如何将序列比对与常见的病原微生物检测标准（例如覆盖深度）一起应用于现代宏基因组样本，从而很容易发现不应该存在的鼠疫耶尔森菌。我还将讨论为什么在解释古代病原体发现时需要谨慎，以及如何使用更好的度量方法（如覆盖范围）来更可靠地检测病原体。
在废弃堆中寻宝
DNA 测序技术的出现有利于许多研究领域，包括个性化医学、人口遗传学和古代DNA (aDNA)。从历史上看，aDNA 研究以人类进化为中心，主要是人类的DNA的提取、测序和分析。然而，最近的研究表明，对以前常被视为副产品的内源性微生物群落中的aDNA进行分析，可以为研究古代流行病、生活方式和过去人口迁移等提供极其宝贵的信息。
有志贺氏痢疾的粪便样本
早在1998年，也许就有第一篇科学文章发表在PNAS杂志上，证明16世纪古代个体的人类牙髓中存在鼠疫耶尔森氏菌（引起鼠疫）。
image.png
Drancourt et al., Proc Natl Acad Sci U S A. 1998 Oct 13;95(21):12637–40. doi: 10.1073/pnas.95.21.12637.
随着第二代测序 (NGS) 技术的发展，从数百个古代样本中筛选微生物群落成为可能。这导致了许多惊人的发现。例如，证明鼠疫耶尔森氏菌在青铜器时代就已经存在，而不仅仅是像以前认为的那样在中世纪存在。
image.png
Rasmussen et al., Cell. 2015 Oct 22;163(3):571–82. doi: 10.1016/j.cell.2015.10.009, PMID: 26496604
从那时起，许多研究报告了在人类和动物的古代样本中发现的微生物病原体。虽然在这里我不是质疑这些发现，因为它们中的大多数都提供了彻底的验证，但我仍然想表达一个谨慎的态度。下面，我将演示在任何随机宏基因组样本中，如果病原体的存在毫无意义，错误地“检测”鼠疫耶尔森菌是多么容易。
覆盖深度可能产生误导
通常，参加科学aDNA研讨会和会议时，我听到以下报告古代样本中病原体发现的常规方式：
检测到鼠疫耶尔森氏菌，比对到~5000 reads，1.2X基因组覆盖率
Yersinia pestis pathogen detected, ~5000 reads mapped, 1.2X genome coverage
这里，“reads”是指通过测序技术读到的其核苷酸序列的小DNA片段。然后，这些“reads”与相关生物体的参考基因组比对。在上述示例中，约5000个小DNA“reads”与鼠疫耶尔森菌参考基因组比对。这相当多，而且似乎是一个很好的证据，证明样本中确实存在鼠疫耶尔森菌。“1.2X基因组覆盖率”通常指鼠疫耶尔森菌参考基因组中的任何随机位置平均覆盖约1.2次“reads”。这个数值越大（更具体地说是覆盖深度），测序和比对的质量越好，样本中存在鼠疫耶尔森菌的可能性也越高。
然而，在处理宏基因组数据分析时，即在对多个微生物进行同时比对时，覆盖深度可能是一个误导性指标。更确切地说，使用来自多个生物体的序列很难确保与特定微生物的参考基因组比对的序列确实来自该物种，而不是来自另一个物种。aDNA的降解和损坏，以及与参考基因组的进化差异，不完整和受污染的数据集，现代污染，跨生物体的 DNA 保守性和水平基因转移，使将序列可靠地分配给特定生物体（包括鼠疫杆菌等病原体）变得复杂。
image.png
具有相同覆盖深度的两种比对方案： a) 在同一个区域比对的reads，b) 均匀分布在参考基因组中的reads
在上图中，我展示了当短序列（蓝色）与参考基因组（黑色）比对时的两种情况。假设序列的长度为 L，为简单起见，参考基因组的长度为 4*L。当使用例如IGV （Integrative Genomics Viewer）工具可视化比对时，经常会遇到所有reads（序列）都在同一个区域上比对的情况。这是一个不幸的场景，表明序列和比对存在问题。另一种情况，当reads在参考基因组中均匀分布时，我们希望看到，因为这是成功比对的证据。但是，如上图所示，对于这两种截然不同的比对场景，计算覆盖深度可能会导致相同的值。覆盖深度度量的更好替代方法是覆盖广度（覆盖均匀度，evenness of of coverage），它简单地计算至少被一个read序列覆盖的参考基因组的一小部分。
随机样本中检测鼠疫耶尔森氏菌
为了演示如何错误地检测鼠疫耶尔森氏菌，让我们使用Broad Institute 的 DIABIMMUNE数据库（3个国家队列）从现代婴儿中随机抽取一份粪便样本。
从Broad Institure的DIABIMMUNE项目中选取的随机宏基因组样本
现在，我们使用Bowtie2序列比对工具将样本中的序列与鼠疫耶尔森氏菌的C092参考基因组进行比对。然而，首先，我们必须为鼠疫耶尔森氏菌的参考基因组建立一个Bowtie2索引。
samtools faidx Y.pestis.fna.gz
bowtie2-build --large-index Y.pestis.fna.gz Y.pestis.fna.gz --threads 4
time bowtie2 --large-index -x Y.pestis.fna.gz --end-to-end --threads 4 \
--very-sensitive -U G69146_pe_1.fastq.gz | samtools view -bS -q 1 -h -@ 4 -\
> G69146_pe_1.fastq.gz.aligned_to_Y.pestis.bam
在这里，我们惊讶地发现19 987条read序列可以比对到鼠疫耶尔森氏菌参考基因组。比对到的read数目很多，但是毫无意义，因为这是一个现代婴儿的随机样本，很难相信鼠疫耶尔森菌会存在于其中。请注意，通过说“19 987个比对到的read”，我们使用覆盖深度作为病原体检测的标准。上述分析中的一个明显的缺陷是，我们仅提供了鼠疫耶尔森氏菌病原体参考基因组用于比对，而这是一个应该包含大量人类DNA的人类粪便样本。因此，许多与鼠疫耶尔森氏菌参考基因组比对的read实际上可能来自人类，并且它们可能已比对到保守区域/直系同源区域。为了解决这个问题，让我们将人类参考基因组添加到鼠疫耶尔森氏菌参考基因组中，然后再次执行比对。
wget ftp://hgdownload.cse.ucsc.edu/goldenPath/hg38/bigZips/hg38.fa.gz
cat Ypestis.fna.gz hg38.fa.gz > Ypestis_plus_hg38.fna.gz
bowtie2-build --large-index Ypestis_plus_hg38.fna.gz Ypestis_plus_hg38.fna.gz --threads 4
bowtie2 --large-index -x Ypestis_plus_hg38.fna.gz --end-to-end --threads 4
--very-sensitive -U G69146_pe_1.fastq.gz | samtools view -bS -q 1 -h -@ 4 - \
> G69146_pe_1.fastq.gz_Aligned_to_Ypestis_plus_hg38_bowtie2.bam
samtools view G69146_pe_1.fastq.gz_Aligned_to_Ypestis_plus_hg38_bowtie2.bam \
NZ_CP009973.1 | wc -l
现在，我们可以看到“只有”8148个read与鼠疫耶尔森菌参考基因组对齐上，而它们可以选择与人类hg38参考基因组进行对齐。随着所用数据库（2个基因组而不是1个基因组）规模的增加，分析变得更加稳健，因为在这里我们不强制read仅与鼠疫耶尔森菌对齐。然而，我们仍然发现许多与鼠疫耶尔森氏菌的明显假阳性比对。这些可能是在分类学上更接近鼠疫耶尔森菌而不是真核生物的微生物序列。由于没有提供（包含在数据库中的）微生物参考基因组，除了鼠疫杆菌之外，许多实际上来自非耶尔森氏菌微生物的序列被迫与鼠疫耶尔森氏菌比对。为了解决这个问题，我们需要下载所有可用的细菌参考基因组，随机选择其中的 10、100、1000 和 10 000 个，与鼠疫杆菌 + 人类参考基因组合并，并与越来越大的数据库进行比对。每次，我们都会检查唯一比对到鼠疫杆菌的read数量，同时记住read可能对其他微生物具有更高的亲和力。
随着更多基因组添加到数据库中，比对到鼠疫杆菌的read数量减少
我们可以看到，随着数据库规模的增长，与鼠疫耶尔森氏菌唯一比对的read数量逐渐减少，当我们将10 000 个随机细菌纳入数据库时，read数量为6。因此，我们得出结论，有限数目的数据库可能会导致微生物的假阳性检测，因为我们得到越来越多的read比对到鼠疫杆菌，这显然是一个假阳性发现，我们提供的竞争性比对的基因组太少。此外，富含病原微生物的数据库可能会导致古代和现代宏基因组样本中病原体的错误发现。
解决方案在哪里？
那么，为什么我们似乎在现代婴儿样本中检测到了鼠疫耶尔森菌，而它不应该存在？显然是因为我们简单地计算了比对到鼠疫杆菌的read数，因此使用覆盖深度作为检测标准。这些read可能来自另一种同源微生物，但无法分配给该参考微生物，因为其参考基因组未包含在数据库中。我们可以过滤掉这些假阳性结果吗？我们可以！通过使用比覆盖深度/read读取数更好的标准。如果我们使用覆盖范围的广度，即read在参考基因组中分布的均匀程度，我们可以很容易地将鼠疫耶尔森氏菌作为假阳性过滤掉。
IGV可视化与鼠疫耶尔森氏菌参考比对的read
当使用IGV工具将鼠疫耶尔森氏菌参考基因组的比对可视化时，假阳性变得明显，如上图所示。read不是沿基因组均匀分布，而是比对上几个区域（最可能是保守的）。因此，对于未来的研究，报导古代样本中真正存在鼠疫耶尔森氏菌的更可靠和更有说服力的方法如下（注意覆盖范围为 85%）：
检测到鼠疫耶尔森氏菌病原体，比对约 5000条read，1.2X genome coverage，覆盖了85%的参考基因组
ersinia pestis pathogen detected, ~5000 reads mapped, 1.2X genome coverage, 85% of reference genome covered
总结
在本文中，我们了解到对古代样本微生物群落的DNA进行测序，可以提供有关古代人类和动物普遍死亡的惊人信息。然而，对古代微生物病原体的可靠检测仍然非常具有挑战性。由于忘记应用多级验证，人们很容易以错误的发现告终，因此需要对每个古代样本进行量身定制的仔细分析，以准确筛选其是否真正存在古代微生物病原体。
附上github代码
