第一个问题：研究最热门的基因是什么

在NCBI的ftp里面关于人的一些基因信息， 在 ：ftp://ftp.ncbi.nlm.nih.gov//gene 下载即可！

其中 gene2pubmed.gz 这个是NCBI的entrez ID号对应着该基因发表过的文章的ID号

链接是：ftp://ftp.ncbi.nlm.nih.gov//gene/DATA/gene2pubmed.gz

读入数据统计，绘制词云图

library(data.table)
library(wordcloud)
library(org.Hs.eg.db)
g2p <- fread('/data/NCBI/gene2pubmed.gz',data.table = F)
head(g2p)tb <- as.data.frame(table(g2p$GeneID))
tb <- tb[order(tb$Freq,decreasing = T),]
colnames(tb)[1]='gene_id'
head(tb)
ids=toTable(org.Hs.egSYMBOL)
head(ids)
tbs=merge(ids,tb,by='gene_id')
wordcloud(words = tbs$symbol, freq = tbs$Freq, min.freq = 1,max.words=200, random.order=FALSE, rot.per=0.35,colors=brewer.pal(8, "Dark2"))

我们发现TP53这个基因研究的最多。

第二个问题：你研究的基因都发表过什么文章？

接下来我们可以对单个基因进行探索。这里TP53这个基因太多了，我以ADORA1这个基因为例。不过这种方式不太适合该基因对应文献过多的情况。

pid <- g2p$PubMed_ID[g2p$GeneID==tbs$gene_id[tbs$symbol=="ADORA1"]]
#### 网络爬虫
url <- "https://pubmed.ncbi.nlm.nih.gov/"
library(tidyr)
library(rvest)
library(dplyr)
library(RCurl)
library(XML)
library(stringr)
pubmedinfo <- lapply(1:length(pid),function(x){id <- pid[x]u <- paste0(url,id)htdata <- read_html(u, encoding = "utf-8")message(x)title <- html_nodes(htdata,"h1")[[1]]%>% html_text() %>% str_trim()Sys.sleep(1)#doi <- html_nodes(htdata,"header div ul li span a")[[1]]%>% html_text() %>% str_trim()info <- html_nodes(htdata,"header div div div.article-source span")[[2]]%>% html_text() %>% str_trim()Sys.sleep(1)year <- gsub("\\D","",unlist(strsplit(info,";"))[1])year <- substr(year,1,4)jur <- html_nodes(htdata,"header div div div.article-source button")[[1]]%>% html_text() %>% str_trim()info <- data.frame(PubMed_ID = id,Journal = jur,Year = year,#DOI = doi,Title = title)outlines = paste0(id,jur,year,title,collapse='\t')writeLines(outlines, con=output)return(info)
})
pubmed_Info <- do.call(rbind,pubmedinfo)

绘制一个柱形图看看

staty = as.data.frame(table(pubmed_Info$Year))
dim(staty)
library(ggplot2)
ggplot(staty,aes(x = Var1,y=Freq))+geom_bar(stat = "identity")

有一个R包RISmed是可以用来探索pubmed数据库的数据的，有时候还是会挂。你可以尝试一下。下面代码，所以每一个循环设置了睡眠3s。可以延长睡眠时间解决。但文献条数太多就不建议使用，因为时间太长，获取太多也会挂。

library(RISmed)
linfo <- lapply(pid, function(uid){message(uid)query <- paste0(uid,"[UID]")Sys.sleep(3)search_query = EUtilsSummary(query = query,retmax=1)records<- EUtilsGet(search_query)pubmed_data <- data.frame(Title = records@ArticleTitle,Year = records@YearPubDate,journal = ISOAbbreviation(records),PubMed_ID = uid)return(pubmed_data)
})
ADORA1_pubmeddata <- do.call(rbind,linfo)

可以对文章标题绘制一个词云图，代码和前面【生信笔记 | 文本挖掘的一般流程

library("tm")
library("SnowballC")
docs <- Corpus(VectorSource(pubmed_Info$Title))
toSpace <- content_transformer(function (x , pattern ) gsub(pattern, " ", x))
docs <- tm_map(docs, toSpace, "/")
docs <- tm_map(docs, toSpace, "@")
docs <- tm_map(docs, toSpace, "\\|")
docs <- tm_map(docs, content_transformer(tolower))
docs <- tm_map(docs, removeNumbers)
docs <- tm_map(docs, removeWords, stopwords("english"))
docs <- tm_map(docs, removeWords, c("characterization", "molecular","comprehensive",'cell','analysis','landscape'))
docs <- tm_map(docs, removePunctuation)
docs <- tm_map(docs, stripWhitespace)
dtm <- TermDocumentMatrix(docs)
m <- as.matrix(dtm)
v <- sort(rowSums(m),decreasing=TRUE)
d <- data.frame(word = names(v),freq=v)
head(d, 10)
wordcloud(words = d$word, freq = d$freq, min.freq = 1,max.words=200, random.order=FALSE, rot.per=0.35,shape = 'pentagon',size=0.7,colors=brewer.pal(8, "Dark2"))

