-added fgsea

Author: pcahan1

NAMESPACE

@@ -5,6 +5,7 @@ export(annSetUp)
 export(butter_classGenes)
 export(butter_classify)
 export(cn_correctZmat)
+export(compileGSEA)
 export(corplot_sub)
 export(corplot_triangle)
 export(cutTreeToNames)
@@ -15,6 +16,7 @@ export(easy_assess)
 export(enrDiff)
 export(expression)
 export(extractImport)
+export(fgsea.wrapper.set)
 export(findVarGenes)
 export(find_genes_byGo)
 export(find_top_pc)

R/annotation.R

@@ -86,6 +86,75 @@ ks.wrapper.set<-function(
     ansList
 }
 
+# fgsea.wrapper.set
+
+#' run fgsea on list of diffExprs
+#'
+#' find gene sets enriched in provided ranked gene list
+#'
+#' @param diffExprList named list of diffExpr results, for example from running gnrAll()
+#' @param geneSets named list of gene sets
+#' @param minSize defaults to 10
+#' @param nPerm defaults to 1e-4
+#'
+#' @return list of  and (2) list of ks.wrapper results
+#' 
+#' @export
+#'
+fgsea.wrapper.set<-function(
+    diffExprList,
+    geneSets,
+    minSize=10,
+    nPerm=1e4)
+{
+
+
+    ansList<-list()
+    qnames<-names(diffExprList)
+    for(qname in qnames){
+        qdiffRes<-diffExprList[[qname]]
+        teststat<-qdiffRes$cval
+        names(teststat)<-rownames(qdiffRes)
+        fRes<-as.data.frame(fgsea(pathways=geneSets, stats=teststat, minSize=minSize, nperm=nPerm))
+        ansList[[qname]]<-fRes
+    }
+    ansList
+}
+
+
+#' compileGSEA
+#' 
+#' make a long table from fgsea.wrapper.set
+#'
+#' make a long table from fgsea.wrapper.set
+#'
+#' @param gseaEnr from rnunning gsea.wrapper.set
+#' @param thresh defaults 0.05
+#'
+#' @return df df
+#' 
+#' @export
+#'
+compileGSEA<-function(
+    gseaEnr,
+    thresh=0.05){
+  
+        longX<-data.frame()
+        for(xname in names(gseaEnr)){
+           tmpDF<-gseaEnr[[xname]]
+           tmpDF<-cbind(tmpDF, cluster=rep(xname, nrow(tmpDF)))
+           tmpDF<-tmpDF[order(tmpDF$NES, decreasing=TRUE),]
+           xi<-which(tmpDF$padj<thresh)
+           sig_gs<-rep(0, nrow(tmpDF))
+           sig_gs[xi]<-1
+           tmpDF<-cbind(tmpDF, isSig=sig_gs)
+           longX<-rbind(longX, tmpDF)
+        }
+        #longX<-longX[longX$padj<thresh,]
+        longX<-cbind(longX, lpval=log10(longX$padj))
+    longX[longX$NES>0,]$lpval<- -1 * longX[longX$NES>0,]$lpval
+    longX
+}
 
 # ks.extract
 
@@ -325,3 +394,5 @@ ks.test.2 <- function (x, y, ..., alternative = c("two.sided", "less", "greater"
 
 
 
+
+

R/wash.R

@@ -74,7 +74,7 @@ prewash<-function
 }
 
 
-
+if(FALSE){
 #' weighted subtraction from mapped reades
 #'
 #' Simulate expression profile of  _total_ mapped reads
@@ -90,12 +90,39 @@ total=1e5){
 
   totalSignal<-sum(vector)
   wAve<-vector/totalSignal
-  resid<-sum(vector)-total #num to subtract from sample
+###  resid<-sum(vector)-total #num to subtract from sample
+  resid<-totalSignal-total #num to subtract from sample
   residW<-wAve*resid # amount to substract from each gene
   ans<-vector-residW
   ans[which(ans<0)]<-0
   ans
 }
+}
+
+#' weighted subtraction from mapped reades
+#'
+#' Simulate expression profile of  _total_ mapped reads
+#' @param vector of total mapped reads per gene/transcript
+#' @param total post transformation sum of read counts
+#'
+#' @return vector of downsampled read mapped to genes/transcripts
+#'
+#' @export
+downSampleW<-function
+(vector,
+ total=1e5,
+ dThresh=0){ 
+
+  totalSignal<-sum(vector)
+  wAve<-vector/totalSignal
+###  resid<-sum(vector)-total #num to subtract from sample
+  resid<-totalSignal-total #num to subtract from sample
+  residW<-wAve*resid # amount to substract from each gene
+  ans<-vector-residW
+  ans[which(ans<dThresh)]<-0
+  ans
+}
+
 
 #' weighted subtraction from mapped reades, applied to all
 #'
@@ -108,9 +135,10 @@ total=1e5){
 #' @export
 weighted_down<-function
 (expRaw,
- total
+ total,
+ dThresh=0
  ){
-    expCountDnW<-apply(expRaw, 2, downSampleW, total)
+    expCountDnW<-apply(expRaw, 2, downSampleW, total=total, dThresh=dThresh)
     #log(1+expCountDnW)
     expCountDnW
   }