R/GC.R

In seqinr: Biological Sequences Retrieval and Analysis

################################
# G+C content
#################################

GC <- function(seq, forceToLower = TRUE, exact = FALSE, NA.GC = NA, oldGC = FALSE,  alphabet = s2c("acgtswmkryvhdb"))
{
    
    
    
    #
    # NA propagation:
    #
    if(length(seq) == 1 && is.na(seq)) return(NA)
    #
    # Check that sequence is a vector of chars:
    #
    if(nchar(seq[1]) > 1) stop("sequence is not a vector of chars")
    #
    # Force to lower-case letters if requested:
    #
    if(forceToLower) seq <- tolower(seq)
    
    seq <- seq[seq %in% alphabet]
    
    #
    # Compute the count of each base:
    #
    nc <- sum( seq == "c" )
    ng <- sum( seq == "g" )
    na <- sum( seq == "a" )
    nt <- sum( seq == "t" )
    
    #
    # oldGC case:
    #
    if(oldGC){
        warning("argument oldGC is deprecated")
        return( (nc + ng)/length(seq) )
    }
    
    #
    # General case:
    #
    if(! exact){
        if(na + nc + ng + nt == 0){
            result <- NA.GC
        } else {
            result <- (nc + ng)/(na + nc + ng + nt)
        }
    } else {
        #
        # We have our first estimate of GC vs. AT base counts:
        #
        ngc <- ng + nc
        nat <- na + nt
        
        #
        # weak and strong bases are 100% informative with respect
        # to the GC content, we just add them:
        #
        # s : Strong (g or c)
        # w : Weak (a or t)
        #
        ngc <- ngc + sum( seq == "s" )
        nat <- nat + sum( seq == "w" )
        
        ##########################
        # Ambiguous base section #
        ##########################
        
        #
        # m : Amino (a or c)
        #
        if(na + nc != 0){
            nm <- sum( seq == "m")
            ngc <- ngc + nm*nc/(na + nc)
            nat <- nat + nm*na/(na + nc)
        }
        
        #
        # k : Keto (g or t)
        #
        if(ng + nt != 0){
            nk <- sum( seq == "k" )
            ngc <- ngc + nk*ng/(ng + nt)
            nat <- nat + nk*nt/(ng + nt)
        }
        
        #
        # r : Purine (a or g)
        #
        if(ng + na != 0){
            nr <- sum( seq == "r" )
            ngc <- ngc + nr*ng/(ng + na)
            nat <- nat + nr*na/(ng + na)
        }
        
        #
        # y : Pyrimidine (c or t)
        #
        if(nc + nt != 0){
            ny <- sum( seq == "y" )
            ngc <- ngc + ny*nc/(nc + nt)
            nat <- nat + ny*nt/(nc + nt)
        }
        
        #
        # v : not t (a, c or g)
        #
        if(na + nc + ng != 0){
            nv <- sum( seq == "v" )
            ngc <- ngc + nv*(nc + ng)/(na + nc + ng)
            nat <- nat + nv*na/(na + nc + ng)
        }
        #
        # h : not g (a, c or t)
        #
        if(na + nc + nt != 0){
            nh <- sum( seq == "h" )
            ngc <- ngc + nh*nc/(na + nc + nt)
            nat <- nat + nh*(na + nt)/(na + nc + nt)
        }
        #
        # d : not c (a, g or t)
        #
        if(na + ng + nt != 0){
            nd <- sum( seq == "d" )
            ngc <- ngc + nd*ng/(na + ng + nt)
            nat <- nat + nd*(na + nt)/(na + ng + nt)
        }
        #
        # b : not a (c, g or t)
        #
        if(nc + ng + nt != 0){
            nb <- sum( seq == "b" )
            ngc <- ngc + nb*(nc + ng)/(nc + ng + nt)
            nat <- nat + nb*nt/(nc + ng + nt)
        }
        #
        # n : any (a, c, g or t) is not informative, so
        # we compute the G+C content as:
        #
        if( ngc + nat == 0){
            result <- NA.GC
        } else {
            result <- ngc/(ngc + nat)
        }
    }
    return(result)
}

######################
# GCpos		     #
######################

GCpos <- function(seq, pos, frame = 0, ...){
    if(nchar(seq[1]) > 1){
        warning("sequence is not a vector of chars, I'm trying to cast it into one")
        seq <- s2c(seq[1])
    }
    #
    # Take frame into account:
    #
    if(frame != 0) seq <- seq[(1 + frame):length(seq)]
    #
    # Return result:
    #
    GC(seq[seq(pos, length(seq), by = 3)], ...)
}

######################
# GC1		     #
######################

GC1 <- function(seq, frame = 0, ...) GCpos(seq = seq, pos = 1, frame = frame, ...)

######################
# GC2		     #
######################

GC2 <- function(seq, frame = 0, ...) GCpos(seq = seq, pos = 2, frame = frame, ...)

######################
# GC3		     #
######################

GC3 <- function(seq, frame = 0, ...) GCpos(seq = seq, pos = 3, frame = frame, ...)