Quickest way to change DNA character string based on quality information (both imported from BAM file) in R

user_15469461 0
Viewed 56

I am attempting to alter the characters in a DNA character string based on the associated quality information. There are a number of tools I have seen for quality filtering (cutadapt, trimmomatic among others), but none of them seem to do prescisely what I want.

I want to convert the basecall to another value (using "Y" in below script for no particular reason) based on a quality threshold. The purpose of this is to exclude that base from downstream analysis, which depends on the presence of A / T / G / C / -.

The quality information is imported from the BAM file (using GenomicAlignments + ScanBam) and I am using Rsubread to convert this to Q-scores.

I have the below for loop working fine, but it takes quite a while on large data sets (~100k reads). Are there ways to maintain a similar output but in a more efficient manner, maybe with apply? I don't know how to turn this into a parallel computation (it seems to only use one core).

Dummy data included for completeness.

#Make dummy sequence data frame
seq1 <- "CCGGGGGATCCGGAGGCACCGGCGGCGGGTCCGGCGGCGGCAGTAATACGACTCACTATAGGGGCACCGGTGGACTGTTCGAGGCCATCGAGGGATTCATCGAAAACGGATGGGAAGGCATGATCGACGGCTGGTACGGCTTTAGGTGACACCAGAACGCCCAGGGCGAGGGCACGGCCGCTTAATAGGCGGCCGCGACTCTAGATCATAATCAGCCATACCACAT"
seq1DF <- data.frame(seq1)
seq2 <- "CCGGGGGATCCGGAGGCACCGGCGGCGGGTCCGGCGGCGGCAGTAATACGACTCACTATAGGGGCACCGGTGGACTGTTCGAGGCCATCGAGGGATTCATCGAAAACGGATGGGAAGGCATGATCGACGGCTGGTACGGCTTTAGGTGACACCAGAACGCCCAGGGCGAGGGCACGGCCGCTTAATAGGCGGCCGCGACTCTAGATCATAATCAGCCATACCACAT"
seq2DF <- data.frame(seq2)

seqDF <- as.data.frame(rbind(seq1,seq2))

#Make dummy quality data frame
quality1 <- data.frame(t(c(replicate(220, 40),30:25)))
quality2 <- data.frame(t(c(replicate(220, 40),35:30)))

qualityDF <- as.data.frame(rbind(quality1,quality2))

#Define quality threshold
threshold <- 30

#number of rows in seqDF changes depending on the number of reads
#some reads have slightly longer lengths, so check for the max nchar and go over all sequences to this extent
#the quality information can contain NA where there is no sequencing data, include !is.na to circumvent this

for (x in c(1:nrow(seqDF))){
  for (n in c(1:max(nchar(seqDF[,1])))) {
    if(qualityDF[x,n] <= threshold && !is.na(qualityDF[x,n])){
      substr(seqDF[x,1],n , n) <- "Y"
    }
  }
}
bioconductor r performance
2 Answers

This should be more efficient. It's unnecessary to loop through the whole string, only the targeted position ( <= threshould ) needed to be replaced by "Y".

# seq to data.frame split as column
df1 <- data.frame(seq1 = strsplit(seq1,"")[[1]])
df2 <- data.frame(seq2 = strsplit(seq2,"")[[1]])

# NA will be skipped 
df1[which(quality1 <= threshold),] <- "Y"
df2[which(quality2 <= threshold),] <- "Y"

# collapse from vector to string
seq1_convert = paste(df1$seq1, sep = '', collapse = '')
seq2_convert = paste(df2$seq2, sep = '', collapse = '')
user_75227640

Maybe the following is a little bit faster.

MNC <- max(nchar(seqDF[,1]))
for(i in c(1:nrow(seqDF))) {
  for(j in which(!is.na(qualityDF[i,]) & qualityDF[i,] <= threshold)) {
    substr(seqDF[i,1], j, j) <- "Y"
  }
  substr(seqDF[i,1], nchar(seqDF[i,1])+1, MNC)
}

Or you convert the data.frames to matrix.

#Convert to matrix
seqM <- strsplit(seqDF[,1], "", TRUE)
seqM <- simplify2array(lapply(tt, "[", seq_len(max(lengths(tt)))))
qualityM <- t(as.matrix(qualityDF))

#overwrite with `Y` where quality is below threshold
seqM[!is.na(qualityM) & qualityM <= threshold] <- "Y"

#In case needed - Convert it back to string
apply(seqM, 2, function(x) paste(x[!is.na(x)], collapse=""))
user_104885040
Related