快速获取功能相关基因

一个Python爬虫，快捷批量爬取与感兴趣功能有关的基因

转录组数据分析中，老师会碰到这样的问题：

（1）这么多差异表达基因中，哪些基因参与了特定的功能，我如何寻找与预期功能有关的基因？

（2）我的测序物种比较新，数据库中信息比较缺乏，难以通过已知数据寻找目标功能基因。

这个时候，如果一个个去设法查询每个基因的功能，工作量可能是非常大的。那么，有没有快捷的方法呢？答案是肯定有，比如可以转换下思路，首先根据预期的功能，在GO、KEGG等数据库中检索基因功能分类，并进而在搜索到的功能条目中细化基因。这样，就可以获得了与预期功能有关的一个基因集。随后，再根据感兴趣的差异基因名称，在该基因集中查找，如果匹配到相似的，就能够明确差异基因的功能了。

当然，在数据库中手动检索功能分类并匹配基因的过程，可能仍然略微繁琐。如果能够将这个过程交给计算机自动实现，那该多么节省时间。为此，我们提供了一个Python爬虫，帮助您在GO数据库中自动搜索并下载与其功能的相关基因。

就从这个思想出发，本文以探讨“氧化应激”相关功能的基因为例，教大家如何快速实现。

1 搜索关键GO功能条目

首先进入GO网站（http://geneontology.org/），输入想知道的功能名称。

例如我想关注和“氧化应激”相关的功能，输入“Oxidative stress”查询。

结果中都是以“Oxidative stress”为关键词找到的GO条目，它们都是涉及“氧化应激”相关的功能。接下来将这些GO条目下载下来，点击右上方“Custom”可获取GO条目列表。下载时，可以视情况过滤一些不重要的，选择要下载的打勾，如果不手动选择将默认下载全部搜索结果。

该列表中的GO名称可手动粘贴到一个文本文档中，例如命名为“Oxidative_stress.txt”。

2 爬取GO条目涉及的基因、物种等

接下来我们用Python编写了一个爬虫，将根据上述的GO条目名称列表爬取一系列关键信息，例如这些GO条目中涉及的基因、物种等。

使用时，只需修改最开始的“data”和“name”中的项即可。“data”中为搜索并下载的GO条目名称列表，可以一次提供多个文件，中间以逗号分隔。“name”则对应了“data”中各文件代表的功能名称，用作输出文件名称。

如下示例，对于上述GO条目名称列表文件，修改“data = ['Oxidative_stress.txt']”和“name = ['Oxidative_stress']”后，在Python中执行即可。

#!/usr/bin/Python3.6
# -*- coding: utf-8 -*-
​
import os
import urllib
import urllib.request
​
#GO列表文件名称，可以一次提供多个文件，依次在列表中添加
#data = ['Reactive Oxygen.txt', 'Gene Expression.txt']
#GO搜索关键词名称，依次在列表中添加，和上述顺序对应
#name = ['Reactive Oxygen', 'Gene Expression']
​
data = ['Oxidative_stress.txt']
name = ['Oxidative_stress']
​
for i in data:
  i1 = name[data.index(i)]
  i2 = '_'.join(i1.split())
  os.system(f'mkdir -p {i2}_go')
  dat_input = open(i, 'r')
  dat_output = open(f'{i2}.go.txt', 'w')
  for line in dat_input:
    line = line.strip().split('\t')
    url = f"http://golr-aux.geneontology.io/solr/select?defType=edismax&qt=standard&indent=on&wt=csv&rows=100000&start=0&fl=bioentity_label,bioentity,bioentity_name,qualifier,annotation_class,annotation_extension_json,assigned_by,taxon,evidence_type,evidence_with,panther_family,type,bioentity_isoform,reference,date&facet=true&facet.mincount=1&facet.sort=count&json.nl=arrarr&facet.limit=25&hl=true&hl.simple.pre=%3Cem%20class=%22hilite%22%3E&hl.snippets=1000&csv.encapsulator=&csv.separator=%09&csv.header=false&csv.mv.separator=%7C&fq=document_category:%22annotation%22&fq=regulates_closure:%22{line[0]}%22&facet.field=aspect&facet.field=taxon_subset_closure_label&facet.field=type&facet.field=evidence_subset_closure_label&facet.field=regulates_closure_label&facet.field=annotation_class_label&facet.field=qualifier&facet.field=annotation_extension_class_closure_label&facet.field=assigned_by&facet.field=panther_family_label&q=*:*"
    html = urllib.request.urlopen(url).read().decode('utf-8').strip()
​
    os.system(f'mkdir -p {i2}_go/{line[0]}')
    tmp_file = open(f'{i2}_go/{line[0]}/go.txt', 'w')
    print(html, file = tmp_file)
    tmp_file.close()
​
    tmp_file = open(f'{i2}_go/{line[0]}/go.txt', 'r')
    for line1 in tmp_file:
      print(f'{i1}\t{line[0]}\t{line[1]}\t{line1.strip()}', file = dat_output)
    tmp_file.close()
​
  dat_input.close()
  dat_output.close()

Python脚本执行后，将自动根据“Oxidative_stress.txt”中提到的GO条目的id，在GO数据库中爬取相关的基因名称等信息。对于上述示例，结果最后默认生成在当前路径下，获得文件“Oxidative_stress_go”和“Oxidative_stress.go.txt”。

文件“Oxidative_stress.go.txt”中，包含了搜索GO条目时用的关键词、检索到的GO条目名称、这些GO条目中的基因、物种、以及子GO项等信息。这个文件是整合后的，包含了一开始搜索下载的“Oxidative_stress.txt”中提到的所有GO条目。如果期望分开查看，可点开另一结果文件夹“Oxidative_stress_go”，按不同的GO条目分开整理的基因、物种、以及子GO项等信息。

3 查询基因

这样，我们一开始期望的，寻找与“氧化应激”相关的功能基因的目的就实现了。

继续分析，只需要根据已知的物种分类，挑选出对应的物种出来就可以了。例如，人类物种在NCBI中的分类登记号为“9606”，因此如果想关注人类物种中与“氧化应激”相关的基因，只需要匹配“NCBITaxon:9606”就可以了。

以及如果做了RNA-seq，找到了一些重要的差异表达基因，也可根据基因名称在该文件中进行匹配，查看哪些重要的差异基因参与了“氧化应激”有关的功能。例如继续，文献中已经报道了蛋白“SOD”与氧化应激有关，RNA-seq的差异基因中也包含它，因此我们在这里查询验证一下“SOD”相关的基因。