1

我正在尝试使用 R 包“rentrez”从 GenBank 中提取一些信息（元数据）以及我在此处找到的示例https://ajrominger.github.io/2018/05/21/gettingDNA.html. 具体来说，对于一组特定的生物，我搜索所有具有地理坐标的记录，然后想要提取有关入藏号、分类单元、测序基因座、国家、纬度和收集日期的数据。作为输出，我想要一个 csv 文件，其中每条记录的数据位于单独的行中。下面的代码似乎可以完成这项工作，但在某些时候，行会被来自不同记录的数据与相邻行重叠。例如，rentrez 从文件中的 NCBI 109 记录中检索到的 157 条记录看起来像我想要实现的，但其余的却一团糟。我将非常感谢有关如何解决此问题的任何建议，因为我是 R 的新手，并且弄清楚每个步骤都需要花费大量时间。

setwd ("C:/R-Works")
library('XML')
library('rentrez')
argasid <- entrez_search(db="nuccore", term = "Argasidae[Organism] AND [lat]", use_history=TRUE, retmax=15000)
x <- entrez_fetch (db="nuccore", id=argasid$ids, rettype= "native", retmode="xml", parse=TRUE)
x <-xmlToList(x)

cleanEntrez <- function(x) {
  basePath <- 'Seq-entry_seq.Bioseq'
  c(
    genbank = as.character(x[paste(basePath, 
                                   'Bioseq_id', 'Seq-id', 'Seq-id_genbank', 
                                   'Textseq-id', 'Textseq-id_accession', 
                                   sep = '.')]),
    taxon = as.character(x[paste(basePath,
                                 'Bioseq_descr', 'Seq-descr', 'Seqdesc', 
                                 'Seqdesc_source', 'BioSource', 'BioSource_org', 
                                 'Org-ref', 'Org-ref_taxname', 
                                 sep = '.')]),
    bseqdesc_title = as.character(x[paste(basePath,
                                          'Bioseq_descr', 'Seq-descr', 'Seqdesc', 
                                          'Seqdesc_title', 
                                          sep = '.')]),
      lat_lon = as.character(x[grep('lat-lon', x) + 1]),
      geo_description = as.character(x[grep('country', x) + 1]),
      coll_date = as.character(x[grep('collection-date', x) + 1])
  )
}
getGenbankMeta  <- function(ids) {
  allRec <- entrez_fetch(db = 'nuccore', id = ids, 
                         rettype = 'native', retmode = 'xml', 
                         parsed = TRUE)
  allRec <- xmlToList(allRec)[[1]]
  
  o <- lapply(allRec, function(x) {
    cleanEntrez(unlist(x))
  })
  
  temp <- array(unlist(o), dim = c(length(o[[1]]), length(ids)))
  seqVec <- temp[nrow(temp), ]
  seqDF <- as.data.frame(t(temp[-nrow(temp), ]))
  names(seqDF) <- names(o[[1]])[-nrow(temp)]
  
  return(list(seq = seqVec, data = seqDF))
} 
write.csv(getGenbankMeta(argasid$ids), 'argasid_georef.csv')
4

0 回答 0