R语言挖掘Pubmed数据库

生信喵 发表于 2025-11-19 10:52:03

<h1>需求描述</h1>
<p>通过API检索NCBI的Pubmed数据库，批量获得文章信息，整理出摘要里基因的词频。</p>
<h1>应用场景</h1>
<ul>
<li>老板交给我的这个基因好陌生，用这套代码跑出表格，一眼望去，发现它的好朋友我都很熟，于是，熟悉的思路、方法都能用上了。</li>
<li>文章做到第三部分卡住了，我的基因到底调控了谁，或者谁调控了它，或者它跟谁是好朋友一同发挥作用？用这套代码跑出表格，一眼望去，发现目的基因跟我熟悉的一类分子一同出现在某篇文章的摘要里，我们实验室有非常完善的研究这类分子的实验体系，接下来就瞄准这个分子了。</li>
</ul>
<p>只要设置好检索词，就可以批量检索并整理文献和摘要，让你的头脑风暴比别人快10倍。</p>
<p>包含三个模块，相对独立，可分别运行。全部运行完将让你对目的基因有个整体的认识。</p>
<p>【模块一】文章的增长趋势</p>
<p>【模块二】发表了哪些文章</p>
<p>【模块三】从摘要找好朋友</p>
<h1>环境设置</h1>
<p>使用国内镜像安装包</p>
<pre><code class="language-r">options("repos"= c(CRAN="https://mirrors.tuna.tsinghua.edu.cn/CRAN/"))
options(BioC_mirror="http://mirrors.ustc.edu.cn/bioc/")
BiocManager::install("rentrez")
BiocManager::install("pubmed.mineR")
install.packages("jiebaR")
install.packages("bibliometrix")
install.packages("DT")
install.packages("htmlwidgets")
</code></pre>
<p>加载包</p>
<pre><code class="language-{r}">library(data.table)
library(stringr)
library(ggplot2)
library(jiebaR)
library(rentrez)
library(pubmed.mineR)
library(bibliometrix)
library(ggrepel)
library(ggthemes)
library(AnnotationDbi)
library(DT)
library(htmlwidgets)

Sys.setenv(LANGUAGE = "en") #显示英文报错信息
options(stringsAsFactors = FALSE) #禁止chr转成factor
</code></pre>
<h1>输入文件</h1>
<h2>常用词</h2>
<p>common_words_new.v2.txt，常用词。【模块三】计算词频时，要去除这些常用词。</p>
<pre><code class="language-{r}">common_words_new <- read.table("common_words_new.v2.txt")
common_words_new <- common_words_new$V1
</code></pre>
<h2>基因名</h2>
<p>每行一个基因。包含基因名、别名、位点，【模块三】将从摘要中统计这些基因名。还可以换成你想统计词频的其他类型的词。</p>
<pre><code class="language-{r}">(load("HGNCdata.rda")) #位于当前文件夹
head(HGNCdata)
</code></pre>
<h1>参数设置</h1>
<p>在这里设置检索词、检索的文章发表时间范围</p>
<pre><code class="language-{r}">targetGene <- "PD-L1" #检索的基因名
year <- 2017:2025 #检索的文章的发表年份范围
</code></pre>
<p>这里用的检索词跟你在网页版Pubmed中用的检索词是一样的，初级用法就是用AND或OR连接你希望检索到的内容，例如基因名、疾病、处理、物种等。医学领域推荐直接用MeSH terme来限定，Medical Subject Headings(MeSH) refer to wiki https://en.wikipedia.org/wiki/Medical_Subject_Headings</p>
<p>如果用基因名作为检索词，建议把基因的别名都放进来，可以通过下面的代码提取：</p>
<pre><code class="language-{r}">targetGene.Symbol <- HGNCdata
targetGene.Synonyms <- HGNCdata
targetGeneInfo <- rbind(targetGene.Symbol, targetGene.Synonyms)
targetGeneInfo

term <- paste(as.character(targetGeneInfo$Approved.Symbol), str_replace_all(targetGeneInfo$Synonyms, ", ", " OR "), sep = " OR ")
term
# "CD274 OR B7-H OR B7H1 OR PD-L1 OR PDL1 OR B7-H1"

# 还可以加上疾病的名字，例如GBM的MeSH term: Glioblastoma
#term <- paste0 ("(", term, ") AND Glioblastoma")
#term

# 或者自己手写基因名作为检索词
#term <- "CD274 OR PD-L1 OR PDL1"
</code></pre>
<h1>【模块一】文章的增长趋势</h1>
<pre><code class="language-{r}">#查询NCBI里的数据库简称
#entrez_dbs()
#pubmed里有哪些可以查询的fields
#entrez_db_searchable(db = "pubmed")

#先写个函数，获取特定年份范围内发表的带有检索词的文章数量
search_year <- function(year, term){
query <- paste(term, "AND (", year, ")")
entrez_search(db="pubmed", term=query, retmax=0)$count
}

