Python 生物信息学秘籍（三）
在这里，我们将介绍 Biopython 的PDB模块，用于处理 PDB。我们将使用代表部分p53蛋白质的三个模型。你可以在 http://www.rcsb.org/pdb/101/motm.do?momID=31和阅读更多关于和p53的文件。银河(galaxyproject.org/tutorials/g101/)是一个开源系统，让非计算用户能够做计算生物学。这是最广泛使用的，用户友好的管道系统。
七、系统发育学
系统发育学是分子测序的应用，用于研究生物之间的进化关系。说明这一过程的典型方法是使用系统进化树。从基因组数据计算这些树是一个活跃的研究领域，具有许多现实世界的应用。
在本书中，我们将把提到的实用方法提升到一个新的水平:这里的大多数食谱都受到了最近一项关于埃博拉病毒的研究的启发，该研究研究了最近在非洲爆发的埃博拉病毒。这项研究名为基因组监测，阐明了 2014 年埃博拉病毒爆发期间的起源和传播，作者 Gire 等人，发表在科学上。在pubmed.ncbi.nlm.nih.gov/25214632/有售。在这里，我们将尝试遵循类似的方法，以达到类似的论文结果。
在这一章中，我们将使用 DendroPy(一个系统发育库)和 Biopython。Docker 映像包括所有必要的软件。
在本章中，我们将介绍以下配方:
- 为系统发育分析准备数据集
- 对齐遗传和基因组数据
- 比较序列
- 重建系统进化树
- 递归地玩树
- 可视化系统发育数据
准备用于系统发育分析的数据集
在这个菜谱中，我们将下载并准备用于我们分析的数据集。该数据集包含埃博拉病毒的完整基因组。我们将使用 DendroPy 来下载和准备数据。
准备就绪
我们将从 GenBank 下载完整的基因组；这些基因组是从各种埃博拉疫情中收集的，包括 2014 年爆发的几次疫情。请注意，有几种病毒会导致埃博拉病毒病；2014 年疫情涉及的物种(EBOV 病毒，正式名称为扎伊尔埃博拉病毒)是最常见的，但这种疾病是由埃博拉病毒属的更多物种引起的；另外四种也以序列形式提供。你可以在en.wikipedia.org/wiki/Ebolavirus了解更多。
如果您已经阅读了前面的章节，您可能会对这里涉及的潜在数据大小感到恐慌；这根本不是问题，因为这些病毒的基因组大小都在 19 kbp 左右。所以，我们大约 100 个基因组实际上很轻。
为了进行这个分析，我们需要安装dendropy。如果您正在使用 Anaconda，请执行以下操作:
conda install –c bioconda dendropy
像往常一样，这个信息可以在相应的 Jupyter 笔记本文件中找到，该文件可以在Chapter07/Exploration.py找到。
怎么做…
看看下面的步骤:
-
首先，让我们从使用树复制指定我们的数据源开始，如下:
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 文件，这些文件要么包括所有基因组，要么只包括 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树类，并计算分类群的数量。然后，我们提取每个序列的所有氨基酸(我们排除了由-确定的缺口)来计算长度并报告最小、平均和最大大小。关于 API 的更多细节，请看一下树图文档。
-
让我们检查 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 森林病毒=TAFV，苏丹病毒=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，在 http://etetoolkit.org/ ?? 有售。
参见
- 美国疾病控制中心在 https://www.cdc.gov/vhf/ebola/history/summaries.xhtml 有一个很好的关于埃博拉病毒疾病的介绍页。
- 系统发育学中的参考应用是乔·费尔森斯坦的 Phylip ，可以在 http://evolution.genetics.washington.edu/phylip.xhtml 的找到。
- 我们将在未来的食谱中使用 Nexus 和 Newick 格式(
evolution . genetics . Washington . edu/phylip/new icktree . XHTML)，但也要查看 PhyloXML 格式(【http://en.wikipedia.org/wiki/PhyloXML】)。
比对遗传和基因组数据
在我们能够进行任何系统发育分析之前，我们需要比对我们的遗传和基因组数据。这里，我们将使用 maft(【http://mafft.cbrc.jp/alignment/software/】)对进行基因组分析。基因分析将使用肌肉进行(www.drive5.com/muscle/)。
准备就绪
要执行基因组比对，您需要安装 MAFFT。此外，为了进行基因比对，将使用肌肉。此外，我们将使用 trimAl(trimal.cgenomics.org/)以自动方式去除虚假序列和排列不良的区域。所有包装均可从 Bioconda 获得:
conda install –c bioconda mafft trimal muscle=3.8
通常，这些信息可以在Chapter07/Alignment.py的相应 Jupyter 笔记本文件中找到。您需要事先运行上一个笔记本，因为它将生成此处需要的文件。在本章中，我们将使用 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)、田岛的 D (taj_d)和沃特森的θ(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 DataFrame 来打印结果，因为它真的是为处理这样的操作而定制的。我们将用一个 NumPy 多维数组初始化我们的数据框架，该数组有四行(基因)和四个统计数据乘以两个物种。
统计数据，如分离位点的数量、核苷酸多样性、Tajima’s D 和 Watterson’s theta，由 DendroPy 计算。注意数组中单个数据点的位置(坐标计算)。
看最后一行:如果你在 Jupyter 中，只需在末尾加上df就能显示数据帧和单元格输出。如果你不在笔记本里，就用print(df)(你也可以在笔记本里执行这个，但是看起来不会那么好看)。
-
现在，让我们提取类似的信息，但是是全基因组的，而不仅仅是全基因的。在这种情况下，我们将使用两次 EBOV 暴发(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 年 EBOV 爆发的数据集很大，我们对其进行了二次抽样，只有 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，一个基于最大似然推理大型系统树的程序，你可以在 http://sco.h-its.org/exelixis/software.xhtml查看。Bioconda 也包括它，但它被命名为raxml。注意，这个二进制文件叫做raxmlHPC。您可以执行以下命令来安装它:
conda install –c bioconda raxml
前面的代码很简单，但是执行起来需要时间，因为它将调用 RAxML(这是计算密集型的)。如果您选择使用树型接口，它也可能会占用大量内存。我们将与 RAxML、DendroPy 和 Biopython 交互，让您选择使用哪个接口；DendroPy 为您提供了一种访问结果的简单方法，而 Biopython 对内存的要求较低。尽管在本章的后面有可视化的方法，我们还是会在这里绘制一个生成的树。
通常，这些信息可在Chapter07/Reconstruction.py的相应笔记本中找到。您将需要上一个食谱的输出来完成这个。
怎么做…
看看下面的步骤:
-
对于 DendroPy，我们将首先加载数据，然后重建 genus 数据集，如下所示:
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 错误。
我们将使用 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()
与 Biopython 相比，DendroPy 有一个更加声明性的接口，所以你可以处理一些额外的事情。您应该指定种子(如果您不这样做，Biopython 将设置一个固定的默认值 10，000)和工作目录。对于 RAxML，工作目录规范需要绝对路径。
-
让我们来看看 Biopython 运行的结果。虽然对于 DendroPy 和 Biopython 来说，RAxML 输出是相同的(除了随机性)，但是 DendroPy 抽象了一些东西。使用 Biopython，您需要自己处理结果。您也可以使用树状视图来执行此操作；但是，在这种情况下，它是可选的:
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/clad istics 5 . 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 的节点表示有一个 level 方法(用于比较)，但这里的重点是引入一个递归算法，所以无论如何我们都会实现。
注意函数是如何工作的；用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加上它整个后代的最大高度。
注意，我们使用 lambda 函数上的 map 来获取当前节点的所有子节点的高度。然后，我们选择最大值(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 的情况具有相同的物种和爆发年份，我们将删除所有子节点、叶子和内部子树节点。
如果它们不属于同一个物种，我们会递归下去，直到发生这种情况。最坏的情况是，当你已经在一个叶节点时，该算法简单地解析为当前节点的种类。
还有更多…
有大量关于树和数据结构的计算机科学文献；如果你想了解更多，维基百科在 http://en.wikipedia.org/wiki/Tree_%28data_structure%29 提供了很棒的介绍。
注意，作为 Python 方言，不鼓励使用lambda函数和map；你可以在 http://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 爆发年份)获得一种颜色。如果它是一个内部节点，并且它下面的所有后代节点都是同一物种，那么它将获得该物种的颜色；如果在那之后有几个物种，它将被涂成灰色。其实颜色功能是可以改的，以后也会改。仅使用边缘颜色(标签将以黑色打印)。
请注意，就清晰的视觉外观而言，梯形化(在前面的树型方法中执行)非常有帮助。
我们也深抄属树来着色一份拷贝；记住前面的方法，一些树遍历函数可以改变状态，在这种情况下，我们希望保留一个没有任何着色的版本。
注意使用 lambda 函数来清除 trimAl 更改的名称，如下图所示:
图 7.6-带有完整埃博拉病毒数据集的分层彩色系统进化树
还有更多…
树和图形可视化是一个复杂的主题；可以说，在这里，树的可视化是严格的，但远不美观。DendroPy 的另一个选择是 ETE(etetoolkit.org/)，它有更多的可视化特性。绘制树和图的一般备选方案有胞景(cytoscape.org/)和格菲(gephi.github.io/)。如果你想了解更多关于渲染树和图的算法，可以查看 http://en.wikipedia.org/wiki/Graph_drawing 的维基百科页面，了解这个有趣话题的介绍。
但是，注意不要用形式来换取实质。例如，这本书的前一版使用图形渲染库对系统进化树进行了漂亮的渲染。虽然这显然是那一章中最美丽的图像，但它在分支长度方面有误导性。
八、使用蛋白质数据库
蛋白质组学是对蛋白质的研究，包括它们的功能和结构。该领域的主要目标之一是表征蛋白质的三维结构。蛋白质组学领域最广为人知的计算资源之一是蛋白质数据库 ( PDB )，这是一个拥有大生物分子结构数据的存储库。当然，许多数据库反而关注蛋白质一级结构；这些有点类似于我们在 第二章 、了解 NumPy、pandas、Arrow 和 Matplotlib 中看到的基因组数据库。
在本章中，我们将主要关注处理来自 PDB 的数据。我们将看看如何解析 PDB 文件，执行一些几何计算，并可视化分子。我们将使用旧的 PDB 文件格式，因为从概念上讲，它允许您在稳定的环境中执行大多数必要的操作。话虽如此，预定取代 PDB 格式的较新的 mmCIF 也将出现在用 Biopython 方法解析 mmCIF 文件的中。我们将使用 Biopython 并引入 PyMOL 进行可视化。我们不会在这里讨论分子对接，因为这可能更适合一本关于化学信息学的书。
在本章中，我们将使用一个蛋白质的经典例子:肿瘤蛋白 p53，一种参与细胞周期调节(例如，凋亡)的蛋白质。这种蛋白质与癌症高度相关。网上有大量关于这种蛋白质的信息。
让我们从你现在应该更熟悉的东西开始:访问数据库，特别是蛋白质的一级结构(如氨基酸序列)。
在本章中，我们将介绍以下配方:
- 在多个数据库中查找蛋白质
- 介绍生物。物理数据库
- 从 PDB 文件中提取更多信息
- 在 PDB 文件上计算分子距离
- 执行几何运算
- 使用 PyMOL 制作动画
- 用 Biopython 解析 mmCIF 文件
在多个数据库中寻找一种蛋白质
在我们开始执行更多的结构生物学之前，我们将看看如何访问现有的蛋白质组数据库，如 UniProt。我们将在 UniProt 中查询我们感兴趣的基因 TP53 ，并从那里获取它。
准备就绪
为了访问数据，我们将使用 Biopython 和 REST API(我们在第五章 、使用基因组中使用了类似的方法)和requests库来访问 web APIs。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
-
深入查看前面的记录会发现很多非常有趣的信息，特别是关于特征、基因本体 ( GO )和数据库
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 条目，以及三个例子。目前，仅这种蛋白质就有 1509 个特征、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')
还有更多
还有更多关于蛋白质信息的数据库——其中一些在前面的记录中提到过。您可以探索其结果，尝试在其他地方找到数据。关于 UniProt 的 REST 接口的详细信息，请参考www.uniprot.org/help/programmatic_access。
介绍简历。物理数据库
在这里，我们将介绍 Biopython 的PDB模块，用于处理 PDB。我们将使用代表部分p53蛋白质的三个模型。你可以在 http://www.rcsb.org/pdb/101/motm.do?momID=31和阅读更多关于和p53的文件。
准备就绪
您应该已经知道模型、链、剩余和原子对象的基本数据模型。在bio python . org/wiki/The _ bio python _ Structural _ bio informatics _ FAQ可以找到关于 Biopython 的结构生物信息学 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.
-
让我们在
1TUP模型中得到除水以外的所有非标准残留物(HETATM)，如下面的代码所示: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笔记本文件里找到这个内容。
怎么做…
我们将按照以下步骤开始:
-
首先我们检索一下
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笔记本文件里找到这个内容。
怎么做…
看看下面的步骤:
-
让我们加载我们的模型，如下:
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') -
我们现在将得到我们的 zincs，稍后我们将进行比较:
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 ng strm:
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]
我们只有三个 zincs，所以计算量大大减少了。然而，想象我们有更多的原子，或者我们正在对集合中的所有原子进行成对比较(记住，在成对的情况下，比较的数量随着原子数量的二次方增长)。虽然我们的案例很小，但是预测用例并不难，而更多的比较会花费很多时间。我们将很快回到这个话题。
-
让我们看看随着距离的增加，我们得到了多少个原子:
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 次，然后以毫秒为单位打印结果。我们将该函数作为一个字符串传递，然后传递另一个带有必要导入的字符串，以使该函数工作。在笔记本上，你可能意识到了%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/prot expl/help _ hyd . htm。
使用 PyMOL 制作动画
这里，我们将制作一个 p53 1TUP模型的视频。为此，我们将使用 PyMOL 可视化库。我们将通过围绕 p53 1TUP模型移动然后放大来开始我们的动画；当我们放大时，我们会改变渲染策略，以便您可以更深入地看到模型。你可以找到你将在odysee.com/@Python:8/protein_video:8制作的视频版本。
准备就绪
这个食谱将以 Python 脚本的形式呈现，而不是以笔记本的形式。这主要是因为输出不是交互式的，而是一组需要进一步后期处理的图像文件。
你需要安装 PyMOL(【http://www.pymol.org】)。在 Debian、Ubuntu 或 Linux 上，可以使用apt-get install pymol命令。如果你使用 Conda，我建议不要使用它，因为依赖性很容易解决——此外，你将安装一个需要许可证的 30 天试用版，而上面的版本是完全开源的。如果你没有使用 Debian 或 Linux，我建议你安装适用于你的操作系统的开源版本。
PyMOL 与其说是一个 Python 库，不如说是一个交互式程序，所以我强烈建议您在继续学习之前先尝试一下。这会很有趣的！这个菜谱的代码作为脚本可以在 GitHub 仓库中找到，还有本章的笔记本文件，在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 界面(带有 GUI)的情况下，才能使用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 -第 200 帧附近的 DNA 分子和场景 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 扩展)。反过来使用它(在没有 GUI 的情况下从 Python 导入)可能会很复杂和令人沮丧。PyMOL 启动一个单独的线程来呈现异步工作的图像。
例如，这意味着您的代码可能位于与呈现器不同的位置。我已经将另一个名为PyMol_Intro.py的脚本放在了 GitHub 库中；您将看到第二个 PNG 调用将在第一个调用完成之前开始。尝试脚本代码，看看您期望它如何运行，以及它实际上是如何运行的。
在www.pymolwiki.org/index.php/MovieSchool有大量从 GUI 角度看 PyMOL 的好文档。如果你想拍电影，这是一个很好的起点，而www.pymolwiki.org是一个信息宝库。
使用 Biopython 解析 mmCIF 文件
mmCIF 文件格式很可能是未来。Biopython 还没有完整的功能来使用它，但我们会看看目前存在的功能。
准备就绪
由于Bio.PDB不能自动下载 mmCIF 文件，您需要获取您的蛋白质文件并将其重命名为1tup.cif。这个可以在1TUP.cif下的github . com/packt publishing/Bioinformatics-with-Python-Cookbook-third-Edition/blob/master/datasets . py找到。
你可以在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解析器有更多的发展。
PDB 的开发者在 http://mmcif . wwpdb . org/docs/SW-examples/Python/html/index . XHTML 提供了一个 Python 库。
九、生物信息管道
管道是任何数据科学环境中的基础。数据处理从来都不是一项单一的任务。许多管道是通过专用脚本实现的。这可以用一种有用的方式来实现，但是在许多情况下，它们不符合几个基本的观点，主要是可再现性、可维护性和可扩展性。
在生物信息学中，你可以找到三种主要类型的管道系统:
- 像 Galaxy(
usegalaxy.org)这样的框架是面向用户的，也就是说，它们公开了易于使用的用户界面，隐藏了大部分底层机制。 - 程序化的工作流程——面向代码接口，虽然是通用的，但源自生物信息学空间。两个例子是 snake make(【https://snakemake.readthedocs.io/】)和 next flow(【https://www.nextflow.io/】)。
- 完全通用的工作流系统，如 Apache air flow(
airflow.incubator.apache.org/)，它们采用了一种不那么以数据为中心的工作流管理方法。
在这一章中，我们将讨论 Galaxy，它对于支持不太倾向于自己编写解决方案的用户的生物信息学家来说尤其重要。虽然您可能不是这些管道系统的典型用户，但您可能仍然需要支持它们。幸运的是，Galaxy 提供了 API，这将是我们的主要关注点。
我们还将讨论 Snakemake 和 Nextflow，它们是起源于生物信息学领域的具有编程接口的通用工作流工具。我们将涵盖这两种工具，因为它们是该领域中最常见的。我们将使用 Snakemake 和 Nextflow 解决一个类似的生物信息学问题。我们将尝试两种框架，并希望能够决定最喜欢的一种。
这些食谱的代码并不是以笔记本的形式出现，而是以 Python 脚本的形式出现在该书知识库的Chapter09目录中。
在本章中，您将找到以下配方:
- 银河服务器简介
- 使用 API 访问 Galaxy
- 用 Snakemake 开发变体分析管道
- 用 netflow 开发变体分析管道
银河服务器简介
银河(galaxyproject.org/tutorials/g101/)是一个开源系统，让非计算用户能够做计算生物学。这是最广泛使用的，用户友好的管道系统。任何用户都可以将 Galaxy 安装在一台服务器上，但是网上也有很多其他的服务器可以公开访问，旗舰服务器是 http://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
这就拉了 bjrn grüning 的惊人的 Docker 星系图像。使用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 -银河 Docker 主页
- 现在用默认的用户名和密码组合:
admin和password登录(见顶栏)。 - 从顶部菜单中选择用户，在里面选择偏好。
- 现在，选择管理 API 密钥。
不要更改 API 密钥。前面练习的目的是让您知道 API 键在哪里。在实际情况下，您必须到此屏幕获取您的密钥。只需注意 API 键:fakekey。顺便说一下，在正常情况下，这将是一个 MD5 散列。
因此，在这个阶段，我们用以下(默认)凭证安装了服务器:用户名为admin，密码为password，API 密匙为fakekey。接入点是localhost:8080。
还有更多
bjrn grüning 的形象在本章中的使用方式非常简单；毕竟这不是一本关于系统管理或者 DevOps 的书，而是一本编程的书。如果你访问 https://github.com/bgruening/docker-galaxy-stable，你会发现有无数种方法来配置镜像，并且都有很好的文档记录。我们这里的简单方法适用于我们的开发目的。
如果你不想在你的本地电脑上安装 Galaxy ，你可以使用一个公共服务器如usegalaxy.org来做下一个菜谱。这不是 100%万无一失的，因为服务会随着时间的推移而变化，但可能会非常接近。采取以下步骤:
- 在公共服务器上创建一个帐户(
usegalaxy.org或其他)。 - 按照前面的说明访问您的 API 密钥。
- 在下一个方法中，您必须替换主机、用户、密码和 API 密钥。
使用 API 访问 Galaxy
虽然 Galaxy 的主要用例是通过一个易于使用的 web 界面，但它也提供了一个用于编程访问的 REST API。有几个语言提供的接口，比如 bio blend(bio blend . readthedocs . io提供 Python 支持。
在这里，我们将开发一个脚本，将一个 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，并指定我们将创建的星系历史的名称(
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' -
最后，我们能够连接到银河服务器:
gi = GalaxyInstance(url=rest_url, key=api_key) gi.verify = False -
我们现在将列出所有可用的
histories: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)
如果你愿意，你可以在网络界面上查看，你会在那里找到新的历史记录。
-
我们现在要上传文件；这需要 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
-
让我们检查一下文件
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 访问了银河。
使用 Snakemake 部署变体分析管道
Galaxy 主要面向不太喜欢编程的用户。知道如何处理它是很重要的，因为它无处不在，即使你更喜欢对程序员友好的环境。令人欣慰的是，有一个 API 可以与 Galaxy 交互。但是如果你想要一个对程序员更友好的管道，有很多选择。在这一章中，我们探索两个广泛使用的程序员友好的管道:snakemake和 Nextflow。在这个食谱中，我们考虑snakemake。
Snakemake 是用 Python 实现的，与它有许多共同的特点。也就是说，它的基本灵感是 Makefile，古老的make构建系统使用的框架。
这里，我们将使用snakemake开发一个迷你变体分析管道。这里的目标不是让科学部分正确——我们将在其他章节讨论——而是看看如何用 ?? 创建管道。我们的迷你管道将下载 HapMap 数据，以 1%的比例对其进行二次采样，进行简单的 PCA，并绘制它。
准备就绪
您需要将 Plink 2 安装在snakemake旁边。为了显示执行策略，您还需要 Graphviz 来绘制执行。我们将定义以下任务:
- 下载数据
- 解压缩它
- 以 1%的比例进行子采样
- 计算 1%子样本的主成分分析
- 绘制 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 脚本。该字符串是一个 Pythonf-string-注意该字符串中对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构造了一个操作的有向无环图 ( DAG )来执行。在任何时候，您都可以要求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-snake make 管道产生的非常粗糙的 PCA
还有更多
前面的配方在一个简单的配置snakemake上运行。在snakemake中有更多的方法来构造规则。
我们没有讨论的最重要的问题是，snakemake可以在许多不同的环境中执行代码，从本地计算机(如我们的例子)、本地集群到云。要求使用本地计算机来尝试snakemake是不合理的，但是不要忘记snakemake可以管理复杂的计算环境。
请记住，snakemake虽然是用 Python 实现的，但在概念上是基于make的。这是一个主观分析，以决定你是否喜欢(蛇)作出的设计。对于另一种设计方法，检查下一个配方，它使用 Nextflow。
使用 Nextflow 部署变体分析管道
生物信息学中的流水线框架空间有两个主要玩家:snakemake和 Nextflow。它们提供流水线功能，同时具有不同的设计方法。Snakemake 基于 Python，但它的语言和哲学来自于用于编译具有依赖关系的复杂程序的make工具。Nextflow 是基于 Java 的(更准确地说，它是用 Groovy 实现的——一种在 Java 虚拟机上工作的语言),并且有自己的领域特定语言 ( DSL )用于实现管道。这个配方(和之前的配方)的主要目的是给你一个 Nextflow 的味道，这样你就可以和snakemake比较，选择一个更适合你需求的。
小费
关于如何评估管道系统，有许多观点。在这里，我们基于用于指定管道的语言呈现一个透视图。然而，在选择管道系统时，你还应该考虑其他因素。例如，它在多大程度上支持您的执行环境(比如本地集群或云)，它是否支持您的工具(或者允许轻松开发扩展来处理新工具)，它是否提供良好的恢复和监控功能？
在这里，我们将使用 Nextflow 开发一个管道，提供与我们使用snakemake实现的相同的功能，从而允许从管道设计的角度进行公平的比较。这里的目标不是让科学部分正确——我们在其他章节中讨论——而是看看如何用snakemake创建管道。我们的迷你管道将下载 HapMap 数据，以 1%对其进行子采样，进行简单的 PCA，并绘制它。
准备就绪
您需要将 Plink 2 与 Nextflow 一起安装。Nextflow 本身需要一些来自 Java 领域的软件:特别是 Java 运行时环境和 Groovy。
我们将定义以下任务:
- 下载数据
- 解压缩它
- 以 1%的比例进行子采样
- 计算 1%子样本的主成分分析
- 绘制 PCA 图
可以在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 } -
我们可以通过在 shell 上运行以下命令来执行工作流:
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:
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
-
最后，让我们绘制 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，一定要比较其他因素。这不仅包括两个管道可以在哪里执行(本地、集群或云中)，还包括框架的设计，因为它们使用完全不同的范例。
DAMO开发者矩阵，由阿里巴巴达摩院和中国互联网协会联合发起，致力于探讨最前沿的技术趋势与应用成果，搭建高质量的交流与分享平台，推动技术创新与产业应用链接，围绕“人工智能与新型计算”构建开放共享的开发者生态。
更多推荐
21
0- 0
已为社区贡献57条内容
所有评论(0)