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

I'm using the package Biostrings




will return


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 2.1 years ago by James W. MacDonald52k • written 2.1 years ago by beausoleilmo0
Answer: How to prevent readDNAStringSet() function from switching masked sequences (lowe
gravatar for James W. MacDonald
2.1 years ago by
United States
James W. MacDonald52k wrote:

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


     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

     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.


 [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 2.1 years ago by James W. MacDonald52k

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 2.1 years 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 2.1 years ago by James W. MacDonald52k

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... 
     verbose = TRUE)

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

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



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


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