如果我们想探索一下什么基因研究的最多，那就是检索pubmed数据库资源。在 NCBI的ftp里面关于人的一些基因信息 ：
代码语言：javascript复制ftp://ftp.ncbi.nlm.nih.gov//gene下载即可！
其中 gene2pubmed.gz 这个是NCBI的entrez ID号对应着该基因发表过的文章的ID号
代码语言：javascript复制ftp://ftp.ncbi.nlm.nih.gov//gene/DATA/gene2pubmed.gz下载后的文件读入，进行词云可视化。
代码语言：javascript复制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这个基因研究的最多。
代码语言：javascript复制tbs <- tbs[order(tbs$Freq,decreasing = T),]
head(tbs)> head(tbs)
gene_id symbol Freq
33878 7157 TP53 10539
17376 1956 EGFR 6094
33847 7124 TNF 6072
23610 3569 IL6 4889
34713 7422 VEGFA 4882
23387 348 APOE 4650接下来，这里来探索一下你研究的基因，发表的文献，可以看看都发表在什么期刊，对题目进行文本挖掘，可以统计每年的发表文献数量等等。。。。
下面是探索ADORA1这个基因的研究情况。
信息是通过网络爬虫的形式获取，中间可能会断，所以下面代码是爬取一个就写入一个到本地文件ADORA1.txt中，如果中断，可以修改一下代码，接着爬，但这还只是适用于数量小的情况。
接着文章【什么基因研究最多？？】的g2p数据框。
代码语言：javascript复制pid <- g2p$PubMed_ID[g2p$GeneID==tbs$gene_id[tbs$symbol=="ADORA1"]]
#### 网络爬虫
url <- "https://pubmed.ncbi.nlm.nih.gov/"
output <- file("data/ADORA1.txt", open="wt")
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)
})
close(output)
pubmed_Info <- do.call(rbind,pubmedinfo)如果没有断，那么我们直接使用pubmedinfo变量
代码语言：javascript复制pubmed_Info <- do.call(rbind,pubmedinfo)
这里就可以看见基本的信息了。
先看看每一年发表的文章。
代码语言：javascript复制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数据库的数据的，有时候还是会挂。你可以尝试一下。下面代码，所以每一个循环设置了睡眠2s。可以延长睡眠时间解决。但文献条数太多就不建议使用。
代码语言：javascript复制library(RISmed)
linfo <- lapply(pid, function(uid){
message(uid)
query <- paste0(uid,"[UID]")
Sys.sleep(2)
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)
save(ADORA1_pubmeddata,file= "data/NCBI/pubmeddata.Rdata")下面就看看这些文章题目的词云图
代码语言：javascript复制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, "\|")
# Convert the text to lower case
docs <- tm_map(docs, content_transformer(tolower))
# Remove numbers
docs <- tm_map(docs, removeNumbers)
# Remove english common stopwords
docs <- tm_map(docs, removeWords, stopwords("english"))
# Remove your own stop word
# specify your stopwords as a character vector
docs <- tm_map(docs, removeWords, c("characterization", "molecular",
"comprehensive",'cell',
'analysis','landscape'))
# Remove punctuations
docs <- tm_map(docs, removePunctuation)
# Eliminate extra white spaces
docs <- tm_map(docs, stripWhitespace)
# Text stemming
# docs <- tm_map(docs, stemDocument)
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"))
我们可以看到，出现最多的就是ADORA1(Adenosine A1 Receptor) 这个基因的名称，理所当然，你可以把这几个词去掉再绘制一个图，再看看。当然你还可以统计一下他的期刊或其他信息。
上面2种方式，获取文献信息，对于文献个数多，总是会挂的，毕竟是基于网络请求的方法。在NCBI上，是有所有文献的信息数据的，可以直接下载，唯一的缺点就是数据量不小，而且是xml的文件，需要全部整理统计出来。下面是下载地址：
代码语言：javascript复制https://ftp.ncbi.nlm.nih.gov/pubmed/后续介绍.....................
参考：
研究最热门的基因是什么
什么基因研究最多？？
