R/gbk2g2.R

In seqinr: Biological Sequences Retrieval and Analysis

#########################################################################
#
#                                 gbk2g2
#
# Conversion of a genbank format into a run-glimmer2 like format
#
#
#########################################################################

gbk2g2 <- function(
        gbkfile =  "https://pbil.univ-lyon1.fr/datasets/seqinr/data/ct.gbk",
        g2.coord = "g2.coord")
{
    input <- readLines(gbkfile)
    
    outfile = file( description = g2.coord, open ="w")
    #
    # Keep lines with CDS flag:
    #
    input <- input[ substring(input,1,8) == "     CDS"]
    
    #
    # Extract boudaries strings
    #
    get.boundaries <- function( line )
    {
        line <- unlist(strsplit(line, split = " "))
        line <- line[nchar(line) > 1]
        return(line[2])
    }
    input <- sapply(input, get.boundaries)
    names(input) <- NULL
    
    #
    # Look for 5' partial genes:
    #
    idx5p <- grep(">", input)
    if( length( idx5p != 0 ) )
        warning("5' partial genes encountered (no output):", idx5p)
    
    #
    # Look for 3' partial genes:
    #
    idx3p <- grep("<", input)
    if( length( idx3p != 0 ) )
        warning("3' partial genes encountered (no output):", idx3p)
    
    #
    # Look for join in features:
    #
    idxjoin <- grep("join", input)
    if( length( idxjoin != 0 ) )
        warning("join encountered (no output):", idxjoin)
    
    #
    # Define partials and join:
    #
    censored <-  c(idx5p, idx3p, idxjoin)
    
    #
    # Extract boundaries:
    #
    for( i in seq(from = 1, to = length(input), by = 1 ) )
    {
        if( i %in% censored )
            next
        
        tmp <- unlist( strsplit(input[i], split="\\.\\.") )
        if( length( grep("complement", input[i]) ) == 1 )
        {
            end <- as.integer( substring(tmp[1], first = 12 ) ) + as.integer(3)
            start <- as.integer( substring(tmp[2], first = 1, last = nchar(tmp[2]) - 1 ))
            line <- sprintf(fmt = "%d %d %d", as.integer(i), start, end)
        }
        else #direct strand
        {
            start <- as.integer(tmp[1])
            end   <- as.integer(tmp[2]) + as.integer(-3)
            line <- sprintf(fmt = "%d %d %d", as.integer(i), start, end)
        }
        writeLines( line, con = outfile )
    }
    close(outfile)
    invisible(input)
}