# 用前面“参数设置”里定义的检索词检索，文章越多，等的越久
papers <- sapply(year, search_year, term=term, USE.NAMES=FALSE)

plot(year, papers, type='b', main="The PD-L1 papers")
</code></pre>
<p><img src="https://tcapi.voiceclouds.cn:8444/uploads/2025/11/19/mi5e8gkx87waiue9zrs.png" alt="" /></p>
<h1>【模块二】发表了哪些文章</h1>
<h2>检索并从结果中提取出PMID、发表日期、文章名、期刊。</h2>
<pre><code class="language-{r}">#我们先把retmax设为0，先不让它返回结果，这步只为了查看一共多少条记录
pre.result <- entrez_search(db="pubmed", term=term, retmode = "xml", retmax = 0)
pre.result$count #一共有42340条记录

#返回所有检索到的结果
#result <- entrez_search(db="pubmed", term=term, retmode = "xml", retmax = pre.result$count)
#网速慢的话，先返回前200条记录
result <- entrez_search(db="pubmed", term=term, retmode = "xml", retmax = 200)

#download all the paper info
#记录太多会被拒绝，因此每次提取200篇
n <- 200 #每次读入的记录数量
res <- c()
for (i in seq(1,length(result$ids),n)) {
multi_summ <- entrez_summary(db="pubmed",id=result$ids)
date_and_cite <- extract_from_esummary(multi_summ, c("uid","pubdate", "authors","title", "fulljournalname","elocationid"))
res1 <- data.frame(t(date_and_cite))
res1 <- data.frame(lapply(res1, as.character), stringsAsFactors=FALSE)
res <- rbind(res,res1)
}

# 整理author name
tmp <- sub('^...............','', res$authors) # 如果想从左侧删除N个字符用
tmp <- gsub("\", \"", ", ", tmp)
tmp <- data.frame(sapply(tmp, function(x) unlist(strsplit(x,'\"'))),stringsAsFactors = F)[,1]
tmp <- data.frame(sapply(tmp, function(x) unlist(strsplit(x,'"),'))),stringsAsFactors = F)[,1]
tmp <- data.frame(sapply(tmp, function(x) unlist(strsplit(x,'",'))),stringsAsFactors = F)[,1]
res$authors <- tmp

#把结果保存到文本文件
write.table(res, "output_paper.txt", quote = F, sep = "\t", row.names = F)
</code></pre>
<p><img src="https://tcapi.voiceclouds.cn:8444/uploads/2025/11/19/mi5ed11rjwqio36k0kg.png" alt="" /></p>
<h2>添加文章的Pubmed链接，保存成网页格式</h2>
<pre><code class="language-{r}">uid <- res$uid
pubdate <- res$pubdate
pmcrefcount <- res$pmcrefcount
authors <- res$authors
title <- res$title
fulljournalname <- res$fulljournalname
elocationid <- res$elocationid

# 写个函数
createLink <- function(base,val) {
sprintf('<a href="%s" class="btn btn-link" target="_blank" >%s</a>',base,val)
}

# 按照网址规律给PMID和题目加上Pubmed链接
res <- data.frame(uid = createLink(paste0("https://www.ncbi.nlm.nih.gov/pubmed/?term=",uid),uid),
               pubdate = pubdate,
               authors = authors,
               title = createLink(paste0("https://www.ncbi.nlm.nih.gov/pubmed/?term=",uid),title),
               fulljournalname = fulljournalname,
               elocationid = elocationid,
               stringsAsFactors = F)
res <- na.omit(res)
y <- DT::datatable(res,escape = F,rownames=F)

#保存到网页格式的文件
DT::saveWidget(y,"output_paper.html")
</code></pre>
<p><img src="https://tcapi.voiceclouds.cn:8444/uploads/2025/11/19/mi5eepgivmt4phbiw7.png" alt="" /></p>

生信喵 发表于 2025-11-19 10:52:46

<h1>【模块三】从摘要找好朋友</h1>
<h2>检索</h2>
<p>这段检索跟前面“【模块二】发表了哪些文章”的开头一样，就不重复运行了。</p>
<pre><code class="language-r">pre.result <- entrez_search(db="pubmed", term=term, retmode = "xml", retmax = 0)
pre.result$count

#返回所有检索到的结果
#result <- entrez_search(db="pubmed", term=term, retmode = "xml", retmax = pre.result$count)
#此处返回前200条记录
result <- entrez_search(db="pubmed", term=term, retmode = "xml", retmax = 200)
</code></pre>
<h2>把xml文件整理成数据框格式</h2>
<p>xml中包含title、authors、year、journal、key_words、doi、pmid、abstract等信息，整理成data.frame格式。</p>
<pre><code class="language-{r}">n <- 200 #每次读入的记录数量
abstract0 <- c()
for (i in seq(1, length(result$ids), n)) {
rec <- parse_pubmed_xml(entrez_fetch(db = "pubmed", id = result$ids, rettype = "xml"))
abstract1 <- as.data.frame(do.call('rbind', rec))
abstract0 <- rbind(abstract0, abstract1)
}
abstract1

