aux <- readRDS('tcgaIsoHer2LumA.rds');
counts <- aux$counts;
conds <- aux$conds;
annot <- aux$annot; rownames(annot) <- annot[,'ENST'];
rm(aux);

transcAnnot <- data.frame(annot[,c('SYMBOL', 'ENST')]);
colnames(transcAnnot) <- c('GeneID', 'IsoID');

table(colnames(counts) == rownames(conds));
table(rownames(counts) == rownames(annot));


allGenes <- sort(unique(annot[,'SYMBOL']));
geneCounts <- do.call(rbind, lapply(allGenes, function(actGene) {
  colSums(counts[rownames(annot)[annot[,'SYMBOL'] == actGene],,drop=F])
}))
rownames(geneCounts) <- allGenes; rm(allGenes);

library('limma');

myDesign <- model.matrix(~PAM50, data.frame(PAM50=as.character(conds$PAM50)));

## transcripts
countsDge <- DGEList(counts, genes=transcAnnot);
keep <- filterByExpr(countsDge, myDesign); table(keep);
countsDge <- countsDge[keep,, keep.lib.sizes=FALSE]
countsDge <- calcNormFactors(countsDge)
countsV <- voom(countsDge, myDesign)
transcFit <- lmFit(countsV, myDesign)
transcEx <- diffSplice(transcFit, geneid='GeneID')
table(transcEx$p.value[,2] < 0.05);
# FALSE  TRUE
# 36982  7871

## genes
geneCountsDge <- rowsum(countsDge[rownames(transcEx$genes),], transcEx$genes$GeneID, reorder=FALSE);
# geneCountsDge <- DGEList(counts=geneCounts);
keep <- filterByExpr(geneCountsDge, myDesign); table(keep);
geneCountsDge <- geneCountsDge[keep,,keep.lib.sizes=FALSE];
geneCountsDge <- calcNormFactors(geneCountsDge);
geneCountsV <- voom(geneCountsDge, design=myDesign);
genesFit <- lmFit(geneCountsV, myDesign);
genesFit <- eBayes(genesFit);
table(genesFit$p.value[,2] < 0.05);
# FALSE  TRUE
# 6908  9758