使用 Python 分析 23AndMe 数据获取遗传起源

发布于 2024-09-23 21:08:37 字数 19411 浏览 10 评论 0

从公开可用的基因数据我们可以学到什么？

你的 DNA 包含了关于你的主线，易患疾病以及复杂特性，包括身高、体重、五官和行为等丰富的信息。使用来自 23andMe，一家直接面向消费者的遗传学公司，的公众可获取数据，我们将展示如何确定在网上找到的来自 23andMe 的一份匿名样本的祖先。

在本教程的最后，你会学到并且能够做到三件事： 1. 如何理解和分析来自两个不同源的基因数据 (23andMe 和千人基因组计划)。 2. 主成分分析的基本知识，一个用于无监督聚类、降维和数据探索的工具。 3. 如何使用主成分分析来预测一个匿名基因数据集的祖先。

首先，你可以下载一组公开可用的 23andMe 数据。这个文件的大小大概是 23Mb。我们会使用 python 的 urllib 包来下载文件。

import urllib

u = urllib.URLopener()
u.retrieve("https://s3-eu-west-1.amazonaws.com/dm-23andme-v3/dm_23andme_v3_110219.txt", "dm_23andme_v3_110219.txt")

我们将使用 pandas 库来加载数据为 dataframe，并看看形状和内容。这个文件有一些以‘#’开头的标题行。我们将跳过这些行。

import pandas as pd
anon = pd.read_table("dm_23andme_v3_110219.txt", sep = "\t", comment = "#", header = None)

检查数据有多少行和列，以及数据是咋样的：

print("The 23andMe datset has {} rows and {} columns.".format(anon.shape[0], anon.shape[1]))
print(anon.head())

The 23andMe datset has 966977 rows and 4 columns.
            0  1       2   3
0   rs4477212  1   72017  AA
1   rs3094315  1  742429  AG
2   rs3131972  1  742584  AG
3  rs12124819  1  766409  --
4  rs11240777  1  788822  AG

在本教程中，你只需要使用这个 dataframe 约 1000 行即可 —— 使用 pandas 来切出前 1000 行：

anon = anon.iloc[0:1000,:]

基因数据是如何组织的

人类基因组是一个由 A / T / G / C 组成的超过 30 亿的字母（称为“碱基”）的序列。这 30 亿个字母序列被分成 23 个不同的段，称为染色体。其中最大的染色体 1 号，有近 2.5 亿个碱基长。我们的基因有两个副本（一个来自母亲，一个来自父亲），因此，我们每个人拥有超过 60 亿个字母的 DNA 序列。

下面的示意图显示了单个染色体，以及百万个碱基对是如何被压缩并打包到染色体中的。

在你刚刚打印的 dataframe 中，第二和第三列对应于染色体及位置。第四列包含在该染色体和位置的两个碱基。第一列称为 rsid，它是基因变异数据库中使用的标识符。

所有这些列的命名信息位于我们加载的文件的头部，但在上面，我们用 comment = "#" 跳过了它。

让我们重命名相应的数据集的列：

anon.columns = ["rsid", "chrom", "pos", "genotype"]
print(anon.head())

rsid chrom     pos genotype
0   rs4477212     1   72017       AA
1   rs3094315     1  742429       AG
2   rs3131972     1  742584       AG
3  rs12124819     1  766409       --
4  rs11240777     1  788822       AG

处理来自三大人群的公开可用基因组数据

接下来，我们将下载来自于三个主要人群的样本，它是千人基因组计划编目世界各地的人类遗传变异序列化的一部分。这三组是 CEU (来自犹他州的北欧人)，YRI (尼日利亚伊巴丹的约鲁巴人) 以及 CHB/JPT (中国汉人和东京的日本人)。

有关处理基因数据的第一个重要教训是，它可能会非常大。序列化涵盖 60 亿左右个位点的单个人类基因组生成的原始数据通常最少 30Gb。通过减少到已知在人类群体内各异的位点子集(23andMe 使用大约 1 百万个位点)，可以减少大小。我们所使用的匿名的 23andMe 数据只有 23Mb。

我们将使用 tabix，这一个用以按顺序索引和查询大的基因组数据集的流行工具，来从千人基因组文件中提取来自我们所感兴趣的位点的数据。tabix 为我们节省了大量的时间，并让大文件更容易访问 —— 更多关于 tabix 的信息，见这里： http://www.htslib.org/doc/tabix.html 。

