最近成为了硕导，克隆基因的基本工作就是获取对应的基因组序列，启动子，内含子等信息。
为了降低学生学习获取核酸序列的难度，准备写个傻瓜脚本用来fetch基因组核酸序列文件。
因为不想用selenium，因为还要开启浏览器和驱动所以用BeautifulSoup和request。
# -*- coding: utf-8 -*-
"""
Created on Tue Dec 21 11:05:15 2021
@author: Bohan
"""
# -*- coding: utf-8 -*-
"""
Created on Mon Dec 20 17:18:07 2021
@author: Bohan
"""
import urllib
import requests
from bs4 import BeautifulSoup
while True:
print('该程序版权属于湖南农业大学农学院水稻生理生化团队所有，需要网络运行！')
print("输入NCBI的基因ID，如LOC4342505，本程序不支持其他数据库ID")
geneid=input("输入代码：") #输入NCBI的基因LOCID
geneidtrans=(''.join(filter(str.isdigit,geneid))) #去除所有非数字信息
genepage= "https://www.ncbi.nlm.nih.gov/gene/"
geneurl=genepage+geneidtrans #合成基因页面url
print("基因网址："+geneurl) #还是要有个Checkpoint
genepage=urllib.request.urlopen(geneurl)
genepagecontents=genepage.read()
soup1=BeautifulSoup(genepagecontents,"lxml")
LOCNM=(soup1.find("span",{"class":"gn"})).string #定位基因名称
LOCID0='LOC'+str((soup1.find("meta",{"name":"ncbi_uidlist"}))).split('content="')[1].split(" name=")[0]
LOCID="LOC"+(''.join(filter(str.isdigit,LOCID0)))
print("基因 ID："+LOCID) #还是要有个Checkpoint
print("基因名称："+LOCNM) #还是要有个Checkpoint
print("物 种："+(soup1.find("em",{"class":"tax"})).string) #还是要有个Checkpoint
print("基因注释："+(soup1.find("h1",{"id":"gene-name"})).text) #还是要有个Checkpoint
urlgenebank="https://www.ncbi.nlm.nih.gov"+soup1.find("a",{"title":"Nucleotide GenBank report"}).attrs['href'] #获取基因信息页面对应的基因组核酸链接url
print("基因组序列genebank链接："+urlgenebank)
ncid1=urlgenebank.split("https://www.ncbi.nlm.nih.gov/nuccore/")[1].split("?report=genbank&from=")[0] #获取关键核酸ID
ncid2=urlgenebank.split("https://www.ncbi.nlm.nih.gov/nuccore/")[1].split("?report=genbank&from=")[1].split("&to=")[0] #获取核酸抬头位置
ncid3=urlgenebank.split("https://www.ncbi.nlm.nih.gov/nuccore/")[1].split("?report=genbank&from=")[1].split("&to=")[1].split("&strand=true")[0] #获取核酸结束位置
ncid2pre=str(int(ncid2)-4000) #计算基因前4000bp的位置
ncid3aft=str(int(ncid3)+4000) #计算基因后4000bp的位置
gfurl="https://www.ncbi.nlm.nih.gov/nuccore/"+ncid1+"?report=genbank&from="+ncid2pre+"&to="+ncid3aft+"&strand=true" #计算基因组前后4000bp链接
print("基因组前后4000bp链接："+gfurl)
gfpage=urllib.request.urlopen(gfurl)
gfcontents=gfpage.read()
soup2=BeautifulSoup(gfcontents,"lxml")
valkeyword0=str(soup2.find("div",{"class":"seq gbff"})).split(" val=")[1].split(" virtualsequence=")[0]
valkeyword=(''.join(filter(str.isdigit,valkeyword0)))
#懒得解释了，大概就是获取NCBI核酸的母本ID
genbankfileurl="https://www.ncbi.nlm.nih.gov/sviewer/viewer.cgi?tool=portal&save=file&log$=seqview&db=nuccore&report=gbwithparts&id="+valkeyword+"&from="+ncid2pre+"&to="+ncid3aft+"&strand=on&withparts=on" #懒得解释了，大概就是获取NCBI核酸下载的真实链接
print("genebank文件链接："+genbankfileurl)
filename="GF"+'-'+LOCID+"-"+LOCNM+'.gb'
locateID=soup1.find("div",{"class":"genomerefseqs_slave"})
SSS=locateID.find_all(name='a')
mRNAID=(str(SSS[-5])).split('>')[1].split('<')[0]
proteinID=(str(SSS[-4])).split('>')[1].split('<')[0]
print("mRNA序列ID:"+mRNAID)
print("蛋白质序列ID:"+proteinID)
mRNAurl="https://www.ncbi.nlm.nih.gov/nuccore/"+mRNAID
proteinurl="https://www.ncbi.nlm.nih.gov/protein/"+proteinID
mrnapage=urllib.request.urlopen(mRNAurl)
mrnacontents=mrnapage.read()
soup3=BeautifulSoup(mrnacontents,"lxml")
locateM=str(soup3.find("meta",{"name":"ncbi_uidlist"}))
mrnacc=(''.join(filter(str.isdigit,locateM)))
proteinpage=urllib.request.urlopen(proteinurl)
proteincontents=proteinpage.read()
soup4=BeautifulSoup(proteincontents,"lxml")
locateP=str(soup4.find("meta",{"name":"ncbi_uidlist"}))
proteinacc=(''.join(filter(str.isdigit,locateP)))
proteinsequrl="https://www.ncbi.nlm.nih.gov/sviewer/viewer.cgi?tool=portal&save=file&log$=seqview&db=protein&report=gpwithparts&id="+proteinacc+"&extrafeat=undefined&withparts=on"
mrnasequrl="https://www.ncbi.nlm.nih.gov/sviewer/viewer.cgi?tool=portal&save=file&log$=seqview&db=nuccore&report=gbwithparts&id="+mrnacc+"&withparts=on"
print("正在下载基因组序列 "+filename+" 文件中")
fileread = requests.get(genbankfileurl,allow_redirects=True) #下载基因组序列，不解释了
with open (filename, 'wb') as f:
f.write(fileread.content)
f.close
print(filename+"已下载完成")
mrnafilename="mRNA"+'-'+LOCID+"-"+mRNAID+'.gb'
profilename="Protein"+'-'+LOCID+"-"+proteinID+'.gp'
print("正在下载mRNA序列 "+mrnafilename+" 文件中")
mfileread = requests.get(mrnasequrl,allow_redirects=True) #下载mRNA序列，不解释了
with open (mrnafilename, 'wb') as fm:
fm.write(mfileread.content)
fm.close
print(mrnafilename+"已下载完成")
print("正在下载蛋白序列 "+profilename+" 文件中")
pfileread = requests.get( proteinsequrl,allow_redirects=True) #下载蛋白序列，不解释了
with open (profilename,'wb') as fp:
fp.write(pfileread.content)
fp.close
print(profilename+"已下载完成")
control=input("是否继续运行？请输入y/n \n")
if control==('n'):
break
运行效果如下：
输入代码：LOC4342505
基因网址：https://www.ncbi.nlm.nih.gov/gene/4342505
基因 ID：LOC4342505
基因名称：GA20OX3
物 种：Oryza sativa Japonica Group
基因注释：GA20OX3 Gibberellin 20 oxidase 3 [ Oryza sativa Japonica Group (Japanese rice)
]
基因组序列genebank链接：https://www.ncbi.nlm.nih.gov/nuccore/NC_029262.1?report=genbank&from=3699596&to=3704922
基因组前后4000bp链接：https://www.ncbi.nlm.nih.gov/nuccore/NC_029262.1?report=genbank&from=3695596&to=3708922&strand=true
genebank文件链接：https://www.ncbi.nlm.nih.gov/sviewer/viewer.cgi?tool=portal&save=file&log$=seqview&db=nuccore&report=gbwithparts&id=996703426&from=3695596&to=3708922&strand=on&withparts=on
mRNA序列ID:XM_026026987.1
蛋白质序列ID:XP_025882772.1
正在下载基因组序列 GF-LOC4342505-GA20OX3.gb 文件中
GF-LOC4342505-GA20OX3.gb已下载完成
正在下载mRNA序列 mRNA-LOC4342505-XM_026026987.1.gb 文件中
mRNA-LOC4342505-XM_026026987.1.gb已下载完成
正在下载蛋白序列 Protein-LOC4342505-XP_025882772.1.gp 文件中
Protein-LOC4342505-XP_025882772.1.gp已下载完成
是否继续运行？请输入y/n
已生成可运行的单个exe文件，请有需要的同学自行下载。
如有转载请注明出处。
链接：https://pan.baidu.com/s/1A8nJ81uq_obqKMEgONLh9Q
提取码：0213
网友评论