1
$\begingroup$

I wish to download a fasta sequence from NCBI using "rentrez" package and then save it in an object, then I want to have this object as a DNAStringSet object to work on the sequence using "Biostrings" capabilities.

I leave my incorrect code here so I would be grateful if you could help me out.

cov2 <- entrez_fetch(db="nucleotide" , rettype = "fasta", id = "NC_045512")
cov2 <- readDNAStringSet(cov2)
...
Error in .Call2("new_input_filexp", filepath, PACKAGE = "XVector") : 
  cannot open file '>NC_045512.2 
$\endgroup$
2
  • $\begingroup$ Please clarify your specific problem or provide additional details to highlight exactly what you need. As it's currently written, it's hard to tell exactly what you're asking. $\endgroup$
    – Community Bot
    Commented Oct 14, 2023 at 11:11
  • 1
    $\begingroup$ Check what the value of cov2 is after the first command. entrez_fetch() returns the fasta record, but readDNAStringSet() requires a file path as input. $\endgroup$
    – Cloudberry
    Commented Oct 15, 2023 at 15:15

2 Answers 2

1
$\begingroup$

Your issue is that entrez_fetch() returns the fasta record(s) as a single string with the definition lines and sequences concatenated. Also, newline characters (\n) appear in place of line breaks, with double newlines (\n\n) separating different records. Here's a function that creates a DNAStringSet given a list of nucleotide identifiers without writing each fasta to a file:

entrez2dss <-
  function(id_list) {
    
    require(rentrez)    # v1.2.3
    require(Biostrings) # v2.66.0
    require(stringr)    # v1.5.0
    
    # fetch all sequences, trim empty elements
    raw_ <-
      entrez_fetch(db = "nucleotide",
                   rettype = "fasta",
                   id = id_list) |>
      str_split_1(pattern = "\n\n") |>
      head(-1)
    
    # get sequence names
    names_ <- 
      gsub(pattern = "\n.*",
           replacement = "",
           x = raw_) |>
      gsub(pattern = "^>",
           replacement  = "")
    
    # get sequences
    seqs_ <-
      sub(pattern = ".*?\n",
          replacement = "",
          x = raw_) |>
      gsub(pattern = "\n",
           replacement = "")
    
    # create named DNAStringSet object
    dss <- DNAStringSet(x = seqs_)
    names(dss) <- names_
    dss
    
  }

seqlist <- c("NC_045512", "OR558592")
cov2 <- entrez2dss(seqlist)
$\endgroup$
1
$\begingroup$

The previous error message was because the entrez_fetch() function returns the sequence as a character string, not a file path. To save the sequence to a file, I can use the writeLines() function. For example, the following code will save the sequence to a file named cov2.fasta:

# Download the FASTA sequence
cov2 <- entrez_fetch(db = "nucleotide", rettype = "fasta", id = c("NC_045512", "OR558592"))
# Save the sequence to a file
writeLines(cov2, "cov2.fasta")
# Read the FASTA sequence into a DNAStringSet object
cov2_dna <- readDNAStringSet("cov2.fasta")
$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.