[v0.7.1.9000] Gene based analysis option

DESCRIPTION

@@ -1,7 +1,7 @@
 Package: coRge
 Type: Package
 Title: An R package for the Examination of Multiple Correction Methodologies in Accurate Genomic Environments.
-Version: 0.6.7.9000
+Version: 0.7.0.9000
 Date: 2016-21-07
 Author: Christopher B. Cole
 Maintainer: Christopher B. Cole <[email protected]>
@@ -12,6 +12,7 @@ Imports: Rcpp (>= 0.11.3),
    foreach,
    devtools,
    lazyeval,
+   futile.logger,
    data.table,
    magrittr,
    dplyr

NAMESPACE

@@ -15,18 +15,16 @@ export(ld_cor)
 export(phen)
 export(rand)
 export(rand0)
-export(rebuild)
 export(returnLD)
 export(sim.gen)
 export(stratify)
 export(sub_job)
 export(th)
 import(devtools)
 import(foreach)
-importFrom(Rcpp,compileAttributes)
+import(futile.logger)
 importFrom(Rcpp,evalCpp)
 importFrom(data.table,fread)
-importFrom(devtools,document)
 importFrom(dplyr,"%>%")
 importFrom(dplyr,filter)
 importFrom(dplyr,filter_)

R/analyze.R

@@ -1,4 +1,4 @@
-#' Correction of Genomes in R
+#' @title Calculate FDR and sFDR for a simulated GWAS
 #'
 #' Software for the Examination of Multiple Correction Methodologies in Accurate Genomic Environments
 #'
@@ -13,15 +13,16 @@
 #' @param pc Proportion of causal SNPs in the second strata
 #' @param pnc Proportion of noncausal SNPs in the second strata
 #' @param nc Number of causal SNPs
-#' @param gen Gen matrix if local 
+#' @param gen Gen matrix if local
 #' @param summary Summary file if preread
 #' @param maf Maf or no
 #' @param maf_range MAF range
 #'
+#' @import futile.logger
 #' @import foreach
 #' @import devtools
 #' @importFrom lazyeval interp
-#' @importFrom data.table fread 
+#' @importFrom data.table fread
 #' @importFrom magrittr %<>%
 #' @importFrom dplyr mutate mutate_ filter filter_ select select_ sample_n %>%
 #' @importFrom Rcpp evalCpp
@@ -33,12 +34,24 @@
 
 
 
