Python-生物信息学秘籍第三版-三-Python 生物信息学秘籍第三版（三）第八章：系统发育学系统发育学是应用分

Python 生物信息学秘籍第三版（三）

原文：annas-archive.org/md5/9694cf42f7d741c69225ff1cf52b0efe

译者：飞龙

协议：CC BY-NC-SA 4.0

第八章：系统发育学

系统发育学是应用分子测序技术来研究生物之间的进化关系。通常用系统发育树来表示这一过程。从基因组数据中计算这些树是一个活跃的研究领域，并具有许多现实世界中的应用。

本书将把之前提到的实际操作方法推向一个新高度：这里的大部分内容都受到了关于埃博拉病毒的最新研究的启发，研究了最近在非洲爆发的埃博拉疫情。该研究名为Genomic surveillance elucidates Ebola virus origin and transmission during the 2014 outbreak，由Gire et al.完成，发表于Science期刊。你可以在pubmed.ncbi.nlm.nih.gov/25214632/阅读到这篇文章。在这里，我们将尝试采用类似的方法，以便得出与该论文相似的结果。

本章我们将使用 DendroPy（一个系统发育学库）和 Biopython。bioinformatics_phylo Docker 镜像包含了所有必需的软件。

本章我们将涵盖以下内容：

为系统发育分析准备数据集
对齐遗传和基因组数据
比较序列
重建系统发育树
递归地玩转树状图
可视化系统发育数据

为系统发育分析准备数据集

在这个步骤中，我们将下载并准备用于分析的数据集。该数据集包含了埃博拉病毒的完整基因组。我们将使用 DendroPy 来下载和准备数据。

准备就绪

我们将从 GenBank 下载完整的基因组，这些基因组来自多个埃博拉疫情爆发，其中包括 2014 年疫情中的多个样本。请注意，导致埃博拉病毒病的病毒种类有很多；2014 年疫情中的主要病毒是 EBOV（即正式命名为扎伊尔埃博拉病毒），这是最常见的一种，但此病还由更多的埃博拉病毒属物种引起；另外四种物种的基因组序列也已被测序。你可以在en.wikipedia.org/wiki/Ebolavirus阅读更多信息。

如果你已经完成了前几章的内容，看到这里可能会对涉及的数据量感到担忧；但这完全不成问题，因为这些是大约 19 kbp 大小的病毒基因组。因此，我们的约 100 个基因组实际上相当小巧。

为了进行这个分析，我们需要安装 dendropy。如果你使用的是 Anaconda，请执行以下操作：

conda install –c bioconda dendropy

像往常一样，这些信息可以在相应的 Jupyter Notebook 文件中找到，文件位于Chapter07/Exploration.py。

如何操作...

请查看以下步骤：

首先，让我们使用 DendroPy 指定我们的数据源，如下所示：

import dendropy
from dendropy.interop import genbank
def get_ebov_2014_sources():
    #EBOV_2014
    #yield 'EBOV_2014', genbank.GenBankDna(id_range=(233036, 233118), prefix='KM')
    yield 'EBOV_2014', genbank.GenBankDna(id_range=(34549, 34563), prefix='KM0')
def get_other_ebov_sources():
    #EBOV other
    yield 'EBOV_1976', genbank.GenBankDna(ids=['AF272001', 'KC242801'])
    yield 'EBOV_1995', genbank.GenBankDna(ids=['KC242796', 'KC242799'])
    yield 'EBOV_2007', genbank.GenBankDna(id_range=(84, 90), prefix='KC2427')
def get_other_ebolavirus_sources():
    #BDBV
    yield 'BDBV', genbank.GenBankDna(id_range=(3, 6), prefix='KC54539')
    yield 'BDBV', genbank.GenBankDna(ids=['FJ217161']) #RESTV
    yield 'RESTV', genbank.GenBankDna(ids=['AB050936', 'JX477165', 'JX477166',  'FJ621583', 'FJ621584', 'FJ621585'])
    #SUDV
    yield 'SUDV', genbank.GenBankDna(ids=['KC242783', 'AY729654', 'EU338380', 'JN638998', 'FJ968794', 'KC589025', 'JN638998'])
    #yield 'SUDV', genbank.GenBankDna(id_range=(89, 92), prefix='KC5453')
    #TAFV
    yield 'TAFV', genbank.GenBankDna(ids=['FJ217162'])

在这里，我们有三个函数：一个用于检索最新 EBOV 爆发的数据，另一个用于检索以前 EBOV 爆发的数据，第三个用于检索其他物种爆发的数据。

请注意，DendroPy 的 GenBank 接口提供了几种不同的方式来指定要检索的记录列表或范围。一些行已被注释掉，包括下载更多基因组的代码。对于我们的目的，我们将下载的子集已足够。

现在，我们将创建一组 FASTA 文件；我们将在这里以及未来的食谱中使用这些文件：

other = open('other.fasta', 'w')
sampled = open('sample.fasta', 'w')
for species, recs in get_other_ebolavirus_sources():
    tn = dendropy.TaxonNamespace()
    char_mat = recs.generate_char_matrix(taxon_namespace=tn,
        gb_to_taxon_fn=lambda gb: tn.require_taxon(label='%s_%s' % (species, gb.accession)))
    char_mat.write_to_stream(other, 'fasta')
    char_mat.write_to_stream(sampled, 'fasta')
other.close()
ebov_2014 = open('ebov_2014.fasta', 'w')
ebov = open('ebov.fasta', 'w')
for species, recs in get_ebov_2014_sources():
    tn = dendropy.TaxonNamespace()
    char_mat = recs.generate_char_matrix(taxon_namespace=tn,
        gb_to_taxon_fn=lambda gb: tn.require_taxon(label='EBOV_2014_%s' % gb.accession))
    char_mat.write_to_stream(ebov_2014, 'fasta')
    char_mat.write_to_stream(sampled, 'fasta')
    char_mat.write_to_stream(ebov, 'fasta')
ebov_2014.close()
ebov_2007 = open('ebov_2007.fasta', 'w')
for species, recs in get_other_ebov_sources():
    tn = dendropy.TaxonNamespace()
    char_mat = recs.generate_char_matrix(taxon_namespace=tn,
        gb_to_taxon_fn=lambda gb: tn.require_taxon(label='%s_%s' % (species, gb.accession)))
    char_mat.write_to_stream(ebov, 'fasta')
    char_mat.write_to_stream(sampled, 'fasta')
    if species == 'EBOV_2007':
        char_mat.write_to_stream(ebov_2007, 'fasta')
ebov.close()
ebov_2007.close()
sampled.close()

我们将生成几个不同的 FASTA 文件，包含所有基因组、仅 EBOV 基因组或仅包含 2014 年爆发的 EBOV 样本的文件。在本章中，我们将主要使用包含所有基因组的sample.fasta文件。

请注意使用dendropy函数创建 FASTA 文件，这些文件是通过转换从 GenBank 记录中提取的。FASTA 文件中每个序列的 ID 是通过一个 lambda 函数生成的，该函数使用物种和年份，以及 GenBank 的登录号。

让我们提取病毒中的四个基因（总共七个），如下所示：

my_genes = ['NP', 'L', 'VP35', 'VP40']
def dump_genes(species, recs, g_dls, p_hdls):
    for rec in recs:
        for feature in rec.feature_table:
            if feature.key == 'CDS':
                gene_name = None
                for qual in feature.qualifiers:
                    if qual.name == 'gene':
                        if qual.value in my_genes:
                            gene_name = qual.value
                    elif qual.name == 'translation':
                        protein_translation = qual.value
                if gene_name is not None:
                    locs = feature.location.split('.')
                    start, end = int(locs[0]), int(locs[-1])
                    g_hdls[gene_name].write('>%s_%s\n' % (species, rec.accession))
                    p_hdls[gene_name].write('>%s_%s\n' % (species, rec.accession))
                    g_hdls[gene_name].write('%s\n' % rec.sequence_text[start - 1 : end])
                    p_hdls[gene_name].write('%s\n' % protein_translation)
g_hdls = {}
p_hdls = {}
for gene in my_genes:
    g_hdls[gene] = open('%s.fasta' % gene, 'w')
    p_hdls[gene] = open('%s_P.fasta' % gene, 'w')
for species, recs in get_other_ebolavirus_sources():
    if species in ['RESTV', 'SUDV']:
        dump_genes(species, recs, g_hdls, p_hdls)
for gene in my_genes:
    g_hdls[gene].close()
    p_hdls[gene].close()

我们首先搜索第一个 GenBank 记录中的所有基因特征（请参阅第三章，下一代测序，或国家生物技术信息中心（NCBI）文档获取更多细节；虽然我们将在这里使用 DendroPy 而不是 Biopython，但概念是相似的），并将数据写入 FASTA 文件，以提取基因。我们将每个基因放入不同的文件中，只取两个病毒物种。我们还获取了转译的蛋白质，这些蛋白质在每个基因的记录中都有提供。

让我们创建一个函数，从比对中获取基本统计信息，如下所示：

def describe_seqs(seqs):
    print('Number of sequences: %d' % len(seqs.taxon_namespace))
    print('First 10 taxon sets: %s' % ' '.join([taxon.label for taxon in seqs.taxon_namespace[:10]]))
    lens = []
    for tax, seq in seqs.items():
        lens.append(len([x for x in seq.symbols_as_list() if x != '-']))
    print('Genome length: min %d, mean %.1f, max %d' % (min(lens), sum(lens) / len(lens), max(lens)))

我们的函数采用DnaCharacterMatrix DendroPy 类，并统计分类单元的数量。然后，我们提取每个序列中的所有氨基酸（排除由-表示的缺失）来计算长度，并报告最小、平均和最大大小。有关 API 的更多细节，请查看 DendroPy 文档。

让我们检查一下 EBOV 基因组的序列并计算基本统计数据，如前所示：

ebov_seqs = dendropy.DnaCharacterMatrix.get_from_path('ebov.fasta', schema='fasta', data_type='dna')
print('EBOV')
describe_seqs(ebov_seqs)
del ebov_seqs

然后我们调用一个函数，得到 25 个序列，最小大小为 18,700，平均大小为 18,925.2，最大大小为 18,959。与真核生物相比，这是一个较小的基因组。

请注意，在最后，内存结构已被删除。这是因为内存占用仍然相当大（DendroPy 是一个纯 Python 库，在速度和内存方面有一些开销）。在加载完整基因组时，要小心内存使用。

现在，让我们检查另一个埃博拉病毒基因组文件，并统计不同物种的数量：

print('ebolavirus sequences')
ebolav_seqs = dendropy.DnaCharacterMatrix.get_from_path('other.fasta', schema='fasta', data_type='dna')
describe_seqs(ebolav_seqs)
from collections import defaultdict
species = defaultdict(int)
for taxon in ebolav_seqs.taxon_namespace:
    toks = taxon.label.split('_')
    my_species = toks[0]
    if my_species == 'EBOV':
        ident = '%s (%s)' % (my_species, toks[1])
    else:
        ident = my_species
    species[ident] += 1
for my_species, cnt in species.items():
    print("%20s: %d" % (my_species, cnt))
del ebolav_seqs

每个分类单元的名称前缀表明了物种，我们利用这一点来填充一个计数字典。

接下来详细介绍物种和 EBOV 的分类（图例中 Bundibugyo 病毒=BDBV，Tai Forest 病毒=TAFV，Sudan 病毒=SUDV，Reston 病毒=RESTV；我们有 1 个 TAFV，6 个 SUDV，6 个 RESTV 和 5 个 BDBV）。

让我们提取病毒中基因的基本统计信息：

gene_length = {}
my_genes = ['NP', 'L', 'VP35', 'VP40']
for name in my_genes:
    gene_name = name.split('.')[0]
    seqs =    
