bioinfoer»论坛 › 生信基础 › 数据库 › R语言挖掘Pubmed数据库

0 收听	12 听众	400 主题

R语言挖掘Pubmed数据库

发表于 3 天前 | 查看: 19| 回复: 2

需求描述

通过API检索NCBI的Pubmed数据库，批量获得文章信息，整理出摘要里基因的词频。

应用场景

老板交给我的这个基因好陌生，用这套代码跑出表格，一眼望去，发现它的好朋友我都很熟，于是，熟悉的思路、方法都能用上了。
文章做到第三部分卡住了，我的基因到底调控了谁，或者谁调控了它，或者它跟谁是好朋友一同发挥作用？用这套代码跑出表格，一眼望去，发现目的基因跟我熟悉的一类分子一同出现在某篇文章的摘要里，我们实验室有非常完善的研究这类分子的实验体系，接下来就瞄准这个分子了。

只要设置好检索词，就可以批量检索并整理文献和摘要，让你的头脑风暴比别人快10倍。

包含三个模块，相对独立，可分别运行。全部运行完将让你对目的基因有个整体的认识。

【模块一】文章的增长趋势

【模块二】发表了哪些文章

【模块三】从摘要找好朋友

环境设置

使用国内镜像安装包

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

加载包

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

输入文件

常用词

common_words_new.v2.txt，常用词。【模块三】计算词频时，要去除这些常用词。

common_words_new <- read.table("common_words_new.v2.txt")
common_words_new <- common_words_new$V1

基因名

每行一个基因。包含基因名、别名、位点，【模块三】将从摘要中统计这些基因名。还可以换成你想统计词频的其他类型的词。

(load("HGNCdata.rda")) #位于当前文件夹
head(HGNCdata)

参数设置

在这里设置检索词、检索的文章发表时间范围

targetGene <- "PD-L1" #检索的基因名
year <- 2017:2025 #检索的文章的发表年份范围

这里用的检索词跟你在网页版Pubmed中用的检索词是一样的，初级用法就是用AND或OR连接你希望检索到的内容，例如基因名、疾病、处理、物种等。医学领域推荐直接用MeSH terme来限定，Medical Subject Headings(MeSH) refer to wiki https://en.wikipedia.org/wiki/Medical_Subject_Headings

如果用基因名作为检索词，建议把基因的别名都放进来，可以通过下面的代码提取：

targetGene.Symbol <- HGNCdata[HGNCdata$Approved.Symbol %like% targetGene,]
targetGene.Synonyms <- HGNCdata[HGNCdata$Synonyms %like% targetGene,]
targetGeneInfo <- rbind(targetGene.Symbol, targetGene.Synonyms)
targetGeneInfo

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

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

# 或者自己手写基因名作为检索词
#term <- "CD274 OR PD-L1 OR PDL1"

【模块一】文章的增长趋势

#查询NCBI里的数据库简称
#entrez_dbs()
#pubmed里有哪些可以查询的fields
#entrez_db_searchable(db = "pubmed")

#先写个函数，获取特定年份范围内发表的带有检索词的文章数量
search_year <- function(year, term){
    query <- paste(term, "AND (", year, "[PDAT])")
    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")

【模块二】发表了哪些文章

检索并从结果中提取出PMID、发表日期、文章名、期刊。

#我们先把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[i:(i+n-1)])
  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,'\"'))[1]),stringsAsFactors = F)[,1]
tmp <- data.frame(sapply(tmp, function(x) unlist(strsplit(x,'"),'))[1]),stringsAsFactors = F)[,1]
tmp <- data.frame(sapply(tmp, function(x) unlist(strsplit(x,'",'))[1]),stringsAsFactors = F)[,1]
res$authors <- tmp

#把结果保存到文本文件
write.table(res, "output_paper.txt", quote = F, sep = "\t", row.names = F)

添加文章的Pubmed链接，保存成网页格式

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

收藏0 回复显示全部楼层道具举报

生信喵

发表于 3 天前

文中所用数据可以关注公众号“生信喵实验柴”

发送关键词“20250428”获取

回复显示全部楼层道具举报

生信喵

发表于 3 天前

【模块三】从摘要找好朋友

检索

这段检索跟前面“【模块二】发表了哪些文章”的开头一样，就不重复运行了。

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)

把xml文件整理成数据框格式

xml中包含title、authors、year、journal、key_words、doi、pmid、abstract等信息，整理成data.frame格式。

n <- 200 #每次读入的记录数量
abstract0 <- c()
for (i in seq(1, length(result$ids), n)) {
  rec <- parse_pubmed_xml(entrez_fetch(db = "pubmed", id = result$ids[i : (i + n - 1)], rettype = "xml"))
  abstract1 <- as.data.frame(do.call('rbind', rec))
  abstract0 <- rbind(abstract0, abstract1)
}
abstract1[1,]

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

#查看其中某一篇文章的摘要
#abstract[30]

提取abstract中的基因名，并统计出现频率。

先把摘要分割成单个的单词，然后根据输入文件common_words_new.v2.txt去掉常用词，再根据HGNCdata.rda挑出基因名。

# 先用标点符号分割出单词
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[tempm,]

# 按照基因名出现频率排序
tempn2<- tempn[order(as.numeric(tempn$Freq), decreasing = T),]
tempo <- unlist(lapply(as.character(tempn2$tempf), function(x){return(which(x == templ))}));
Synonyms <- as.character(HGNCdata$Synonyms[tempo]);
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)

给基因添加GeneCards链接，保存成网页格式

# 跟【模块二】一样的函数
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")

功能扩展

【模块一】的文章数量可以借鉴时间数据的展示方式，画出酷炫的图；

【模块二】的结果可以添加影响因子；

【模块三】的结果可以画成词云。

参考资料

NCBI里的各种数据库，都可以用rentrez这个R包通过API访问，达到批量检索、处理的目的。

参考资料：https://bioconnector.github.io/workshops/r-ncbi.html

回复显示全部楼层道具举报

返回列表

		自动登录	找回密码
密码			立即注册