基于PubMed数据库挖掘研究最多的基因与以及有关某基因发表了多少篇文献？这些文献有什么特点？？？

DoubleHelix

发布于 2022-11-24 10:36:02

7270

发布于 2022-11-24 10:36:02

如果我们想探索一下什么基因研究的最多，那就是检索pubmed数据库资源。在 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这个基因研究的最多。

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数据框。

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变量

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数据库的数据的，有时候还是会挂。你可以尝试一下。下面代码，所以每一个循环设置了睡眠2s。可以延长睡眠时间解决。但文献条数太多就不建议使用。

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")

下面就看看这些文章题目的词云图

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"))