Committing current progress towards implementing LOCO for GEMMA, not done yet

4 files changed, 81 insertions, 21 deletions

diff --git a/wqflask/wqflask/marker_regression/gemma_mapping.py b/wqflask/wqflask/marker_regression/gemma_mapping.py
index 31e94266..2c99d1fc 100644
--- a/wqflask/wqflask/marker_regression/gemma_mapping.py
+++ b/wqflask/wqflask/marker_regression/gemma_mapping.py
@@ -3,12 +3,12 @@ import os, math
 from base import webqtlConfig
 from base.trait import GeneralTrait
 from base.data_set import create_dataset
-from utility.tools import flat_files, GEMMA_COMMAND
+from utility.tools import flat_files, GEMMA_COMMAND, GEMMA_WRAPPER_COMMAND
 
 import utility.logger
 logger = utility.logger.getLogger(__name__ )
 
-def run_gemma(this_dataset, samples, vals, covariates, method):
+def run_gemma(this_dataset, samples, vals, covariates, method, use_loco):
     """Generates p-values for each marker using GEMMA"""
 
     print("INSIDE GEMMA_MAPPING")
@@ -18,7 +18,13 @@ def run_gemma(this_dataset, samples, vals, covariates, method):
     if not os.path.isfile("{}{}_output.assoc.txt".format(webqtlConfig.GENERATED_IMAGE_DIR, this_dataset.group.name)):
         open("{}{}_output.assoc.txt".format(webqtlConfig.GENERATED_IMAGE_DIR, this_dataset.group.name), "w+")
 
-    logger.debug("COVARIATES_GEMMA:", covariates)
+    this_chromosomes = this_dataset.species.chromosomes.chromosomes
+    chr_list_string = ""
+    for i in range(this_chromosomes):
+        if i < (len(this_chromosomes) - 1):
+            chr_list_string += this_chromosomes[i+1].name + ","
+        else:
+            chr_list_string += this_chromosomes[i+1].name
 
     if covariates != "":
         gen_covariates_file(this_dataset, covariates)
@@ -35,19 +41,38 @@ def run_gemma(this_dataset, samples, vals, covariates, method):
             # use GEMMA_RUN in the next one, create a unique temp file
         else:
             logger.debug("FLAT FILES:", flat_files('mapping'))
-            #gemma_command = GEMMA_COMMAND + ' -g %s/%s_geno.txt -p %s/%s_pheno.txt -a %s/%s_snps.txt -k %s/%s.sXX.txt -lmm 1 -maf 0.1 -c %s/%s_covariates.txt -o %s_output' % (flat_files('genotype/bimbam'),
-            gemma_command = GEMMA_COMMAND + ' -g %s/%s_geno.txt -p %s/%s_pheno.txt -a %s/%s_snps.txt -k %s/%s.cXX.txt -lmm 1 -maf 0.1 -c %s/%s_covariates.txt -outdir %s -debug -o %s_output' % (flat_files('genotype/bimbam'),
-                                                                                            this_dataset.group.name,
-                                                                                            flat_files('genotype/bimbam'),
-                                                                                            this_dataset.group.name,
-                                                                                            flat_files('genotype/bimbam'),
-                                                                                            this_dataset.group.name,
-                                                                                            flat_files('genotype/bimbam'),
-                                                                                            this_dataset.group.name,
-                                                                                            flat_files('mapping'),
-                                                                                            this_dataset.group.name,
-                                                                                            webqtlConfig.GENERATED_IMAGE_DIR,
-                                                                                            this_dataset.group.name)
+            if use_loco:
+                generate_k_command = GEMMA_WRAPPER_COMMAND + ' --json --loco ' + chr_list_string + ' -- -g %s/%s_geno.txt -p %s/%s_pheno.txt -a %s/%s_snps.txt -gk -debug > /home/zas1024/tmp/gn2/K.json' % (flat_files('genotype/bimbam'),
+                                                                                                this_dataset.group.name,
+                                                                                                flat_files('genotype/bimbam'),
+                                                                                                this_dataset.group.name,
+                                                                                                flat_files('genotype/bimbam'),
+                                                                                                this_dataset.group.name)
+                logger.debug("k_command:" + generate_k_command)
+                os.system(generate_k_command)
+
+                run_loco_command = GEMMA_WRAPPER_COMMAND + ' --json --loco --input /home/zas1024/tmp/gn2/K.json -- -g %s/%s_geno.txt -p %s/%s_pheno.txt -c %s/%s_covariates.txt -a %s/%s_snps.txt -lmm 1 -maf 0.1 -debug > /home/zas1024/tmp/gn2/GWA.json' % (flat_files('genotype/bimbam'),
+                                                                                                this_dataset.group.name,
+                                                                                                flat_files('genotype/bimbam'),
+                                                                                                this_dataset.group.name,
+                                                                                                flat_files('mapping'),
+                                                                                                this_dataset.group.name,
+                                                                                                flat_files('genotype/bimbam'),
+                                                                                                this_dataset.group.name)
+                os.system(run_loco_command)
+            else:
+                gemma_command = GEMMA_COMMAND + ' -g %s/%s_geno.txt -p %s/%s_pheno.txt -a %s/%s_snps.txt -k %s/%s.cXX.txt -lmm 1 -maf 0.1 -c %s/%s_covariates.txt -outdir %s -debug -o %s_output' % (flat_files('genotype/bimbam'),
+                                                                                                this_dataset.group.name,
+                                                                                                flat_files('genotype/bimbam'),
+                                                                                                this_dataset.group.name,
+                                                                                                flat_files('genotype/bimbam'),
+                                                                                                this_dataset.group.name,
+                                                                                                flat_files('genotype/bimbam'),
+                                                                                                this_dataset.group.name,
+                                                                                                flat_files('mapping'),
+                                                                                                this_dataset.group.name,
+                                                                                                webqtlConfig.GENERATED_IMAGE_DIR,
+                                                                                                this_dataset.group.name)
     else:
         if method == "gemma":
             #gemma_command = GEMMA_COMMAND + ' -bfile %s/%s -k %s/%s.sXX.txt -lmm 1 -maf 0.1 -o %s_output' % (flat_files('mapping'),
