R/cpmCountsEP.R

In harshsharma-cb/FASE: Analysis of RNA-Sequencing data using FASE (Finding Alternative Splicing Events).

#' cpmCountsEP
#'
#' @description This function requires ranking of cassette exon/intron retention events as generated by \code{\link{getPvaluesByContrast}} to generate raw read counts and log2cpm expression for the ranking of events in a given contrast only.
#'
#' @param filename file in which ExonPointer/IntronPointer ranking of a contrast is saved. \cr *(Filename should be in csv format.)
#' @param designM design matrix.
#' @param Groups list of sample groups. \cr
#'        Example: If there are two sample groups with three samples each, 'Groups' should be formed as: c(1, 1, 1, 2, 2, 2).
#' @param Gcount list; contains gene-wise  matrix of meta-features read counts times samples, generated by \code{\link{countMatrixGenes}}.
#'
#' @import edgeR
#'
#' @return Read counts and log2cpm expression of ranked cassette exon/intron retention events for the given contrast.
#'
#' @references \enumerate{
#' \item Robinson, M. D., McCarthy, D. J. & Smyth, G. K. edgeR: A Bioconductor package for differential expression analysis of digital gene expression data. Bioinformatics 26, 139–140 (2009)
#'  }
#' @export
#'
cpmCountsEP<- function(filename, designM, Groups, Gcount)
{
        #Function to write two more file for ouput file of Exon or intronPointer
        #one for counts data another for log2cpm data
        test <- read.csv(filename)
        filename<- strsplit(filename,'.',fixed=TRUE)[[1]][1]

        eventName<- as.vector(test$event)
        geneName<- as.vector(test$gene)
        #load('Gcount.Rdata')
        index <- match(unique(geneName), names(Gcount))
        Gcount <- Gcount[index]
        Gcount <- do.call('rbind',Gcount)
        rnames<- unlist(lapply(rownames(Gcount),function(x) strsplit(x,'.',fixed=TRUE)[[1]][[2]]))
        rownames(Gcount) <- rnames
        index <- match(eventName,rnames)
        outCounts <- Gcount[index,,drop=FALSE]

        countfilename <- paste(filename,'_count.csv',sep="",collapse="")
        write.csv(outCounts,file=countfilename)

        counts <- DGEList(counts=outCounts,group=Groups)
        counts <- calcNormFactors(counts)
        fit<- voom(counts,designM)

        E<- ((outCounts !=0)*1) * fit$E
        l_filename <- paste(filename,'_log2cpm.csv',sep="",collapse="")
        write.csv(E,file=l_filename)


}