Search
Question: How to prevent readDNAStringSet() function from switching masked sequences (lowercase) to uppercase sequence?
0
gravatar for beausoleilmo
13 days ago by
beausoleilmo0 wrote:

I'm using the package Biostrings

 

library("Biostrings")

DNAString("TATCAAATACTCAAGCACtaaggaaacaggaaaatct")

will return

37-letter "DNAString" instance seq: TATCAAATACTCAAGCACTAAGGAAACAGGAAAATCT

Why is aaggaaacaggaaaatct not staying in lowercase?

Is there a way to prevent the transformation to uppercase? I need to keep the lowercase letters. 

ADD COMMENTlink modified 12 days ago by James W. MacDonald45k • written 13 days ago by beausoleilmo0
0
gravatar for James W. MacDonald
12 days ago by
United States
James W. MacDonald45k wrote:

The simple answer is 'Because that is how it is designed to work'. From ?DNAString, which you should always check first:

Details:

     The DNAString class is a direct XString subclass (with no
     additional slot).  Therefore all functions and methods described
     in the XString man page also work with a DNAString object
     (inheritance).

     Unlike the BString container that allows storage of any single
     string (based on a single-byte character set) the DNAString
     container can only store a string based on the DNA alphabet (see
     below).  In addition, the letters stored in a DNAString object are
     encoded in a way that optimizes fast search algorithms.

And

> DNA_ALPHABET
 [1] "A" "C" "G" "T" "M" "R" "W" "S" "Y" "K" "V" "H" "D" "B" "N" "-" "+" "."

Since there are no lower case letters in the DNA alphabet, you can't have them in a DNAString. But do note that:

> BString("TATCAAATACTCAAGCACtaaggaaacaggaaaatct")
  37-letter "BString" instance
seq: TATCAAATACTCAAGCACtaaggaaacaggaaaatct

Where the BString class inherits from XString, just like DNAString, so it should probably suffice?

ADD COMMENTlink written 12 days ago by James W. MacDonald45k

The thing is that I want to read a big fasta file, not just individual sequences. 

readDNAStringSet(filepath = "~/Desktop/UCSC/my.genome.fa.gz")

But this one is not taking the masked sequences... is there a readDNAStringSet that can keep the masked sequences?

ADD REPLYlink written 12 days ago by beausoleilmo0

The first place you should be looking is the help page for the function you are trying to use! I already tried to nudge you in that direction, so let me be more direct. Reading the help page first will save you lots of time, and that is what you should always do. If it seems to people who might help that you are using the support site as a stand in for doing your own homework, you will quickly find any help you might get dwindling down to nothing.

So do you see anything in the help page for readDNAStringSet that might be applicable to your situation, particularly in relation to what I have already told you about BStrings?
 

ADD REPLYlink written 12 days ago by James W. MacDonald45k

I actually look at help pages. I haven't found anything related to that. I also found somebody else, on another forum, asking for the same question and they said to ask the question here. Now, it seems that there are no easy way to do want I'm trying to aim. 

This is my solution:

I used the raw genome that is masked, and then I needed to transform it in a file that was readable by R. 

# Create one big line
sed -e ':a' -e 'N' -e '$!ba' -e 's/\n//g' ~/Desktop/UCSC/my.fasta1.fa > ~/Desktop/UCSC/my.fasta2.fa
# brew install gnu-sed # Because I'm on a mac 
gsed  's/>[A-Z]\{2\}[0-9]\{6\}/\n&\n/g;s/^\n//' ~/Desktop/UCSC/my.fasta2.fa > ~/Desktop/UCSC/my.fasta3.fa

Then from this, I imported my file in R and did whatever I wanted. I actually merged that file with the one from the readDNAStringSet function to benefit from both worlds. Thanks for your help. 

# I actually found a way to print the masked sequences, so here they are... 
load("~/Desktop/UCSC/whateverIdidtothefastafileinbetween.RData",
     verbose = TRUE)

df = merge(df,
      ucsc.genome.masked, 
      by.x = "seq_name", 
      by.y = "scaffolds.names")

If you have another solution, I'd be happy to see it. 

 

 

ADD REPLYlink modified 9 days ago • written 9 days ago by beausoleilmo0
Please log in to add an answer.

Help
Access

Use of this site constitutes acceptance of our User Agreement and Privacy Policy.
Powered by Biostar version 2.2.0
Traffic: 205 users visited in the last hour