我们会下载两个文件 - 一个是基因数据本身，另一个是由 tabix 生成的索引，它可以让查询更快更省内存。这个文件被压缩，大小约为 1Gb (969Mb)，因此，取决于你的互联网连接，可能要花几分钟来下载它。你不需要解压缩它 —— 如果你决定解压缩，看看里面有啥，那么解压缩后大小大概为 7Gb，因此，不要尝试一次性看完整个文件！

你需要运行命令'pip install -user pytabix'，才能在 python 中使用 tabix。

u = urllib.URLopener()
u.retrieve("ftp://ftp.1000genomes.ebi.ac.uk/vol1/ftp/pilot_data/paper_data_sets/a_map_of_human_variation/low_coverage/snps/YRI.low_coverage.2010_09.genotypes.vcf.gz", "YRI.low_coverage.2010_09.genotypes.vcf.gz")
u.retrieve("ftp://ftp.1000genomes.ebi.ac.uk/vol1/ftp/pilot_data/paper_data_sets/a_map_of_human_variation/low_coverage/snps/YRI.low_coverage.2010_09.genotypes.vcf.gz.tbi", "YRI.low_coverage.2010_09.genotypes.vcf.gz.tbi")

import tabix

YRI_file = "YRI.low_coverage.2010_09.genotypes.vcf.gz"
yri = tabix.open(YRI_file)

基因坐标通常是基于 0 的，意味着一个序列中的第一个碱基将位于位置 0 和位置 1“之间”。下面，我们将通过查询 1 号染色体位置 742428 和 742429（在序列中，将会是第 742429 个碱基）之间的位点，检查 YRI 数据集，并将其赋给 rs3094315 变量。下面是显式 0 碱基 vs 1 碱基坐标之间差异的简单的图形，来自这篇 biostars 文章： https://www.biostars.org/p/84686/

rs3094315 是我们的 23andMe 数据集中的第二个变量，在 YRI 数据集中，它也被序列化：

rs3094315 = yri.query("1", 742428, 742429)  # an iterator object
print(rs3094315.next())

['1', '742429', 'rs3094315', 'G', 'A', '.', 'PASS', 'AA=g;AC=46;AN=118;DP=237;HM2;GP=1:752566;BN=103',
'GT:DP:CB', '0|1:3:SM', '1|1:4:MB', '1|0:5:SMB', '0|1:2:SMB', '1|0:6:SMB', '1|0:7:SMB', '0|1:4:SMB',
'0|0:4:SMB', '1|1:0:SMB', '1|1:12:SMB', '0|1:4:SMB', '0|1:2:SMB', '1|0:4:MB', '0|0:7:SMB',
'1|1:4:SMB', '0|0:4:SMB', '1|1:6:SMB', '0|0:5:SMB', '0|1:4:SMB', '1|1:5:MB', '0|0:6:SMB',
'0|0:5:SMB', '0|1:1:SMB', '1|1:2:SMB', '0|0:9:SMB', '0|0:1:SMB', '0|0:10:SMB', '0|1:9:SMB',
'1|0:9:SMB', '0|1:2:SMB', '0|1:8:SMB', '1|1:4:SMB', '0|1:9:SMB', '0|0:2:SMB', '1|0:5:SMB',
'0|1:2:SMB', '0|0:3:SMB', '0|0:0:SMB', '0|0:4:SMB', '0|1:7:SMB', '1|0:3:SM', '0|0:2:SMB',
'0|0:0:SMB', '0|1:9:SMB', '0|1:4:SMB', '0|0:1:SMB', '0|0:1:SMB', '0|0:1:SMB', '0|0:3:SMB',
'1|1:2:SMB', '0|0:2:SMB', '1|0:4:SMB', '0|0:2:SMB', '0|0:2:SMB', '1|0:2:SMB', '0|0:0:SMB',
'1|0:2:SMB', '1|1:3:SMB', '1|0:4:SMB']

我们要如何阅读上面的输出呢？正如你可能已经注意到的，我们的 23andMe 数据组织的方式相当不同。

千人基因组文件保存在 VCF (variant call format) 文件中，这与 23andMe 提供的格式不同，因此，我们将需要进行一些数据转换，以使得这两个可以一起工作。