@@ -58,8 +83,28 @@ def run_gemma(this_dataset, samples, vals, covariates, method):
                                                                                             webqtlConfig.GENERATED_IMAGE_DIR,
                                                                                             this_dataset.group.name)
         else:
-            #gemma_command = GEMMA_COMMAND + ' -g %s/%s_geno.txt -p %s/%s_pheno.txt -a %s/%s_snps.txt -k %s/%s.sXX.txt -lmm 1 -maf 0.1 -o %s_output' % (flat_files('genotype/bimbam'),
-            gemma_command = GEMMA_COMMAND + ' -g %s/%s_geno.txt -p %s/%s_pheno.txt -a %s/%s_snps.txt -k %s/%s.cXX.txt -lmm 1 -maf 0.1 -outdir %s -debug -o %s_output' % (flat_files('genotype/bimbam'),
+            if use_loco:
+                generate_k_command = GEMMA_WRAPPER_COMMAND + ' --json --loco ' + chr_list_string + ' -- -g %s/%s_geno.txt -p %s/%s_pheno.txt -a %s/%s_snps.txt -gk -debug > /home/zas1024/tmp/gn2/K.json' % (flat_files('genotype/bimbam'),
+                                                                                                this_dataset.group.name,
+                                                                                                flat_files('genotype/bimbam'),
+                                                                                                this_dataset.group.name,
+                                                                                                flat_files('genotype/bimbam'),
+                                                                                                this_dataset.group.name)
+                logger.debug("k_command:" + generate_k_command)
+                os.system(generate_k_command)
+
+                run_loco_command = GEMMA_WRAPPER_COMMAND + ' --json --loco --input /home/zas1024/tmp/gn2/K.json -- -g %s/%s_geno.txt -p %s/%s_pheno.txt -c %s/%s_covariates.txt -a %s/%s_snps.txt -lmm 1 -maf 0.1 -debug > /home/zas1024/tmp/gn2/GWA.json' % (flat_files('genotype/bimbam'),
+                                                                                                this_dataset.group.name,
+                                                                                                flat_files('genotype/bimbam'),
+                                                                                                this_dataset.group.name,
+                                                                                                flat_files('mapping'),
+                                                                                                this_dataset.group.name,
+                                                                                                flat_files('genotype/bimbam'),
+                                                                                                this_dataset.group.name)
+                os.system(run_loco_command)
+            else:
+                #gemma_command = GEMMA_COMMAND + ' -g %s/%s_geno.txt -p %s/%s_pheno.txt -a %s/%s_snps.txt -k %s/%s.sXX.txt -lmm 1 -maf 0.1 -o %s_output' % (flat_files('genotype/bimbam'),
+                gemma_command = GEMMA_COMMAND + ' -g %s/%s_geno.txt -p %s/%s_pheno.txt -a %s/%s_snps.txt -k %s/%s.cXX.txt -lmm 1 -maf 0.1 -outdir %s -debug -o %s_output' % (flat_files('genotype/bimbam'),
                                                                                             this_dataset.group.name,
                                                                                             flat_files('genotype/bimbam'),
                                                                                             this_dataset.group.name,
diff --git a/wqflask/wqflask/marker_regression/marker_regression.py b/wqflask/wqflask/marker_regression/marker_regression.py
index 55bbacac..24292c35 100644
--- a/wqflask/wqflask/marker_regression/marker_regression.py
+++ b/wqflask/wqflask/marker_regression/marker_regression.py
@@ -84,6 +84,7 @@ class MarkerRegression(object):
             self.manhattan_plot = False
 
         self.maf = start_vars['maf'] # Minor allele frequency
+        self.use_loco = start_vars['use_loco']
         self.suggestive = ""
         self.significant = ""
         self.pair_scan = False # Initializing this since it is checked in views to determine which template to use
@@ -158,7 +159,7 @@ class MarkerRegression(object):
             self.score_type = "-log(p)"
             self.manhattan_plot = True
             with Bench("Running GEMMA"):
-                marker_obs = gemma_mapping.run_gemma(self.dataset, self.samples, self.vals, self.covariates, self.mapping_method)
+                marker_obs = gemma_mapping.run_gemma(self.dataset, self.samples, self.vals, self.covariates, self.mapping_method, self.use_loco)
             results = marker_obs
         elif self.mapping_method == "rqtl_plink":
             results = self.run_rqtl_plink()
diff --git a/wqflask/wqflask/templates/show_trait_mapping_tools.html b/wqflask/wqflask/templates/show_trait_mapping_tools.html
index c24a5853..eadc7481 100644
--- a/wqflask/wqflask/templates/show_trait_mapping_tools.html
+++ b/wqflask/wqflask/templates/show_trait_mapping_tools.html
@@ -302,8 +302,20 @@
                                 <input name="maf_gemma" value="0.01" type="text" class="form-control">
                             </div>
                         </div>
+                        <div class="mapping_method_fields form-group">
+                            <label class="col-xs-4 control-label">Use LOCO</label>
+                            <div style="margin-left: 20px;" class="col-xs-4 controls">
+                                <label class="radio-inline">
+                                    <input type="radio" name="use_loco" value="True">
+                                    Yes
+                                </label>
+                                <label class="radio-inline">
+                                    <input type="radio" name="use_loco" value="False" checked="">
+                                    No
+                               </label>
+                            </div>
+                        </div>
                     </div>
-
                     <div style="padding-top: 5px; padding-bottom: 5px; padding-left: 20px;" class="form-horizontal">
                         <div class="mapping_method_fields form-group">
                             <button type="button" id="select_covariates" class="btn btn-default">
@@ -395,7 +407,7 @@
             <dt>GEMMA</dt>
             <dd>
                 GEMMA is software implementing the Genome-wide Efficient Mixed Model Association algorithm for a standard linear mixed model for genome-wide association studies (GWAS).<br>
-                <font style="color: red;">Note: There currently exists an issue where GEMMA can't be run on traits from the same group simultaneously. If you recieve an error, please wait a few minutes and try again.</font>
+                Note: There currently exists an issue where GEMMA can't be run on traits from the same group simultaneously. If you receive an error, please wait a few minutes and try again.
             </dd>
             <dt>PLINK</dt>
             <dd>PLINK is a free, open-source whole genome association analysis toolset.</dd>
diff --git a/wqflask/wqflask/views.py b/wqflask/wqflask/views.py
index 2d4fd0f2..f340c954 100644
--- a/wqflask/wqflask/views.py
+++ b/wqflask/wqflask/views.py
@@ -516,6 +516,7 @@ def loading_page():
         'LRSCheck',
         'covariates',
         'maf',
+        'use_loco',
         'manhattan_plot',
         'control_marker',
         'control_marker_db',
@@ -572,6 +573,7 @@ def marker_regression_page():
         'LRSCheck',
         'covariates',
         'maf',
+        'use_loco',
         'manhattan_plot',
         'control_marker',
         'control_marker_db',