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

看到九月份学徒在群里提问，写爬虫批量循环抓取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如果不翻墙就很卡，甚至链接中断（使用了vpn后的确会有好转，但是不知道为什么vpn隔几分钟就中断，放弃了，太折腾了）

解决方案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)

最后生成的文件如下：