投稿-批量基因annotation

看到九月份学徒在群里提问,写爬虫批量循环抓取NBCI数据库的基因信息,但是经常掉线,还有可能被封,求助!

我简单指点了他去找基因数据库文件即可,随便邀请他总结投稿如下:

一大早师姐给了个小任务,让我帮忙给注释下一批基因,格式类似如下:

1571816785359

问了具体后原来是ncbi上的信息,相当于在ncbi上在gene库中查找,然后爬取目标信息。如下:

1571817018744

解决方案1:

我的第一反映就是用python爬虫去爬,想倒是挺好想的,但是太久没用python了,语法都忘得差不多了,于是就考虑使用R语言来做:

代码如下(如果可以翻墙这种方法也可行):

#install.packages("RCurl")
#install.packages("XML")
library(RCurl)
library(stringr)
library(XML)
library(clusterProfiler)

rm(list=ls())
options(stringsAsFactors = F)
# 读入基因列表:
genes <- read.table("HSC_MPP1_BMvsoldBM.csv",sep = ",",header = T)[,1]

# 将gene symbol转为entrze ID:
genes <- bitr(genes, fromType="SYMBOL", toType="ENTREZID", OrgDb="org.Hs.eg.db")

# 网址数据框:
genes$NCBI_url <- paste("https://www.ncbi.nlm.nih.gov/gene/",genes$ENTREZID,sep="")
head(genes)

# 根据xpath获取节点内容:
getNodesTxt <- function(html_txt1,xpath_p){
 els1 = getNodeSet(html_txt1, xpath_p)
 # 获得Node的内容,并且去除空字符:
 els1_txt <- sapply(els1,xmlValue)[!(sapply(els1,xmlValue)=="")]
 # 去除\n:
 str_replace_all(els1_txt,"(\\n )+","")
}

# 处理节点格式,为character且长度为0的赋值为NA:
dealNodeTxt <- function(NodeTxt){
 ifelse(is.character(NodeTxt)==T && length(NodeTxt)!=0 , NodeTxt , NA)
}

# 简单使用xpath来获取:
for(i in 1:nrow(genes)){ 
 # 获得网址: 
 doc <- getURL(genes[i,"NCBI_url"]) 
 cat("成功获得网页!\t") 
 # 获得网页内容 
 html_txt1 = htmlParse(doc, asText = TRUE) 
 # 获得Full Name: 
 genes[i,"FullName"] <- dealNodeTxt(getNodesTxt(html_txt1,'//*[@id="summaryDl"]/dd[2]/text()')) 
 cat("写入基因\t") 
 # 获得HGNC ID: 
 genes[i,"HGNC_ID"] <- str_replace_all(dealNodeTxt(getNodesTxt(html_txt1,'//*[@id="summaryDl"]/dd[3]/a')),"HGNC|:","") 
 cat("写入HGNC_ID\t") 
 # 获得Gene type: 
 genes[i,"GeneType"] <- dealNodeTxt(getNodesTxt(html_txt1,'//*[@id="summaryDl"]/dd[5]/text()')) 
 cat("写入GeneType\t") 
 # 获得summary: 
 genes[i,"Summary"] <- ifelse(length(getNodesTxt(html_txt1,'//*[@id="summaryDl"]/dd[10]/text()'))!=0,getNodesTxt(html_txt1,'//*[@id="summaryDl"]/dd[10]/text()'),NA) 
 cat("写入Summary\n") 
 print(paste("完成第",i,"个了!")) 
}

可是,还没爬几条数据就出现了问题:因为网络问题,一是过快访问,存在被封ip的危险,二是访问ncbi如果不翻墙就很卡,甚至链接中断(使用了vpn后的确会有好转,但是不知道为什么vpn隔几分钟就中断,放弃了,太折腾了)

解决方案2

于是咨询了下Jimmy老师,Jimmy说有一种文件就可以直接有这些信息,于是乎:

发朋友圈求助各路大神!!!

于是还是有大神来帮助的:

1571817322401

特别感谢!

于是,思路转变。

进入 http://www.ensembl.org/index.html 官网

找到目的界面:

可以看到里面有很多选项可以自己设置

1571817693474

根据自己的需求自行设置吧:

1571817376724

点击go下载文件,打开后如下:

1571817613453

基本上有点像了。于是导入R中进行后续处理,后面就很简单了,大家可以根据自己的需求自己去设置了。

还是贴上我写的垃圾代码,大家看看就好,别评论,丢不起这个人。

rm(list=ls())
options(stringsAsFactors = F)
# 读入基因列表:
genes <- read.table("HSC_MPP1_BMvsoldBM.csv",sep = ",",header = T)[,1]
tmpdb <- read.delim("mart_export.txt",sep = ",",header = T)
tmp=tmpdb[,c(1:3,5,6)]

#去除重复
tmp=tmp[!duplicated(tmpdb[,1]),]

target=tmp[match(genes,tmp[,1]),]
all <- read.table("HSC_MPP1_BMvsoldBM.csv",sep = ",",header = T)
all=all[,c(1,2,6,3)]
colnames(all)
colnames(tmp)

target=merge(tmp,all,by.y = "X",by.x = "Gene.name")
head(target)
target=target[,c(1,6:8,3:5)]
target=target[order(target$p_val_adj),]

write.table(target,file = "gene_annotation.xls",sep = "\t",row.names = F,col.names = T)

最后生成的文件如下:

1571819075050

可以,很棒,就这样了,滚去复习六级了。

Comments are closed.