第三问题：统计近十年的文章发展状况。

但是这个方式只适合文献数量少的情况，文献数量多，就不适用。在NCBI上有所有文献的信息文件。可以下载后整理分析。

下载地址：https://ftp.ncbi.nlm.nih.gov/pubmed/baseline/

library(tidyr)
library(rvest)
library(dplyr)
library(RCurl)
library(XML)
library(stringr)
library(RISmed)###
ncbi <- "https://ftp.ncbi.nlm.nih.gov/pubmed/baseline/"
hn <- read_html(ncbi, encoding = "utf-8")
dlurl <- html_nodes(hn,"pre a")%>% html_text() %>% str_trim()
dlurl <- dlurl[-c(1:2)]
dlurl <- dlurl[!grepl(".md5$",dlurl)]
dlurl[1:3]
length(dlurl)
lapply(dlurl[1:length(dlurl)],function(u){au <- paste0(ncbi,u)file_destination <- paste0("literatureInfo/",u)download.file(au,destfile =file_destination)
})

下载后处理。

下载的文件都是xml格式，但是不是每个文件里面的标签信息都一致。下面也仅仅是获取文章发表年限，期刊以及文章标题的信息。文件里面没有时间的我用0替代了。

pubmedinfo <- lapply(38:length(xmlfiles), function(len){xml = xmlfiles[len]rootnode <- xmlParse(file = xml) %>% xmlRoot()file = str_split(xml,"/")[[1]][4]sink(paste0("xmlinfo/",file,".txt"))lapply(c(1:xmlSize(rootnode)),function(x){#xmlSize(rootnode)pmid = 0journal = NAArticleDate = 0title = NAxmldata <- xmlToList(rootnode[[x]][[1]])nm = names(xmldata)if("PMID" %in% nm){pmid <- xmldata[["PMID"]][["text"]]if(!is.character(pmid)){message(paste0(xml,"pmid:",x))}}if("Article" %in% nm){journal = xmldata[["Article"]][["Journal"]][["Title"]]if(!is.character(journal)){message(paste0(xml,"journal:",x))}ArticleDate = xmldata[["Article"]][["ArticleDate"]][["Year"]]if(is.null(ArticleDate)){ArticleDate = xmldata[["Article"]][["Journal"]][["JournalIssue"]][["PubDate"]][["Year"]]if(is.null(ArticleDate)){ArticleDate = xmldata[["DateCompleted"]][["Year"]]if(is.null(ArticleDate)){ArticleDate = 0}}}title = xmldata[["Article"]][["ArticleTitle"]]if(is.null(title)){title <- xmldata[["Article"]][["VernacularTitle"]]}if(length(title)>1){v <- unlist(title) %>% as.vector()t <- NULLfor(y in 1:length(v)){t <- paste(t,v[y])}title = t}}cat( paste(c(pmid,journal,ArticleDate,title),collapse ='\t'))cat('\n')rm(xmldata)})sink()rm(rootnode)if(is.integer(len/15)){rm(rootnode)gc()Sys.sleep(300)}onexmlinfo <- do.call(rbind,onexmlinfo)
})

在本地，要把所有文件的运行结束，是需要很长时间的。我这里只跑了部分文件【1008-1062】

xmlinfodir <- dir("xmlinfo/",".txt$",full.names = T)
library(data.table)
literinfo <- data.frame()
for(x in xmlinfodir){onelit <- fread(x,header = F)literinfo = rbind(literinfo,onelit)
}
head(literinfo)
table(literinfo$V3)

有好多文献没有获取到时间。

> literinfo$PubMed_ID[literinfo$ArticleDate == 0] %>% length()
[1] 12621

后面我发现，有点文献出版时间是在Article标签中，但很多是空值，有的xml文件里面有DateCompleted标签可以获取时间，有的又没有。可以通过RISmed包检索获取时间。

整理好所有信息，再绘图。

参考：

研究最热门的基因是什么