library(Rsamtools)
library(GenomicFeatures)
library(rtracklayer)
library(GenomicAlignments)
library(componentSkeleton)

execute <- function(cf) {

    type <- get.parameter(cf, 'entity', type = 'string')
    ref  <- get.parameter(cf, 'reference', type = 'string')
    inputDB    <- get.parameter(cf, 'inputDB', type = 'string')
    inputTable <- get.parameter(cf, 'inputTable', type = 'string')
    entityList <- get.parameter(cf, 'entityList', type = 'string')
    unstranded <- get.parameter(cf, 'ignoreStrands', type = 'boolean')

    if (!(type %in% c('gene','transcript','exon'))){
        write.error(cf, paste('Parameter type must be one of gene, exon or transcript. Encountered ', type))
        return(INVALID_INPUT)
    }
    if (type == 'transcript'){ type <- 'tx' }
    if (entityList == ''){ 
        entityList <- NULL
    } else {
        entityList <- unlist(strsplit(entityList, ','))
    }

    if (ref != ""){
        txdb <- loadFeatures(ref)
    } else {
        #txdb <- makeTranscriptDbFromUCSC(genome         = inputDB,
        #                                 tablename      = inputTable,
        #                                 transcript_ids = entityList)
        txdb <- makeTxDbFromUCSC(genome         = inputDB,
                                 tablename      = inputTable,
                                 transcript_ids = entityList)
    }
    exonRanges <- exonsBy(txdb, type)
    counts     <- data.frame()

    array <- Array.read(cf, 'in')
    for (i in 1:Array.size(array)) {
        file   <- Array.getFile(array, i)
        #aligns <- readBamGappedAlignments(file)
        aligns <- readGAlignments(file)
        if (unstranded){
            exonRanges <- unlist(exonRanges)
            strand(exonRanges) <- "*"
            geneNames <- sub("\\..*$", "", names(exonRanges))
            exonRanges <- split(exonRanges, geneNames)
        }
        if (i == 1) {
            counts <- countOverlaps(exonRanges, aligns)
        } else {
            counts <- cbind(counts, countOverlaps(exonRanges, aligns))
        }
    }
    counts <- data.frame(ID=names(exonRanges), counts)
    CSV.write(get.output(cf, 'out'), counts)
}

main(execute)