从左到右，前五个列编码染色体、位置、rsid、参考碱基 (这是我们大部分时间希望看到的字母)、替代碱基 (这是“变量”字母，它存在于人群中的一些成员中)。更多关于 VCF 中包含的信息的信息，可以在这里找到： http://www.1000genomes.org/wiki/Analysis/vcf4.0/ 。这是遗传学中最常见的文件格式之一，而除了基因类型本身外的许多额外数据可以用于过滤数据，从而提高数据质量。要看看 VCF 中的数据，你可以创建一份你所下载的压缩的 VCF 文件的拷贝，解压它，然后在命令行上使用“less”来看一看。一旦解压，文件大小约为 7Gb。

VCF 的最后一列和（这个列表的最后一条) 显示在这个最终文件中包含的 59 个 YRI 样本的基因类型。列表中的最后 59 条都遵循相似的格式来报告用冒号分隔的基因类型。现在，我们只对第一块信息感兴趣，它是基因类型 (它将是 0|0, 1|0, 0|1, or 1|1)。

正如你可以在上面看到的，在这个位点上，我们的匿名样本基因类型是 AG，它对应于上面列出的格式的中的 1|0。0 是参考(G)，而 1 是替代(A)。有几个 YRI 样本是 1|1，这意味着它们在这个位点是 AA，而 0|0 的 YRI 样本在这个位点上是 GG。在我们更进一步之前，通过添加一列展示在千人基因组 VCF 文件中使用的格式为 0 和 1 的基因类型，我们想要将我们的 23andMe 数据转换成更紧密匹配我们刚刚下载的人群参考格式。(这比 23andMe 使用的格式更常用)。

def convert_anon_genotype(chrom, pos, genotype, vcf_tabix):
    site = vcf_tabix.query(chrom, pos - 1, pos)
    try:
        row = site.next() # this will throw an error (which is caught by 'except' on the next line) if the site we queried is not in the tabix file
    except StopIteration:
        return None # put None in the dataframe if we are missing this genotype in 1000 Genomes
    ref = row[3]
    alt = row[4]
    if genotype == ref+ref:
        return("0|0")
    elif (genotype == ref+alt) | (genotype == alt+ref):
        return("0|1")
    elif genotype == alt+alt:
        return("1|1")
    else: # missing genotype, or incorrect annotation, we assume ref/ref
        return("0|0")

genotypes_1kg_format = []
for chrom, pos, genotype in zip(anon['chrom'], anon['pos'], anon['genotype']):
    genotypes_1kg_format.append(convert_anon_genotype(str(chrom), pos, genotype, yri))

anon['genotype_1kg_format'] = genotypes_1kg_format
print(anon.head())
print(anon.shape)

rsid chrom     pos genotype genotype_1kg_format
0   rs4477212     1   72017       AA                None
1   rs3094315     1  742429       AG                 0|1
2   rs3131972     1  742584       AG                 0|1
3  rs12124819     1  766409       --                None
4  rs11240777     1  788822       AG                None
(1000, 5)

从基因类型创建一个特征空间

记住，我们的首要目标是预测祖先，因此想要将基因数据表示为一个特征空间。我们将构建一个 dataframe，其中，行是样本（人），而我们将每一条染色体和每一个位置上的基因类型当成一个单独的特征。由于我们将匿名 23andMe 样本限制为前一千个位点，因此我们的 dataframe 将有两个列，描述每个样本的人群和名字，加上另外 1000 个列 (每个特性一个)。

# make a data frame with one row for each of the YRI samples
yri_genotypes = pd.DataFrame({"sample": ["YRI" + str(i) for i in range(1, 60)], "population": "YRI"})
print(yri_genotypes.head())


  population sample
0        YRI   YRI1
1        YRI   YRI2
2        YRI   YRI3
3        YRI   YRI4
4        YRI   YRI5

要获取每个人的基因类型，我们要写一个使用 pytabix 的函数：

# extract genotype information for a set of sites
def extract_genotype(chrom, pos, vcf_tabix):
    site = vcf_tabix.query(chrom, pos - 1, pos)
    try:
        g = site.next()[9:]
    except StopIteration:
        return None # put None in the dataframe if we are missing this genotype in 1000 Genomes
    g = [i.split(":")[0] for i in g]  # if present in 1000 genomes, get the genotypes
    return(g)

