R/rearranged.oriloc.R

In seqinr: Biological Sequences Retrieval and Analysis

#############################################################################
#                        Rearranged oriloc                                  #  
#                                                                           #
# Detection of replication-associated effects on base composition asymmetry #
#                                                                           #
#############################################################################

rearranged.oriloc <- function(
    seq.fasta = system.file("sequences/ct.fasta.gz", package = "seqinr") ,
    g2.coord = system.file("sequences/ct.predict", package = "seqinr") )
  {
    
    seq.fasta <- read.fasta(seq.fasta)[[1]]
    g2 <- readLines(g2.coord)
    g2 <- lapply(g2, function(x) unlist(strsplit(x, split=" ")))
    g2 <- lapply(g2, function(x) x[which(x!="")])
        
    start <-as.numeric(unlist(lapply(g2, function(x) x[2]))) 
    end <- as.numeric(unlist(lapply(g2, function(x) x[3])))
    strand <- rep("forward",length(start))
    strand[which(end<start)] <- "reverse"


    meancoord <- (start+end)/2

    ncds <- length(start)

    gcskew <- numeric(ncds)
    atskew <- numeric(ncds)
    
    for(i in seq_len(ncds)){

      cds=seq.fasta[start[i]:end[i]]
      
      g3=sum(cds[seq(from=3,to=length(cds),by=3)]%in%c("g","G"))
      c3=sum(cds[seq(from=3,to=length(cds),by=3)]%in%c("c","C"))
      t3=sum(cds[seq(from=3,to=length(cds),by=3)]%in%c("t","T"))
      a3=sum(cds[seq(from=3,to=length(cds),by=3)]%in%c("a","A"))

      gcskew[i]=(g3-c3)/(g3+c3)
      atskew[i]=(a3-t3)/(a3+t3)
      if(is.na(gcskew[i])){
        gcskew[i]=0
      }
      if(is.na(atskew[i])){
        atskew[i]=0
      }
      
    }

    neworder=c(which(strand=="forward"),which(strand=="reverse"))

    meancoord.rear <- seq_len(ncds)
    gcskew.rear=gcskew[neworder]
    atskew.rear=atskew[neworder]
    strand.rear=strand[neworder]
    meancoord.real=meancoord[neworder]
    data=data.frame(meancoord.rear,gcskew.rear,atskew.rear,strand.rear,neworder,meancoord.real)

    colnames(data)=c("meancoord.rearr","gcskew.rearr","atskew.rearr","strand.rearr","order","meancoord.real")

    return(data)
    
    
  }