about summary refs log tree commit diff
diff options
context:
space:
mode:
authorAlexander_Kabui2024-11-08 12:19:07 +0300
committerAlexander_Kabui2024-11-08 12:19:07 +0300
commit5a996dfe02abf9fa25563844662372f003dbd67e (patch)
tree8496b3e1e78de05c20d6482b948634882b16c83d
parent8c87de4d70795c07098608342d9fee1acc571346 (diff)
downloadgenenetwork3-5a996dfe02abf9fa25563844662372f003dbd67e.tar.gz
feat: init add functionality to estimate qtl effect.
-rw-r--r--scripts/rqtl2_wrapper.R45
1 files changed, 45 insertions, 0 deletions
diff --git a/scripts/rqtl2_wrapper.R b/scripts/rqtl2_wrapper.R
index 966d3c6..b68c82f 100644
--- a/scripts/rqtl2_wrapper.R
+++ b/scripts/rqtl2_wrapper.R
@@ -430,8 +430,50 @@ lod_peaks = find_peaks(
 
 lod_peaks
 cat("Generating the ouput data as  vector\n")
+
+# get the estimated qtl effects
+cat("Getting the estimated qtl effects\n")
+
+get_qtl_effect <- function(chromosome,geno_prob,pheno,covar=NULL,LOCO= NULL){
+     cat("Finding the qtl effect\n")
+     chr_Pr <- geno_prob[,chromosome]
+     if (!is.null(chr_Pr)){
+     str_glue("Find qtl effect for chromosome {chromosome} and pheno {pheno}")
+     if (!is.null(LOCO)) {
+          str_glue("Find qtl effect for chromosome {chromosome} and pheno {pheno} and LOCO {chromosome}")
+     kinship <- calc_kinship(chr_Pr, "loco")[[chromosome]]
+     return(scan1coef(chr_Pr, pheno, kinship, addcovar=covar))
+     }
+     else {
+      return(scan1coef(chr_Pr, pheno, addcovar=covar))
+     }
+     }
+    return(NULL)
+}
+
+
+
+# take the first phenotype in the dataset
+
+# grab phenotypes and covariates; ensure that covariates have names attribute
+# in this case we are using the first chromosome
+pheno <- dataset$pheno[,1]
+if (!is.null(dataset$covar) && !is.null(dataset$covar$sex)){
+ covar <- match(dataset$covar$sex, c("f", "m")) # make numeric
+ names(covar) <- rownames(dataset$covar)
+} else {
+covar  <- NULL
+}
+covar
+
+
+# TODO fix this hardcoded chromosome
+# TODO get all chromosomes and iterate coeff_results for a given chromosomes
+coeff_results  <- get_qtl_effect("5", Pr, pheno)
+
 output = list(lod_peaks = lod_peaks,
              scan_results =scan_results,
+	     genetic_probabilities = Pr,
 	     lod_significance = lod_significance,
 	     permutation_results = perm,
 	     lod_peaks = lod_peaks,
@@ -445,3 +487,6 @@ str_glue("The output file path is  {output_file_path}")
 cat("Writing to the output file\n")
 write(output_json_data, file=output_file_path)
 
+
+
+