for rsid, chrom, pos in zip(anon['rsid'], anon['chrom'], anon['pos']):
    g = extract_genotype(str(chrom), pos, yri)
    yri_genotypes[rsid] = g

print("The dataframe including all of the samples from the YRI population has {}
samples and {} genotypes.".format(yri_genotypes.shape[0], yri_genotypes.shape[1] - 2))

The dataframe including all of the samples from the YRI population has 59 samples and 1000 genotypes.

所有带有'None'的列都是由 23andMe 测序，但在我们的 YRI 数据集中没有的位点。稍后，我们将丢弃所有在我们任意人群中缺失数据的位点。

print(yri_genotypes.iloc[0:10, 0:7])

population sample rs4477212 rs3094315 rs3131972 rs12124819 rs11240777
0        YRI   YRI1      None       0|1       0|1       None       None
1        YRI   YRI2      None       1|1       1|0       None       None
2        YRI   YRI3      None       1|0       1|0       None       None
3        YRI   YRI4      None       0|1       0|1       None       None
4        YRI   YRI5      None       1|0       1|0       None       None
5        YRI   YRI6      None       1|0       0|0       None       None
6        YRI   YRI7      None       0|1       0|1       None       None
7        YRI   YRI8      None       0|0       0|0       None       None
8        YRI   YRI9      None       1|1       1|1       None       None
9        YRI  YRI10      None       1|1       1|0       None       None

现在，我们要下载数据集，然后为 CEU (犹他州的欧洲人) 和 CHB/JPT (中国汉族/日本人) 创建数据帧。再次，应该在命令行使用 curl 或者 wget 来下载这些文件和 tabix 索引：

u = urllib.URLopener()
u.retrieve("ftp://ftp.1000genomes.ebi.ac.uk/vol1/ftp/pilot_data/paper_data_sets/a_map_of_human_variation/low_coverage/snps/CEU.low_coverage.2010_09.genotypes.vcf.gz", "CEU.low_coverage.2010_09.genotypes.vcf.gz")
u.retrieve("ftp://ftp.1000genomes.ebi.ac.uk/vol1/ftp/pilot_data/paper_data_sets/a_map_of_human_variation/low_coverage/snps/CEU.low_coverage.2010_09.genotypes.vcf.gz.tbi", "CEU.low_coverage.2010_09.genotypes.vcf.gz.tbi")
u.retrieve("ftp://ftp.1000genomes.ebi.ac.uk/vol1/ftp/pilot_data/paper_data_sets/a_map_of_human_variation/low_coverage/snps/CHBJPT.low_coverage.2010_09.genotypes.vcf.gz", "CHBJPT.low_coverage.2010_09.genotypes.vcf.gz")
u.retrieve("ftp://ftp.1000genomes.ebi.ac.uk/vol1/ftp/pilot_data/paper_data_sets/a_map_of_human_variation/low_coverage/snps/CHBJPT.low_coverage.2010_09.genotypes.vcf.gz.tbi", "CHBJPT.low_coverage.2010_09.genotypes.vcf.gz.tbi")

# europeans from utah
CEU_file = "CEU.low_coverage.2010_09.genotypes.vcf.gz"
ceu = tabix.open(CEU_file)

number_ceu_samples = len(ceu.query("1", 742428, 742429).next()[9:])

ceu_genotypes = pd.DataFrame({"sample": ["CEU" + str(i) for i in range(1, number_ceu_samples + 1)], "population": "CEU"})

# Han chinese and Japanese
CHBJPT_file = "CHBJPT.low_coverage.2010_09.genotypes.vcf.gz"
chbjpt = tabix.open(CHBJPT_file)

number_chbjpt_samples = len(chbjpt.query("1", 742428, 742429).next()[9:])

chbjpt_genotypes = pd.DataFrame({"sample": ["CHBJPT" + str(i) for i in range(1, number_chbjpt_samples + 1)], "population": "CHBJPT"})

for rsid, chrom, pos in zip(anon['rsid'], anon['chrom'], anon['pos']):
    yri_genotypes[rsid] =  extract_genotype(str(chrom), pos, yri)
    ceu_genotypes[rsid] =  extract_genotype(str(chrom), pos, ceu)
    chbjpt_genotypes[rsid] =  extract_genotype(str(chrom), pos, chbjpt)

此时，我们有了来自我们感兴趣的三个人群的三个独立的 dataframe，包含它们在 1000 个不同位置的基因类型。我们的目标是看看忽略原有的人群标签，仅仅基于其遗传信息，我们是否能够将这些样本分离到各个人群中。

genotypes = yri_genotypes.copy()
genotypes = genotypes.append(ceu_genotypes, ignore_index=True)
genotypes = genotypes.append(chbjpt_genotypes, ignore_index=True)

print("Now the genotypes data frame has {} samples and {} genotypes").format(genotypes.shape[0], genotypes.shape[1]-2)

Now the genotypes data frame has 179 samples and 1000 genotypes

使用主成分分析的无监督聚类

要尝试将我们的数据分离到人群，我们将使用主成分分析 (PCA)，这是一个无监督方法，它对于将来自于大的特征空间的数据整合成一个较小集合的特性，以捕获样本间的大部分差异非常有用。

PCA 需要一个高维的特征空间 (在本例中，是 10 个特性)，并使用所有特性的线性组合创建一个捕获样本间差异的“主成分”。我们的数据可以被想象成高维空间中的点云。第一个主成分是通过这个高维空间的一条线。第二个主成分必须正交（成 90 度角）于第一个主成分，并构造以捕获第二大的变异量。

另一种想象主成分的方法是，作为高维数据的“投影”或“影子”到低维空间。下面的图显示了一个在二维空间中可视化思维特性空间的例子。

这张图来自于一个 scikit-learn 教程，它使用非常流行的 Iris 数据集的数据： http://scikit-learn.org/stable/auto_examples/decomposition/plot_pca_vs_lda.html

现在，我们将在我们的基因数据集上使用 PCA，尝试可视化我们的样本间的差异。我们会使用 scikit-learn 库来进行 PCA。下面，我们将基因类型数据转换成连续值。有参考字母的两份拷贝的样本设为 0，一份参考拷贝和一份替代拷贝的设为 0.5，有两份替代拷贝的设为 1.0。另外，我们会丢弃千人基因组的人群面板中存在'None'的任何变量，因为它们并不含信息。

from sklearn.decomposition import PCA
pca = PCA(n_components = 2)

genotypes_only = genotypes.copy().iloc[:, 2:]  # we make a copy here, otherwise pandas will gripe at us!
genotypes_only[genotypes_only == "1|1"] = 1
genotypes_only[genotypes_only == "0|1"] = 0.5
genotypes_only[genotypes_only == "0/1"] = 0.5
genotypes_only[genotypes_only == "1|0"] = 0.5
genotypes_only[genotypes_only == "0|0"] = 0.0

# remove variants with None
genotypes_only = genotypes_only.dropna(axis=1)

import matplotlib.pyplot as plt
%matplotlib inline

pca.fit(genotypes_only)
pc = pca.transform(genotypes_only)

plt.figure(figsize=(10,6))
plt.scatter(pc[:, 0], pc[:, 1])
plt.title('PCA of 1000 23andMe SNPs')
plt.xlabel('PC1')
plt.ylabel('PC2')
plt.show()

使用来自我们的 23andMe 数据集的少于一千个位点（记住，这只是 23andMe 数据集中总共近 1 百万个位点的一小部分)，我们看到，我们的数据被非常清晰地分到三个集群中。如果我们把人群标签加回去，那么可以看到，这些集群代表了我们开始的三个人群。对于这个特定的分析，PCA 一个最强大的部分是，一开始，我们无需指定寻找的集群数。

