计算分子生物学的 Python 工具
计算分子生物学的 Python 工具
BioGeography 是一个正在开发中的模块，由 Nick Matzke 为 Google Summer of Code 2009 项目开发。它通过 NESCENT 的 Phyloinformatics Summer of Code 2009 运行。请参阅项目提案：用于 BioPython 的生物地理系统发育。指导老师是 Stephen Smith（主要）、Brad Chapman 和 David Kidd。源代码位于 GitHub 上 nmatzke 分支的 Geography 分支 的 Bio/Geography 目录中，您可以在 此处 查看模块正在进行的开发的时间线和其他信息。新模块正在 BioPython wiki 上作为 BioGeography 进行文档记录。
摘要：创建一个 BioPython 模块，使用户能够自动访问和解析来自在线生物多样性数据库的物种地点记录；将这些记录链接到用户指定的系统发育树；计算基本的 alpha 和 beta 系统发育多样性汇总统计数据；生成用于输入各种推断算法的输入文件，这些算法可用于推断历史生物地理；将这些程序的输出转换为适合映射的文件，例如 Google Earth（KML 文件）。
所有类和函数都使用标准文档字符串进行了文档记录。代码可在最新的 GitHub 提交中找到，地址：http://github.com/nmatzke/biopython/commits/Geography
Bio.Geography 是一个用于收集和处理生物地理数据的模块。该模块的主要动机是协助进化生物地理学的分析。有各种推断算法可用于此类分析，例如 DIVA 和 lagrange。此类程序的输入通常是 (a) 系统发育树和 (b) 系统发育树顶端物种所居住的区域。收集了特定群体数据的研究人员可能可以直接访问物种位置数据，但许多大型分析可能需要收集大量出现数据。自动收集/处理出现数据还有其他各种应用，包括物种映射、生态位建模、博物馆记录的错误检查以及监测范围变化。
出现数据主要来自博物馆收藏。此类数据的主要来源是 全球生物多样性信息设施 (GBIF)。GBIF 提供了全球数百个博物馆记录的出现数据。GBIF 出现数据可以 手动搜索，并以各种格式下载结果（请参阅 GBIF 网站上的示例）：电子表格、Google Earth KML 或 XML DarwinCore 格式。
GBIF 也可以通过 API 访问。Bio.Geography 可以处理手动下载的 DarwinCore 结果，或直接访问 GBIF。
了解 Bio.Geography 中类的总体组织很有用。GbifXml 模块中有四个类
对于 GBIF 的一次性使用，您可能会发现只需以电子表格格式（用于分析）或 KML（用于映射）下载出现数据最容易。但是，对于许多群体的分析，或者随着 GBIF 中添加新数据而重复更新分析，自动化是可取的。
在 GBIF 网站上进行的手动搜索可以以符合 DarwinCore 数据标准的 XML 文件形式返回结果。示例文件可以在 Biopython 的 Tests/Geography 目录中找到，名为utric_search_v2.xml。该文件包含 1000 多个关于Utricularia（一种食肉植物属）的出现记录。
将 utric_search_v2.xml 文件保存在您的工作目录中（或从 GBIF 下载类似文件）。以下是在使用 Bio.Geography 的 GbifXml 模块解析文件的建议步骤。首先，导入必要的类和函数，并指定输入文件的名称。
from Bio.Geography.GbifXml import GbifXmlTree, GbifSearchResults
from Bio.Geography.GeneralUtils import fix_ASCII_file
xml_fn = "utric_search_v2.xml"
其次，为了在 python 中将结果显示到屏幕上，我们需要将文件转换为纯 ASCII（GBIF 结果包含来自不同语言的所有许多不寻常字符，以及斜线引号等的规范化；这会导致尝试在 python 或 ipython 中打印到屏幕时崩溃）。
xml_fn_new = fix_ASCII_file(xml_fn)
这会创建一个新文件，文件名中添加了字符串“_fixed.xml”。
接下来，我们将 XML 文件解析为 ElementTree（一个 Python 对象，它将 XML 文件中的数据作为嵌套的列表和字典系列包含在内）。
from xml.etree import ElementTree as ET
xmltree = ET.parse(xml_fn_new)
然后，我们可以将元素树存储为类 GbifXmlTree 的对象
gbif_recs_xmltree = GbifXmlTree(xmltree)
然后，在存储了 xmltree 之后，我们将它解析为单个记录（存储在类 GbifObservationRecord 的单个对象中），然后将它们作为一组存储在类 GbifSearchResults 的对象中。
recs = GbifSearchResults(gbif_recs_xmltree)
recs.extract_occurrences_from_gbif_xmltree(recs.gbif_recs_xmltree)
可以在 recs.obs_recs_list 中访问单个观察记录的列表。这将显示对前五个记录的引用
recs.obs_recs_list[0:4]
获取第一条记录的数据
rec = recs.obs_recs_list[0]
dir(rec)
rec.lat 将返回纬度，rec.long 返回经度等。某些数据属性在所有 GBIF 记录中都找不到；如果它们丢失，则该字段将包含“None”。
要以制表符分隔的表格格式打印所有记录
recs.print_records()
在我们花时间下载数千条记录之前，我们可能想先知道 GBIF 中有多少条记录。用户必须设置一个字典，其中包含字段和搜索词作为键和项。即，
from GbifXml import GbifXmlTree, GbifSearchResults
params = {"format": "darwin", "scientificname": "Genlisea*"}
“‘format’: ‘darwin’” 指定 GBIF 应以 DarwinCore 格式返回结果。
‘scientificname’ 指定要搜索的属名。在名称后添加一个 ‘*’ 将返回以“Genlisea”开头的任何内容。
完整的搜索词列表可以在 GBIF 的 出现记录数据服务 上找到，该服务链接自 使用 GBIF 门户网站中的数据。
指定搜索参数后，启动一个新的 GbifSearchResults 对象并运行 get_numhits 获取命中次数
params = {"format": "darwin", "scientificname": "Genlisea*"}
recs = GbifSearchResults()
numhits = recs.get_numhits(params)
截至 2009 年 8 月，GBIF 中存在 169 条与“Genlisea*”匹配的记录
为了对比，不使用星号 (‘*’) 运行相同的搜索
params = {"format": "darwin", "scientificname": "Genlisea"}
numhits = recs.get_numhits(params)
我们只得到了大约 10 个结果——可能是只识别到属而不是更深层次的标本记录。
可以通过键下载单个记录。要下载单个记录
rec = recs.obs_recs_list[0]
key = rec.gbifkey
# (or manually)
# key = 175067484
xmlrec = recs.get_record(key)
print(xmlrec)
如果要打印 xmlrec ElementTree 对象，请将 xmlrec 存储在 GbifXmlTree 对象中并运行 print_xmltree
GbifXmlTree(xmlrec).print_xmltree()
生物地理区域通常以 alpha 和 beta 多样性统计数据为特征：基本上，这些是区域内或区域之间物种数量的指标。给定一个区域中生物体的系统发育树，可以计算系统发育 alpha 和 beta 多样性统计数据。这在 Webb 等人的 phylocom 包 中得到了充分的实现，但出于某些目的，直接在 python 中计算统计数据很有用。
这里，我们需要从 Newick 树字符串开始
trstr2 = "(((t9:0.385832, (t8:0.445135,t4:0.41401)C:0.024032)B:0.041436, t6:0.392496)A:0.0291131, t2:0.497673, ((t0:0.301171, t7:0.482152)E:0.0268148, ((t5:0.0984167,t3:0.488578)G:0.0349662, t1:0.130208)F:0.0318288)D:0.0273876);"
to2 = Tree(trstr2)
然后，我们创建一个树汇总对象
ts = TreeSum(to2)
函数 test_Tree 将运行指标（MPD = 平均系统发育距离、NRI = 净相关性指数、MNPD = 平均最近邻居系统发育距离、NTI = 最近分类单元指数、PD = 总系统发育距离）并将输出到屏幕
ts.test_Tree()
通过将树子集化到仅存在于一个区域中的分类单元，可以按区域计算统计数据。
GBIF 每次只允许下载最多 1000 条观察记录（KML 记录为 10,000 条）。要获得更多记录，我们需要分阶段下载和处理它们。
我们再次设置参数字典，以及一个“inc”变量来指定每次服务器请求下载的记录数量。
params = {"format": "darwin", "scientificname": "Genlisea*"}
inc = 100
recs3 = GbifSearchResults()
gbif_xmltree_list = recs3.get_all_records_by_increment(params, inc)
与 biopython 与 NCBI 服务器的交互一样，GbifSearchResults 模块会跟踪上次 GBIF 请求的时间，并且在进行新的请求之前需要等待 3 秒。
每次服务器请求都会返回一个 XML 字符串；这些字符串被解析为 GbifXmlTree 对象，并且返回的 GbifXmlTree 对象列表将返回到 gbif_xmltree_list。单个记录也已被解析
recs3.print_records()
生物地理学分析通常需要您确定分类单元生活在哪些区域。区域并不总是明显地划分，分析师可能希望尝试几种不同的可能区域集，并观察这如何影响他们的分析。
在下面，我们设置了一个包含北半球纬度/经度坐标的多边形，然后将每个匹配记录的“area”属性设置为“NorthernHemisphere”
ul = (-180, 90)
ur = (180, 90)
ll = (-180, 0)
lr = (180, 0)
poly = [ul, ur, ll, lr]
polyname = "NorthernHemisphere"
recs3.print_records()
此过程可以对所有感兴趣的多边形重复，直到所有 GBIF 记录都被分类（除了 GBIF 记录本身就缺乏纬度/经度数据，这种情况有时会发生）。
GeogUtils 还包含用于处理 shapefile/dbf 文件的开放访问库——这些是标准的 GIS 文件格式，各种公开可用的 shapefile 可能用作多边形的来源。
警告：如果您的多边形跨越国际日期变更线，则点在多边形内的操作将彻底失败。在这种情况下，最好的解决方案是将跨越日期变更线的任何多边形拆分为两个多边形，一个位于线的每一侧。
GBIF 搜索结果通常包含非 ASCII 字符（例如国际地名）和其他令人困惑的项目，例如尖括号中的网络链接，如果将 GBIF 搜索结果读入 ASCII，然后尝试对其进行解析，则可能会被误解为不匹配的 XML 标记。
一般来说，如果结果在后台处理，Geography 模块会很好地处理它们；但是要将结果打印到屏幕上，会使用 GeneralUtils 中的一系列函数将字符串转换为纯 ASCII。这样可以避免在将数据打印到屏幕时发生崩溃。因此，这些打印到屏幕的结果可能会略微改变原始搜索结果的内容。