6

I am trying to convert DNAbin files to fasta format. Rationale: I want to use the fasta file to calculate non-synonymous/synonymous mutation rate.

my_dnabin1 is a DNAbin file of 55 samples and I am using the following code to convert it into a fasta file.

dnabin_to_fasta <- lapply(my_dnabin1, function(x) as.character(x[1:length(x)])) This generates a list of 55 samples which looks like:

$SS.11.01

[1] "t" "t" "a" "c" "c" "t" "a" "a" "a" "a" "a" "g" "c" "c" "g" "c" "t" "t" "c" "c" "c" "t" "c" "c" "a" "a" [27] "c" "c" "c" "t" "a" "g" "a" "a" "g" "c" "a" "a" "a" "c" "c" "t" "t" "t" "c" "a" "a" "c" "c" "c" "c" "a"

$SS.11.02

[1] "t" "t" "a" "c" "c" "t" "a" "a" "a" "a" "a" "g" "c" "c" "g" "c" "t" "t" "c" "c" "c" "t" "c" "c" "a" "a" [27] "c" "c" "c" "t" "a" "g" "a" "a" "g" "c" "a" "a" "a" "c" "c" "t" "t" "t" "c" "a" "a" "c" "c" "c" "c" "a"

and so on...

However, I want a fasta formatted file as the output that may look something like:

>SS.11.01 ttacctga

>SS.11.02 ttacctga

Nikita
  • 63
  • 4

1 Answers1

5

The following should work:

cat(file="out.fa", paste(paste0(">",names(dnabin_to_fasta)),
          sapply(dnabin_to_fasta, paste, collapse=""), sep="\n"), sep="\n");

A bit of explanation, starting from the middle outwards:

  1. Collapse the letter arrays inside dnabin_to_fasta into single strings
  2. Append ">" to the front of the names of each sequence
  3. Concatenate the results of steps 1 & 2 together, separating with a line break
  4. Write each concatenated string to a file, separating by a line break
gringer
  • 14,012
  • 5
  • 23
  • 79
  • Thank you so much! It solved a major part of my problem. When I print it to the screen, I get the exact output I want. However, it does not produce any output files. Could you please help me with this, too? – Nikita Aug 30 '17 at 14:08
  • Note that R will save files in the "working directory", which is not necessarily the directory that the last file was loaded from. If you're using cat with the file option, it's probably saving somewhere. You can find out where by using the getwd() command. – gringer Aug 30 '17 at 14:18
  • Thank you so much! It solved my problem! I was stuck in it since a long time. – Nikita Aug 30 '17 at 14:19