import numpy as np

plt.figure(figsize=(10,6))

for c, pop in zip("rby", ["YRI", "CEU", "CHBJPT"]):
    plt.scatter(pc[np.where(genotypes['population'] == pop), 0], pc[np.where(genotypes['population'] == pop), 1], c = c, label = pop)
plt.title('PCA of 1000 23andMe SNPs')
plt.xlabel('PC1')
plt.ylabel('PC2')
plt.legend(loc = 'upper left')
plt.show()

分析还告诉我们的最重要的事情之一是，如果我们对从我们的基因数据中学到除了人们的祖先来自哪里的其他东西感兴趣，那么需要控制人群结构 (例如，利用人群/祖先作为协变量)。当我们尝试了解对于常见疾病（例如，糖尿病，肥胖，心脏疾病，遗传性癌症) 的易感性，基因所扮演的角色，那么我们必须确保我们的观察组（具有某疾病的人) 和控制组（不具有某疾病的人) 的祖先并不会有很大差异。

对 PCA 图另一个有趣的观察室，在第一个主成分中，约鲁巴集群（红色）从中国人/日本人和欧洲人集群中分离开来。这是在非洲的现代人类起源和随后迁移出非洲并进入欧洲/亚洲的反映。然后，第二个主成分反映了后续的迁移导致了中国人/日本人集群和欧洲人集群之间的基因差异。