length(abstract1$abstract) # how many abstract
abstract <- unlist(abstract0$abstract)

#查看其中某一篇文章的摘要
#abstract
</code></pre>
<h2>提取abstract中的基因名，并统计出现频率。</h2>
<p>先把摘要分割成单个的单词，然后根据输入文件common_words_new.v2.txt去掉常用词，再根据HGNCdata.rda挑出基因名。</p>
<pre><code class="language-{r}"># 先用标点符号分割出单词
tempa <- unlist(strsplit(abstract, ",",fixed = T));
tempb <- unlist(strsplit(tempa, ":",fixed = T));
tempc <- unlist(strsplit(tempb, ";",fixed = T));
tempd <- unlist(strsplit(tempc, "'",fixed = T));
tempe <- unlist(strsplit(tempd, " ",fixed = T));
tempf <- unlist(strsplit(tempe, "/",fixed = T));
tempf <- unlist(strsplit(tempe, "\\|",fixed = T));
tempf <- unlist(strsplit(tempf, "(",fixed = T));
tempf <- unlist(strsplit(tempf, ")",fixed = T));
tempf <- unlist(strsplit(tempf, "-",fixed = T));
tempf <- unlist(strsplit(tempf, ".",fixed = T));

# 这里都转成大写字母，便于识别基因名，或者根据你自己的需要调整
tempf <- toupper(tempf)
tempi <- as.data.frame(table(tempf));

# 判断是不是常用词，如果是常用词，就删掉
tempj <- unlist(lapply(toupper(common_words_new), function(x){tempoo = which(as.character(tempi[,1]) == x); if (length(tempoo) != 0) return(tempoo)}));
tempk <- tempi[-tempj,];

# 根据Approved.Symbol挑出基因名
templ <- as.character(HGNCdata$Approved.Symbol);
head(templ)
# 有些物种的gene symbol里有小写字母，需要先转成大写，就运行下面这行
#templ <- toupper(templ)
tempm <- unlist(lapply(templ,function(x){return(which(x == as.character(tempk$tempf)))}));
tempn <- tempk

# 按照基因名出现频率排序
tempn2<- tempn
tempo <- unlist(lapply(as.character(tempn2$tempf), function(x){return(which(x == templ))}));
Synonyms <- as.character(HGNCdata$Synonyms);
data_table <- cbind.data.frame(Symbol=as.character(tempn2$tempf), Synonyms, Freq=as.numeric(tempn2$Freq));
colnames(data_table) <- c("Gene_symbol", "Synonyms", "Freq");

# 把基因名的出现频率保存到文件
write.table(data_table, "output_friends.txt", sep = "\t", quote = F, row.names = F)
</code></pre>
<p><img src="https://tcapi.voiceclouds.cn:8444/uploads/2025/11/19/mi5el79en0h8qcvu88.png" alt="" /></p>
<h2>给基因添加GeneCards链接，保存成网页格式</h2>
<pre><code class="language-{r}"># 跟【模块二】一样的函数
createLink <- function(base,val) {
sprintf('<a href="%s" class="btn btn-link" target="_blank" >%s</a>',base,val)
}

# 根据GeneCards的网址规律，给基因名添加链接
res <- data.frame(symbols=createLink(paste0("https://www.genecards.org/cgi-bin/carddisp.pl?gene=",as.character(tempn2$tempf)),as.character(tempn2$tempf)),
               Synonyms,
               as.numeric(tempn2$Freq),
               stringsAsFactors = F)
res <- na.omit(res)
y <- DT::datatable(res,escape = F,rownames=F)

#保存到网页格式的文件
DT::saveWidget(y,"output_friends.html")
</code></pre>
<p><img src="https://tcapi.voiceclouds.cn:8444/uploads/2025/11/19/mi5emvbye32qb2kll2.png" alt="" /></p>
<h1>功能扩展</h1>
<p>【模块一】的文章数量可以借鉴时间数据的展示方式，画出酷炫的图；</p>
<p>【模块二】的结果可以添加影响因子；</p>
<p>【模块三】的结果可以画成词云。</p>
<h1>参考资料</h1>
<p>NCBI里的各种数据库，都可以用rentrez这个R包通过API访问，达到批量检索、处理的目的。</p>
<p>参考资料：https://bioconnector.github.io/workshops/r-ncbi.html</p>

生信喵 发表于 2025-11-19 16:07:11

<p>文中所用数据可以关注公众号“生信喵实验柴”<br />
<img src="https://tcapi.voiceclouds.cn:8444/uploads/2025/09/07/mf9nu72s90ofyseu9hl.png" alt="" /><br />
发送关键词“20250428”获取</p>

页: [1]

bioinfoer's Archiver

R语言挖掘Pubmed数据库