Data now included with package

.Rbuildignore

@@ -0,0 +1,3 @@
+^.*\.Rproj$
+^\.Rproj\.user$
+^data-raw$

.gitignore

@@ -1 +1,4 @@
-/Data/*
+.Rproj.user
+.RDS
+*.pdf
+*.pptx

R/FGEM_Func.R

@@ -47,15 +47,7 @@ glist2df <- function(Genes,feat.name){
   return(retdf)
 }
 
-GO2df <-function(Genes,terms="BP"){
-  require(reshape2)
-  require(topGO)
-  require(dplyr)
-  mdf <- melt(topGO::annFUN.org(terms,feasibleGenes=Genes,mapping="org.Hs.eg.db",ID="symbol"))
-  mdf <- as_data_frame(mdf)
-  mdf <- dplyr::select(mdf,Gene=value,feature=L1) %>% mutate(value=1,class=paste0("GO_",terms))
-  return(mdf)
-}
+
 
 
 cfeat_mat <- function(annodf,datadf){
@@ -171,7 +163,7 @@ sem_df <-function(full_feat,scale=F,prior_mean=0.02){
     feat_mat <- feat_mat[,!is.na(fBeta)[-1]]
     fBeta <- coefficients(glm(tmu~feat_mat+0,family=quasibinomial(link="logit")))
   }
-  retdf <- FGEM(fbeta = fBeta,feat_mat=feat_mat,BF=BF)
+  retdf <- FGEM(fBeta = fBeta,feat_mat=feat_mat,BF=BF)
 
   return(retdf)
 }

analyses/Cluster_FGEM.R

@@ -3,47 +3,24 @@ library(FGEM)
 library(dplyr)
 library(BBmisc)
 
-# chunki <- 1
-# # chunki <- as.integer(commandArgs(trailingOnly = F))
-# sub_all_anno_df<- readRDS("~/Dropbox/BayesianDA/FGEM_Data/sub_annotations.RDS")
-# sub_all_anno_df <- filter(sub_all_anno_df,grepl("GO",feature))
-# datadf <- readRDS("~/Dropbox/BayesianDA/FGEM_Data/TADA.RDS")
-#
-# featureinds <- unique(sub_all_anno_df$feat_ind)
-# tfi <- filter(sub_all_anno_df, feat_ind==2)
-# mchunk <- chunk(featureinds,chunk.size = 50)
-# for(chunki in 90:length(mchunk)){
-#   mchunkl <- mchunk[[chunki]]
-#   miniter <- min(mchunkl)
-#   maxiter <- max(mchunkl)
-#   s_anno_df <- filter(sub_all_anno_df,between(feat_ind,miniter,maxiter))
-#   f_results <- group_by(s_anno_df,feature) %>% do(cfeat_df(.,datadf,impute=F))
-#   %>% do(fisher_comp(.))
-#   #results <- group_by(s_anno_df,feature) %>%  do(cfeat_df(.,datadf,impute=F)) %>%  do(sem_df(.)) %>% ungroup()
-#   saveRDS(f_results,paste0("~/Dropbox/BayesianDA/FGEM_Data/FGEM_chisq_",miniter,"_",maxiter,".RDS"))
-#   gc()
-# }
-
 #chunki <- 1
 
-
 chunki <- as.integer(commandArgs(trailingOnly = T))
 cat(chunki,"\n")
-sub_all_anno_df<- readRDS("/scratch/midway/nwknoblauch/FGEM/remaining_features.RDS")
-datadf <- readRDS("/scratch/midway/nwknoblauch/FGEM/TADA.RDS")
-
-featureinds <- unique(sub_all_anno_df$feat_ind)
+GOdf <- distinct(BPGOdf,feature) %>% mutate(feat_ind=1:n()) %>% inner_join(BPGOdf)
+datadf <- BFdata
+featureinds <- unique(GOdf$feat_ind)
 cat("There are ",length(featureinds),"features\n")
 
-mchunk <- chunk(featureinds,chunk.size = 100)
+mchunk <- chunk(featureinds,chunk.size = 3)
 cat("There are ",length(mchunk)," chunks of size 100\n")
 mchunk <- mchunk[[chunki]]
 miniter <- min(mchunk)
 maxiter <- max(mchunk)
 
-s_anno_df <- filter(sub_all_anno_df,between(feat_ind,miniter,maxiter))
+s_anno_df <- filter(GOdf,between(feat_ind,miniter,maxiter))
 
-results <- group_by(s_anno_df,feature) %>%  do(cfeat_df(.,datadf,impute=F)) %>%  do(sem_df(.)) %>% ungroup()
+results <- select(s_anno_df,-feat_ind) %>% group_by(feature) %>%  do(sem_df(cfeat_df(.,datadf,impute=F))) %>% ungroup()
 
 saveRDS(results,paste0("/scratch/midway/nwknoblauch/FGEM/FGEM_Results/rem_FGEM_",miniter,"_",maxiter,".RDS"))
 

data-raw/BFdata.R

@@ -0,0 +1,8 @@
+library(FGEM)
+library(readr)
+library(dplyr)
+datalink <- "https://www.dropbox.com/s/7k62nly5jzj2j00/TADA_ASC_SSC_results_Dec23.csv?dl=1"
+datadf <- read_delim(datalink,delim=",",col_names=T)
+datadf <- dplyr::select(datadf,Gene,BF,qvalue,pval.TADA)
+BFdata <- datadf
+save(datadf,"data/BFdata.rdata")

data-raw/ExAC.R

@@ -0,0 +1,9 @@
+library(dplyr)
+require(tidyr)
+tmp <- tempfile()
+exacurl <- "ftp://ftp.broadinstitute.org/pub/ExAC_release/release0.3.1/functional_gene_constraint/fordist_cleaned_exac_r03_march16_z_pli_rec_null_data.txt"
+exacf <- download.file(exacurl,tmp)
+exacdf <- read.table(tmp,sep="\t",header=T,stringsAsFactors = F)
+exacdf <- dplyr::select(exacdf,-chr,-transcript)  %>% dplyr::rename(Gene=gene)
+exacdf <- gather(exacdf,feature,value,-Gene)  %>% mutate(class="ExAC")
+save(exacdf,"data/ExAC.rdata")

data-raw/GO_BP.R

@@ -0,0 +1,11 @@
+library(readr)
+require(reshape2)
+require(topGO)
+library(dplyr)
+datalink <- "https://www.dropbox.com/s/7k62nly5jzj2j00/TADA_ASC_SSC_results_Dec23.csv?dl=1"
+datadf <- read_delim(datalink,delim=",",col_names=T)
+genes <- unique(datadf$Gene)
+mdf <- melt(topGO::annFUN.org("BP",feasibleGenes=genes,mapping="org.Hs.eg.db",ID="symbol"))
+mdf <- dplyr::as_data_frame(mdf)
+BPGOdf <- dplyr::select(mdf,Gene=value,feature=L1) %>% dplyr::mutate(value=1,class="GO_BP")
+save(BPGOdf,file="data/BP_GO.rdata")