dendropy.DnaCharacterMatrix.get_from_path('%s.fasta' % name, schema='fasta', data_type='dna')
    gene_length[gene_name] = []
    for tax, seq in seqs.items():
        gene_length[gene_name].append(len([x for x in  seq.symbols_as_list() if x != '-'])
for gene, lens in gene_length.items():
    print ('%6s: %d' % (gene, sum(lens) / len(lens)))

这允许你概览基本的基因信息（即名称和平均大小），如下所示：

NP: 2218
L: 6636
VP35: 990
VP40: 988

还有更多...

这里的大部分工作可能可以通过 Biopython 完成，但 DendroPy 具有更多额外的功能，将在后续的步骤中进行探索。此外，正如你将发现的，它在某些任务（如文件解析）上更为强大。更重要的是，还有另一个 Python 库用于执行系统发育学分析，你应该考虑一下。它叫做 ETE，可以在etetoolkit.org/找到。

另见

美国疾病控制中心（CDC）在www.cdc.gov/vhf/ebola/history/summaries.xhtml上有一个关于埃博拉病毒病的很好的简介页面。
系统发育学中的参考应用是 Joe Felsenstein 的Phylip，可以在evolution.genetics.washington.edu/phylip.xhtml找到。
我们将在未来的步骤中使用 Nexus 和 Newick 格式（evolution.genetics.washington.edu/phylip/newicktree.xhtml），但也可以查看 PhyloXML 格式（en.wikipedia.org/wiki/PhyloXML）。

对齐基因和基因组数据

在进行任何系统发育分析之前，我们需要对基因和基因组数据进行对齐。在这里，我们将使用 MAFFT（mafft.cbrc.jp/alignment/software/）进行基因组分析。基因分析将使用 MUSCLE（www.drive5.com/muscle/）进行。

准备就绪

要执行基因组对齐，你需要安装 MAFFT。此外，为了进行基因对齐，将使用 MUSCLE。我们还将使用 trimAl（trimal.cgenomics.org/）以自动化方式去除虚假序列和对齐不良的区域。所有包都可以从 Bioconda 获取：

conda install –c bioconda mafft trimal muscle=3.8

如常，这些信息可以在相应的 Jupyter Notebook 文件Chapter07/Alignment.py中找到。你需要先运行之前的 Notebook，因为它会生成这里所需的文件。在本章中，我们将使用 Biopython。

如何操作...

查看以下步骤：

现在，我们将运行 MAFFT 来对齐基因组，如下面的代码所示。这个任务是 CPU 密集型和内存密集型的，且将花费相当长的时间：

from Bio.Align.Applications import MafftCommandline
mafft_cline = MafftCommandline(input='sample.fasta', ep=0.123, reorder=True, maxiterate=1000, localpair=True)
print(mafft_cline)
stdout, stderr = mafft_cline()
with open('align.fasta', 'w') as w:
    w.write(stdout)

之前的参数与论文附录中指定的参数相同。我们将使用 Biopython 接口调用 MAFFT。

让我们使用 trimAl 来修剪序列，如下所示：

os.system('trimal -automated1 -in align.fasta -out trim.fasta -fasta')

在这里，我们只是通过os.system调用应用程序。-automated1参数来自补充材料。

此外，我们可以运行MUSCLE来对蛋白质进行比对：

from Bio.Align.Applications import MuscleCommandline
my_genes = ['NP', 'L', 'VP35', 'VP40']
for gene in my_genes:
    muscle_cline = MuscleCommandline(input='%s_P.fasta' % gene)
    print(muscle_cline)
    stdout, stderr = muscle_cline()
    with open('%s_P_align.fasta' % gene, 'w') as w:
    w.write(stdout)

我们使用 Biopython 来调用外部应用程序。在这里，我们将对一组蛋白质进行比对。

请注意，为了进行分子进化分析，我们必须比较对齐后的基因，而不是蛋白质（例如，比较同义突变和非同义突变）。然而，我们只对齐了蛋白质。因此，我们必须将对齐数据转换为基因序列形式。

让我们通过找到三个对应于每个氨基酸的核苷酸来对基因进行比对：

from Bio import SeqIO
from Bio.Seq import Seq
from Bio.SeqRecord import SeqRecord
for gene in my_genes:
    gene_seqs = {}
    unal_gene = SeqIO.parse('%s.fasta' % gene, 'fasta')
    for rec in unal_gene:
        gene_seqs[rec.id] = rec.seq
    al_prot = SeqIO.parse('%s_P_align.fasta' % gene, 'fasta')
    al_genes = []
    for protein in al_prot:
        my_id = protein.id
        seq = ''
        pos = 0
        for c in protein.seq:
            if c == '-':
                seq += '---'
            else:
                seq += str(gene_seqs[my_id][pos:pos + 3])
                pos += 3
        al_genes.append(SeqRecord(Seq(seq), id=my_id))
    SeqIO.write(al_genes, '%s_align.fasta' % gene, 'fasta')

该代码获取蛋白质和基因编码。如果在蛋白质中发现空缺，则写入三个空缺；如果发现氨基酸，则写入相应的基因核苷酸。

比较序列

在这里，我们将比较在上一配方中比对的序列。我们将进行基因范围和基因组范围的比较。

准备开始

我们将使用 DendroPy，并且需要前两个配方的结果。像往常一样，这些信息可以在对应的笔记本Chapter07/Comparison.py中找到。

如何操作……

看一下接下来的步骤：

让我们开始分析基因数据。为简便起见，我们将仅使用来自伊波拉病毒属的另外两种物种的数据，这些数据已包含在扩展数据集中，即雷斯顿病毒（RESTV）和苏丹病毒（SUDV）：

import os
from collections import OrderedDict
import dendropy
from dendropy.calculate import popgenstat
genes_species = OrderedDict()
my_species = ['RESTV', 'SUDV']
my_genes = ['NP', 'L', 'VP35', 'VP40']
for name in my_genes:
    gene_name = name.split('.')[0]
    char_mat = dendropy.DnaCharacterMatrix.get_from_path('%s_align.fasta' % name, 'fasta')
    genes_species[gene_name] = {}

    for species in my_species:
        genes_species[gene_name][species] = dendropy.DnaCharacterMatrix()
    for taxon, char_map in char_mat.items():
        species = taxon.label.split('_')[0]
        if species in my_species:
            genes_species[gene_name][species].taxon_namespace.add_taxon(taxon)
            genes_species[gene_name][species][taxon] = char_map

我们得到在第一步中存储的四个基因，并在第二步中对其进行了比对。

我们加载所有文件（格式为 FASTA），并创建一个包含所有基因的字典。每个条目本身也是一个字典，包含 RESTV 或 SUDV 物种，及所有读取的数据。这些数据量不大，只有少量基因。

让我们打印所有四个基因的一些基本信息，如分离位点数（seg_sites）、核苷酸多样性（nuc_div）、Tajima’s D（taj_d）和 Waterson’s theta（wat_theta）（请查看本配方中*更多...*部分的链接，了解这些统计数据的相关信息）：

import numpy as np
import pandas as pd
summary = np.ndarray(shape=(len(genes_species), 4 * len(my_species)))
stats = ['seg_sites', 'nuc_div', 'taj_d', 'wat_theta']
for row, (gene, species_data) in enumerate(genes_species.items()):
    for col_base, species in enumerate(my_species):
        summary[row, col_base * 4] = popgenstat.num_segregating_sites(species_data[species])
        summary[row, col_base * 4 + 1] = popgenstat.nucleotide_diversity(species_data[species])
        summary[row, col_base * 4 + 2] = popgenstat.tajimas_d(species_data[species])
        summary[row, col_base * 4 + 3] = popgenstat.wattersons_theta(species_data[species])
columns = []
for species in my_species:
    columns.extend(['%s (%s)' % (stat, species) for stat in stats])
df = pd.DataFrame(summary, index=genes_species.keys(), columns=columns)
df # vs print(df)

首先，我们来看一下输出，然后再解释如何构建它：

图 7.1 – 病毒数据集的数据框

我使用pandas数据框打印结果，因为它非常适合处理这样的操作。我们将用一个 NumPy 多维数组初始化数据框，数组包含四行（基因）和四个统计数据乘以两个物种。

这些统计数据，如分离位点数、核苷酸多样性、Tajima’s D 和 Watterson’s theta，都是通过 DendroPy 计算的。请注意单个数据点在数组中的位置（坐标计算）。

看看最后一行：如果你在 Jupyter 中，只需要将 df 放在末尾，它会渲染出 DataFrame 和单元格输出。如果你不在笔记本中，使用 print(df)（你也可以在笔记本中执行这个操作，但显示效果可能不如直接在 Jupyter 中漂亮）。

现在，让我们提取类似的信息，但这次是基因组范围的数据，而不仅仅是基因范围。在这种情况下，我们将使用两个埃博拉爆发（2007 年和 2014 年）的子样本。我们将执行一个函数来显示基本统计信息，如下所示：

def do_basic_popgen(seqs):
    num_seg_sites = popgenstat.num_segregating_sites(seqs)
    avg_pair = popgenstat.average_number_of_pairwise_differences(seqs)
    nuc_div = popgenstat.nucleotide_diversity(seqs)
    print('Segregating sites: %d, Avg pairwise diffs: %.2f, Nucleotide diversity %.6f' % (num_seg_sites, avg_pair, nuc_div))
    print("Watterson's theta: %s" % popgenstat.wattersons_theta(seqs))
    print("Tajima's D: %s" % popgenstat.tajimas_d(seqs))

到现在为止，鉴于之前的例子，这个函数应该很容易理解。

现在，让我们正确地提取数据的子样本，并输出统计信息：

ebov_seqs = dendropy.DnaCharacterMatrix.get_from_path(
    'trim.fasta', schema='fasta', data_type='dna')
sl_2014 = []
drc_2007 = []
ebov2007_set = dendropy.DnaCharacterMatrix()
ebov2014_set = dendropy.DnaCharacterMatrix()
for taxon, char_map in ebov_seqs.items():
    print(taxon.label)
    if taxon.label.startswith('EBOV_2014') and len(sl_2014) < 8:
        sl_2014.append(char_map)
        ebov2014_set.taxon_namespace.add_taxon(taxon)
        ebov2014_set[taxon] = char_map
    elif taxon.label.startswith('EBOV_2007'):
        drc_2007.append(char_map)
        ebov2007_set.taxon_namespace.add_taxon(taxon)
        ebov2007_set[taxon] = char_map
        #ebov2007_set.extend_map({taxon: char_map})
del ebov_seqs
print('2007 outbreak:')
print('Number of individuals: %s' % len(ebov2007_set.taxon_set))
do_basic_popgen(ebov2007_set)
print('\n2014 outbreak:')
print('Number of individuals: %s' % len(ebov2014_set.taxon_set))
do_basic_popgen(ebov2014_set)

在这里，我们将构建两个数据集的两个版本：2014 年爆发和 2007 年爆发。我们将生成一个版本作为 DnaCharacterMatrix，另一个作为列表。我们将在这个食谱的最后使用这个列表版本。

由于 2014 年埃博拉爆发的数据集很大，我们仅用 8 个个体进行子样本抽取，这与 2007 年爆发的数据集的样本大小相当。

再次，我们删除 ebov_seqs 数据结构以节省内存（这些是基因组，而不仅仅是基因）。

如果你对 GenBank 上可用的 2014 年爆发的完整数据集（99 个样本）执行此分析，请准备好等待相当长的时间。

输出如下所示：

2007 outbreak:
Number of individuals: 7
Segregating sites: 25, Avg pairwise diffs: 7.71, Nucleotide diversity 0.000412
Watterson's theta: 10.204081632653063
Tajima's D: -1.383114157484101
2014 outbreak:
Number of individuals: 8
Segregating sites: 6, Avg pairwise diffs: 2.79, Nucleotide diversity 0.000149
Watterson's theta: 2.31404958677686
Tajima's D: 0.9501208027581887

最后，我们对 2007 年和 2014 年的两个子集进行一些统计分析，如下所示：

pair_stats = popgenstat.PopulationPairSummaryStatistics(sl_2014, drc_2007)
print('Average number of pairwise differences irrespective of population: %.2f' % pair_stats.average_number_of_pairwise_differences)
print('Average number of pairwise differences between populations: %.2f' % pair_stats.average_number_of_pairwise_differences_between)
print('Average number of pairwise differences within populations: %.2f' % pair_stats.average_number_of_pairwise_differences_within)
print('Average number of net pairwise differences : %.2f' % pair_stats.average_number_of_pairwise_differences_net)
print('Number of segregating sites: %d' % pair_stats.num_segregating_sites)
print("Watterson's theta: %.2f" % pair_stats.wattersons_theta)
print("Wakeley's Psi: %.3f" % pair_stats.wakeleys_psi)
print("Tajima's D: %.2f" % pair_stats.tajimas_d)

请注意，我们这里执行的操作稍有不同；我们将要求 DendroPy（popgenstat.PopulationPairSummaryStatistics）直接比较两个种群，以便获得以下结果：

Average number of pairwise differences irrespective of population: 284.46
Average number of pairwise differences between populations: 535.82
Average number of pairwise differences within populations: 10.50
Average number of net pairwise differences : 525.32
Number of segregating sites: 549
Watterson's theta: 168.84
Wakeley's Psi: 0.308
Tajima's D: 3.05

现在，分化位点的数量要大得多，因为我们正在处理来自两个合理分化的不同种群的数据。种群间的平均成对差异数值非常大。如预期的那样，这个数值远大于种群内部的平均数值，无论是否有种群信息。

还有更多……

如果你想获取更多的系统发育学和群体遗传学公式，包括这里使用的那些，我强烈推荐你获取 Arlequin 软件套件的手册（cmpg.unibe.ch/software/arlequin35/）。如果你不使用 Arlequin 进行数据分析，那么它的手册可能是实现公式的最佳参考。这本免费的文档可能包含比任何我能记得的书籍更相关的公式实现细节。

重建系统发育树

在这里，我们将为所有埃博拉物种的对齐数据集构建系统发育树。我们将遵循与论文中使用的过程非常相似的步骤。

准备工作

这个食谱需要 RAxML，这是一个用于最大似然推断大型系统发育树的程序，你可以在 sco.h-its.org/exelixis/software.xhtml 上查看它。Bioconda 也包含它，但它的名称是 raxml。请注意，二进制文件叫做 raxmlHPC。你可以执行以下命令来安装它：

conda install –c bioconda raxml

上面的代码很简单，但执行起来需要时间，因为它会调用 RAxML（这是一个计算密集型过程）。如果你选择使用 DendroPy 接口，它也可能会变得内存密集。我们将与 RAxML、DendroPy 和 Biopython 进行交互，给你选择使用哪个接口的自由；DendroPy 给你提供了一个简单的方式来访问结果，而 Biopython 则是内存占用较少的选择。尽管本章后面有一个可视化的食谱，但我们仍然会在这里绘制我们生成的树。

和往常一样，这些信息可以在相应的笔记本 Chapter07/Reconstruction.py 中找到。你需要前一个食谱的输出才能完成这个食谱。

如何操作…

看一下以下步骤：

对于 DendroPy，我们将首先加载数据，然后重建属数据集，如下所示：

import os
import shutil
import dendropy
from dendropy.interop import raxml
ebola_data = dendropy.DnaCharacterMatrix.get_from_path('trim.fasta', 'fasta')
rx = raxml.RaxmlRunner()
ebola_tree = rx.estimate_tree(ebola_data, ['-m', 'GTRGAMMA', '-N', '10'])
print('RAxML temporary directory %s:' % rx.working_dir_path)
del ebola_data

请记住，这个数据结构的大小相当大；因此，确保你有足够的内存来加载它（至少 10 GB）。

要准备好等待一段时间。根据你的计算机，这可能需要超过一个小时。如果时间过长，考虑重新启动进程，因为有时 RAxML 可能会出现 bug。

我们将使用 GTRΓ 核苷酸替代模型运行 RAxML，如论文中所述。我们只进行 10 次重复以加快结果速度，但你可能应该做更多的重复，比如 100 次。在过程结束时，我们将从内存中删除基因组数据，因为它占用了大量内存。

ebola_data 变量将包含最佳的 RAxML 树，并且包括距离信息。RaxmlRunner 对象将可以访问 RAxML 生成的其他信息。让我们打印出 DendroPy 将执行 RAxML 的目录。如果你检查这个目录，你会发现很多文件。由于 RAxML 返回的是最佳树，你可能会忽略所有这些文件，但我们将在替代的 Biopython 步骤中稍作讨论。

我们将保存树用于未来的分析；在我们的案例中，它将是一个可视化，如下代码所示：
```
ebola_tree.write_to_path('my_ebola.nex', 'nexus')
```

我们将把序列写入 NEXUS 文件，因为我们需要存储拓扑信息。FASTA 在这里不够用。

让我们可视化我们的属树，如下所示：

import matplotlib.pyplot as plt
from Bio import Phylo
my_ebola_tree = Phylo.read('my_ebola.nex', 'nexus')
my_ebola_tree.name = 'Our Ebolavirus tree'
fig = plt.figure(figsize=(16, 18))
ax = fig.add_subplot(1, 1, 1)
Phylo.draw(my_ebola_tree, axes=ax)

我们将在稍后介绍合适的步骤时再解释这段代码，但如果你查看下图并将其与论文中的结果进行比较，你会清楚地看到它像是朝着正确方向迈进的一步。例如，所有同一物种的个体被聚集在一起。

你会注意到，trimAl 改变了其序列的名称，例如，通过添加它们的大小。这很容易解决；我们将在可视化系统发育数据这一部分中处理这个问题：

图 7.2 – 我们使用 RAxML 生成的所有埃博拉病毒的系统发育树

让我们通过 Biopython 使用 RAxML 重建系统发育树。Biopython 的接口比 DendroPy 更少声明式，但在内存效率上要高得多。因此，在运行完后，你需要负责处理输出，而 DendroPy 会自动返回最优的树，如下代码所示：

import random
import shutil
from Bio.Phylo.Applications import RaxmlCommandline
raxml_cline = RaxmlCommandline(sequences='trim.fasta', model='GTRGAMMA', name='biopython', num_replicates='10', parsimony_seed=random.randint(0, sys.maxsize), working_dir=os.getcwd() + os.sep + 'bp_rx')
print(raxml_cline)
try:
    os.mkdir('bp_rx')
except OSError:
    shutil.rmtree('bp_rx')
    os.mkdir('bp_rx')
out, err = raxml_cline()

DendroPy 比 Biopython 具有更具声明性的接口，因此你可以处理一些额外的事项。你应该指定种子（如果不指定，Biopython 会自动使用 10,000 作为默认值）以及工作目录。使用 RAxML 时，工作目录的指定要求使用绝对路径。

让我们检查一下 Biopython 运行的结果。虽然 RAxML 的输出（除了随机性）对于 DendroPy 和 Biopython 是相同的，但 DendroPy 隐藏了几个细节。使用 Biopython，你需要自己处理结果。你也可以使用 DendroPy 来执行这个操作；不过在这种情况下，它是可选的：
```
from Bio import Phylo
biopython_tree = Phylo.read('bp_rx/RAxML_bestTree.biopython', 'newick')
```

上面的代码将读取 RAxML 运行中最优的树。文件名后附加了你在前一步中指定的项目名称（在本例中为biopython）。

看看bp_rx目录的内容；在这里，你将找到来自 RAxML 的所有输出，包括所有 10 个备选树。

还有更多……

尽管本书的目的是不教授系统发育分析，但了解为什么我们不检查树拓扑中的共识和支持信息还是很重要的。你应该在自己的数据集里研究这一点。更多信息，请参考www.geol.umd.edu/~tholtz/G331/lectures/cladistics5.pdf。

递归操作树

这不是一本关于 Python 编程的书，因为这个话题非常广泛。话虽如此，入门级 Python 书籍通常不会详细讨论递归编程。通常，递归编程技术非常适合处理树结构。而且，它还是功能性编程方言中的一种必要编程策略，这在进行并发处理时非常有用。这在处理非常大的数据集时是常见的。

系统发育树的概念与计算机科学中的树有所不同。系统发育树可以是有根的（如果是，它们就是普通的树数据结构），也可以是无根的，使其成为无向非循环图。此外，系统发育树的边上可以有权重。因此，在阅读文档时要注意这一点；如果文献是由系统发育学家编写的，你可以期待该树是有根或无根的，而大多数其他文档则会使用无向非循环图来表示无根树。在这个配方中，我们假设所有树都是有根的。

最后，请注意，虽然本配方主要旨在帮助你理解递归算法和树形结构，但最后一部分实际上非常实用，并且对下一个配方的实现至关重要。

准备工作

你需要准备好前一配方中的文件。像往常一样，你可以在 Chapter07/Trees.py 笔记本文件中找到这些内容。在这里，我们将使用 DendroPy 的树表示。请注意，比较其他树表示和库（无论是系统发育的还是非系统发育的）而言，大部分代码都是易于通用的。

如何操作...

看一下以下步骤：

首先，让我们加载由 RAxML 生成的所有埃博拉病毒树，如下所示：

import dendropy
ebola_raxml = dendropy.Tree.get_from_path('my_ebola.nex', 'nexus')

接着，我们需要计算每个节点的层级（到根节点的距离）：

def compute_level(node, level=0):
    for child in node.child_nodes():
        compute_level(child, level + 1)
    if node.taxon is not None:
        print("%s: %d %d" % (node.taxon, node.level(), level))
compute_level(ebola_raxml.seed_node)

DendroPy 的节点表示有一个层级方法（用于比较），但这里的重点是介绍递归算法，所以我们还是会实现它。

注意这个函数的工作原理；它是用 seed_node 调用的（这是根节点，因为代码假设我们处理的是有根树）。根节点的默认层级是 0。然后，该函数会递归调用所有子节点，将层级加一。对于每个非叶子节点（即，它是树的内部节点），这个调用会被重复，直到我们到达叶子节点。

对于叶子节点，我们接着打印出层级（我们也可以对内部节点执行相同的操作），并展示由 DendroPy 内部函数计算出的相同信息。

现在，让我们计算每个节点的高度。节点的高度是从该节点开始的最大向下路径（通向叶子）的边数，如下所示：

def compute_height(node):
    children = node.child_nodes()
    if len(children) == 0:
        height = 0
    else:
    height = 1 + max(map(lambda x: compute_height(x), children))
    desc = node.taxon or 'Internal'
    print("%s: %d %d" % (desc, height, node.level()))
    return height
compute_height(ebola_raxml.seed_node)

在这里，我们将使用相同的递归策略，但每个节点将把它的高度返回给其父节点。如果该节点是叶子节点，则高度为 0；如果不是，则高度为 1 加上其所有后代的最大高度。

请注意，我们使用 map 配合 lambda 函数来获取当前节点所有子节点的高度。然后，我们选择最大值（max 函数在这里执行一个 reduce 操作，因为它总结了所有报告的值）。如果你将其与 MapReduce 框架联系起来，你是正确的；它们的灵感来自于像这样的函数式编程方言。

现在，让我们计算每个节点的子孙数量。到现在为止，这应该很容易理解：

def compute_nofs(node):
    children = node.child_nodes()
    nofs = len(children)
    map(lambda x: compute_nofs(x), children)
    desc = node.taxon or 'Internal'
    print("%s: %d %d" % (desc, nofs, node.level()))
compute_nofs(ebola_raxml.seed_node)

现在，我们将打印所有叶节点（显然，这是微不足道的）：

def print_nodes(node):
    for child in node.child_nodes():
        print_nodes(child)
    if node.taxon is not None:
        print('%s (%d)' % (node.taxon, node.level()))
print_nodes(ebola_raxml.seed_node)

请注意，到目前为止，我们开发的所有函数都对树施加了非常明确的遍历模式。它首先调用它的第一个子节点，然后该子节点会调用它的子节点，依此类推；只有这样，函数才能按深度优先的模式调用下一个子节点。然而，我们也可以采取不同的方式。

现在，让我们以广度优先的方式打印叶节点，也就是说，我们会首先打印离根节点最近的叶节点，按如下方式：

from collections import deque
def print_breadth(tree):
    queue = deque()
    queue.append(tree.seed_node)
    while len(queue) > 0:
        process_node = queue.popleft()
        if process_node.taxon is not None:
            print('%s (%d)' % (process_node.taxon, process_node.level()))
        else:
            for child in process_node.child_nodes():
                queue.append(child)
print_breadth(ebola_raxml)

在我们解释这个算法之前，先看看这个运行的结果与上一次运行的结果有何不同。首先，看看下面的图示。如果你按深度优先顺序打印节点，你会得到 Y、A、X、B 和 C。但是如果你执行广度优先遍历，你会得到 X、B、C、Y 和 A。树的遍历会影响节点的访问顺序；这一点往往非常重要。

关于前面的代码，在这里我们将采用完全不同的方法，因为我们将执行一个迭代算法。我们将使用先进先出（FIFO）队列来帮助我们排序节点。请注意，Python 的 deque 可以像 FIFO 一样高效地使用，也可以用于后进先出（LIFO）。这是因为它在两端操作时实现了高效的数据结构。

算法从将根节点放入队列开始。当队列不为空时，我们将取出队列中的第一个节点。如果是内部节点，我们将把它的所有子节点放入队列。

我们将继续执行前一步骤，直到队列为空。我鼓励你拿起笔和纸，通过执行下图所示的示例来看看这个过程是如何运作的。代码虽小，但并不简单：

图 7.3 – 遍历树；第一个数字表示在深度优先遍历中访问该节点的顺序，第二个数字则表示广度优先遍历的顺序

让我们回到实际的数据集。由于我们有太多数据无法完全可视化，我们将生成一个精简版，去除包含单一物种的子树（以 EBOV 为例，它们有相同的爆发）。我们还将进行树的阶梯化，即按子节点的数量对子节点进行排序：

from copy import deepcopy
simple_ebola = deepcopy(ebola_raxml)
def simplify_tree(node):
    prefs = set()
    for leaf in node.leaf_nodes():
        my_toks = leaf.taxon.label.split(' ')
        if my_toks[0] == 'EBOV':
            prefs.add('EBOV' + my_toks[1])
        else:
            prefs.add(my_toks[0])
    if len(prefs) == 1:
        print(prefs, len(node.leaf_nodes()))
        node.taxon = dendropy.Taxon(label=list(prefs)[0])
        node.set_child_nodes([])
    else:
        for child in node.child_nodes():
            simplify_tree(child)
simplify_tree(simple_ebola.seed_node)
simple_ebola.ladderize()
simple_ebola.write_to_path('ebola_simple.nex', 'nexus')

我们将对树结构进行深拷贝。由于我们的函数和阶梯化过程是破坏性的（它们会改变树结构），我们需要保持原始树结构不变。

DendroPy 能够列举出所有叶节点（在这个阶段，一个好的练习是编写一个函数来执行这个操作）。通过这个功能，我们可以获取某个节点的所有叶子。如果它们与 EBOV 的情况一样，具有相同的物种和爆发年份，我们将移除所有子节点、叶节点和内部子树节点。

如果它们不属于相同物种，我们将递归下去，直到满足条件。最坏的情况是当你已经处于叶节点时，算法会简单地解析为当前节点的物种。

还有更多...

关于树和数据结构的话题，有大量的计算机科学文献；如果你想阅读更多内容，维基百科提供了一个很好的介绍，网址是en.wikipedia.org/wiki/Tree_%28data_structure%29。

请注意，lambda函数和map的使用并不被推荐作为 Python 方言；你可以阅读 Guido van Rossum 关于这个主题的一些（较旧的）观点，网址是www.artima.com/weblogs/viewpost.jsp?thread=98196。我在这里呈现它是因为它是函数式和递归编程中的一种非常常见的方言。更常见的方言将基于列表推导式。

无论如何，基于使用map和reduce操作的函数式方言是 MapReduce 框架的概念基础，你可以使用像 Hadoop、Disco 或 Spark 这样的框架来执行高性能的生物信息学计算。

可视化系统发育数据

在这个配方中，我们将讨论如何可视化系统发育树。DendroPy 仅具有基于绘制文本 ASCII 树的简单可视化机制，但 Biopython 拥有非常丰富的基础设施，我们将在这里利用它。

准备工作

这将要求你完成所有前面的配方。记住，我们拥有完整的埃博拉病毒属的文件，包括 RAxML 树。此外，简化版属版本将在前面的配方中生成。如常，你可以在Chapter07/Visualization.py笔记本文件中找到这些内容。

如何实现...

请看下面的步骤：

让我们加载所有系统发育数据：

from copy import deepcopy
from Bio import Phylo
ebola_tree = Phylo.read('my_ebola.nex', 'nexus')
ebola_tree.name = 'Ebolavirus tree'
ebola_simple_tree = Phylo.read('ebola_simple.nex', 'nexus')
ebola_simple_tree.name = 'Ebolavirus simplified tree'

对于我们读取的所有树，我们将更改树的名称，因为稍后会打印出该名称。

现在，我们可以绘制树的 ASCII 表示：

Phylo.draw_ascii(ebola_simple_tree)
Phylo.draw_ascii(ebola_tree)

简化版属树的 ASCII 表示如下所示。在这里，我们不会打印完整版本，因为它将占用好几页。但如果你运行前面的代码，你将能够看到它实际上是非常易于阅读的：

图 7.4 – 简化版埃博拉病毒数据集的 ASCII 表示

Bio.Phylo通过使用matplotlib作为后端来实现树的图形表示：

import matplotlib.pyplot as plt
fig = plt.figure(figsize=(16, 22))
ax = fig.add_subplot(111)
Phylo.draw(ebola_simple_tree, branch_labels=lambda c: c.branch_length if c.branch_length > 0.02 else None, axes=ax)

在这种情况下，我们将在边缘处打印分支长度，但会去除所有小于 0.02 的长度，以避免杂乱。这样做的结果如下图所示：

图 7.5 – 一个基于 matplotlib 的简化数据集版本，并添加了分支长度

现在我们将绘制完整的数据集，但每个树的部分将使用不同的颜色。如果一个子树只有单一的病毒种类，它将拥有自己独特的颜色。埃博拉病毒（EBOV）将有两种颜色，也就是说，一种用于 2014 年的疫情，另一种用于其他年份，如下所示：

fig = plt.figure(figsize=(16, 22))
ax = fig.add_subplot(111)
from collections import OrderedDict
my_colors = OrderedDict({
'EBOV_2014': 'red',
'EBOV': 'magenta',
'BDBV': 'cyan',
'SUDV': 'blue',
'RESTV' : 'green',
'TAFV' : 'yellow'
})
def get_color(name):
    for pref, color in my_colors.items():
        if name.find(pref) > -1:
            return color
    return 'grey'
def color_tree(node, fun_color=get_color):
    if node.is_terminal():
        node.color = fun_color(node.name)
    else:
        my_children = set()
        for child in node.clades:
            color_tree(child, fun_color)
            my_children.add(child.color.to_hex())
        if len(my_children) == 1:
            node.color = child.color
        else:
            node.color = 'grey'
ebola_color_tree = deepcopy(ebola_tree)
color_tree(ebola_color_tree.root)
Phylo.draw(ebola_color_tree, axes=ax, label_func=lambda x: x.name.split(' ')[0][1:] if x.name is not None else None)

这是一个树遍历算法，类似于前面例子中呈现的算法。作为一个递归算法，它的工作方式如下。如果节点是叶子，它将根据其种类（或 EBOV 疫情年份）来着色。如果它是一个内部节点，并且所有下面的后代节点都是同一物种，它将采用该物种的颜色；如果后代节点包含多个物种，它将被着色为灰色。实际上，颜色函数可以更改，并且稍后会进行更改。只有边缘颜色会被使用（标签会以黑色打印）。

注意，阶梯化（在前面的例子中使用 DendroPy 完成）对于清晰的视觉效果帮助很大。

我们还会对属树进行深拷贝，以便对副本进行着色；请记住，在前面的例子中提到的某些树遍历函数可能会改变状态，而在这种情况下，我们希望保留一个没有任何着色的版本。

注意使用了 lambda 函数来清理由 trimAl 修改的名称，如下图所示：

图 7.6 – 一个带有完整埃博拉病毒数据集的阶梯化和着色的系统发育树

还有更多...

树和图的可视化是一个复杂的话题；可以说，这里的树可视化是严谨的，但远不漂亮。DendroPy 的一个替代方案是 ETE（etetoolkit.org/），它具有更多的可视化功能。绘制树和图的常见替代方案包括 Cytoscape（cytoscape.org/）和 Gephi（gephi.github.io/）。如果你想了解更多关于渲染树和图的算法，可以查看 Wikipedia 页面：en.wikipedia.org/wiki/Graph_drawing，了解这个迷人的话题。

不过，要小心不要以牺牲实质内容为代价追求风格。例如，本书的上一版使用图形渲染库绘制了一个美观的系统发育树。虽然它显然是该章节中最漂亮的图像，但在分支长度上却具有误导性。

第九章：使用蛋白质数据银行

蛋白质组学是研究蛋白质的学科，包括蛋白质的功能和结构。该领域的主要目标之一是表征蛋白质的三维结构。在蛋白质组学领域，最广为人知的计算资源之一是蛋白质数据银行（PDB），这是一个包含大分子生物体结构数据的数据库。当然，也有许多数据库专注于蛋白质的初级结构；这些数据库与我们在第二章中看到的基因组数据库有些相似，了解 NumPy、pandas、Arrow 和 Matplotlib。

在本章中，我们将主要关注如何处理来自 PDB 的数据。我们将学习如何解析 PDB 文件，执行一些几何计算，并可视化分子。我们将使用旧的 PDB 文件格式，因为从概念上讲，它允许你在一个稳定的环境中执行大多数必要的操作。话虽如此，新的 mmCIF 格式计划取代 PDB 格式，在使用 Biopython 解析 mmCIF 文件的食谱中也会介绍它。我们将使用 Biopython 并介绍 PyMOL 用于可视化。我们不会在这里讨论分子对接，因为那更适合一本关于化学信息学的书。

在本章中，我们将使用一个经典的蛋白质例子：肿瘤蛋白 p53，它参与细胞周期的调节（例如，凋亡）。该蛋白质与癌症关系密切。网上有大量关于该蛋白质的信息。

让我们从你现在应该更熟悉的内容开始：访问数据库，特别是获取蛋白质的初级结构（即氨基酸序列）。

在本章中，我们将介绍以下内容：

在多个数据库中查找蛋白质
介绍 Bio.PDB
从 PDB 文件中提取更多信息
在 PDB 文件中计算分子距离
执行几何操作
使用 PyMOL 进行动画制作
使用 Biopython 解析 mmCIF 文件

在多个数据库中查找蛋白质

在我们开始进行更多的结构生物学工作之前，我们将看看如何访问现有的蛋白质组学数据库，比如 UniProt。我们将查询 UniProt 以查找我们感兴趣的基因，TP53，并从那里开始。

做好准备

为了访问数据，我们将使用 Biopython 和 REST API（我们在第五章中使用了类似的方法，基因组学工作）以及requests库来访问 Web API。requests API 是一个易于使用的 Web 请求封装库，可以通过标准的 Python 机制（例如，pip和conda）安装。你可以在Chapter08/Intro.py笔记本文件中找到这部分内容。

如何实现...

请查看以下步骤：

首先，我们来定义一个函数来执行对 UniProt 的 REST 查询，代码如下：

import requests
server = 'http://www.uniprot.org/uniprot'
def do_request(server, ID='', **kwargs):
    params = ''
    req = requests.get('%s/%s%s' % (server, ID, params), params=kwargs)
    if not req.ok:
        req.raise_for_status()
    return req

现在，我们可以查询所有已审阅的p53基因：

req = do_request(server, query='gene:p53 AND reviewed:yes', format='tab',
 columns='id,entry name,length,organism,organism-id,database(PDB),database(HGNC)',
 limit='50')

我们将查询p53基因，并请求查看所有已审核的条目（即，手动校对过的）。输出将以表格格式显示。我们将请求最多 50 个结果，并指定所需的列。

我们本可以将输出限制为仅包含人类数据，但为了这个示例，我们将包括所有可用的物种。

让我们查看结果，内容如下：

import pandas as pd
import io
uniprot_list = pd.read_table(io.StringIO(req.text))
uniprot_list.rename(columns={'Organism ID': 'ID'}, inplace=True)
print(uniprot_list)

我们使用pandas来轻松处理制表符分隔的列表并进行美观打印。笔记本的简化输出如下：

图 8.1 - 一个简化的 TP53 蛋白物种列表

现在，我们可以获取人类p53基因 ID，并使用 Biopython 检索并解析SwissProt记录：

from Bio import ExPASy, SwissProt
p53_human = uniprot_list[
    (uniprot_list.ID == 9606) &
    (uniprot_list['Entry name'].str.contains('P53'))]['Entry'].iloc[0] 
handle = ExPASy.get_sprot_raw(p53_human)
sp_rec = SwissProt.read(handle)

然后，我们使用 Biopython 的SwissProt模块来解析记录。9606是人类的 NCBI 分类代码。

和往常一样，如果网络服务出现错误，可能是网络或服务器问题。如果是这样，请稍后再试。

让我们来看一下p53记录，内容如下：

print(sp_rec.entry_name, sp_rec.sequence_length, sp_rec.gene_name)
print(sp_rec.description)
print(sp_rec.organism, sp_rec.seqinfo)
print(sp_rec.sequence)
print(sp_rec.comments)
print(sp_rec.keywords)

输出如下：

P53_HUMAN 393 Name=TP53; Synonyms=P53;
 RecName: Full=Cellular tumor antigen p53; AltName: Full=Antigen NY-CO-13; AltName: Full=Phosphoprotein p53; AltName: Full=Tumor suppressor p53;
 Homo sapiens (Human). (393, 43653, 'AD5C149FD8106131')
 MEEPQSDPSVEPPLSQETFSDLWKLLPENNVLSPLPSQAMDDLMLSPDDIEQWFTED PGPDEAPRMPEAAPPVAPAPAAPTPAAPAPAPSWPLSSSVPSQKTYQGSYGFRLGF LHSGTAKSVTCTYSPALNKMFCQLAKTCPVQLWVDSTPPPGTRVRAMAIYKQSQHM TEVVRRCPHHERCSDSDGLAPPQHLIRVEGNLRVEYLDDRNTFRHSVVVPYEPPEVG SDCTTIHYNYMCNSSCMGGMNRRPILTIITLEDSSGNLLGRNSFEVRVCACPGRDRR TEEENLRKKGEPHHELPPGSTKRALPNNTSSSPQPKKKPLDGEYFTLQIRGRERFEM FRELNEALELKDAQAGKEPGGSRAHSSHLKSKKGQSTSRHKKLMFKTEGPDSD

更深入地查看前面的记录会发现许多非常有趣的信息，特别是在特征、cross_references方面：

from collections import defaultdict
done_features = set()
print(len(sp_rec.features))
for feature in sp_rec.features:
    if feature[0] in done_features:
        continue
    else:
        done_features.add(feature[0])
        print(feature)
print(len(sp_rec.cross_references))
per_source = defaultdict(list)
for xref in sp_rec.cross_references:
    source = xref[0]
    per_source[source].append(xref[1:])
print(per_source.keys())
done_GOs = set()
print(len(per_source['GO']))
for annot in per_source['GO']:
    if annot[1][0] in done_GOs:
        continue
    else:
        done_GOs.add(annot[1][0])
        print(annot)

请注意，我们这里并没有打印所有的信息，仅仅是一个摘要。我们打印了序列的多个特征，每种类型展示一个例子，还有一些外部数据库的引用，以及提到的数据库，还有一些 GO 条目，并附上了三个例子。目前，仅这个蛋白质就有 1,509 个特征、923 个外部引用和 173 个 GO 术语。以下是输出的高度简化版本：

Total features: 1509
type: CHAIN
location: [0:393]
id: PRO_0000185703
qualifiers:
    Key: note, Value: Cellular tumor antigen p53
type: DNA_BIND
location: [101:292]
qualifiers:
type: REGION
location: [0:320]
qualifiers:
    Key: evidence, Value: ECO:0000269|PubMed:25732823
    Key: note, Value: Interaction with CCAR2
[...]
Cross references:  923
dict_keys(['EMBL', 'CCDS', 'PIR', 'RefSeq', 'PDB', 'PDBsum', 'BMRB', 'SMR', 'BioGRID', 'ComplexPortal', 'CORUM', 'DIP', 'ELM', 'IntAct', 'MINT', 'STRING', 'BindingDB', 'ChEMBL', 'DrugBank', 'MoonDB', 'TCDB', 'GlyGen', 'iPTMnet', 'MetOSite', 'PhosphoSitePlus', 'BioMuta', 'DMDM', 'SWISS-2DPAGE', 'CPTAC', 'EPD', 'jPOST', 'MassIVE', 'MaxQB', 'PaxDb', 'PeptideAtlas', 'PRIDE', 'ProteomicsDB', 'ABCD', 'Antibodypedia', 'CPTC', 'DNASU', 'Ensembl', 'GeneID', 'KEGG', 'MANE-Select', 'UCSC', 'CTD', 'DisGeNET', 'GeneCards', 'GeneReviews', 'HGNC', 'HPA', 'MalaCards', 'MIM', 'neXtProt', 'OpenTargets', 'Orphanet', 'PharmGKB', 'VEuPathDB', 'eggNOG', 'GeneTree', 'InParanoid', 'OMA', 'OrthoDB', 'PhylomeDB', 'TreeFam', 'PathwayCommons', 'Reactome', 'SABIO-RK', 'SignaLink', 'SIGNOR', 'BioGRID-ORCS', 'ChiTaRS', 'EvolutionaryTrace', 'GeneWiki', 'GenomeRNAi', 'Pharos', 'PRO', 'Proteomes', 'RNAct', 'Bgee', 'ExpressionAtlas', 'Genevisible', 'GO', 'CDD', 'DisProt', 'Gene3D', 'IDEAL', 'InterPro', 'PANTHER', 'Pfam', 'PRINTS', 'SUPFAM', 'PROSITE'])
Annotation SOURCES: 173
('GO:0005813', 'C:centrosome', 'IDA:UniProtKB')
('GO:0036310', 'F:ATP-dependent DNA/DNA annealing activity', 'IDA:UniProtKB')
('GO:0006914', 'P:autophagy', 'IMP:CAFA')

介绍 Bio.PDB

在这里，我们将介绍 Biopython 的PDB模块，用于处理 PDB 文件。我们将使用三个模型，这些模型代表p53蛋白的部分结构。你可以在www.rcsb.org/pdb/101/motm.do?momID=31了解更多关于这些文件和p53的信息。

准备工作

你应该已经了解了基本的PDB数据模型，包括模型、链、残基和原子对象。关于Biopython 的结构生物信息学 FAQ的详细解释可以在biopython.org/wiki/The_Biopython_Structural_Bioinformatics_FAQ找到。

你可以在Chapter08/PDB.py笔记本文件中找到这些内容。

在我们将要下载的三个模型中，1TUP模型将用于接下来的所有示例。花点时间研究这个模型，它将在后续帮助你。

如何操作……

请看以下步骤：

首先，让我们检索我们感兴趣的模型，如下所示：

from Bio import PDB
repository = PDB.PDBList()
repository.retrieve_pdb_file('1TUP', pdir='.', file_format='pdb')
repository.retrieve_pdb_file('1OLG', pdir='.', file_format='pdb')
repository.retrieve_pdb_file('1YCQ', pdir='.', file_format='pdb')

请注意，Bio.PDB会为您下载文件。此外，只有在没有本地副本的情况下才会进行这些下载。

让我们解析我们的记录，如下所示的代码：

parser = PDB.PDBParser()
p53_1tup = parser.get_structure('P 53 - DNA Binding', 'pdb1tup.ent')
p53_1olg = parser.get_structure('P 53 - Tetramerization', 'pdb1olg.ent')
p53_1ycq = parser.get_structure('P 53 - Transactivation', 'pdb1ycq.ent')

您可能会收到有关文件内容的一些警告。这些通常不会有问题。

让我们检查我们的头文件，如下所示：

def print_pdb_headers(headers, indent=0):
   ind_text = ' ' * indent
   for header, content in headers.items():
       if type(content) == dict:
          print('\n%s%20s:' % (ind_text, header))
          print_pdb_headers(content, indent + 4)
          print()
       elif type(content) == list:
          print('%s%20s:' % (ind_text, header))
          for elem in content:
              print('%s%21s %s' % (ind_text, '->', elem))
      else:
          print('%s%20s: %s' % (ind_text, header, content))
print_pdb_headers(p53_1tup.header)

头文件被解析为字典的字典。因此，我们将使用递归函数来解析它们。此函数将增加缩进以便于阅读，并使用->前缀注释元素列表。有关递归函数的示例，请参见前一章，第七章，系统发生学。有关 Python 中递归的高级讨论，请转到最后一章，第十二章，生物信息学的函数式编程。简化后的输出如下所示：

                name: tumor suppressor p53 complexed with dna
                head: antitumor protein/dna
              idcode: 1TUP
     deposition_date: 1995-07-11
        release_date: 1995-07-11
    structure_method: x-ray diffraction
          resolution: 2.2
 structure_reference:
-> n.p.pavletich,k.a.chambers,c.o.pabo the dna-binding domain of p53 contains the four conserved regions and the major mutation hot spots genes dev. v. 7 2556 1993 issn 0890-9369 
              author: Y.Cho,S.Gorina,P.D.Jeffrey,N.P.Pavletich
            compound:
                       2:
misc: 
molecule: dna (5'-d(*ap*tp*ap*ap*tp*tp*gp*gp*gp*cp*ap*ap*gp*tp*cp*tp*a p*gp*gp*ap*a)-3') 
                       chain: f
                  engineered: yes
has_missing_residues: True
    missing_residues:
                   -> {'model': None, 'res_name': 'ARG', 'chain': 'A', 'ssseq': 290, 'insertion': None}
keywords: antigen p53, antitumor protein/dna complex
             journal: AUTH   Y.CHO,S.GORINA,P.D.JEFFREY,N.P.PAVLETICHTITL   CRYSTAL STRUCTURE OF A P53 TUMOR SUPPRESSOR-DNATITL 2 COMPLEX: UNDERSTANDING TUMORIGENIC MUTATIONS.REF    SCIENCE57

我们想要知道这些文件中每条链的内容；为此，让我们看一下COMPND记录：

print(p53_1tup.header['compound'])
print(p53_1olg.header['compound'])
print(p53_1ycq.header['compound'])

这将打印出在前面的代码中列出的所有化合物头文件。不幸的是，这不是获取链信息的最佳方式。另一种方法是获取DBREF记录，但 Biopython 的解析器目前无法访问这些记录。话虽如此，使用诸如grep之类的工具将轻松提取这些信息。

注意，对于1TUP模型，链A、B和C来自蛋白质，而链E和F来自 DNA。这些信息将在未来很有用。

让我们对每个PDB文件进行自上而下的分析。现在，让我们只获取所有的链、残基数和每条链中的原子数，如下所示：

def describe_model(name, pdb):
print()
for model in pdb:
    for chain in model:
        print('%s - Chain: %s. Number of residues: %d. Number of atoms: %d.' %
              (name, chain.id, len(chain), len(list(chain.get_atoms()))))
describe_model('1TUP', p53_1tup)
describe_model('1OLG', p53_1olg)
describe_model('1YCQ', p53_1ycq)

在稍后的配方中，我们将采用自下而上的方法。以下是1TUP的输出：

1TUP - Chain: E. Number of residues: 43\. Number of atoms: 442.
1TUP - Chain: F. Number of residues: 35\. Number of atoms: 449.
1TUP - Chain: A. Number of residues: 395\. Number of atoms: 1734.
1TUP - Chain: B. Number of residues: 265\. Number of atoms: 1593.
1TUP - Chain: C. Number of residues: 276\. Number of atoms: 1610.

1OLG - Chain: A. Number of residues: 42\. Number of atoms: 698.
1OLG - Chain: B. Number of residues: 42\. Number of atoms: 698.
1OLG - Chain: C. Number of residues: 42\. Number of atoms: 698.
1OLG - Chain: D. Number of residues: 42\. Number of atoms: 698.

1YCQ - Chain: A. Number of residues: 123\. Number of atoms: 741.
1YCQ - Chain: B. Number of residues: 16\. Number of atoms: 100.

让我们获取所有非标准残基（HETATM），除了水，在1TUP模型中，如下所示的代码：

for residue in p53_1tup.get_residues():
    if residue.id[0] in [' ', 'W']:
        continue
print(residue.id)

我们有三个锌原子，每个蛋白链一个。

让我们来看一个残基：

res = next(p53_1tup[0]['A'].get_residues())
print(res)
for atom in res:
    print(atom, atom.serial_number, atom.element)
p53_1tup[0]['A'][94]['CA']

这将打印出某个残基中的所有原子：

<Residue SER het=  resseq=94 icode= >
 <Atom N> 858 N
 <Atom CA> 859 C
 <Atom C> 860 C
 <Atom O> 861 O
 <Atom CB> 862 C
 <Atom OG> 863 O
 <Atom CA>

注意最后一句话。它只是为了向您展示，您可以通过解析模型、链、残基和最终原子来直接访问一个原子。

最后，让我们将蛋白质片段导出到一个 FASTA 文件中，如下所示：

from Bio.SeqIO import PdbIO, FastaIO
def get_fasta(pdb_file, fasta_file, transfer_ids=None):
    fasta_writer = FastaIO.FastaWriter(fasta_file)
    fasta_writer.write_header()
    for rec in PdbIO.PdbSeqresIterator(pdb_file):
        if len(rec.seq) == 0:
            continue
        if transfer_ids is not None and rec.id not in transfer_ids:
            continue
        print(rec.id, rec.seq, len(rec.seq))
        fasta_writer.write_record(rec)

get_fasta(open('pdb1tup.ent'), open('1tup.fasta', 'w'), transfer_ids=['1TUP:B'])
get_fasta(open('pdb1olg.ent'), open('1olg.fasta', 'w'), transfer_ids=['1OLG:B'])
get_fasta(open('pdb1ycq.ent'), open('1ycq.fasta', 'w'), transfer_ids=['1YCQ:B'])

如果您检查蛋白质链，您会发现它们在每个模型中都是相等的，因此我们导出一个单独的链。在1YCQ的情况下，我们导出最小的一个，因为最大的一个与p53无关。正如您在这里看到的，我们使用的是Bio.SeqIO，而不是Bio.PDB。

还有更多

PDB 解析器不完整。很可能不会很快见到完整的解析器，因为社区正在迁移到 mmCIF 格式。

尽管未来是 mmCIF 格式（mmcif.wwpdb.org/），PDB 文件仍然存在。从概念上讲，文件解析后的许多操作是类似的。

从 PDB 文件中提取更多信息

在这里，我们将继续探索Bio.PDB从 PDB 文件生成的记录结构。

准备工作

有关我们正在使用的 PDB 模型的详细信息，请参见前面的章节。

你可以在Chapter08/Stats.py Notebook 文件中找到这些内容。

如何操作…

我们将通过以下步骤开始：

首先，让我们提取1TUP，如下所示：

from Bio import PDB
repository = PDB.PDBList()
parser = PDB.PDBParser()
repository.retrieve_pdb_file('1TUP', pdir='.', file_format='pdb') p53_1tup = parser.get_structure('P 53', 'pdb1tup.ent')

然后，提取一些与原子相关的统计信息：

from collections import defaultdict
atom_cnt = defaultdict(int)
atom_chain = defaultdict(int)
atom_res_types = defaultdict(int)
for atom in p53_1tup.get_atoms():
    my_residue = atom.parent
    my_chain = my_residue.parent
    atom_chain[my_chain.id] += 1
    if my_residue.resname != 'HOH':
        atom_cnt[atom.element] += 1
    atom_res_types[my_residue.resname] += 1
print(dict(atom_res_types))
print(dict(atom_chain))
print(dict(atom_cnt))

这将打印出原子残基类型、每条链的原子数量和每种元素的数量，如下所示：

{' DT': 257, ' DC': 152, ' DA': 270, ' DG': 176, 'HOH': 384, 'SER': 323, 'VAL': 315, 'PRO': 294, 'GLN': 189, 'LYS': 135, 'THR': 294, 'TYR': 288, 'GLY': 156, 'PHE': 165, 'ARG': 561, 'LEU': 336, 'HIS': 210, 'ALA': 105, 'CYS': 180, 'ASN': 216, 'MET': 144, 'TRP': 42, 'ASP': 192, 'ILE': 144, 'GLU': 297, ' ZN': 3}
 {'E': 442, 'F': 449, 'A': 1734, 'B': 1593, 'C': 1610}
 {'O': 1114, 'C': 3238, 'N': 1001, 'P': 40, 'S': 48, 'ZN': 3}

请注意，前面提到的残基数不是正确的残基数，而是某个残基类型被引用的次数（它加起来等于原子数，而不是残基数）。

注意水（W）、核苷酸（DA，DC，DG，DT）和锌（ZN）残基，它们与氨基酸残基一起出现。

现在，让我们统计每个残基的实例数量和每条链的残基数量：

res_types = defaultdict(int)
res_per_chain = defaultdict(int)
for residue in p53_1tup.get_residues():
res_types[residue.resname] += 1
res_per_chain[residue.parent.id] +=1
print(dict(res_types))
print(dict(res_per_chain))

以下是输出结果：

{' DT': 13, ' DC': 8, ' DA': 13, ' DG': 8, 'HOH': 384, 'SER': 54, 'VAL': 45, 'PRO': 42, 'GLN': 21, 'LYS': 15, 'THR': 42, 'TYR': 24, 'GLY': 39, 'PHE': 15, 'ARG': 51, 'LEU': 42, 'HIS': 21, 'ALA': 21, 'CYS': 30, 'ASN': 27, 'MET': 18, 'TRP': 3, 'ASP': 24, 'ILE': 18, 'GLU': 33, ' ZN': 3}
 {'E': 43, 'F': 35, 'A': 395, 'B': 265, 'C': 276}

我们还可以获取一组原子的边界：

import sys
def get_bounds(my_atoms):
    my_min = [sys.maxsize] * 3
    my_max = [-sys.maxsize] * 3
    for atom in my_atoms:
        for i, coord in enumerate(atom.coord):
            if coord < my_min[i]:
                my_min[i] = coord
            if coord > my_max[i]:
                my_max[i] = coord
    return my_min, my_max
chain_bounds = {}
for chain in p53_1tup.get_chains():
    print(chain.id, get_bounds(chain.get_atoms()))
    chain_bounds[chain.id] = get_bounds(chain.get_atoms())
print(get_bounds(p53_1tup.get_atoms()))

一组原子可以是整个模型、一条链、一种残基或任何你感兴趣的子集。在这种情况下，我们将打印所有链条和整个模型的边界。数字表达不太直观，因此我们将采用更具图形化的方式。

为了了解每条链的大小，绘图可能比前面代码中的数字更具信息量：

import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D
fig = plt.figure(figsize=(16, 9))
ax3d = fig.add_subplot(111, projection='3d')
ax_xy = fig.add_subplot(331)
ax_xy.set_title('X/Y')
ax_xz = fig.add_subplot(334)
ax_xz.set_title('X/Z')
ax_zy = fig.add_subplot(337)
ax_zy.set_title('Z/Y')
color = {'A': 'r', 'B': 'g', 'C': 'b', 'E': '0.5', 'F': '0.75'}
zx, zy, zz = [], [], []
for chain in p53_1tup.get_chains():
    xs, ys, zs = [], [], []
    for residue in chain.get_residues():
        ref_atom = next(residue.get_iterator())
        x, y, z = ref_atom.coord
        if ref_atom.element == 'ZN':
            zx.append(x)
            zy.append(y)
            zz.append(z)
            continue
        xs.append(x)
        ys.append(y)
        zs.append(z)
    ax3d.scatter(xs, ys, zs, color=color[chain.id])
    ax_xy.scatter(xs, ys, marker='.', color=color[chain.id])
    ax_xz.scatter(xs, zs, marker='.', color=color[chain.id])
    ax_zy.scatter(zs, ys, marker='.', color=color[chain.id])
ax3d.set_xlabel('X')
ax3d.set_ylabel('Y')
ax3d.set_zlabel('Z')
ax3d.scatter(zx, zy, zz, color='k', marker='v', s=300)
ax_xy.scatter(zx, zy, color='k', marker='v', s=80)
ax_xz.scatter(zx, zz, color='k', marker='v', s=80)
ax_zy.scatter(zz, zy, color='k', marker='v', s=80)
for ax in [ax_xy, ax_xz, ax_zy]:
    ax.get_yaxis().set_visible(False)
    ax.get_xaxis().set_visible(False)

目前有很多分子可视化工具。实际上，我们稍后会讨论 PyMOL。不过，matplotlib对于简单的可视化已经足够了。关于matplotlib最重要的一点是它非常稳定，且容易集成到可靠的生产代码中。

在下图中，我们对链条进行了三维绘制，DNA 部分为灰色，蛋白质链条用不同颜色表示。我们还在下图的左侧绘制了平面投影（X/Y，X/Z，和Z/Y）：

图 8.2 - 蛋白质链的空间分布——主图是一个 3D 图，左侧子图是平面视图（X/Y，X/Z，和 Z/Y）

计算 PDB 文件中的分子距离

在这里，我们将找到与1TUP模型中三个锌原子接近的原子。我们将考虑这些锌原子之间的几种距离，并借此机会讨论算法的性能。

准备工作

你可以在Chapter08/Distance.py Notebook 文件中找到这些内容。

如何操作…

请查看以下步骤：

让我们加载我们的模型，如下所示：

from Bio import PDB
repository = PDB.PDBList()
parser = PDB.PDBParser()
repository.retrieve_pdb_file('1TUP', pdir='.', file_format='pdb')
p53_1tup = parser.get_structure('P 53', 'pdb1tup.ent')

现在，我们来提取锌原子，后续我们将以这些原子为基准进行比较：

zns = []for atom in p53_1tup.get_atoms():
if atom.element == 'ZN':
zns.append(atom)
for zn in zns:
    print(zn, zn.coord)

你应该能看到三个锌原子。

现在，让我们定义一个函数来获取一个原子与一组其他原子之间的距离，如下所示：

import math
def get_closest_atoms(pdb_struct, ref_atom, distance):
    atoms = {}
    rx, ry, rz = ref_atom.coord
    for atom in pdb_struct.get_atoms():
        if atom == ref_atom:
            continue
        x, y, z = atom.coord
        my_dist = math.sqrt((x - rx)**2 + (y - ry)**2 + (z - rz)**2)
        if my_dist < distance:
            atoms[atom] = my_dist
    return atoms

我们获取参考原子的坐标，然后遍历我们希望比较的原子列表。如果某个原子足够接近，它会被添加到return列表中。

现在我们计算接近锌原子的原子，距离在我们的模型中可以达到 4 埃（Å）：

for zn in zns:
    print()
    print(zn.coord)
    atoms = get_closest_atoms(p53_1tup, zn, 4)
    for atom, distance in atoms.items():
        print(atom.element, distance, atom.coord)

这里，我们展示了第一个锌原子的结果，包括元素、距离和坐标：

[58.108 23.242 57.424]
 C 3.4080117696286854 [57.77  21.214 60.142]
 S 2.3262243799594877 [57.065 21.452 58.482]
 C 3.4566537492335123 [58.886 20.867 55.036]
 C 3.064120559761192 [58.047 22.038 54.607]
 N 1.9918273537290707 [57.755 23.073 55.471]
 C 2.9243719601324525 [56.993 23.943 54.813]
 C 3.857729198122736 [61.148 25.061 55.897]
 C 3.62725094648044 [61.61  24.087 57.001]
 S 2.2789209624943494 [60.317 23.318 57.979]
 C 3.087214470667822 [57.205 25.099 59.719]
 S 2.2253158446520818 [56.914 25.054 57.917]

我们只有三个锌原子，因此计算量大大减少。然而，假设我们有更多原子，或者我们在整个原子集之间进行成对比较（记住，在成对比较的情况下，比较次数是随着原子数量的平方增长的）。虽然我们的案例较小，但预测使用案例并不难，更多的比较会消耗大量时间。我们很快会回到这个问题。

让我们看看随着距离增加，我们会得到多少原子：

for distance in [1, 2, 4, 8, 16, 32, 64, 128]:
    my_atoms = []
    for zn in zns:
        atoms = get_closest_atoms(p53_1tup, zn, distance)
        my_atoms.append(len(atoms))
    print(distance, my_atoms)

结果如下：

1 [0, 0, 0]
2 [1, 0, 0]
4 [11, 11, 12]
8 [109, 113, 106]
16 [523, 721, 487]
32 [2381, 3493, 2053]
64 [5800, 5827, 5501]
128 [5827, 5827, 5827]

如我们之前所见，这个特定的情况并不太昂贵，但我们还是计时看看：

import timeit
nexecs = 10
print(timeit.timeit('get_closest_atoms(p53_1tup, zns[0], 4.0)',
      'from __main__ import get_closest_atoms, p53_1tup, zns',
      number=nexecs) / nexecs * 1000)

在这里，我们将使用timeit模块执行这个函数 10 次，然后以毫秒为单位打印结果。我们将函数作为字符串传递，并传递另一个包含必要导入的字符串以使函数正常工作。在 Notebook 中，你可能知道%timeit魔法命令，它能让你的生活变得更加轻松。在测试代码的机器上，这大约需要 40 毫秒。显然，在你的电脑上，你会得到稍有不同的结果。

我们能做得更好吗？让我们考虑一个不同的distance函数，如下面的代码所示：

def get_closest_alternative(pdb_struct, ref_atom, distance):
    atoms = {}
    rx, ry, rz = ref_atom.coord
    for atom in pdb_struct.get_atoms():
        if atom == ref_atom:
            continue
        x, y, z = atom.coord
        if abs(x - rx) > distance or abs(y - ry) > distance or abs(z - rz) > distance:
            continue
        my_dist = math.sqrt((x - rx)**2 + (y - ry)**2 + (z - rz)**2)
        if my_dist < distance:
            atoms[atom] = my_dist
    return atoms

所以，我们对原始函数进行修改，加入一个非常简单的if条件来处理距离。这样做的理由是，平方根计算和可能的浮点幂运算非常昂贵，因此我们将尽量避免它。然而，对于任何维度中距离小于目标距离的原子，这个函数会变得更加昂贵。

现在，让我们来计时：

print(timeit.timeit('get_closest_alternative(p53_1tup, zns[0], 4.0)',
      'from __main__ import get_closest_alternative, p53_1tup, zns',
      number=nexecs) / nexecs * 1000)

在我们之前使用的同一台机器上，它需要 16 毫秒，这意味着它大约快了三倍。

然而，这总是更好吗？让我们比较不同距离下的成本，如下所示：

print('Standard')
for distance in [1, 4, 16, 64, 128]:
    print(timeit.timeit('get_closest_atoms(p53_1tup, zns[0], distance)',
          'from __main__ import get_closest_atoms, p53_1tup, zns, distance',
          number=nexecs) / nexecs * 1000)
print('Optimized')
for distance in [1, 4, 16, 64, 128]:
    print(timeit.timeit('get_closest_alternative(p53_1tup, zns[0], distance)',
          'from __main__ import get_closest_alternative, p53_1tup, zns, distance',
          number=nexecs) / nexecs * 1000)

结果显示在以下输出中：

Standard
 85.08649739999328
 86.50681579999855
 86.79630599999655
 96.95437099999253
 96.21982420001132
 Optimized
 30.253444099980698
 32.69531210000878
 52.965772600009586
 142.53310030001103
 141.26269519999823

请注意，标准版本的成本大致是恒定的，而优化版本的成本则取决于最近原子的距离；距离越大，使用额外if和平方根进行计算的情况就越多，从而使得函数变得更昂贵。

这里的关键点是，你可能可以通过聪明的计算捷径编写更高效的函数，但复杂度成本可能会发生质的变化。在前面的例子中，我建议第二个函数在所有现实且有趣的情况下更高效，特别是在你尝试找到最接近的原子时。然而，在设计你自己的优化算法时，你必须小心。

执行几何操作

现在我们将使用几何信息进行计算，包括计算链条和整个模型的质心。

准备工作

你可以在 Chapter08/Mass.py 笔记本文件中找到这些内容。

如何操作...

让我们看一下以下步骤：

首先，让我们获取数据：

from Bio import PDB
repository = PDB.PDBList()
parser = PDB.PDBParser()
repository.retrieve_pdb_file('1TUP', pdir='.', file_format='pdb')
p53_1tup = parser.get_structure('P 53', 'pdb1tup.ent')

然后，使用以下代码回顾我们拥有的残基类型：

my_residues = set()
for residue in p53_1tup.get_residues():
    my_residues.add(residue.id[0])
print(my_residues)

所以，我们有 H_ ZN（锌）和 W（水），它们是 HETATM 类型；绝大多数是标准 PDB 原子。

使用以下代码计算所有链条、锌和水的质量：

def get_mass(atoms, accept_fun=lambda atom: atom.parent.id[0] != 'W'):
    return sum([atom.mass for atom in atoms if accept_fun(atom)])
chain_names = [chain.id for chain in p53_1tup.get_chains()]
my_mass = np.ndarray((len(chain_names), 3))
for i, chain in enumerate(p53_1tup.get_chains()):
    my_mass[i, 0] = get_mass(chain.get_atoms())
    my_mass[i, 1] = get_mass(chain.get_atoms(),
        accept_fun=lambda atom: atom.parent.id[0] not in [' ', 'W'])
    my_mass[i, 2] = get_mass(chain.get_atoms(),
        accept_fun=lambda atom: atom.parent.id[0] == 'W')
masses = pd.DataFrame(my_mass, index=chain_names, columns=['No Water','Zincs', 'Water'])
print(masses)

get_mass 函数返回通过接收标准函数筛选的列表中所有原子的质量。这里，默认的接收标准是排除水分子残基。

然后，我们计算所有链条的质量。我们有三个版本：只有氨基酸、锌和水。锌仅在此模型中检测每条链中的单个原子。输出结果如下：

图 8.3 - 所有蛋白质链的质量

让我们计算模型的几何中心和质心，如下所示：

def get_center(atoms,
    weight_fun=lambda atom: 1 if atom.parent.id[0] != 'W' else 0):
    xsum = ysum = zsum = 0.0
    acum = 0.0
    for atom in atoms:
        x, y, z = atom.coord
        weight = weight_fun(atom)
        acum += weight
        xsum += weight * x
        ysum += weight * y
        zsum += weight * z
    return xsum / acum, ysum / acum, zsum / acum
print(get_center(p53_1tup.get_atoms()))
print(get_center(p53_1tup.get_atoms(),
    weight_fun=lambda atom: atom.mass if atom.parent.id[0] != 'W' else 0))

首先，我们定义一个加权函数来获取质心的坐标。默认的函数会将所有原子视为相等，只要它们不是水分子残基。

然后，我们通过重新定义 weight 函数，将每个原子的值设置为其质量，来计算几何中心和质心。几何中心的计算不考虑分子量。

例如，你可能想要计算没有 DNA 链的蛋白质的质心。

让我们计算每个链条的质心和几何中心，如下所示：

my_center = np.ndarray((len(chain_names), 6))
for i, chain in enumerate(p53_1tup.get_chains()):
    x, y, z = get_center(chain.get_atoms())
    my_center[i, 0] = x
    my_center[i, 1] = y
    my_center[i, 2] = z
    x, y, z = get_center(chain.get_atoms(),
        weight_fun=lambda atom: atom.mass if atom.parent.id[0] != 'W' else 0)
    my_center[i, 3] = x
    my_center[i, 4] = y
    my_center[i, 5] = z
weights = pd.DataFrame(my_center, index=chain_names,
    columns=['X', 'Y', 'Z', 'X (Mass)', 'Y (Mass)', 'Z (Mass)'])
print(weights)

结果如图所示：

图 8.4 - 每个蛋白质链的质心和几何中心

还有更多

虽然这本书不是基于蛋白质结构解析技术的，但需要记住的是，X 射线晶体学方法无法检测氢原子，因此残基的质量计算可能基于非常不准确的模型；有关更多信息，请参阅 www.umass.edu/microbio/chime/pe_beta/pe/protexpl/help_hyd.htm。

使用 PyMOL 动画

在这里，我们将创建一个关于 p53 1TUP 模型的视频。为此，我们将使用 PyMOL 可视化库。我们将通过移动 p53 1TUP 模型开始动画，然后进行缩放；随着缩放，我们改变渲染策略，以便可以更深入地观察模型。你可以在 odysee.com/@Python:8/protein_video:8 找到你将生成的视频版本。

准备工作

这个配方将以 Python 脚本的形式呈现，而不是以 Notebook 的形式。这主要是因为输出不是交互式的，而是一组需要进一步后期处理的图像文件。

你需要安装 PyMOL（www.pymol.org）。在 Debian、Ubuntu 或 Linux 系统中，你可以使用apt-get install pymol命令。如果你使用的是 Conda，我建议不要使用它，因为依赖项会很容易解决——而且你将安装一个仅限 30 天试用的版本，需要许可证，而上述版本是完全开源的。如果你不是使用 Debian 或 Linux，我建议你安装适用于你操作系统的开源版本。

PyMOL 更多是一个交互式程序，而非 Python 库，因此我强烈建议你在继续操作前先进行一些探索。这可以是很有趣的！这个配方的代码可以在 GitHub 仓库中找到，作为脚本文件以及本章的 Notebook 文件，位于Chapter08。我们将在这个配方中使用PyMol_Movie.py文件。

如何操作...

请查看以下步骤：

让我们初始化并获取我们的 PDB 模型，并准备渲染，如下所示：

import pymol
from pymol import cmd
#pymol.pymol_argv = ['pymol', '-qc'] # Quiet / no GUI
pymol.finish_launching()
cmd.fetch('1TUP', async=False)
cmd.disable('all')
cmd.enable('1TUP')
cmd.hide('all')
cmd.show('sphere', 'name zn')

请注意，pymol_argv这一行会使代码保持静默。在第一次执行时，你可能想要注释掉这行代码，看看用户界面。

对于电影渲染，这将非常有用（如我们将很快看到的）。作为一个库，PyMOL 的使用相当复杂。例如，导入后，你需要调用finish_launching。接着，我们获取我们的 PDB 文件。

接下来是一些 PyMOL 命令。许多关于交互式使用的网页指南对于理解发生的事情非常有帮助。在这里，我们将启用所有模型以供查看，隐藏所有模型（因为默认视图是线条表示，这样不够好），然后将锌渲染为球体。

在这个阶段，除锌外，其他都不可见。

为了渲染我们的模型，我们将使用以下三个场景：

cmd.show('surface', 'chain A+B+C')
cmd.show('cartoon', 'chain E+F')
cmd.scene('S0', action='store', view=0, frame=0, animate=-1)
cmd.show('cartoon')
cmd.hide('surface')
cmd.scene('S1', action='store', view=0, frame=0, animate=-1)
cmd.hide('cartoon', 'chain A+B+C')
cmd.show('mesh', 'chain A')
cmd.show('sticks', 'chain A+B+C')
cmd.scene('S2', action='store', view=0, frame=0, animate=-1)

我们需要定义两个场景。一个场景对应于我们围绕蛋白质移动（基于表面，因此是透明的），另一个场景则对应于我们深入观察（基于卡通式）。DNA 始终以卡通形式渲染。

我们还定义了第三个场景，当我们在最后进行缩小时使用。蛋白质将被渲染为棒状，并且我们将给 A 链添加一个网格，以便更清楚地展示它与 DNA 的关系。

让我们定义视频的基本参数，如下所示：
```
cmd.set('ray_trace_frames', 0)
cmd.mset(1, 500)
```

我们定义了默认的光线追踪算法。这一行并非必需，但请尝试将数字增加到1、2或3，并准备好等待很长时间。

如果你打开了 OpenGL 接口（即图形界面），那么只能使用0，因此对于这个快速版本，你需要打开 GUI（pymol_argv应该保持注释状态）。

然后，我们通知 PyMOL 我们将使用 500 帧。

在前 150 帧中，我们使用初始场景围绕蛋白质移动。我们稍微环绕模型，然后使用以下代码接近 DNA。

cmd.frame(0)
cmd.scene('S0')
cmd.mview()
cmd.frame(60)
cmd.set_view((-0.175534308,   -0.331560850,   -0.926960170,
             0.541812420,     0.753615797,   -0.372158051,
             0.821965039,    -0.567564785,    0.047358301,
             0.000000000,     0.000000000, -249.619018555,
             58.625568390,   15.602619171,   77.781631470,
             196.801528931, 302.436492920,  -20.000000000))
cmd.mview()
cmd.frame(90)
cmd.set_view((-0.175534308,   -0.331560850,   -0.926960170,
              0.541812420,    0.753615797,   -0.372158051,
              0.821965039,   -0.567564785,    0.047358301,
              -0.000067875,    0.000017881, -249.615447998,
              54.029174805,   26.956727982,   77.124832153,
             196.801528931,  302.436492920,  -20.000000000))
cmd.mview()
cmd.frame(150)
cmd.set_view((-0.175534308,   -0.331560850,   -0.926960170,
              0.541812420,    0.753615797,   -0.372158051,
              0.821965039,   -0.567564785,    0.047358301,
              -0.000067875,    0.000017881,  -55.406421661,
              54.029174805,   26.956727982,   77.124832153,
              2.592475891,  108.227416992,  -20.000000000))
cmd.mview()

我们定义了三个点；前两个点与 DNA 对齐，最后一个点进入其中。我们通过在交互模式下使用 PyMOL，使用鼠标和键盘导航，并使用 get_view 命令来获取坐标（所有这些数字），然后可以剪切并粘贴。

第一个帧如下所示：

图 8.5 - 第 0 帧和场景 DS0

我们现在更改场景，为进入蛋白质内部做准备：
```
cmd.frame(200)
cmd.scene('S1')
cmd.mview()
```

以下截图显示了当前的位置：

图 8.6 - DNA 分子附近的第 200 帧和场景 S1

我们进入蛋白质内部，并在结束时使用以下代码更改场景：

cmd.frame(350)
cmd.scene('S1')
cmd.set_view((0.395763457,   -0.173441306,    0.901825786,
              0.915456235,    0.152441502,   -0.372427106,
             -0.072881661,    0.972972929,    0.219108686,
              0.000070953,    0.000013039,  -37.689743042,
             57.748500824,   14.325904846,   77.241867065,
             -15.123448372,   90.511535645,  -20.000000000))
cmd.mview()
cmd.frame(351)
cmd.scene('S2')
cmd.mview()

我们现在完全进入其中，如以下截图所示：

图 8.7 - 第 350 帧 - 场景 S1 即将更改为 S2

最后，我们让 PyMOL 返回到原始位置，然后播放、保存并退出：

cmd.frame(500)
cmd.scene('S2')
cmd.mview()
cmd.mplay()
cmd.mpng('p53_1tup')
cmd.quit()

这将生成 500 个以 p53_1tup 为前缀的 PNG 文件。

这是一个接近结束的帧（450）：

图 8.8 - 第 450 帧和场景 S2

还有更多

该 YouTube 视频是在 Linux 上使用 ffmpeg 以每秒 15 帧的速度生成的，如下所示：

ffmpeg -r 15 -f image2 -start_number 1 -i "p53_1tup%04d.png" example.mp4

有很多应用程序可以用来从图像生成视频。PyMOL 可以生成 MPEG 格式的视频，但需要安装额外的库。

PyMOL 是为从其控制台交互式使用而创建的（可以在 Python 中扩展）。反向使用（从 Python 导入且没有图形界面）可能会变得复杂且令人沮丧。PyMOL 启动一个单独的线程来渲染图像，这些图像异步工作。

例如，这意味着你的代码可能与渲染器的位置不同。我已经在 GitHub 仓库中放了另一个名为 PyMol_Intro.py 的脚本；你会看到第二个 PNG 调用将在第一个还没有完成之前就开始。试试这个脚本代码，看看它应该如何运行，以及它实际是如何运行的。

从 GUI 角度来看，PyMOL 有很多优秀的文档，访问 www.pymolwiki.org/index.php/MovieSchool 可以获得。如果你想制作电影，这是一个很好的起点，www.pymolwiki.org 是一个信息的宝库。

使用 Biopython 解析 mmCIF 文件

mmCIF 文件格式可能是未来的趋势。Biopython 目前还没有完全支持它的功能，但我们将看看当前有哪些功能。

正在准备

由于Bio.PDB不能自动下载 mmCIF 文件，你需要获取你的蛋白质文件并将其重命名为1tup.cif。它可以在github.com/PacktPublishing/Bioinformatics-with-Python-Cookbook-third-Edition/blob/master/Datasets.py中的1TUP.cif找到。

你可以在Chapter08/mmCIF.py笔记本文件中找到这些内容。

如何操作...

看一下以下步骤：

让我们解析文件。我们只需使用 MMCIF 解析器，而不是 PDB 解析器：

from Bio import PDB
parser = PDB.MMCIFParser()
p53_1tup = parser.get_structure('P53', '1tup.cif')

让我们检查以下链条：

def describe_model(name, pdb):
    print()
    for model in p53_1tup:
        for chain in model:
            print('%s - Chain: %s. Number of residues: %d. Number of atoms: %d.' %
                  (name, chain.id, len(chain), len(list(chain.get_atoms()))))
describe_model('1TUP', p53_1tup)

输出将如下所示：

1TUP - Chain: E. Number of residues: 43\. Number of atoms: 442.
1TUP - Chain: F. Number of residues: 35\. Number of atoms: 449.
1TUP - Chain: A. Number of residues: 395\. Number of atoms: 1734.
1TUP - Chain: B. Number of residues: 265\. Number of atoms: 1593.
1TUP - Chain: C. Number of residues: 276\. Number of atoms: 1610.

许多字段在解析的结构中不可用，但可以通过使用更低级别的字典来检索这些字段，如下所示：

mmcif_dict = PDB.MMCIF2Dict.MMCIF2Dict('1tup.cif')
for k, v in mmcif_dict.items():
    print(k, v)
    print()

不幸的是，这个列表很大，需要一些后处理才能理解它，但它是可以获得的。

还有更多内容

你仍然可以获取来自 Biopython 提供的 mmCIF 文件的所有模型信息，因此解析器仍然非常有用。我们可以预期mmCIF解析器会有更多的开发，而不是PDB解析器。

这有一个 Python 库，由 PDB 的开发者提供，网址为mmcif.wwpdb.org/docs/sw-examples/python/html/index.xhtml。

第十章：生物信息学管道

管道在任何数据科学环境中都至关重要。数据处理从来不是一项单一的任务。许多管道是通过临时脚本实现的。虽然这种方式可以有效使用，但在许多情况下，它们未能满足几个基本视角，主要是可重复性、可维护性和可扩展性。

在生物信息学中，您可以找到三种主要类型的管道系统：

像 Galaxy (usegalaxy.org) 这样的框架，专为用户设计，即它们提供了易于使用的用户界面，并隐藏了大部分底层机制。
编程工作流 — 针对代码接口的工作流，虽然是通用的，但来源于生物信息学领域。两个例子是 Snakemake (snakemake.readthedocs.io/) 和 Nextflow (www.nextflow.io/)。
完全通用的工作流系统，例如 Apache Airflow (airflow.incubator.apache.org/)，它采用一种不太以数据为中心的工作流管理方式。

在本章中，我们将讨论 Galaxy，尤其是对于那些支持不太倾向于编写代码解决方案的用户的生物信息学专家而言，Galaxy 尤为重要。虽然您可能不是这些管道系统的典型用户，但您仍然可能需要为它们提供支持。幸运的是，Galaxy 提供了 API，这是我们主要关注的内容。

我们还将讨论 Snakemake 和 Nextflow 作为通用工作流工具，这些工具有编程接口，最初来源于生物信息学领域。我们将涵盖这两种工具，因为它们是该领域中最常见的工具。我们将用 Snakemake 和 Nextflow 解决一个类似的生物信息学问题。我们将体验这两个框架，并希望能够决定哪个是我们的最爱。

这些配方的代码不是以笔记本的形式呈现，而是作为 Python 脚本，存放在本书仓库的Chapter09目录下。

本章中，您将找到以下配方：

介绍 Galaxy 服务器
使用 API 访问 Galaxy
使用 Snakemake 开发变异分析管道
使用 Nextflow 开发变异分析管道

介绍 Galaxy 服务器

Galaxy (galaxyproject.org/tutorials/g101/) 是一个开源系统，旨在帮助非计算用户进行计算生物学研究。它是目前最广泛使用且用户友好的管道系统。任何用户都可以在服务器上安装 Galaxy，但网络上也有许多公开访问的其他服务器，其中最著名的是 usegalaxy.org。

接下来食谱的重点将是 Galaxy 的编程方面：使用 Galaxy API 进行接口连接，以及开发 Galaxy 工具以扩展其功能。在开始之前，强烈建议你先作为用户使用 Galaxy。你可以通过在usegalaxy.org创建一个免费账户并稍作体验来实现这一点。建议你了解工作流的基本知识。

准备就绪

在此食谱中，我们将使用 Docker 在本地安装 Galaxy 服务器。因此，需要一个本地 Docker 安装。不同操作系统的复杂度不同：在 Linux 上容易，在 macOS 上中等，在 Windows 上中到难。

此安装推荐用于接下来的两个食谱，但你也可以使用现有的公共服务器。请注意，公共服务器的接口可能会随时间变化，因此今天有效的内容明天可能会失效。如何使用公共服务器来执行接下来的两个食谱的说明，详见*更多信息...*部分。

如何操作…

看一下以下步骤。这些步骤假设你已经启用了 Docker 命令行：

首先，我们使用以下命令拉取 Galaxy Docker 镜像：
```
docker pull bgruening/galaxy-stable:20.09
```

这将拉取 Björn Grüning 的精彩 Docker Galaxy 镜像。请使用20.09标签，如前面命令所示；任何更新的版本可能会破坏此食谱和下一个食谱。

在系统上创建一个目录。此目录将保存 Docker 容器在多次运行中的持久输出。

注意

Docker 容器对于磁盘空间是临时的。这意味着当你停止容器时，所有磁盘上的更改都会丢失。可以通过在 Docker 中挂载来自主机的卷来解决此问题，如下一个步骤所示。挂载卷中的所有内容将会持久化。

现在我们可以通过以下命令运行镜像：

docker run -d -v YOUR_DIRECTORY:/export -p 8080:80 -p 8021:21 bgruening/galaxy-stable:20.09

将YOUR_DIRECTORY替换为你在步骤 2中创建的目录的完整路径。如果前面的命令失败，请确保你有权限运行 Docker。不同操作系统的权限设置可能不同。

检查YOUR_DIRECTORY中的内容。第一次运行镜像时，它会创建所有需要的文件，以便在 Docker 运行之间持久化执行。这意味着会保持用户数据库、数据集和工作流。

将浏览器指向http://localhost:8080。如果遇到任何错误，稍等几秒钟。你应该能看到如下屏幕：

图 9.1 - Galaxy Docker 首页

现在使用默认的用户名和密码组合：admin和password，登录（参见顶部栏）。
从顶部菜单中选择用户，然后选择首选项。
现在，选择管理 API 密钥。

不要更改 API 密钥。前面的练习目的是让你知道 API 密钥的位置。在实际场景中，你需要进入这个页面获取你的密钥。请记下 API 密钥：fakekey。顺便说一句，正常情况下，这将是一个 MD5 哈希。

到目前为止，我们已经在服务器上安装了以下（默认）凭据：用户为admin，密码为password，API 密钥为fakekey。访问点是localhost:8080。

还有更多内容

Björn Grüning 的镜像将在本章中使用的方式相当简单；毕竟，这不是一本关于系统管理或 DevOps 的书，而是一本关于编程的书。如果你访问github.com/bgruening/docker-galaxy-stable，你会看到有无数种配置镜像的方式，而且都有详细的文档。在这里，我们简单的方法足够满足我们的开发需求。

如果你不想在本地计算机上安装 Galaxy，可以使用公共服务器，例如usegalaxy.org来进行下一步操作。这个方法不是 100%万无一失的，因为服务会随时间变化，但它可能会非常接近。请按以下步骤进行：

在公共服务器上创建一个帐户（usegalaxy.org或其他服务器）。
请按照之前的说明获取你的 API 密钥。
在下一个步骤中，你需要替换主机、用户、密码和 API 密钥。

使用 API 访问 Galaxy

虽然 Galaxy 的主要用法是通过易于使用的 Web 界面，但它也提供了一个 REST API 供程序化访问。它提供了多种语言的接口，例如，Python 支持可以通过 BioBlend 获得(bioblend.readthedocs.io)。

在这里，我们将开发一个脚本，加载一个 BED 文件到 Galaxy，并调用一个工具将其转换为 GFF 格式。我们将通过 Galaxy 的 FTP 服务器加载文件。

准备工作

如果你没有完成前面的步骤，请阅读相应的*更多内容...*部分。代码已在前面步骤中准备的本地服务器上进行了测试，因此如果你在公共服务器上运行，可能需要进行一些调整。

我们的代码将需要在 Galaxy 服务器上进行身份验证，以执行必要的操作。由于安全性是一个重要问题，本教程在这方面不会完全天真。我们的脚本将通过 YAML 文件进行配置，例如：

rest_protocol: http
server: localhost
rest_port: 8080
sftp_port: 8022
user: admin
password: password
api_key: fakekey

我们的脚本不会接受这个文件作为纯文本，而是要求它被加密。也就是说，我们的安全计划中存在一个大漏洞：我们将使用 HTTP（而不是 HTTPS），这意味着密码将在网络上传输时以明文形式发送。显然，这不是一个好的解决方案，但考虑到空间的限制，我们只能做到这一点（特别是在前面的步骤中）。真正安全的解决方案需要使用 HTTPS。

我们将需要一个脚本，该脚本接受一个 YAML 文件并生成加密版本：

import base64
import getpass
from io import StringIO
import os
from ruamel.yaml import YAML
from cryptography.fernet import Fernet
from cryptography.hazmat.backends import default_backend
from cryptography.hazmat.primitives import hashes
from cryptography.hazmat.primitives.kdf.pbkdf2 import PBKDF2HMAC
password = getpass.getpass('Please enter the password:').encode()
salt = os.urandom(16)
kdf = PBKDF2HMAC(algorithm=hashes.SHA256(), length=32, salt=salt,
                 iterations=100000, backend=default_backend())
key = base64.urlsafe_b64encode(kdf.derive(password))
fernet = Fernet(key)
with open('salt', 'wb') as w:
    w.write(salt)
yaml = YAML()
content = yaml.load(open('galaxy.yaml', 'rt', encoding='utf-8'))
print(type(content), content)
output = StringIO()
yaml.dump(content, output)
print ('Encrypting:\n%s' % output.getvalue())
enc_output = fernet.encrypt(output.getvalue().encode())
with open('galaxy.yaml.enc', 'wb') as w:
    w.write(enc_output)

前面的文件可以在 GitHub 仓库中的Chapter09/pipelines/galaxy/encrypt.py找到。

你将需要输入一个密码来进行加密。

前面的代码与 Galaxy 无关：它读取一个YAML文件，并使用用户提供的密码进行加密。它使用cryptography模块进行加密，并用ruaml.yaml处理YAML文件。输出两个文件：加密后的YAML文件和加密用的salt文件。出于安全考虑，salt文件不应公开。

这种保护凭证的方法远非复杂；它主要是为了说明在处理认证令牌时，你必须小心代码中硬编码的安全凭证。在网络上有很多硬编码安全凭证的实例。

如何操作…

请查看以下步骤，它们可以在Chapter09/pipelines/galaxy/api.py中找到：

我们从解密我们的配置文件开始。我们需要提供一个密码：

import base64
from collections import defaultdict
import getpass
import pprint
import warnings
from ruamel.yaml import YAML
from cryptography.fernet import Fernet
from cryptography.hazmat.backends import default_backend
from cryptography.hazmat.primitives import hashes
from cryptography.hazmat.primitives.kdf.pbkdf2 import PBKDF2HMAC
import pandas as pd
Import paramiko
from bioblend.galaxy import GalaxyInstance
pp = pprint.PrettyPrinter()
warnings.filterwarnings('ignore')
# explain above, and warn
with open('galaxy.yaml.enc', 'rb') as f:
    enc_conf = f.read()
password = getpass.getpass('Please enter the password:').encode()
with open('salt', 'rb') as f:
    salt = f.read()
kdf = PBKDF2HMAC(algorithm=hashes.SHA256(), length=32, salt=salt,
                 iterations=100000, backend=default_backend())
key = base64.urlsafe_b64encode(kdf.derive(password))
fernet = Fernet(key)
yaml = YAML()
conf = yaml.load(fernet.decrypt(enc_conf).decode())

最后一行总结了所有内容：YAML模块将从解密后的文件中加载配置。请注意，我们还读取了salt，以便能够解密文件。

现在我们将获取所有配置变量，准备服务器 URL，并指定我们将要创建的 Galaxy 历史记录的名称（bioinf_example）：

server = conf['server']
rest_protocol = conf['rest_protocol']
rest_port = conf['rest_port']
user = conf['user']
password = conf['password']
ftp_port = int(conf['ftp_port'])
api_key = conf['api_key']
rest_url = '%s://%s:%d' % (rest_protocol, server, rest_port)
history_name = 'bioinf_example'

最后，我们能够连接到 Galaxy 服务器：

gi = GalaxyInstance(url=rest_url, key=api_key)
gi.verify = False

我们将列出所有可用的历史记录：

histories = gi.histories
print('Existing histories:')
for history in histories.get_histories():
    if history['name'] == history_name:
        histories.delete_history(history['id'])
    print('  - ' + history['name'])
print()

在第一次执行时，你将获得一个未命名的历史记录，但在之后的执行中，你还将获得bioinf_example，我们将在此阶段删除它，以便从干净的状态开始。

之后，我们创建了bioinf_example历史记录：

ds_history = histories.create_history(history_name)

如果你愿意，你可以在 Web 界面上检查，你会在那里看到新的历史记录。

现在我们将上传文件；这需要一个 SFTP 连接。文件通过以下代码提供：

print('Uploading file')
transport = paramiko.Transport((server, sftp_port))
transport.connect(None, user, password)
sftp = paramiko.SFTPClient.from_transport(transport)
sftp.put('LCT.bed', 'LCT.bed')
sftp.close()
transport.close()

现在我们将告诉 Galaxy 将 FTP 服务器上的文件加载到其内部数据库中：
```
gi.tools.upload_from_ftp('LCT.bed', ds_history['id'])
```

让我们总结一下我们历史记录的内容：

def summarize_contents(contents):
 summary = defaultdict(list)
 for item in contents:
 summary['íd'].append(item['id'])
 summary['híd'].append(item['hid'])
 summary['name'].append(item['name'])
 summary['type'].append(item['type'])
 summary['extension'].append(item['extension'])
 return pd.DataFrame.from_dict(summary)
print('History contents:')
pd_contents = summarize_contents(contents)
print(pd_contents)
print()

我们只有一个条目：

                 íd  híd     name  type extension
0  f2db41e1fa331b3e    1  LCT.bed  file      auto

让我们检查一下BED文件的元数据：

print('Metadata for LCT.bed')
bed_ds = contents[0]
pp.pprint(bed_ds)
print()

结果包括以下内容：

{'create_time': '2018-11-28T21:27:28.952118',
 'dataset_id': 'f2db41e1fa331b3e',
 'deleted': False,
 'extension': 'auto',
 'hid': 1,
 'history_content_type': 'dataset',
 'history_id': 'f2db41e1fa331b3e',
 'id': 'f2db41e1fa331b3e',
 'name': 'LCT.bed',
 'purged': False,
 'state': 'queued',
 'tags': [],
 'type': 'file',
 'type_id': 'dataset-f2db41e1fa331b3e',
 'update_time': '2018-11-28T21:27:29.149933',
 'url': '/api/histories/f2db41e1fa331b3e/contents/f2db41e1fa331b3e',
 'visible': True}

让我们将注意力转向服务器上现有的工具，并获取有关它们的元数据：

print('Metadata about all tools')
all_tools = gi.tools.get_tools()
pp.pprint(all_tools)
print()

这将打印出一个长长的工具列表。

现在让我们获取一些关于我们工具的信息：

bed2gff = gi.tools.get_tools(name='Convert BED to GFF')[0]
print("Converter metadata:")
pp.pprint(gi.tools.show_tool(bed2gff['id'], io_details=True, link_details=True))
print()

工具的名称在前面的步骤中可用。请注意，我们获取的是列表中的第一个元素，因为理论上可能安装了多个版本的工具。简化后的输出如下：

{'config_file': '/galaxy-central/lib/galaxy/datatypes/converters/bed_to_gff_converter.xml',
 'id': 'CONVERTER_bed_to_gff_0',
 'inputs': [{'argument': None,
             'edam': {'edam_data': ['data_3002'],
                      'edam_formats': ['format_3003']},
             'extensions': ['bed'],
             'label': 'Choose BED file',
             'multiple': False,
             'name': 'input1',
             'optional': False,
             'type': 'data',
             'value': None}],
 'labels': [],
 'link': '/tool_runner?tool_id=CONVERTER_bed_to_gff_0',
 'min_width': -1,
 'model_class': 'Tool',
 'name': 'Convert BED to GFF',
 'outputs': [{'edam_data': 'data_1255',
              'edam_format': 'format_2305',
              'format': 'gff',
              'hidden': False,
              'model_class': 'ToolOutput',
              'name': 'output1'}],
 'panel_section_id': None,
 'panel_section_name': None,
 'target': 'galaxy_main',
 'version': '2.0.0'}

最后，让我们运行一个工具，将我们的BED文件转换为GFF：

def dataset_to_param(dataset):
    return dict(src='hda', id=dataset['id'])
tool_inputs = {
    'input1': dataset_to_param(bed_ds)
    }
gi.tools.run_tool(ds_history['id'], bed2gff['id'], tool_inputs=tool_inputs)

可以在前面的步骤中检查工具的参数。如果你进入 Web 界面，你会看到类似以下的内容：

图 9.2 - 通过 Galaxy 的 Web 界面检查我们脚本的结果

因此，我们通过其 REST API 访问了 Galaxy。

使用 Snakemake 部署变异分析管道

Galaxy 主要面向那些不太倾向于编程的用户。即使你更喜欢编程友好的环境，了解如何使用 Galaxy 仍然很重要，因为它的广泛应用。幸运的是，Galaxy 提供了一个 API 可以进行交互。不过，如果你需要一个更适合编程的流程，也有很多替代方案。在本章中，我们将探讨两个广泛使用的编程友好型流程：snakemake 和 Nextflow。在本配方中，我们考虑使用 snakemake。

Snakemake 是用 Python 实现的，并且与 Python 共享许多特性。话虽如此，它的基本灵感来源于 Makefile，这是久负盛名的 make 构建系统所使用的框架。

在这里，我们将使用 snakemake 开发一个迷你变异分析流程。这里的目标不是做出正确的科学部分—我们在其他章节中会讨论这个—而是展示如何使用 snakemake 创建流程。我们的迷你流程将下载 HapMap 数据，对其进行 1% 的子采样，做一个简单的 PCA，并绘制图形。

准备工作

你需要在安装了 snakemake 的同时安装 Plink 2。为了展示执行策略，你还需要 Graphviz 来绘制执行过程。我们将定义以下任务：

下载数据
解压它
对其进行 1% 的子采样
在 1% 子样本上计算 PCA
绘制 PCA 图

我们的流程配方将包含两个部分：用 snakemake 编写的实际流程和一个 Python 支持脚本。

这段 snakemake 代码可以在 Chapter09/snakemake/Snakefile 中找到，而 Python 支持脚本则在 Chapter09/snakemake/plot_pca.py 中。

如何执行……

第一个任务是下载数据：

from snakemake.remote.HTTP import RemoteProvider as HTTPRemoteProvider 
HTTP = HTTPRemoteProvider()
download_root = "https://ftp.ncbi.nlm.nih.gov/hapmap/genotypes/hapmap3_r3"
remote_hapmap_map = f"{download_root}/plink_format/hapmap3_r3_b36_fwd.consensus.qc.poly.map.gz"
remote_hapmap_ped = f"{download_root}/plink_format/hapmap3_r3_b36_fwd.consensus.qc.poly.ped.gz"
remote_hapmap_rel = f"{download_root}/relationships_w_pops_041510.txt"

rule plink_download:
    input:
        map=HTTP.remote(remote_hapmap_map, keep_local=True),
        ped=HTTP.remote(remote_hapmap_ped, keep_local=True),
        rel=HTTP.remote(remote_hapmap_rel, keep_local=True)

    output:
        map="scratch/hapmap.map.gz",
        ped="scratch/hapmap.ped.gz",
        rel="data/relationships.txt"

    shell:
        "mv {input.map} {output.map};"
        "mv {input.ped} {output.ped};"
        "mv {input.rel} {output.rel}"

Snakemake 的语言依赖于 Python，正如你从最前面的几行代码可以看到的那样，这些代码从 Python 的角度来看应该非常容易理解。核心部分是规则。它有一组输入流，在我们的案例中通过 HTTP.remote 进行渲染，因为我们处理的是远程文件，接着是输出。我们将两个文件放在 scratch 目录（那些尚未解压的文件）中，另一个文件放在 data 目录。最后，我们的流程代码是一个简单的 shell 脚本，它将下载的 HTTP 文件移到最终位置。注意 shell 脚本如何引用输入和输出。

使用这个脚本，下载文件变得非常简单。只需在命令行中运行以下命令：
```
snakemake -c1 data/relationships.txt
```

这告诉 snakemake 你希望生成 data/relationships.txt 文件。我们将使用单核 -c1。由于这是 plink_download 规则的输出，规则将会被执行（除非该文件已经存在—如果是这样，snakemake 就什么都不做）。以下是简化版的输出：

Building DAG of jobs...
Using shell: /usr/bin/bash
Provided cores: 1 (use --cores to define parallelism)
Rules claiming more threads will be scaled down.
Job stats:
job               count    min threads    max threads
--------------  -------  -------------  -------------
plink_download        1              1              1
total                 1              1              1

Select jobs to execute...

[Mon Jun 13 18:54:26 2022]
rule plink_download:
    input: ftp.ncbi.nlm.nih.gov/hapmap/ge [...]
    output: [..], data/relationships.txt
    jobid: 0
    reason: Missing output files: data/relationships.txt
    resources: tmpdir=/tmp

Downloading from remote: [...]relationships_w_pops_041510.txt
Finished download.
[...]
Finished job 0.
1 of 1 steps (100%) done

Snakemake 会给你一些关于将要执行的任务的信息，并开始运行这些任务。

现在我们有了数据，接下来看看解压文件的规则：

PLINKEXTS = ['ped', 'map']
rule uncompress_plink:
    input:
        "scratch/hapmap.{plinkext}.gz"
    output:
        "data/hapmap.{plinkext}"
    shell:
        "gzip -dc {input} > {output}"

这里最有趣的特点是我们可以指定多个文件进行下载。注意PLINKEXTS列表是如何在代码中转换成离散的plinkext元素的。你可以通过请求规则的输出执行它。

现在，让我们将数据下采样到 1%：

rule subsample_1p:
    input:
        "data/hapmap.ped",
        "data/hapmap.map"

    output:
        "data/hapmap1.ped",
        "data/hapmap1.map"

    run:
        shell(f"plink2 --pedmap {input[0][:-4]} --out {output[0][:-4]} --thin 0.01 --geno 0.1 --export ped")

新的内容在最后两行：我们没有使用script，而是使用run。这告诉snakemake执行是基于 Python 的，并且有一些额外的功能可用。在这里，我们看到的是 shell 函数，它执行一个 shell 脚本。这个字符串是一个 Python 的f-字符串——注意字符串中对snakemake的input和output变量的引用。你可以在这里放入更复杂的 Python 代码——例如，你可以遍历输入数据。

提示

在这里，我们假设 Plink 是可用的，因为我们已经预安装了它，但snakemake确实提供了一些功能来处理依赖关系。更具体地说，snakemake规则可以通过一个指向conda依赖的YAML文件进行注释。

现在我们已经对数据进行了下采样，让我们计算 PCA。在这种情况下，我们将使用 Plink 的内部 PCA 框架来进行计算：

rule plink_pca:
    input:
        "data/hapmap1.ped",
        "data/hapmap1.map"
    output:
        "data/hapmap1.eigenvec",
        "data/hapmap1.eigenval"
    shell:
        "plink2 --pca --file data/hapmap1 -out data/hapmap1"

和大多数管道系统一样，snakemake构建了一个snakemake来为你展示执行的 DAG，以生成你的请求。例如，为了生成 PCA，使用以下命令：
```
snakemake --dag data/hapmap1.eigenvec | dot -Tsvg > bio.svg
```

这将生成以下图像：

图 9.3 - 计算 PCA 的 DAG

最后，让我们为 PCA 生成plot规则：

rule plot_pca:
    input:
        "data/hapmap1.eigenvec",
        "data/hapmap1.eigenval"

    output:
        "pca.png"

    script:
        "./plot_pca.py"

plot规则引入了一种新的执行类型，script。在这种情况下，调用了一个外部的 Python 脚本来处理规则。

我们用来生成图表的 Python 脚本如下：

import pandas as pd

eigen_fname = snakemake.input[0] if snakemake.input[0].endswith('eigenvec') else snakemake.input[1]
pca_df = pd.read_csv(eigen_fname, sep='\t') 
ax = pca_df.plot.scatter(x=2, y=3, figsize=(16, 9))
ax.figure.savefig(snakemake.output[0])

Python 脚本可以访问snakemake对象。这个对象暴露了规则的内容：注意我们如何使用input来获取 PCA 数据，使用output来生成图像。

最后，生成粗略图表的代码如下：

图 9.4 - Snakemake 管道生成的非常粗略的 PCA

还有更多

上面的食谱是为了在简单配置的snakemake上运行制作的。snakemake中还有许多其他构建规则的方式。

我们没有讨论的最重要问题是，snakemake可以在多种不同环境中执行代码，从本地计算机（如我们的案例）、本地集群，到云端。要求除了使用本地计算机运行snakemake外更多的功能是不合理的，但不要忘了snakemake可以管理复杂的计算环境。

记住，虽然snakemake是用 Python 实现的，但其概念上基于make。是否喜欢这种（蛇形）设计是一个主观的分析。如果想了解一种替代的设计方法，查看下一个食谱，它使用 Nextflow。

使用 Nextflow 部署变异分析管道

在生物信息学的管道框架领域，有两个主要的参与者：snakemake和 Nextflow。它们提供管道功能，但设计方法不同。Snakemake 基于 Python，但它的语言和理念来源于用于编译有依赖关系的复杂程序的make工具。Nextflow 是基于 Java 的（更准确地说，它是用 Groovy 实现的——一种运行在 Java 虚拟机上的语言），并且有自己的snakemake，可以选择适合你需求的那个。

提示

评估管道系统有许多视角。这里，我们呈现了一个基于用来指定管道的语言的视角。然而，选择管道系统时，你还应该考虑其他方面。例如，它是否能很好地支持你的执行环境（如本地集群或云环境），是否支持你的工具（或者是否允许轻松开发扩展以应对新工具），以及它是否提供良好的恢复和监控功能？

在这里，我们将开发一个 Nextflow 管道，提供与我们使用snakemake实现的相同功能，从而允许从管道设计角度进行公平比较。这里的目标不是搞清楚科学部分——我们将在其他章节中讨论这个——而是展示如何使用snakemake创建管道。我们的迷你管道将下载 HapMap 数据，进行 1%的子抽样，做一个简单的 PCA，并绘制出来。

做好准备

你需要安装 Plink 2 以及 Nextflow。Nextflow 本身需要一些来自 Java 领域的软件：特别是 Java 运行时环境和 Groovy。

我们将定义以下任务：

下载数据
解压它
对其进行 1%的子抽样
对 1%子抽样数据进行 PCA 计算
绘制 PCA 图

相关的 Nextflow 代码可以在Chapter09/nextflow/pipeline.nf找到。

如何做到这一点…

第一个任务是下载数据：

nextflow.enable.dsl=2
download_root = "https://ftp.ncbi.nlm.nih.gov/hapmap/genotypes/hapmap3_r3"
 process plink_download {
  output:
  path 'hapmap.map.gz'
  path 'hapmap.ped.gz'
  script:
  """
  wget $download_root/plink_format/hapmap3_r3_b36_fwd.consensus.qc.poly.map.gz -O hapmap.map.gz
  wget $download_root/plink_format/hapmap3_r3_b36_fwd.consensus.qc.poly.ped.gz -O hapmap.ped.gz
   """
}

请记住，管道的基础语言不是 Python，而是 Groovy，因此语法会有些不同，比如使用大括号表示代码块或忽略缩进。

我们创建了一个名为plink_download的过程（Nextflow 中的管道构建块），用于下载 Plink 文件。它只指定了输出。第一个输出将是hapmap.map.gz文件，第二个输出将是hapmap.ped.gz。该过程将有两个输出通道（Nextflow 中的另一个概念，类似于流），可以被另一个过程消费。

该过程的代码默认是一个 bash 脚本。值得注意的是，脚本输出的文件名与输出部分是同步的。同时，也要注意我们是如何引用管道中定义的变量的（例如我们的例子中的download_root）。

现在我们来定义一个过程，使用 HapMap 文件并解压它们：

process uncompress_plink {
  publishDir 'data', glob: '*', mode: 'copy'

  input:
  path mapgz
  path pedgz

  output:
  path 'hapmap.map'
  path 'hapmap.ped'

  script:
  """
  gzip -dc $mapgz > hapmap.map
  gzip -dc $pedgz > hapmap.ped
  """
}

在这个过程中有三个需要注意的问题：我们现在有了一些输入（记住我们从上一个过程中得到了一些输出）。我们的脚本现在也引用了输入变量（$mapgz和$pedgz）。最后，我们通过使用publishDir发布输出。因此，任何未发布的文件将只会被暂时存储。

让我们指定一个下载并解压文件的工作流的第一个版本：
```
workflow {
    plink_download | uncompress_plink
}
```
我们可以通过在命令行中运行以下命令来执行工作流：
```
nextflow run pipeline.nf -resume
```

最后的resume标志将确保管道从已完成的步骤继续执行。步骤会在本地执行时存储在work目录中。

如果我们删除work目录，我们不希望下载 HapMap 文件，如果它们已经被发布。由于这些文件不在work目录中，因此不直接跟踪，我们需要修改工作流，以便在已发布的目录中跟踪数据：
```
workflow {
    ped_file = file('data/hapmap.ped')
    map_file = file('data/hapmap.map')
    if (!ped_file.exists() | !map_file.exists()) {
        plink_download | uncompress_plink
    }
}
```

还有其他方法可以做到这一点，但我想介绍一些 Groovy 代码，因为你可能有时需要在 Groovy 中编写代码。正如你很快会看到的，也有方法可以使用 Python 代码。

现在，我们需要对数据进行亚采样：

process subsample_1p {
  input:
  path 'hapmap.map'
  path 'hapmap.ped'

  output:
  path 'hapmap1.map'
  path 'hapmap1.ped'

  script:
  """
  plink2 --pedmap hapmap --out hapmap1 --thin 0.01 --geno 0.1 --export ped
  """
}

现在，让我们使用 Plink 计算 PCA：

process plink_pca {
  input:
  path 'hapmap.map'
  path 'hapmap.ped'
  output:
  path 'hapmap.eigenvec'
  path 'hapmap.eigenval'
   script:
  """
  plink2 --pca --pedmap hapmap -out hapmap
  """
}

最后，让我们绘制 PCA 图：

process plot_pca {
  publishDir '.', glob: '*', mode: 'copy'

  input:
  path 'hapmap.eigenvec'
  path 'hapmap.eigenval'

  output:
  path 'pca.png'

  script:
  """
  #!/usr/bin/env python
  import pandas as pd

  pca_df = pd.read_csv('hapmap.eigenvec', sep='\t') 
  ax = pca_df.plot.scatter(x=2, y=3, figsize=(16, 9))
  ax.figure.savefig('pca.png')
  """
}

该代码的新特性是我们使用 shebang（#!）操作符指定了 bash 脚本，这使我们能够调用外部脚本语言来处理数据。

这是我们的最终工作流：

workflow {
    ped_file = file('data/hapmap.ped')
    map_file = file('data/hapmap.map')
    if (!ped_file.exists() | !map_file.exists()) {
        plink_download | uncompress_plink | subsample_1p | plink_pca | plot_pca
    }
    else {
        subsample_1p(
            Channel.fromPath('data/hapmap.map'),
            Channel.fromPath('data/hapmap.ped')) | plink_pca | plot_pca
    }
}

我们要么下载数据，要么使用已经下载的数据。

虽然有其他方言用于设计完整的工作流，但我想让你注意我们如何在文件可用时使用subsample_1p；我们可以显式地将两个通道传递给一个过程。

我们可以运行管道并请求执行 HTML 报告：

nextflow run pipeline.nf -with-report report/report.xhtml

报告非常详细，能够帮助你从不同角度了解管道中哪些部分是开销较大的，无论是与时间、内存、CPU 消耗还是 I/O 相关。

还有更多内容

这是一个简单的 Nextflow 入门示例，希望能让你对这个框架有所了解，尤其是让你能够与snakemake进行比较。Nextflow 有更多的功能，建议你查看它的网站。

与snakemake一样，我们没有讨论的最重要问题是 Nextflow 可以在许多不同的环境中执行代码，从本地计算机、现场集群到云。请查阅 Nextflow 的文档，了解当前支持哪些计算环境。

无论底层语言多么重要，Groovy 与 Nextflow，Python 与snakemake一样，都必须比较其他因素。这不仅包括两个管道可以在哪些地方执行（本地、集群或云），还包括框架的设计，因为它们使用的是截然不同的范式。