使用 R 从 NCBI 基因数据库中获取数据

Posted

技术标签:

【中文标题】使用 R 从 NCBI 基因数据库中获取数据【英文标题】：Obtaining data from NCBI gene database with R 【发布时间】：2022-01-15 22:31:01 【问题描述】：

Rentrez 包

根据this manual，我在 Linux（Ubuntu 20.04.2）的实验室计算机上发现了 RStudio（版本 1.1.442）中的rentrez 包。 但是，后来当我想在笔记本电脑上的 Windows 8 Pro (RStudio 2021.09.0) 上运行相同的代码时

library (rentrez)
entrez_dbs()
entrez_db_searchable("gene")
#res <- entrez_search (db = "gene", term = "(Vibrio[Organism] OR vibrio[All Fields]) AND (16s[All Fields]) AND (rna[All Fields]) AND (owensii[All Fields] OR navarrensis[All Fields])", retmax = 500, use_history = TRUE)

即使关闭会话或重新安装rentrez 包后，我也无法摆脱此错误

curl::curl_fetch_memory(url, handle = handle) 中的错误：schannel: 下一个 InitializeSecurityContext 失败：SEC_E_ILLEGAL_MESSAGE (0x80090326) - 当致命的 SSL/TLS 警报出现时，通常会发生此错误 收到（例如握手失败）。

这是我面临的主要问题。

RSelenium 包

后来我决定处理包含有关genes 及其sequences in FASTA format 修改我以前使用的代码的详细信息的页面。它使用rvest 和rselenium 包，结果很完美。

# Specifying a webpage

url <- "https://www.ncbi.nlm.nih.gov/gene/66940694"    # the last 9 numbers is gene id

library(rvest)
library(RSelenium)

# Opening a browser

driver <- rsDriver(browser = c("firefox"))
remDr <- driver[["client"]]
remDr$errorDetails
remDr$navigate(url)

# Clicked outside in an empty space next to the FASTA button and copied a full xPath (redirecting to a FASTA data containing webpage)

remDr$findElement(using = "xpath", value = '/html/body/div[1]/div[1]/form/div[1]/div[5]/div/div[6]/div[2]/div[3]/div/div/div[3]/div/p/a[2]')$clickElement()

webElem <- remDr$findElement("css", "body")

#scrolling to the end of a webpage: left it from the old code for the case of a long gene

for (i in 1:5)
  Sys.sleep(2)
  webElem$sendKeysToElement(list(key = "end"))

# Let's get gene FASTA, for example

page <- read_html(remDr$getPageSource()[[1]])
fasta <- page %>%
html_nodes('pre') %>%
  html_text()
print(fasta)

输出：">NZ_QKKR01000022.1:c3037-151 副霍乱弧菌菌株 2016V-1111 2016V-1111_ori_contig_18，全基因组霰弹枪 序列\nGGT...

该代码可以很好地获取有关基因的其他详细信息，例如其登录号、位置、生物体等。

循环处理多个基因 ID

后来我尝试更改代码以同时获取 the explanations 之后的几个基因 ID 的相同信息，我来到这里是为了我的另一个项目。

# Specifying a list of gene IDs

res_id <- c('57838769','61919208','66940694')
dt <- res_id    # <lapply> looping function refused to work if an argument had a different name rather than <dt>
  
driver <- rsDriver(browser = c("firefox"))  
remDr <- driver[["client"]]

## Writing a function of GET_FASTA dependent on GENE_ID (x)

get_fasta <- function(x)
  link = paste0('https://www.ncbi.nlm.nih.gov/gene/',x)
  remDr$navigate(link)
  remDr$findElement(using = "xpath", value = '/html/body/div[1]/div[1]/form/div[1]/div[5]/div/div[6]/div[2]/div[3]/div/div/div[3]/div/p/a[2]')$clickElement()

...下面有一个继续但是这里出现了错误，说找不到之前成功使用的同一个xPath。

错误：摘要：NoSuchElement 详细信息：找不到元素 在页面上使用给定的搜索参数。班级： org.openqa.selenium.NoSuchElementException 更多细节：运行 errorDetails 方法

我尝试在 xPath 末尾删除/a[2] 以获取/html/.../p，因为它在初始代码中工作，但稍后再次出现错误。

  webElem <- remDr$findElement("css", "body")

  for (i in 1:5)
    Sys.sleep(2)
    webElem$sendKeysToElement(list(key = "end"))
   

 # Addressing selectors of FASTA on the website
  
  fasta <- remDr$getPageSource()[[1]] %>% 
    read_html() %>%
    html_nodes('pre') %>%
    html_text()
  fasta
  return(fasta)


## Writing a function of GET_ACC_NUM dependent on GENE_ID (x)