-analyze <- function(i = double(), j = double(), mode = "default", path.base = "/scratch/hpc2862/CAMH/perm_container/container_", summary.file = "/scratch/hpc2862/CAMH/perm_container/snp_summary2.out", output = "~/repos/coR-ge/data/test_run2.txt", test = TRUE, safe = TRUE, local = FALSE, h2 = 0.45, pc = 0.5, pnc = 0.5, nc = 1000, gen = NULL, summary = summary, maf = FALSE, maf_range = NULL){
-
-
-#message(paste0("coR-ge v", packageVersion("coRge"), " \t \t http://github.io/Chris1221/coR-ge \n \n \t Parameters in use: \n \t \t i: ",i, "\n \t \t j: ", j, "\n \t \t path.base: ", path.base,"\n \t \t summary.file: ", summary.file,"\n \t \t output: ", output,"\n \t \t th2: ", h2,"\n \t \t pc: ", pc,"\n \t \t pnc: ", pnc,"\n \t \t nc: ", nc, "\n \n \t Local options: \n \t \t local: ", local,"\n \t \t gen ", gen,"\n \t \t summary: ", summary, " \n \n \t Testing options: \n \t \t test: ", test,"\n \t \t safe: ", safe, "\n \n LOG: \n"))
-
-
+analyze <- function(i = NA,
+		    j = NA,
+		    mode = "default",
+		    path.base = "/scratch/hpc2862/CAMH/perm_container/container_",
+		    summary.file = "/scratch/hpc2862/CAMH/perm_container/snp_summary2.out",
+		    output = "~/repos/coR-ge/data/test_run2.txt",
+		    test = TRUE,
+		    safe = TRUE,
+		    local = FALSE,
+		    h2 = 0.45,
+		    pc = 0.5,
+		    pnc = 0.5,
+		    nc = 1000,
+		    gen = NULL,
+		    summary = summary,
+		    maf = FALSE,
+		    maf_range = NULL,
+		    gene_kb = 0){
 
 	if(local) {
 
@@ -49,7 +62,7 @@ analyze <- function(i = double(), j = double(), mode = "default", path.base = "/
 
 			message("Error Checking")
 
-		if(any(is.null(c(i,j,path.base, summary.file)))) stop("Please complete all input arguemnets")
+		#if(any(is.null(c(path.base, summary.file)))) stop("Please complete all input arguemnets")
 
 		path <- paste0(path.base,i,"_",j,"/")
 		setwd(path)
@@ -132,7 +145,7 @@ analyze <- function(i = double(), j = double(), mode = "default", path.base = "/
 
 	out <- correct(strata=strata, n_strata = n_strata, assoc = P_list, group = FALSE)
 
-	
+
 
 
 #            ___                                       _
@@ -241,24 +254,83 @@ analyze <- function(i = double(), j = double(), mode = "default", path.base = "/
 	# Initiate it with 0 and then sequentially overwrite.
 	strata$ld <- 0
 
-	# Calculate the LD from ld.cpp 
+	# Calculate the LD from ld.cpp
 	LdList <- ld_cor(snps, gen)
-	
+
 	# Loop through each of the LD TP thresholds.
 	for(th in c(0.2, 0.4, 0.6, 0.8, 0.9, 1)){
-		
+
 		# Create a list of SNPs which have higher than the threshold
 		# LD.
 		snp_b <- LdList$rsid[LdList$ld > th]
-		
+
 		# Assign the threshold value to each of these SNPs.
 		strata$ld[strata$rsid %in% snp_b] <- th
 	}
 
 	 out <- correct(strata=strata, n_strata = n_strata, assoc = P_list, group = FALSE, mode = "ld")
 
 
-  }
+  } else if(mode == "genes"){
+
+
+    		message("Selecting Causal SNPs")
+
+    	snps <- causal.snps(summary, nc = nc, maf = maf)
+    	colnames(snps)[3] <- "V3"
+
+
+    		message("Merging together and converting from Oxford to R format...")
+
+    	comb <- as.data.frame(merge(gen, snps, by = "V3"))
+
+    		comb$rsid <- NULL
+    		comb$chromosome <- NULL
+    		comb$all_maf <- NULL
+    		comb$k <- NULL
+    		comb$chromosomes <- NULL
+
+    		WAS <- calculate_was(gen = comb, snps = snps, h2 = h2)
+
+    	samp$Z <- as.character(foreach(q = 1:length(WAS), .combine = 'c') %do% WAS[q] + rnorm(1, 0, sd = sqrt(1-h2)))
+
+
+    	var <- data.frame(0, 0, 0, "P")
+    	samp$Z <- as.character(samp$Z)
+    	colnames(var) <- colnames(samp)
+    	samp <- rbind(var, samp)
+
+	P_list <- assoc_wrapper(gen, samp)
+
+    	message("Performing correction")
+
+	snps %>% select(rsid) %>% as.vector -> snp_list
+
+	#this is a major assumption so leave it
+	n_strata <- 2
+	strata <- stratify(snp_list = snp_list, summary = summary, n_strata = n_strata, pc = pc, pnc = pnc, mode = "genes", gene_kb = gene_kb)
+
+	# Initialize a column for the LD to sit in
+	# Initiate it with 0 and then sequentially overwrite.
+	strata$ld <- 0
+
+	# Calculate the LD from ld.cpp
+	LdList <- ld_cor(snps, gen)
+
+	# Loop through each of the LD TP thresholds.
+	for(th in c(0.2, 0.4, 0.6, 0.8, 0.9, 1)){
+
+		# Create a list of SNPs which have higher than the threshold
+		# LD.
+		snp_b <- LdList$rsid[LdList$ld > th]
+
+		# Assign the threshold value to each of these SNPs.
+		strata$ld[strata$rsid %in% snp_b] <- th
+	}
+	out <- correct(strata=strata, n_strata = n_strata, assoc = P_list, group = FALSE, mode="ld")
+
+
+}
 
 
   out$h2 <- h2
@@ -270,10 +342,10 @@ analyze <- function(i = double(), j = double(), mode = "default", path.base = "/
   out$maf_l <- maf_range[1]
   out$maf_u <- maf_range[2]
 
+ write.table(out, file = output, append = T, quote = F, row.names = F, col.names = F)
+ # if(!file.exists(output)) suppressWarnings(write.table(out, output, row.names = F, col.names = TRUE, quote = F, append = T)) else if(file.exists(output)) write.table(out, output, row.names = F, col.names = F, quote = F, append = T)
 
-  if(!file.exists(output)) suppressWarnings(write.table(out, output, row.names = F, col.names = TRUE, quote = F, append = T)) else if(file.exists(output)) write.table(out, output, row.names = F, col.names = F, quote = F, append = T)
-
-
+ message("This has passed the writting")
   return(0)
 
 }

R/onAttach.R

@@ -3,7 +3,7 @@
   if (interactive()) {
     v = packageVersion("coRge")
 
-    packageStartupMessage(paste0("coRge v", v, ". This version is unstable and should not be used in development. More information and guide at https://github.com/chris1221/cor-ge."))
+    packageStartupMessage(paste0("coRge v", v, ".\nGuides, issues, and FAQ: https://github.com/Chris1221/coRge \nPrepublication:\t\t arXiv: XXXXX \nCitation:\t\t citation(\"coRge\")"))
 
 
   }

R/stratify.R

@@ -13,23 +13,82 @@
 #'
 #' @export
 
-stratify <- function(snp_list = NULL, summary = NULL, pnc = NULL, pc = NULL, n_strata = NULL){
+stratify <- function(snp_list = NULL, summary = NULL, pnc = NULL, pc = NULL, n_strata = NULL, mode = "default", gene_kb = 0){
+
+	flog.debug("Entering into the stratify() function block")
 
 	if(any(c(is.null(snp_list), is.null(summary), is.null(n_strata)))) stop("Missing a required input arguement")
 
-	h0_summary <- summary %>% filter(!(rsid %in% unlist(snp_list))) %>% mutate(h1 = F)
-	h1_summary <- summary %>% filter(rsid %in% unlist(snp_list)) %>% mutate(h1 = T)
-
-	h1_s1 <- h1_summary %>% sample_n(floor(nrow(h1_summary)*pc))
-	h1_s2 <- h1_summary %>% filter(!(rsid %in% h1_s1$rsid)) 
+	if(mode == "genes"){
+
+		flog.debug("Gene mode selected")
+
+		h0_summary <- summary %>%
+			filter(!(rsid %in% unlist(snp_list))) %>%
+			mutate(h1 = F)
+		h1_summary <- summary %>%
+			filter(rsid %in% unlist(snp_list)) %>%
+			mutate(h1 = T)
+
+		flog.trace(paste0("There are ", nrow(h1_summary), " causal SNPs in the current analysis and ", nrow(h0_summary), " non causal SNPs."))
+
+		# Run through each SNP in the h1 summary and give it a gene
+		# identifyier to be used in the gene column.
+
+		# Randomly put floor(pc*nrow(h1_summary)) causal into the h1
+		# 	and the rest into the second strata.
+		h1_s1 <- h1_summary %>%
+			sample_n(floor(nrow(h1_summary)*pc))
+		h1_s2 <- h1_summary %>%
+			filter(!(rsid %in% h1_s1$rsid))
+
+		s1 <- h0_summary %>%
+			filter(position > h1_s1$position[1] - gene_kb*1000 & position < h1_s1$position[1] + gene_kb*1000) %>% nrow
+
+		for(i in 2:nrow(h1_s1)){
+			s1 <- h0_summary %>%
+			filter(position > h1_s1$position[1] - gene_kb*1000 & position < h1_s1$position[1] + gene_kb*1000) %>%
+			rbind(s1)
+		}
+
+		s2 <- h0_summary[!(h0_summary$rsid %in% s1$rsid),]
+
+		s1 %<>% rbind(h1_s1) %>%
+			mutate(s = 1)
+		s2 %<>% rbind(h1_s2) %>%
+			mutate(s = 2)
+
+		strata <- rbind(s1, s2)
+
+
+
+	} else {
+
+		h0_summary <- summary %>%
+			filter(!(rsid %in% unlist(snp_list))) %>%
+			mutate(h1 = F)
+		h1_summary <- summary %>%
+			filter(rsid %in% unlist(snp_list)) %>%
+			mutate(h1 = T)
+
+		h1_s1 <- h1_summary %>%
+			sample_n(floor(nrow(h1_summary)*pc))
+		h1_s2 <- h1_summary %>%
+			filter(!(rsid %in% h1_s1$rsid))
+
+		s1 <- h0_summary %>%
+			sample_n(floor(nrow(h0_summary)*pnc) - length(snp_list))
+		s2 <- h0_summary %>%
+			filter(!(rsid %in% s1$rsid))
 
-	s1 <- h0_summary %>% sample_n(floor(nrow(h0_summary)*pnc) - length(snp_list))
-	s2 <- h0_summary %>% filter(!(rsid %in% s1$rsid))
+		s1 %<>% rbind(h1_s1) %>%
+			mutate(s = 1)
+		s2 %<>% rbind(h1_s2) %>%
+			mutate(s = 2)
 
-	s1 %<>% rbind(h1_s1) %>% mutate(s = 1)
-	s2 %<>% rbind(h1_s2) %>% mutate(s = 2)
+		strata <- rbind(s1, s2)
+	}
 
-	strata <- rbind(s1, s2)
 	return(strata)
 
 }

README.md

@@ -5,16 +5,16 @@
 
 | Branch | Travis-CI | Appveyor | Coverage | CRAN | Downloads | Publication |
 | :--- | :---: | :---: | :--: | :---: | :---: | :---: |
-| `master` | ![Build Status](https://travis-ci.org/Chris1221/coR-ge.svg?branch=master) | ![Build status](https://ci.appveyor.com/api/projects/status/v64oe85q29btxln9?svg=true) | [![codecov.io](https://codecov.io/github/Chris1221/coR-ge/coverage.svg?branch=master)](https://codecov.io/github/Chris1221/coR-ge?branch=master) | ![CRAN_Status_Badge](http://www.r-pkg.org/badges/version/coRge) | ![](http://cranlogs.r-pkg.org/badges/coRge) | GitXiv |
-| `devel` |![Build Status](https://travis-ci.org/Chris1221/coR-ge.svg?branch=devel) | [![Build status](https://ci.appveyor.com/api/projects/status/v64oe85q29btxln9?svg=true)](https://ci.appveyor.com/project/Chris1221/miner) | [![codecov.io](https://codecov.io/github/Chris1221/coR-ge/coverage.svg?branch=devel)](https://codecov.io/github/Chris1221/coR-ge?branch=devel) | ![CRAN_Status_Badge](http://www.r-pkg.org/badges/version/coRge) | ![](http://cranlogs.r-pkg.org/badges/coRge) | GitXiv | 
+| `master` | ![Build Status](https://travis-ci.org/Chris1221/coRge.svg?branch=master) | ![Build status](https://ci.appveyor.com/api/projects/status/v64oe85q29btxln9?svg=true) | [![codecov.io](https://codecov.io/github/Chris1221/coRge/coverage.svg?branch=master)](https://codecov.io/github/Chris1221/coRge?branch=master) | ![CRAN_Status_Badge](http://www.r-pkg.org/badges/version/coRge) | ![](http://cranlogs.r-pkg.org/badges/coRge) | GitXiv |
+| `devel` |![Build Status](https://travis-ci.org/Chris1221/coRge.svg?branch=devel) | [![Build status](https://ci.appveyor.com/api/projects/status/v64oe85q29btxln9?svg=true)](https://ci.appveyor.com/project/Chris1221/miner) | [![codecov.io](https://codecov.io/github/Chris1221/coRge/coverage.svg?branch=devel)](https://codecov.io/github/Chris1221/coRge?branch=devel) | ![CRAN_Status_Badge](http://www.r-pkg.org/badges/version/coRge) | ![](http://cranlogs.r-pkg.org/badges/coRge) | GitXiv | 
 
 -----------------------------------
 
 To install:
 
 ```R
 if(!require(devtools)) install.packages("devtools")
-devtools::install_github("Chris1221/coR-ge")
+devtools::install_github("Chris1221/coRge")
 ```
 
 To run:

data/.gitignore

@@ -0,0 +1,19 @@
+# History files
+.Rhistory
+.Rapp.history
+
+# Session Data files
+.RData
+
+# Example code in package build process
+*-Ex.R
+
+# RStudio files
+.Rproj.user/
+
+# produced vignettes
+vignettes/*.html
+vignettes/*.pdf
+
+# OAuth2 token, see https://github.com/hadley/httr/releases/tag/v0.3
+.httr-oauth

data/R/test_run.R

@@ -0,0 +1,26 @@
+#!/usr/bin/Rscript
+
+zz <- file("~/repos/coR-ge/data/test_run.sink", open = "wt")
+
+sink(zz)
+sink(zz, type = "message")
+
+if(!require(devtools)) install.packages("devtools")
+devtools::install_github("Chris1221/coR-ge", ref = "devel")
+library(coRge)
+
+library(dplyr)
+library(foreach)
+library(magrittr)
+library(data.table)
+
+for(h2 in seq(0.1, 0.9, by = 0.1)){
+
+	coRge::analyze(i = 3, j = 3, h2 = h2)
+
+	}
+
+sink(type = "message")
+sink()
+
+system("touch ~/repos/coR-ge/data/done")

data/makefile

@@ -0,0 +1,10 @@
+# Makefile for analysis 
+# Going from analysis script to output file
+#
+# Will have to change this probably
+
+all: test_run.txt
+
+test_run.txt: pipeline2.sh pipeline2.R
+	bash pipeline2.sh
+