Excel与基因名的故事（续）
8号的时候，刷知乎遇到一个问题。如何评价科学家重命名了多个人类基因，以避免被 Excel 自动纠正？没有别的方法吗？
虽然不知道知乎从何种途径观察到我在研究这个问题然后推了相关问题给我。但是这个问题确实推荐的不错。
从问题以及相关搜索中可以知道，7号，HGNC，也就是人类基因组组织基因命名委员会，发了一篇NG的Comment，介绍了他们新的基因命名指南。
Guidelines for human gene nomenclature 一文中提到：
Symbols that affect data handling and retrieval.
For example, all symbols that autoconverted to dates in Microsoft Excel have been changed (for example, SEPT1 is now SEPTIN1; MARCH1 is now MARCHF1); tRNA synthetase symbols that were also common words have been changed (for example, WARS is now WARS1; CARS is now CARS1).
也就是说，HGNC决定改名了！
啊，真是喜大普奔。虽然微软在该事中被诟病，但是，生物学家还是屈服了。怪不得我之前搜索的时候显示Sept2的官方名为SEPTIN2，原来他已经改完名了。
那为什么之前转换成功的是Sept2呢，一方面是数据库没有及时更新，还有就是，之前是鼠的基因，而HGNC管理的主要是人的基因，人的基因名中除了orf之外的字母都是全大写的。
本篇文章均为人的基因名。
续集
那为什么还有续集呢？
因为，我想测试一下代码的兼容性，替换sep和mar够不够，还需不需要替换其他的月份。那么，我就需要一个被Excel影响的基因名清单。
找遍全网居然没有？
哦不对，有一个。在公司微信群里，有人提到了一个R包HGNChelper。它的核心功能就是将不合法的基因名转变为合法的基因名（也就是HGNC官方认可的）。
稍微介绍一下：
library(HGNChelper)
human = c("FN1", "tp53", "UNKNOWNGENE","7-Sep", "9/7", "1-Mar", "Oct4", "4-Oct",
"OCT4-PG4", "C19ORF71", "C19orf71")
checkGeneSymbols(human)
out:
Human gene symbols should be all upper-case except for the 'orf' in open reading frames. The case of some letters was corrected.x contains non-approved gene symbols x Approved Suggested.Symbol
1 FN1 TRUE FN1
2 tp53 FALSE TP53
3 UNKNOWNGENE FALSE <NA>
4 7-Sep FALSE SEPT7
5 9/7 FALSE SEPT7
6 1-Mar FALSE MARC1 /// MARCH1
7 Oct4 FALSE POU5F1
8 4-Oct FALSE POU5F1
9 OCT4-PG4 FALSE POU5F1P4
10 C19ORF71 FALSE C19orf71
11 C19orf71 TRUE C19orf71
本例来源于：拯救那些被EXCEL篡改的基因名
看起来很厉害。但是，考虑到python中Pandas读取错误的基因名读进来后是日期类型，而R读进来是数字类型。虽然也不排除其他方式可能能够读取Excel中原始的字符2-Sep，但是对我而言，这个包并没有什么用。我暂时也不知道他们开发这个包是用于处理何种来源的原始数据。
不过，这个包的转换方式采用的是存储数据直接映射，那么，是不是就可以从中提取可能被影响的基因列表呢。
从他的GitHub下载了rda文件。
hgnc = new.env()
load('hgnc.table.rda',envir = hgnc)
print(unique(hgnc$hgnc.table$Approved.Symbol)[1:10]) # 数据太多，取前10
out:
[1] "NFYC-AS1" "IFITM2" "IFITM3" "PRDX6" "SEC24B-AS1"
[6] "ALDH1L1" "KNOP1" "CYB561D2" "KIR2DL4" "ARHGAP9"
NFYC-AS1也会被影响的吗？
Alias symbols是0808y08y，但是这个在Excel中不会被改变啊。不是很懂。
放弃此路！
自己动手，丰衣足食
找不到清单，那我就搞一个清单。
基因名，就从最官方的HGNC爬。
本来打算直接上爬虫的，但是仔细一看，HGNC贴心的提供了HGNC REST web-service。
那就容易了啊！
使用到的主要是两个接口。
| 示例 | 类型 | 功能 |
|---|---|---|
| http://rest.genenames.org/fetch/symbol/ZNF3 | fetch | 返回基因名对应的信息 |
| http://rest.genenames.org/search/symbol/ZNF3 | search | 搜索基因名，返回基因列表 |
第二个接口接受类UNIX的匹配符，如ZNF?、ZNF*，接口不区分大小写。
所有易导致错误的基因名都应当是以月份的简称开头的，所以，搜索sep2*就可以获得sep开头的基因列表，然后拿到每个基因的Previous symbols和Alias symbols，就可以获得基因所有的名字了。
首先编写一个兼容两个接口的查询函数。
import requests
from urllib.parse import urljoin # 用于拼接url
from collections import namedtuple
request_result = namedtuple('request_result',['request_type', 'content']) # 命名元组可以通过属性访问值
class RequestTypeError(Exception): # 自定义异常
pass
def request_gene(request_content,request_type = 'search',request_content_type = 'symbol'):
if request_type.lower().startswith('f'):
request_type = 'fetch'
elif request_type.lower().startswith('s'):
request_type = 'search'
else:
raise RequestTypeError # 当request_type得小写不是以f或s开头时就抛出异常
headers = {'Accept': 'application/json'}
domin = 'http://rest.genenames.org/'
try:
url = urljoin(domin,"/".join([request_type,request_content_type,request_content]))
res = requests.get(url,headers = headers)
content = res.json()
except Exception: # 兼容可能出现的所有错误
return None
num_found = content['response']['numFound'] #该内容给出结果的数量
if num_found:
return request_result(request_type,content['response']['docs']) #docs内存储了具体的内容，为列表
else:
return None
然后，写循环执行操作。
def main():
months = ['Jan','Feb','Mar','Apr','May','Jun','Jul','Aug','Sep','Oct','Nov','Dec']
results = []
for month in months:
month_result = request_gene(month+'*',request_type='search') #使用*匹配月份名开头的基因名
if month_result is None:
continue
for i in month_result.content: #content有多个结果
gene_symbol = i['symbol']
gene_result = request_gene(gene_symbol,request_type='fetch')
if gene_result is None:
continue
gene_result = gene_result.content[0] # fetch只有一个结果
if 'alias_symbol' in gene_result:
results.extend([gene_symbol,'alias_symbol',alias_symbol] for alias_symbol in gene_result['alias_symbol'])
if 'prev_symbol' in gene_result:
results.extend([gene_symbol,'prev_symbol',prev_symbol] for prev_symbol in gene_result['prev_symbol'])
return results
获得结果，使用Pandas生成Excel。
results = main()
import pandas as pd
genes = pd.DataFrame(results,columns=['gene_symbol','other_symbol_type','other_symbol'])
genes.to_excel('genes.xlsx',index=False)
打开Excel。
？？？
说好的自动转换呢。看来Pandas存储的Excel在打开时还是能保存原样的。不过双击单元格再Enter就变过去了，这Excel有点意思。
那怎么办呢，生成csv再打开应该就可以了。
genes.to_csv('genes.csv',index=False)
打开。
果然，可以了。然后另存为genes_after_open.xlsx
after_open= pd.read_excel('genes_after_open.xlsx')
after_open.head()
整合，然后筛选出所有被改变的基因名。
from datetime import datetime
genes['gene_symbol_after_open'] = after_open['gene_symbol']
genes['other_symbol_after_open'] = after_open['other_symbol']
genes[genes.gene_symbol_after_open.apply(lambda x:isinstance(x,datetime))]
官方名字没有被错误更改的（看来改的没有遗漏）。
再看下other_symbol。
genes[genes.other_symbol_after_open.apply(lambda x:isinstance(x,datetime))]
out:
看了下都是SEP基因和MARCH开头的基因。咦，不是说还有其他月份的吗？而且，不是说有一个MARC1和MARCH1被共同转换成1-Mar了吗？
查了下，原来MARC1基因的更改后名字居然是MTR开头的。
那，看来只进行基因名的匹配搜索已经不够了啊。
看了下接口，官方也提供了全文搜索的选项，接口类似http://rest.genenames.org/search/ZNF3。
改一下request_gene函数和main中的调用方式。
import requests
from urllib.parse import urljoin
from collections import namedtuple
request_result = namedtuple('request_result',['request_type', 'content'])
class RequestTypeError(Exception):
pass
def request_gene(request_content,request_type = 'search',request_content_type = None):
if request_type.startswith('f'):
request_type = 'fetch'
elif request_type.startswith('s'):
request_type = 'search'
else:
raise RequestTypeError
headers = {'Accept': 'application/json'}
domin = 'http://rest.genenames.org/'
try:
if request_content_type is None:
url = urljoin(domin,"/".join([request_type,request_content]))
else:
url = urljoin(domin,"/".join([request_type,request_content_type,request_content]))
res = requests.get(url,headers = headers)
content = res.json()
except Exception:
return None
num_found = content['response']['numFound']
if num_found:
return request_result(request_type,content['response']['docs'])
else:
return None
def main():
months = ['Jan','Feb','Mar','Apr','May','Jun','Jul','Aug','Sep','Oct','Nov','Dec']
results = []
for month in months:
month_result = request_gene(month+'*',request_type='search')
if month_result is None:
continue
for i in month_result.content:
gene_symbol = i['symbol']
gene_result = request_gene(gene_symbol,request_type='fetch',request_content_type='symbol')
if gene_result is None:
continue
gene_result = gene_result.content[0]
if 'alias_symbol' in gene_result:
results.extend([gene_symbol,'alias_symbol',alias_symbol] for alias_symbol in gene_result['alias_symbol'])
if 'prev_symbol' in gene_result:
results.extend([gene_symbol,'prev_symbol',prev_symbol] for prev_symbol in gene_result['prev_symbol'])
return results
results = main()
genes = pd.DataFrame(results,columns=['gene_symbol','other_symbol_type','other_symbol'])
genes.to_csv('genes.csv',index=False)
重新跑，重复之前的操作保存生成Excel。
after_open= pd.read_excel('genes_after_open.xlsx')
genes['gene_symbol_after_open'] = after_open['gene_symbol']
genes['other_symbol_after_open'] = after_open['other_symbol']
from datetime import datetime
pd.set_option('display.max_rows', None) #显示所有行
genes[genes.other_symbol_after_open.apply(lambda x:isinstance(x,datetime))]
看了下，确实包含了之前没找到的MTAR1基因。
统计一下都有哪些开头。
import re
incorrect_symbols = genes[genes.other_symbol_after_open.apply(lambda x:isinstance(x,datetime))].other_symbol
incorrect_symbols.apply(lambda x:re.match(r'([A-Za-z]*)\d*',x).group(1)).value_counts()
out:
SEPT 16
MARCH 11
OCT 9
Mar 9
SEP 6
APR 4
DEC 3
MARC 2
FEB 2
Oct 1
NOV 1
看来确实只有sep和mar需要变换。Mar开头的就离谱。看一下。
genes[genes.other_symbol_after_open.apply(lambda x:isinstance(x,datetime))][incorrect_symbols.apply(lambda x:x.lower().startswith('mar')).values]
out:
有一些基因的alias_symbol是Mar开头的。现在通用的基因名基本都是大写，这些应该不会再使用了吧，不管他们。MARC和MARCH的问题确实解决不了。
最后转换一下，添加一个大写。
genes['other_symbol_after_transform'] = genes.other_symbol_after_open.apply(lambda x: datetime.strftime(x,'%b%#d').replace('Sep','SEPTIN').replace('Mar','MARCHF').upper() if isinstance(x,datetime) else x)
看一下prev_symbol的复原情况。
genes[(genes.other_symbol_after_open.apply(lambda x:isinstance(x,datetime))) & (genes.other_symbol_type == 'prev_symbol')]
out:
除了可恶的MARC之外都可以。这个表共31行，应该包含了所有27个刚被改名的基因。
再看一下alias_symbol。
genes[(genes.other_symbol_after_open.apply(lambda x:isinstance(x,datetime))) & (genes.other_symbol_type == 'alias_symbol')]
out:
这个转换就比较不行了，主要是因为各种牛鬼蛇神都有，SEPTIN6、SEPTIN8居然还有SEP2的alias_symbol，只能在一些symbol上试图性的恢复一下。看来有一些基因的alias_symbol也不可靠。
总结
本文主要是通过找到所有被影响的基因名来验证之前的替换是否足够完备。从结果来看，之前的替换已经基本够用，就是需要在更换物种的时候进行相应的大小写切换和替换更改。
