Question: How to prevent readDNAStringSet() function from switching masked sequences (lowercase) to uppercase sequence?
0
23 months 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.

modified 23 months ago by James W. MacDonald51k • written 23 months ago by beausoleilmo0
0
23 months ago by
United States
James W. MacDonald51k 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?

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?

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?

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.