根据NCBI中的蛋白名在uniprot中查找对应蛋白登录号及基因名

想法来源

首先,我要感谢我从事蛋白质组分析工作中遇到的形(xi)形(qi)色(gu)色(guai)的客户,我还是非常热衷于为客户服务的。最后完成的脚本算是我为了完成客户需求一点点摸索出来的成果。

第一个客户

第一次涉及到这个需求的客户指出,在我们完成蛋白质组结果中(客户指定使用NCBI数据库)没有提供蛋白对应的基因名,需要我们在提供的结果中添加一下对应的基因名。

???
???
来看看数据库,除了特殊的自建库及客户提供的转录组数据单框翻译外,我们常用的蛋白数据库基本就是Uniprot和NCBI,就以这两个数据的.fasta文件对比。

>ACD58760.1 hypothetical protein PXO_00731 [Xanthomonas oryzae pv. oryzae PXO99A]
MGILIYLVPASALWALIATALAFVRGRQLRAESSELASTQDRLGSYQAALSQQKARAAATTLELESLQRS
YAVLKQSLEQQEQNASEQQEQNASEQQAATAGQVIPMVLVQRLDIASEIGTLFGHVARVARILRRYSAYS
RGHNAPEPATARYDLHWLADCLHSFDQLGHALVHGNVAALITACQDLLSMYEHYLKDGSGYNSRDTFQRL
SHDVPLSEATDAIRSIIVKATLAQDVQDAVQNDEVTANVG
上面这个是NCBI的,下面这个是uniprot
>tr|A0A0K0GK93|A0A0K0GK93_XANOP Uncharacterized protein OS=Xanthomonas oryzae pv. oryzae (strain PXO99A) OX=360094 GN=PXO_00731 PE=4 SV=1
MGILIYLVPASALWALIATALAFVRGRQLRAESSELASTQDRLGSYQAALSQQKARAAATTLELESLQRS
YAVLKQSLEQQEQNASEQQEQNASEQQAATAGQVIPMVLVQRLDIASEIGTLFGHVARVARILRRYSAYS
RGHNAPEPATARYDLHWLADCLHSFDQLGHALVHGNVAALITACQDLLSMYEHYLKDGSGYNSRDTFQRL
SHDVPLSEATDAIRSIIVKATLAQDVQDAVQNDEVTANVG

