diff --git a/R/RcppExports.R b/R/RcppExports.R index dab9c08..81e2ac6 100644 --- a/R/RcppExports.R +++ b/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 diff --git a/R/correct.R b/R/correct.R index d644a95..15b5b29 100644 --- a/R/correct.R +++ b/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) diff --git a/R/gen_phen_df.R b/R/gen_phen_df.R index a06fee1..d13534c 100644 --- a/R/gen_phen_df.R +++ b/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)) diff --git a/TODO.md b/TODO.md new file mode 100644 index 0000000..38837f0 --- /dev/null +++ b/TODO.md @@ -0,0 +1,4 @@ +## TODO + +- Simulate one and look at MAF and other things +- Change the referene panel. diff --git a/inst/analysis/new_ref/makefile b/inst/analysis/new_ref/makefile new file mode 100644 index 0000000..7960d89 --- /dev/null +++ b/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 diff --git a/inst/analysis/new_ref/sub.sh b/inst/analysis/new_ref/sub.sh new file mode 100644 index 0000000..f61a9ae --- /dev/null +++ b/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 diff --git a/inst/bash/test.txt b/inst/bash/test.txt new file mode 100644 index 0000000..bc9006f --- /dev/null +++ b/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 diff --git a/inst/bash/test2.txt b/inst/bash/test2.txt new file mode 100644 index 0000000..e69de29 diff --git a/inst/data b/inst/data index 26da8db..4dfd1c7 160000 --- a/inst/data +++ b/inst/data @@ -1 +1 @@ -Subproject commit 26da8db5dfc54ac8168d37704d8a9cc39461c6b3 +Subproject commit 4dfd1c72f2fc29beb4364d39834498ff65f2aabe diff --git a/src/RcppExports.cpp b/src/RcppExports.cpp index 26ad132..dc40ea9 100644 --- a/src/RcppExports.cpp +++ b/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 @@ -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 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 >::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 } diff --git a/src/lm.cpp b/src/lm.cpp index 2cb1600..9dc2f3a 100644 --- a/src/lm.cpp +++ b/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;