在这张地图中，我们的匿名 23andMe 样本落在哪里呢？23andMe 数据格式与 VCF 不同，因此，我们会写一个快速的函数来检查 VCF (使用 tabix)，以映射 23andMe 基因类型，这将把字母对应到我们现在习惯了的 0|0，0|1 和 1|1。

# keep only the genotypes used in our PCA above
anon = anon.loc[anon['rsid'].isin(genotypes_only.columns.values), :]  # only keep the 23andMe data where we have no missing data in 1000 genomes

anon_genotypes = anon.copy()["genotype_1kg_format"]

anon_genotypes[anon_genotypes == "1|1"] = 1
anon_genotypes[anon_genotypes == "0|1"] = 0.5
anon_genotypes[anon_genotypes == "1|0"] = 0.5
anon_genotypes[anon_genotypes == "0|0"] = 0.0
#anon_genotypes[anon_genotypes == None] = 0.0
anon_genotypes = anon_genotypes.reshape(1,-1) # reshape, otherwise sci-kit learn will throw a deprecation warning

# assume any missing data in our 23andme sample is ref/ref
#anon_genotypes[anon_genotypes is None] = "0|0"

anon_pca = pca.transform(anon_genotypes)  # pca was fit on the 1000 genomes data and we use it to transform the anonymous genotypes

在上面的第一行中，我们丢弃了那些不属于千人基因组位点的 23andMe 位点，因为为了在主成分空间中计算匿名样本的值，我们必须有两个数据集（匿名 23andMe 和千人基因组）之间的相同的基础数据。

print(anon.head())

rsid chrom     pos genotype genotype_1kg_format
1  rs3094315     1  742429       AG                 0|1
2  rs3131972     1  742584       AG                 0|1
6  rs4970383     1  828418       CC                 0|0
7  rs4475691     1  836671       CC                 0|0
8  rs7537756     1  844113       AA                 0|0

plt.figure(figsize=(10,6))

for c, pop in zip("rgb", ["YRI", "CEU", "CHBJPT"]):
    plt.scatter(pc[np.where(genotypes['population'] == pop), 0], pc[np.where(genotypes['population'] == pop), 1], c = c, label = pop)

# take the code above and add in the anonymous sample

plt.scatter(anon_pca[0,0], anon_pca[0,1], c = "yellow", label = "Anonymous 23andMe Sample", marker = (5,1,0), s = 200)

plt.title('PCA of 1000 23andMe SNPs')
plt.xlabel('PC1')
plt.ylabel('PC2')
plt.legend(loc = 'upper left')
plt.show()

正如我们在上面图中所看到的那样，我们的匿名样本（黄色的星星）正对欧洲人集群。在实践中，一个人群参考板比这里使用的三个具有多得多的子人群。本教程展示了如何用几个 Python 段进行简单的祖先分析，但我们可以更进一步，使用一个更复杂的模型来预测祖先。

千人基因组计划的数据集的最终版本拥有来自 26 个不同人群的数据 ( http://www.1000genomes.org/category/population/) 。在这种情况下，将来自千人基因组计划的基因数据当成训练数据，例如使用高斯混合模型或 K-近邻，你就可以构建一个分类器来预测比例祖先 (例如，1/2 的爱尔兰，1/2 中国汉人)。

除了祖先，基因数据可以用于预测疾病风险，药物副作用，甚至构建脸部模型 ( https://www.newscientist.com/article/mg22129613-600-genetic-mugshot-recreates-faces-from-nothing-but-dna/) 。

回来看看关于使用 Python 分析基因数据的更多免费教程，或者注册我们的数据科学课程！

分享到QQ

分享到微博