虽然在NCBI的数据库中蛋白描述部分确实有PXO_00731,但是这是这个蛋白缺少验证数据的原因产生的特殊情况,在NCBI的数据库文件中很少会在蛋白描述部分放入蛋白的基因名。(其实这个PXO_00731也不是基因名,只是还没系统命名导致的一个简单命名,暂代了基因名),而在uniprot数据库中基因名是有专门的位置标注了“GN=PXO_00731”。
这也就是客户想要的结果,但是咋弄啊,一个是NCBI的数据库,一个是uniprot的数据库。(你是不是在心里说了句一个个查,一个个粘?
再来看看NCBI中搜索这个蛋白的内容NCBI中该蛋白的详细信息
在这段信息中该蛋白作为CDS区翻译的结果对应的基因序列名是有标注的。
我看这段内容看着挺简单的呀,要是把每个蛋白的这部分信息用python爬虫爬出来,再用perl处理一下是不是就能拿到想要的结果了?(也许你会很奇怪为啥不直接用python处理,因为我的python真的很L J等于没学,perl相对来说用的强点)
于是我参考了论坛里大佬分享的一个爬虫教学Python爬虫超详细讲解(零基础入门,老年人都看的懂)
但是,但是啊,爬不出来,真的爬不出来。
我把完整的网页源码加载出来放在.txt文件里查找要的信息,发现NCBI对蛋白描述那部分的内容并不在源码中,加上我对前端,js,css啥的也是基本等于0,也不知道是不是我水平不到家,还是NCBI有意防着我。捣鼓了半天也没个结果,我选择了放弃,跟我这位可爱的客户说了“臣妾做不到啊!”

第二位客户

第二位客户,跟第一位客户的背景是一样的,但是第二位客户客户给我的售后申请中附加了一个他期望的结果样式图(就是一个表格的截图),我就不放原图了,害怕客户看到我用他的图片写文章会给我一个爱(fen)的(nu)拥(tou)抱(su)。

蛋白Group主要蛋白蛋白描述基因名
来自uniprot的机密机密机密机密

我本能的又想说一句“臣妾。。。好像。。。又行了”。
这位客户在第一列中加入的蛋白名并不是我们给他使用的NCBI数据库中的名字,而是这个蛋白在uniprot中的蛋白名。对呀,这几个数据库其实是有互相保留对方的联系。。呸,蛋白名的。如果我拿NCBI中的蛋白名去uniprot中搜索会不会。。。。。,嘿嘿,嘿嘿嘿。

看看情况如何
在uniprot中搜索NCBI蛋白名获得的结果
基因名有了!!!uniprot中的蛋白名都有了。
那如果我再爬uniprot源码中的这两个结果是不是就可以获得想要的结果了?
于是我又去看那个老年人教程,遇事不决,再来一遍。
阔德行了!这次真的阔德行了。谢谢大佬,谢谢我亲爱的客户。

代码分享

原本思路是python和perl结合的,python爬网页代码,perl从所有代码中摘取数据。
但是突然想在论坛里发篇文章还是好(拉)好(高)改(逼)善(格)一下。所以正则部分还是循着大佬的思路用python进行。

import requests
import urllib.request, urllib.error
import re
import os
#from bs4 import BeautifulSoup

# 判断文件夹是否存在,不存在则创建一个
#if not os.path.exists("./temp"):#这个文件夹是用来存放网页源码的.txt文件
#   os.mkdir("./temp")

with open("protein.txt", "r") as f:
    data = f.readlines()
   # print (len(data))
   # data = data.strip('\n')
   # print (data[0])

    #print (data[i])
baseurl = "https://www.uniprot.org/uniprot/?query="
size = "&sort=score"  #用来拼接和protein.txt中的蛋白名拼接成链接baseurl+“name”+size
findGn = re.compile(r'<span class="shortName">(.*?)</span>', re.S)
findUn = re.compile(r'id=\"checkbox_(.*?)\" type="checkbox\"/',re.S)#其实应该要摘取的不是这个正则的,但是要摘的位置的正则有些麻烦,取了个巧,用这个正则可以获取相同的信息。

def askURL(url):
    head = {   
    "User-Agent": "Mozilla / 5.0(Windows NT 10.0; Win64; x64) AppleWebKit / 537.36(KHTML, like Gecko) Chrome / 80.0.3987.122  Safari / 537.36"
    }
    # 用户代理,表示告诉豆瓣服务器,我们是什么类型的机器、浏览器(本质上是告诉浏览器,我们可以接收什么水平的文件内容)

    request = urllib.request.Request(url, headers=head)
    html = ""
    try:
        response = urllib.request.urlopen(request)
        html = response.read().decode("utf-8")
    except urllib.error.URLError as e:
        if hasattr(e, "code"):
            print(e.code)
        if hasattr(e, "reason"):
            print(e.reason)
    return html

a = ("./temp.txt")
with open (a,"w",encoding='utf-8') as first:
    first.write("name\tuniprot name\tdes\n")
for i in range(0, len(data)):  # 准备的蛋白名有多少个文件有多少行
    name = data[i].strip('\n')
    url = baseurl + name + size
    html = askURL(url)  # 保存获取到的网页源码
    m = ("./temp/"+name+".txt")
    #with open (m,"w",encoding="utf-8") as htlm:  #注释掉的部分是将网页源码以.txt格式存放入temp文件夹中
    #   htlm.write(html)
    gn = re.findall(findGn, html)[0]
    un = re.findall(findUn, html)[0]
    with open (a,"a",encoding='utf-8') as out:
        out.write(name + ("\t")+ un + ("\t")+ gn +("\n"))
    print(name)

b = ("./result.txt")
with open (b,"w",encoding="utf-8") as refirst:
    refirst.write ("name\tuniprot name\tGene Name\n")

with open("temp.txt", "r") as temp:
    teda = temp.readlines()

for j in range(1,len(teda)):
    res1 = teda[j].replace("\t ","\t")
    res2 = res1.replace(" \n","\n")
    res3 = res2.replace("<strong>","")
    res4 = res3.replace(" ",",")
    res5 = res4.replace("</strong>","")#删除不要的标签,将uniprot中出现的gene name用,拼接放在结果文件中
    #我承认这个循环。。。。很low
    with open (b,"a",encoding="utf-8") as reinf:
        reinf.write(res5)

总的思路还是比较简单,事先准备一个protein.txt文件夹把要爬取的蛋白名(NCBI)粘进去,换行符分隔。temp中存放的是爬取的网页源码(代码中该部分被注释掉了),temp.txt文件中存放的是为处理过的正则爬取结果,result.txt中就是结果啦。

反思

其实这个思路还是存在一些潜在问题的:
1.通常NCBI中获取的物种蛋白库会比uniprot大很多,这是因为NCBI在处理冗余序列时并没有那么细,而uniprot则会对那些高重复的冗余序列进行归并,所以两个网站提供的蛋白名可能会出现不一一对应的情况(至少我目前测的两个物种没遇到,感觉会有这种情况)。
2.另外两个网站都是会周期性维护更新的,会不会某一方对一个蛋白进行了信息更新,而另一方没有跟上脚步呢?还是有可能的,这应该也会导致二者对应不上吧,我感觉。。。。
3.还有如果客户是很早之前做的分析,很晚才来做的售后,这个数据库早就更新好几代了,那以前的蛋白名(NCBI)在现在的uniprot官网中搜还能搜到吗?我不知道。。。。

Logo

为开发者提供学习成长、分享交流、生态实践、资源工具等服务,帮助开发者快速成长。

更多推荐