爬虫那么危险，干嘛不直接基因数据库下载文件呢？

生信技能树

发布于 2019-10-25 03:30:26

2K0

发布于 2019-10-25 03:30:26

看到九月份学徒在群里提问，写爬虫批量循环抓取NBCI数据库的基因信息，但是经常掉线，还有可能被封，求助！我简单指点了他去找基因数据库文件即可，随便邀请他总结投稿如下： 分割线 一大早师姐给了个小任务，让我帮忙给注释下一批基因，格式类似如下：

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

解决方案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如果不访问外国网站就很卡，甚至链接中断（使用了V**后的确会有好转，但是不知道为什么V**隔几分钟就中断，放弃了，太折腾了）

解决方案2

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

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

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

特别感谢！

于是，思路转变。

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

找到目的界面：

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

于是导入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)

最后生成的文件如下：