get_acc_num <- function(x)
  link = paste0( 'https://www.ncbi.nlm.nih.gov/gene/', x)
  remDr$navigate(link)
  remDr$findElement(using = "xpath", value = '/html/body/div[1]/div[1]/form/div[1]/div[5]/div/div[6]/div[2]/div[3]/div/div/div[3]/div/p')$clickElement()
  webElem <- remDr$findElement("css", "body")

  for (i in 1:5)
    Sys.sleep(2)
    webElem$sendKeysToElement(list(key = "end"))
   
  
  # Addressing selectors of ACC_NUM on the website
  
  acc_num <- remDr$getPageSource()[[1]] %>% 
    read_html() %>%
    html_nodes('.itemid') %>%
    html_text() %>%
    str_sub(start= -17)
  acc_num
  return(acc_num)


## Collecting all FUNCTION into tibble

get_data_table <- function(x)

      # Extract the Basic information from the HTML
      fasta <- get_fasta(x)
      acc_num <- get_acc_num(x)

      # Combine into a tibble
      combined_data <- tibble( Acc_Number = acc_num,
                               FASTA = fasta)


## Running FUNCTION for all x

df <- lapply(dt, get_data_table)

head(df)

我也试过写代码

仅限rvest， 用for (i in res_id) 编写循环， 使用if () else 引入两个以/html/.../p/a[2] 或.../p 结尾的不同xPath

但结果更令人困惑。

我在从事此类任务的同时正在学习 R 编码，因此欢迎任何建议和批评。

【问题讨论】：

fasta 返回为空。如果存在这样的节点，它应该是.pre 而不是pre。您能否提供您尝试用于单基因的fasta 的输出。和错误w.r.t。 rentrez package 我看到很多人问类似的问题，也许可以在那里找到解决方案？包含r 标签，以便您可以覆盖更多人。 @NadPat 你是对的。完整的 xPath 应为 /html/.../p/a[2]，否则 fasta 为空。已经更正了上面的问题。我检查了 fasta 选择器，它是pre，没有点。在此处检查：(ncbi.nlm.nih.gov/nuccore/…) 具有一个基因 id 的代码的输出如下：[1] ">NZ_ADAF01000001.1:257558-260444 模拟弧菌 MB451 菌株 MB-451 Contig43，全基因组鸟枪法序列\nGGT.. .\nCCGTGAGGCTTAACCAT\n\n" 后来我用 xPath 进行了实验，它说没有这样的元素。 @NadPat 我检查了 curl 错误的解释，与 rentrez 包没有连接。但他们超越了我。通常，“rentrez”是专门为从 NCBI 数据库中获取信息而设计的，最初是用 Python 编写的。但是直接使用选择器会使工作更短，不包括在 .gb 格式的一组文件中使用正则表达式的步骤 【参考方案1】：

节点pre 无效。我们必须在 class 或 'id` 等内部寻找价值。

webElem$sendKeysToElement(list(key = "end")你不需要这个命令，因为没有必要滚动页面。

下面是获取基因序列的代码。

首先我们必须通过rvest获取基因序列的链接

library(rvest)
library(dplyr)
res_id <- c('57838769','61919208','66940694')

link = vector()
for(i in res_id)
  url = paste0('https://www.ncbi.nlm.nih.gov/gene/', i)
  df = url %>%
    read_html() %>% 
    html_node('.note-link') 
  
  link1 = xml_attrs(xml_child(df, 3))[["href"]]
  link1 = paste0('https://www.ncbi.nlm.nih.gov', link1)
  link = rbind(link, link1)


link1 "https://www.ncbi.nlm.nih.gov/nuccore/NZ_ADAF01000001.1?report=fasta&from=257558&to=260444"       
link1 "https://www.ncbi.nlm.nih.gov/nuccore/NZ_VARQ01000103.1?report=fasta&from=64&to=2616&strand=true" 
link1 "https://www.ncbi.nlm.nih.gov/nuccore/NZ_QKKR01000022.1?report=fasta&from=151&to=3037&strand=true"

在获得链接后，我们将得到我们通过RSelenium做的基因序列。我尝试使用rvest 进行操作，但无法获得序列。

启动浏览器

library(RSelenium)
driver = rsDriver(browser = c("firefox"))
remDr <- driver[["client"]]

获取序列的函数

get_seq = function(link)
  remDr$navigate(link)
  Sys.sleep(5)
  df = remDr$getPageSource()[[1]] %>% 
    read_html() %>% 
    html_nodes(xpath = '//*[@id="viewercontent1"]') %>% 
    html_text()
  return(df)


df = lapply(link, get_seq)

现在我们有列表 df 包含所有信息。

【讨论】：

考虑到我们工作的目标，代码在我们需要的一些调整下完美运行。代码在开始时出错，直到我安装了“xml2”包。关于rentrez这个包，我们电脑上更新了Ubuntu之后就开始工作了，所以之前的错误原因一直很神秘。

以上是关于使用 R 从 NCBI 基因数据库中获取数据的主要内容，如果未能解决你的问题，请参考以下文章

如何从NCBI下载基因组数据

如何在genbank查找某个细菌的基因序列？

NCBI批量下载基因组

人类基因组在三大基因数据库中的不同版本

基因序列碱基数测算

如何计算每个基因的覆盖度与深度