R/gbk2g2.euk.R

In seqinr: Biological Sequences Retrieval and Analysis

gbk2g2.euk <- function(
        gbkfile = system.file("sequences/ame1.gbk", package = "seqinr"), g2.coord = "g2.coord")
{
    input <- readLines(gbkfile)
    outfile <- file( description = g2.coord, open ="w+")
    
    
    #
    # Keep lines with annotation flag:
    #
    
    cds <- which(substring(input,1,8) == "     CDS")
    
    features <- which(substr(input,6,6)!=" ")
    
    features <- c(features,length(input))
    
    genes <- which(substr(input,22,26)=="/gene")
    
    genes.cds <- character(length(cds))
    
    for(i in seq_len(length(cds))){
        print(i)
        this.gene=genes[which(genes>cds[i])[1]]
        
        nextfeat=features[which(features>cds[i])[1]]
        
        if(this.gene<nextfeat){
            genes.cds[i]=unlist(strsplit(input[this.gene],split="\""))[2]
            
        }
        else{
            print(paste("no id for cds #",i))
        }
        
    }
    
    
    #
    # Extract boundaries strings
    #
    get.boundaries <- function( index.line )
    {
        
        join <- grep("join",input[index.line])
        
        if(length(join)>0){ ## there are introns !!
            
            
            end.par <- grep("\\)",input[index.line])
            
            
            if(length(end.par)==0){ ## the exons are written on more than one line
                
                next.end.par <-grep("\\)",input)
                next.end.par <-next.end.par[which(next.end.par>index.line)[1]] ## we take the first ending parenthesis
                exons <- paste(input[index.line],paste(substr(input[(index.line+1):next.end.par],22,nchar(input[(index.line+1):next.end.par])),sep="",collapse=""),collapse="",sep="")
            }
            else{
                exons <- input[index.line]
            }
            
            
            
            
            complement <- grep("complement",exons)
            
            if(length(complement)==0){
                
                exons=unlist(strsplit(exons, split="\\("))[2]
                exons=unlist(strsplit(exons, split="\\)"))[1]
                
                exons=unlist(strsplit(exons,split=","))
                exons.begin=unlist(lapply(exons, function(x) unlist(strsplit(x, split="\\.\\."))[1]))
                exons.end=unlist(lapply(exons, function(x) unlist(strsplit(x, split="\\.\\."))[2]))
                
                
                
            }
            
            else{
                exons=unlist(strsplit(exons, split="\\("))[3]
                exons=unlist(strsplit(exons, split="\\)"))[1]
                
                exons=unlist(strsplit(exons,split=","))
                exons.end=unlist(lapply(exons, function(x) unlist(strsplit(x, split="\\.\\."))[1]))
                exons.begin=unlist(lapply(exons, function(x) unlist(strsplit(x, split="\\.\\."))[2]))
                
                
                
            }
        }
        else{
            complement <- grep("complement",input[index.line])
            
            exons=unlist(strsplit(input[index.line],split=" "))
            exons=exons[exons!=""]
            exons=exons[2]
            
            if(length(complement)==0){
                exons.begin=unlist(strsplit(exons, split="\\.\\."))[1]
                exons.end=unlist(strsplit(exons, split="\\.\\."))[2]
                
            }
            else{
                exons=unlist(strsplit(exons, split="\\("))[2]
                exons=unlist(strsplit(exons, split="\\)"))[1]
                exons.begin=unlist(strsplit(exons, split="\\.\\."))[2]
                exons.end=unlist(strsplit(exons, split="\\.\\."))[1]
            }
            
        }
        
        
        return(list(exons.begin,exons.end))
        
        
    }
    
    
    
    
    
    
    already=character(0)
    
    if(length(cds)>0){
        for(i in seq_len(length(cds))){
            
            boundaries=get.boundaries(cds[i])
            exons.begin=boundaries[[1]]
            exons.end=boundaries[[2]]
            
            
            
            for(j in seq_len(length(exons.begin))){
                phrase=paste(genes.cds[i],exons.begin[j],exons.end[j])
                
                if(!phrase%in%already){
                    writeLines(paste(genes.cds[i],exons.begin[j],exons.end[j],sep=" "),outfile)
                    already=c(already,phrase)
                }
            }
        }
        
    } 
    
    
    
    close(outfile)
    
}