(v)matchPattern DNAStringSetList of Codons to Reference DNAString

user_15077940 0
Viewed 96

I am assessing the impact of hotspot single nucleotide polymorphism (SNPs) from a next generation sequencing (NGS) experiment on the protein sequence of a virus. I have the reference DNA sequence and a list of hotspots. I need to first figure out the reading frame of where these hotspots are seen. To do this, I generated a DNAStringSetList with all human codons and want to use a vmatchpattern or matchpattern from the Biostrings package to figure out where the hotspots land in the codon reading frame.

I often struggle with lapply and other apply functions, so I tend to utilize for loops instead. I am trying to improve in this area, so welcome a apply solution should one be available.

Here is the code for the list of codons:

alanine <- DNAStringSet("GCN")
arginine <- DNAStringSet(c("CGN", "AGR", "CGY", "MGR"))  
asparginine <- DNAStringSet("AAY") 
aspartic_acid <- DNAStringSet("GAY")
asparagine_or_aspartic_acid <- DNAStringSet("RAY")
cysteine <- DNAStringSet("TGY")
glutamine <- DNAStringSet("CAR")
glutamic_acid <- DNAStringSet("GAR")
glutamine_or_glutamic_acid <- DNAStringSet("SAR")
glycine <- DNAStringSet("GGN")
histidine <- DNAStringSet("CAY")
start <- DNAStringSet("ATG")
isoleucine <- DNAStringSet("ATH")
leucine <- DNAStringSet(c("CTN", "TTR", "CTY", "YTR"))
lysine <- DNAStringSet("AAR") 
methionine <- DNAStringSet("ATG") 
phenylalanine <- DNAStringSet("TTY") 
proline <- DNAStringSet("CCN")
serine <- DNAStringSet(c("TCN", "AGY"))
threonine <- DNAStringSet("ACN")
tyrosine <- DNAStringSet("TGG")
tryptophan <- DNAStringSet("TAY")
valine <- DNAStringSet("GTN")
stop <- DNAStringSet(c("TRA", "TAR"))

codons <- DNAStringSetList(list(alanine, arginine, asparginine, aspartic_acid, asparagine_or_aspartic_acid,
                           cysteine, glutamine, glutamic_acid, glutamine_or_glutamic_acid, glycine,
                           histidine, start, isoleucine, leucine, lysine, methionine, phenylalanine,
                           proline, serine, threonine, tyrosine, tryptophan, valine, stop))

Current for loop code:

reference_stringset <-  DNAStringSet(covid)

codon_locations <- list()

for (i in 1:length(codons)) {
  pattern <- codons[[i]]
  codon_locations[i] <- vmatchPattern(pattern, reference_stringset)
}

Current error code. I am filtering the codon DNAStringSetList so that it is a DNAStringSet.

Error in normargPattern(pattern, subject) : 'pattern' must be a single string or an XString object

I can't give out the exact nucleotide sequence, but here is the COVID genome (link: https://www.ncbi.nlm.nih.gov/nuccore/NC_045512.2?report=fasta) to use as a reprex:

#for those not used to using .fasta files, first copy and past genome into notepad and save as a .fasta file
#use readDNAStringSet from Biostrings package to read in the .fasta file
filepath = #insert file path
covid <- readDNAStringSet(filepath)
bioconductor r
1 Answers

For the current code, change the way the codons is formed. Currently the output of codons looks like this:

DNAStringSetList of length 24
[[1]] GCN
[[2]] CGN AGR CGY MGR
[[3]] AAY
[[4]] GAY
[[5]] RAY
[[6]] TGY
[[7]] CAR
[[8]] GAR
[[9]] SAR
[[10]] GGN
...
<14 more elements>

Change it from DNAStringSetList to a conglomerate DNAStringSet of the amino acids.

codons <- DNAStringSet(c(alanine, arginine, asparginine, aspartic_acid, asparagine_or_aspartic_acid,
                                cysteine, glutamine, glutamic_acid, glutamine_or_glutamic_acid, glycine,
                                histidine, start, isoleucine, leucine, lysine, methionine, phenylalanine,
                                proline, serine, threonine, tyrosine, tryptophan, valine, stop))

codons
DNAStringSet object of length 32:
     width seq
 [1]     3 GCN
 [2]     3 CGN
 [3]     3 AGR
 [4]     3 CGY
 [5]     3 MGR
 ...   ... ...
[28]     3 TGG
[29]     3 TAY
[30]     3 GTN
[31]     3 TRA
[32]     3 TAR

When I run the script I get the following output with the SARS-CoV-2 isolate listed for the example (I'm showing a small slice)

codon_locations[27:28]
[[1]]
MIndex object of length 1
$`NC_045512.2 Severe acute respiratory syndrome coronavirus 2 isolate Wuhan-Hu-1, complete genome`
IRanges object with 0 ranges and 0 metadata columns:
       start       end     width
   <integer> <integer> <integer>


[[2]]
MIndex object of length 1
$`NC_045512.2 Severe acute respiratory syndrome coronavirus 2 isolate Wuhan-Hu-1, complete genome`
IRanges object with 554 ranges and 0 metadata columns:
            start       end     width
        <integer> <integer> <integer>
    [1]        89        91         3
    [2]       267       269         3
    [3]       283       285         3
    [4]       352       354         3
    [5]       358       360         3
    ...       ...       ...       ...
  [550]     29261     29263         3
  [551]     29289     29291         3
  [552]     29472     29474         3
  [553]     29559     29561         3
  [554]     29793     29795         3

Looking at the ones that had an output, only those with the standard nucleotides ("ATCG", no wobbles) found matches. Those will need to be changed as well to search.

If you're on twitter, I suggest linking the question using the #rstats, #bioconductor, and #bioinformatics hashtags to generate some more traction, I've noticed that bioinformatic specific questions on SO don't generate as much buzz.

user_137098260
Related