diff options
Diffstat (limited to 'wqflask')
| -rw-r--r-- | wqflask/wqflask/marker_regression/gemma_mapping.py | 141 | ||||
| -rw-r--r-- | wqflask/wqflask/marker_regression/marker_regression.py | 6 | ||||
| -rw-r--r-- | wqflask/wqflask/marker_regression/marker_regression_gn1.py | 4 | ||||
| -rw-r--r-- | wqflask/wqflask/templates/marker_regression_gn1.html | 1 | ||||
| -rw-r--r-- | wqflask/wqflask/templates/show_trait_mapping_tools.html | 16 | ||||
| -rw-r--r-- | wqflask/wqflask/views.py | 2 |
6 files changed, 131 insertions, 39 deletions
diff --git a/wqflask/wqflask/marker_regression/gemma_mapping.py b/wqflask/wqflask/marker_regression/gemma_mapping.py index 31e94266..f23c390e 100644 --- a/wqflask/wqflask/marker_regression/gemma_mapping.py +++ b/wqflask/wqflask/marker_regression/gemma_mapping.py @@ -1,14 +1,14 @@ -import os, math +import os, math, string, random, json 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,62 +18,88 @@ 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(len(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) - if method == "gemma": - #gemma_command = GEMMA_COMMAND + ' -bfile %s/%s -k %s/%s.sXX.txt -lmm 1 -maf 0.1 -c %s/%s_covariates.txt -o %s_output' % (flat_files('mapping'), - gemma_command = GEMMA_COMMAND + ' -bfile %s/%s -k %s/%s.cXX.txt -lmm 1 -maf 0.1 -c %s/%s_covariates.txt -outdir %s -o %s_output' % (flat_files('mapping'), - this_dataset.group.name, - flat_files('mapping'), - this_dataset.group.name, - flat_files('mapping'), - this_dataset.group.name, - webqtlConfig.GENERATED_IMAGE_DIR, - this_dataset.group.name) - # use GEMMA_RUN in the next one, create a unique temp file + + if method == "gemma": + gemma_command = GEMMA_COMMAND + ' -bfile %s/%s -k %s/%s.cXX.txt -lmm 1 -maf 0.1' % (flat_files('mapping'), + this_dataset.group.name, + flat_files('mapping'), + this_dataset.group.name) + if covariates != "": + gemma_command += ' -c %s/%s_covariates.txt -outdir %s -o %s_output' % (flat_files('mapping'), + this_dataset.group.name, + webqtlConfig.GENERATED_IMAGE_DIR, + this_dataset.group.name) 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'), + #gemma_command = GEMMA_COMMAND + ' -bfile %s/%s -k %s/%s.sXX.txt -lmm 1 -maf 0.1 -o %s_output' % (flat_files('mapping'), + gemma_command += ' -outdir %s -o %s_output' % (webqtlConfig.GENERATED_IMAGE_DIR, + this_dataset.group.name) + else: + if use_loco == "True": + k_output_filename = this_dataset.group.name + "_K_" + ''.join(random.choice(string.ascii_uppercase + string.digits) for _ in range(6)) + 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/%s.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, + k_output_filename) + logger.debug("k_command:" + generate_k_command) + os.system(generate_k_command) + + gemma_command = GEMMA_WRAPPER_COMMAND + ' --json --loco --input /home/zas1024/tmp/gn2/%s.json -- -g %s/%s_geno.txt -p %s/%s_pheno.txt' % (k_output_filename, 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'), - gemma_command = GEMMA_COMMAND + ' -bfile %s/%s -k %s/%s.cXX.txt -lmm 1 -maf 0.1 -outdir %s -o %s_output' % (flat_files('mapping'), - this_dataset.group.name, - flat_files('mapping'), - this_dataset.group.name, - webqtlConfig.GENERATED_IMAGE_DIR, + flat_files('genotype/bimbam'), this_dataset.group.name) + + gwa_output_filename = this_dataset.group.name + "_GWA_" + ''.join(random.choice(string.ascii_uppercase + string.digits) for _ in range(6)) + if covariates != "": + gemma_command += ' -c %s/%s_covariates.txt -a %s/%s_snps.txt -lmm 1 -maf 0.1 -debug > /home/zas1024/tmp/gn2/%s.json' % (flat_files('mapping'), + this_dataset.group.name, + flat_files('genotype/bimbam'), + this_dataset.group.name, + gwa_output_filename) + else: + gemma_command += ' -a %s/%s_snps.txt -lmm 1 -maf 0.1 -debug > /home/zas1024/tmp/gn2/%s.json' % (flat_files('genotype/bimbam'), + this_dataset.group.name, + gwa_output_filename) + 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'), + 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' % (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, - webqtlConfig.GENERATED_IMAGE_DIR, this_dataset.group.name) - logger.debug("gemma_command:" + gemma_command) + if covariates != "": + gemma_command += ' -c %s/%s_covariates.txt -outdir %s -debug -o %s_output' % (flat_files('mapping'), + this_dataset.group.name, + webqtlConfig.GENERATED_IMAGE_DIR, + this_dataset.group.name) + else: + gemma_command += ' -outdir %s -debug -o %s_output' % (webqtlConfig.GENERATED_IMAGE_DIR, + this_dataset.group.name) + + logger.debug("gemma_command:" + gemma_command) os.system(gemma_command) - marker_obs = parse_gemma_output(this_dataset) + if use_loco == "True": + marker_obs = parse_loco_output(this_dataset, gwa_output_filename) + else: + marker_obs = parse_gemma_output(this_dataset) return marker_obs @@ -172,3 +198,46 @@ def parse_gemma_output(this_dataset): p_values.append(float(line.split("\t")[10])) return marker_obs + +def parse_loco_output(this_dataset, gwa_output_filename): + + output_filelist = [] + with open("/home/zas1024/tmp/gn2/" + gwa_output_filename + ".json") as data_file: + data = json.load(data_file) + + files = data['files'] + for file in files: + output_filelist.append(file[2]) + + included_markers = [] + p_values = [] + marker_obs = [] + previous_chr = 0 + + for this_file in output_filelist: + #with open("/home/zas1024/gene/wqflask/output/{}_output.assoc.txt".format(this_dataset.group.name)) as output_file: + with open(this_file) as output_file: + for line in output_file: + if line.startswith("chr"): + continue + else: + marker = {} + marker['name'] = line.split("\t")[1] + if line.split("\t")[0] != "X" and line.split("\t")[0] != "X/Y": + marker['chr'] = int(line.split("\t")[0]) + else: + marker['chr'] = line.split("\t")[0] + marker['Mb'] = float(line.split("\t")[2]) / 1000000 + marker['p_value'] = float(line.split("\t")[10]) + if math.isnan(marker['p_value']) or (marker['p_value'] <= 0): + marker['lod_score'] = 0 + #marker['lrs_value'] = 0 + else: + marker['lod_score'] = -math.log10(marker['p_value']) + #marker['lrs_value'] = -math.log10(marker['p_value']) * 4.61 + marker_obs.append(marker) + + included_markers.append(line.split("\t")[1]) + p_values.append(float(line.split("\t")[10])) + + return marker_obs \ No newline at end of file diff --git a/wqflask/wqflask/marker_regression/marker_regression.py b/wqflask/wqflask/marker_regression/marker_regression.py index 55bbacac..9f2f06f3 100644 --- a/wqflask/wqflask/marker_regression/marker_regression.py +++ b/wqflask/wqflask/marker_regression/marker_regression.py @@ -84,6 +84,10 @@ class MarkerRegression(object): self.manhattan_plot = False self.maf = start_vars['maf'] # Minor allele frequency + if "use_loco" in start_vars: + self.use_loco = start_vars['use_loco'] + else: + self.use_loco = None self.suggestive = "" self.significant = "" self.pair_scan = False # Initializing this since it is checked in views to determine which template to use @@ -158,7 +162,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/marker_regression/marker_regression_gn1.py b/wqflask/wqflask/marker_regression/marker_regression_gn1.py index 93d75a03..687cfeb5 100644 --- a/wqflask/wqflask/marker_regression/marker_regression_gn1.py +++ b/wqflask/wqflask/marker_regression/marker_regression_gn1.py @@ -257,6 +257,10 @@ class MarkerRegression(object): self.controlLocus = "" if 'covariates' in start_vars.keys(): self.covariates = start_vars['covariates'] + if 'maf' in start_vars.keys(): + self.maf = start_vars['maf'] + if 'use_loco' in start_vars.keys(): + self.use_loco = start_vars['use_loco'] #try: self.selectedChr = int(start_vars['selected_chr']) diff --git a/wqflask/wqflask/templates/marker_regression_gn1.html b/wqflask/wqflask/templates/marker_regression_gn1.html index c6c6bc23..5c457275 100644 --- a/wqflask/wqflask/templates/marker_regression_gn1.html +++ b/wqflask/wqflask/templates/marker_regression_gn1.html @@ -20,6 +20,7 @@ <input type="hidden" name="value:{{ sample }}" value="{{ vals[loop.index - 1] }}"> {% endfor %} <input type="hidden" name="maf"> + <input type="hidden" name="use_loco" value="{{ use_loco }}"> <input type="hidden" name="selected_chr" value="{{ selectedChr }}"> <input type="hidden" name="manhattan_plot" value="{{ manhattan_plot }}"> <input type="hidden" name="num_perm" value="{{ nperm }}"> 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', |