Integrate ODP correction

R/RcppExports.R

@@ -1,4 +1,4 @@
-# This file was generated by Rcpp::compileAttributes
+# Generated by using Rcpp::compileAttributes() -> do not edit by hand
 # Generator token: 10BE3573-1514-4C36-9D1C-5A225CD40393
 
 #' @export

R/correct.R

@@ -1,9 +1,9 @@
-#' Correct strata by group
+#' @title Correct strata by group
 #'
-#' Note that group name must be a factor
+#' @description Note that group name must be a factor
 #'
-#' @param strata Strata
-#' @param n_strata Number of strata
+#' @param strata Strata object returned from \code{stratify}.
+#' @param n_strata Number of strata you wish to consider. Note: Only 2 currently supported.
 #' @param assoc Assoc file path
 #'
 #' @importFrom data.table fread
@@ -15,6 +15,15 @@
 
 correct <- function(strata = NULL, n_strata = NULL, assoc = NULL, group = FALSE, group_name = NULL, mode = "default"){
 
+	if(mode == "odp"){
+
+		flog.info("Using ODP for correction.")
+
+
+
+	}
+
+
   if(mode == "default"){
 
   strata <- as.data.frame(strata)

R/gen_phen_df.R

@@ -10,10 +10,12 @@ assoc_wrapper <- function(gen, samp){
 
 	phen <- samp$Z[-1] %>% as.vector %>% as.double
 
-	t <- assoc(as.matrix(gen[,-c(1:5)]), phen)
+	out <- assoc(as.matrix(gen[,-c(1:5)]), phen)
+	t <- out[,1] / out[,2]
+
 	P <- pt(t, df = length(phen) -2, lower = FALSE)*2
 
-	return(data.frame(rsid = gen[,2], P = P))
+	return(data.frame(rsid = gen[,2], out[,1], out[,2], P = P))
 
 
 

TODO.md

@@ -0,0 +1,4 @@
+## TODO 
+
+- Simulate one and look at MAF and other things
+- Change the referene panel.

inst/analysis/new_ref/makefile

@@ -0,0 +1,45 @@
+## This analysis is centered around using a new reference panel with more SNPs
+#	hopefully attenuating the increase in false positives that we have been 
+#	observing in the previous analyses.
+#
+#	This is going hopefully towards publication eventually. 
+
+# Set up some variables first.
+home=/scratch/hpc2862/CAMH/new_corge
+hapgen2=$(home)/bin/hapgen2
+
+# The default rule
+# 	For now just update the git repo.
+.all: $(home)/sim/.done
+
+# Submission
+#
+sub:
+	qsub -N "new_run" sub.sh
+
+# Git add rule
+# 	
+git:
+	cd ../../.. && git add -A
+	git commit -am "update coRge for new analysis see changelog" 
+
+# Download and unpack the reference files.
+$(home)/ref/.done: 
+	cd $(home)/ref && wget https://mathgen.stats.ox.ac.uk/impute/1000GP_Phase3.tgz 
+	cd $(home)/ref && tar -xvzf 1000GP_Phase3.tgz
+	cd $(home)/ref/1000GP_Phase3 && mv * ../
+	cd $(home)/ref && rm -rf 1000GP_Phase3.tgz 1000GP_Phase3
+	touch $(home)/ref/.done # Done
+
+# This is the simulation protocol
+#	Note that the target is a placeholder generated at the end of the simulation
+#	as opposed to any of the actual files for the sake of simplicity.
+$(home)/sim/.done: $(home)/ref/.done 
+	mkdir -p $(home)/sim
+	cd $(home)/ref && $(hapgen2) -h 1000GP_Phase3_chr1.hap \
+		-l 1000GP_Phase3_chr1.legend \
+		-m genetic_map_chr1_combined_b37.txt \
+		-dl 723891 1 1.0 1.0 \
+		-n 1000 \
+		-int 0 500000000 \
+		-o  ../sim/chr1_sim1

inst/analysis/new_ref/sub.sh

@@ -0,0 +1,12 @@
+#!/bin/bash
+#$ -S /bin/bash
+#$ -q abaqus.q
+#$ -l qname=abaqus.q
+#$ -cwd
+#$ -V
+#$ -l hostname=sw0050
+#$ -j y
+#$ -o /home/hpc2862/repos/coR-ge/inst/logs/$JOB_NAME.txt
+
+cd /home/hpc2862/repos/coR-ge/inst/analysis/new_ref
+make .all

inst/bash/test.txt

@@ -0,0 +1,12 @@
+job-ID  prior   name       user         state submit/start at     queue                          slots ja-task-ID 
+-----------------------------------------------------------------------------------------------------------------
+1841411 0.50500 hi         hpc2862      r     10/26/2016 15:42:22 abaqus.q@sw0050.hpcvl.org          1 2
+1841411 0.50500 hi         hpc2862      t     10/26/2016 15:42:22 abaqus.q@sw0050.hpcvl.org          1 3
+1841411 0.50500 hi         hpc2862      t     10/26/2016 15:42:23 abaqus.q@sw0050.hpcvl.org          1 4
+1841411 0.50500 hi         hpc2862      t     10/26/2016 15:42:23 abaqus.q@sw0050.hpcvl.org          1 5
+1841411 0.50500 hi         hpc2862      t     10/26/2016 15:42:23 abaqus.q@sw0050.hpcvl.org          1 6
+1841411 0.50500 hi         hpc2862      t     10/26/2016 15:42:23 abaqus.q@sw0050.hpcvl.org          1 7
+1841411 0.50500 hi         hpc2862      t     10/26/2016 15:42:23 abaqus.q@sw0050.hpcvl.org          1 8
+1841411 0.50500 hi         hpc2862      t     10/26/2016 15:42:24 abaqus.q@sw0050.hpcvl.org          1 9
+1841411 0.50500 hi         hpc2862      t     10/26/2016 15:42:24 abaqus.q@sw0050.hpcvl.org          1 10
+1841411 0.50500 hi         hpc2862      r     10/26/2016 15:42:22 abaqus.q@sw0053.hpcvl.org          1 1

src/RcppExports.cpp

@@ -1,4 +1,4 @@
-// This file was generated by Rcpp::compileAttributes
+// Generated by using Rcpp::compileAttributes() -> do not edit by hand
 // Generator token: 10BE3573-1514-4C36-9D1C-5A225CD40393
 
 #include <RcppArmadillo.h>
@@ -10,61 +10,61 @@ using namespace Rcpp;
 arma::mat gen_cor(arma::vec causal, arma::mat all);
 RcppExport SEXP coRge_gen_cor(SEXP causalSEXP, SEXP allSEXP) {
 BEGIN_RCPP
-    Rcpp::RObject __result;
-    Rcpp::RNGScope __rngScope;
+    Rcpp::RObject rcpp_result_gen;
+    Rcpp::RNGScope rcpp_rngScope_gen;
     Rcpp::traits::input_parameter< arma::vec >::type causal(causalSEXP);
     Rcpp::traits::input_parameter< arma::mat >::type all(allSEXP);
-    __result = Rcpp::wrap(gen_cor(causal, all));
-    return __result;
+    rcpp_result_gen = Rcpp::wrap(gen_cor(causal, all));
+    return rcpp_result_gen;
 END_RCPP
 }
 // returnLD
 Rcpp::List returnLD(arma::uvec cIndex, arma::mat gen, arma::vec bpVec);
 RcppExport SEXP coRge_returnLD(SEXP cIndexSEXP, SEXP genSEXP, SEXP bpVecSEXP) {
 BEGIN_RCPP
-    Rcpp::RObject __result;
-    Rcpp::RNGScope __rngScope;
+    Rcpp::RObject rcpp_result_gen;
+    Rcpp::RNGScope rcpp_rngScope_gen;
     Rcpp::traits::input_parameter< arma::uvec >::type cIndex(cIndexSEXP);
     Rcpp::traits::input_parameter< arma::mat >::type gen(genSEXP);
     Rcpp::traits::input_parameter< arma::vec >::type bpVec(bpVecSEXP);
-    __result = Rcpp::wrap(returnLD(cIndex, gen, bpVec));
-    return __result;
+    rcpp_result_gen = Rcpp::wrap(returnLD(cIndex, gen, bpVec));
+    return rcpp_result_gen;
 END_RCPP
 }
 // assoc
-arma::vec assoc(arma::mat gen, arma::colvec y);
+arma::mat assoc(arma::mat gen, arma::colvec y);
 RcppExport SEXP coRge_assoc(SEXP genSEXP, SEXP ySEXP) {
 BEGIN_RCPP
-    Rcpp::RObject __result;
-    Rcpp::RNGScope __rngScope;
+    Rcpp::RObject rcpp_result_gen;
+    Rcpp::RNGScope rcpp_rngScope_gen;
     Rcpp::traits::input_parameter< arma::mat >::type gen(genSEXP);
     Rcpp::traits::input_parameter< arma::colvec >::type y(ySEXP);
-    __result = Rcpp::wrap(assoc(gen, y));
-    return __result;
+    rcpp_result_gen = Rcpp::wrap(assoc(gen, y));
+    return rcpp_result_gen;
 END_RCPP
 }
 // fun
 arma::uvec fun(arma::field<std::string> input_field, std::string id);
 RcppExport SEXP coRge_fun(SEXP input_fieldSEXP, SEXP idSEXP) {
 BEGIN_RCPP
-    Rcpp::RObject __result;
-    Rcpp::RNGScope __rngScope;
+    Rcpp::RObject rcpp_result_gen;
+    Rcpp::RNGScope rcpp_rngScope_gen;
     Rcpp::traits::input_parameter< arma::field<std::string> >::type input_field(input_fieldSEXP);
     Rcpp::traits::input_parameter< std::string >::type id(idSEXP);
-    __result = Rcpp::wrap(fun(input_field, id));
-    return __result;
+    rcpp_result_gen = Rcpp::wrap(fun(input_field, id));
+    return rcpp_result_gen;
 END_RCPP
 }
 // th
 arma::vec th(Rcpp::StringVector strata_rsid, Rcpp::StringVector rsid, Rcpp::NumericVector r2);
 RcppExport SEXP coRge_th(SEXP strata_rsidSEXP, SEXP rsidSEXP, SEXP r2SEXP) {
 BEGIN_RCPP
-    Rcpp::RObject __result;
-    Rcpp::RNGScope __rngScope;
+    Rcpp::RObject rcpp_result_gen;
+    Rcpp::RNGScope rcpp_rngScope_gen;
     Rcpp::traits::input_parameter< Rcpp::StringVector >::type strata_rsid(strata_rsidSEXP);
     Rcpp::traits::input_parameter< Rcpp::StringVector >::type rsid(rsidSEXP);
     Rcpp::traits::input_parameter< Rcpp::NumericVector >::type r2(r2SEXP);
-    __result = Rcpp::wrap(th(strata_rsid, rsid, r2));
-    return __result;
+    rcpp_result_gen = Rcpp::wrap(th(strata_rsid, rsid, r2));
+    return rcpp_result_gen;
 END_RCPP
 }

src/lm.cpp

@@ -24,10 +24,10 @@ using namespace arma;
 
 //' @export
 // [[Rcpp::export]]
-arma::vec assoc(arma::mat gen, arma::colvec y){
+arma::mat assoc(arma::mat gen, arma::colvec y){
 
 	// Initialize output as a vector.
-	arma::vec outputList(gen.n_rows);
+	arma::mat outputList(gen.n_rows, 2);
 
 	// Split regression into each row of the gen file
 	// meaning each SNP is regressed seperately.
@@ -87,7 +87,8 @@ arma::vec assoc(arma::mat gen, arma::colvec y){
 
 		std::string name = patch::to_string(i);
 
-		outputList[i] = t;
+		outputList[i,0] = coef;
+		outputList[i,1] = stderrest;
 	}
 
 	return outputList;