aboutsummaryrefslogtreecommitdiff
diff options
context:
space:
mode:
-rw-r--r--wqflask/wqflask/marker_regression/gemma_mapping.py65
1 files changed, 54 insertions, 11 deletions
diff --git a/wqflask/wqflask/marker_regression/gemma_mapping.py b/wqflask/wqflask/marker_regression/gemma_mapping.py
index c17f21aa..ffe26190 100644
--- a/wqflask/wqflask/marker_regression/gemma_mapping.py
+++ b/wqflask/wqflask/marker_regression/gemma_mapping.py
@@ -36,8 +36,9 @@ def run_gemma(this_dataset, samples, vals, covariates, use_loco, maf=0.01):
if covariates != "":
gen_covariates_file(this_dataset, covariates)
+ k_output_filename = this_dataset.group.name + "_K_" + ''.join(random.choice(string.ascii_uppercase + string.digits) for _ in range(6))
+ gwa_output_filename = this_dataset.group.name + "_GWA_" + ''.join(random.choice(string.ascii_uppercase + string.digits) for _ in range(6))
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 + ' -- ' + GEMMAOPTS + ' -g %s/%s_geno.txt -p %s/%s_pheno.txt -a %s/%s_snps.txt -gk > %s/gn2/%s.json' % (flat_files('genotype/bimbam'),
genofile_name,
flat_files('genotype/bimbam'),
@@ -53,32 +54,58 @@ def run_gemma(this_dataset, samples, vals, covariates, use_loco, maf=0.01):
genofile_name,
flat_files('genotype/bimbam'),
genofile_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 2 -maf 0.1 > %s/gn2/%s.json' % (flat_files('mapping'),
- this_dataset.group.name,
- flat_files('genotype/bimbam'),
- genofile_name,
- TEMPDIR,
- gwa_output_filename)
+ gemma_command += ' -c %s/%s_covariates.txt -a %s/%s_snps.txt -lmm 2 -maf %s > %s/gn2/%s.json' % (flat_files('mapping'),
+ this_dataset.group.name,
+ flat_files('genotype/bimbam'),
+ genofile_name,
+ maf,
+ TEMPDIR,
+ gwa_output_filename)
else:
- gemma_command += ' -a %s/%s_snps.txt -lmm 2 -maf 0.1 > %s/gn2/%s.json' % (flat_files('genotype/bimbam'),
+ gemma_command += ' -a %s/%s_snps.txt -lmm 2 -maf %s > %s/gn2/%s.json' % (flat_files('genotype/bimbam'),
genofile_name,
+ maf,
TEMPDIR,
gwa_output_filename)
else:
- gemma_command = GEMMA_COMMAND + ' ' + GEMMAOPTS + ' -g %s/%s_geno.txt -p %s/%s_pheno.txt -a %s/%s_snps.txt -k %s/%s.cXX.txt -lmm 2 -maf %s' % (flat_files('genotype/bimbam'),
+ generate_k_command = GEMMA_COMMAND + ' ' + GEMMAOPTS + ' -g %s/%s_geno.txt -p %s/%s_pheno.txt -a %s/%s_snps.txt -gk -outdir %s/gn2/ -o %s' % (flat_files('genotype/bimbam'),
+ genofile_name,
+ flat_files('genotype/bimbam'),
genofile_name,
flat_files('genotype/bimbam'),
genofile_name,
+ TEMPDIR,
+ k_output_filename)
+ #generate_k_command = GEMMA_WRAPPER_COMMAND + ' --json -- ' + GEMMAOPTS + ' -g %s/%s_geno.txt -p %s/%s_pheno.txt -a %s/%s_snps.txt -gk > %s/gn2/%s.json' % (flat_files('genotype/bimbam'),
+ # genofile_name,
+ # flat_files('genotype/bimbam'),
+ # genofile_name,
+ # flat_files('genotype/bimbam'),
+ # genofile_name,
+ # TEMPDIR,
+ # k_output_filename)
+
+ logger.debug("k_command:" + generate_k_command)
+ os.system(generate_k_command)
+
+ gemma_command = GEMMA_COMMAND + ' ' + GEMMAOPTS + ' -g %s/%s_geno.txt -p %s/%s_pheno.txt -a %s/%s_snps.txt -k %s/gn2/%s.cXX.txt -lmm 2 -maf %s' % (flat_files('genotype/bimbam'),
+ genofile_name,
flat_files('genotype/bimbam'),
genofile_name,
flat_files('genotype/bimbam'),
genofile_name,
+ TEMPDIR,
+ k_output_filename,
maf)
+ #gemma_command = GEMMA_WRAPPER_COMMAND + ' --json --input %s/gn2/%s.json -- ' % (TEMPDIR, k_output_filename) + GEMMAOPTS + ' -g %s/%s_geno.txt -p %s/%s_pheno.txt' % (flat_files('genotype/bimbam'),
+ # genofile_name,
+ # flat_files('genotype/bimbam'),
+ # genofile_name)
+
+
if covariates != "":
gemma_command += ' -c %s/%s_covariates.txt -outdir %s -o %s_output' % (flat_files('mapping'),
this_dataset.group.name,
@@ -88,12 +115,28 @@ def run_gemma(this_dataset, samples, vals, covariates, use_loco, maf=0.01):
gemma_command += ' -outdir %s -o %s_output' % (webqtlConfig.GENERATED_IMAGE_DIR,
genofile_name)
+ #if covariates != "":
+ # gemma_command += ' -c %s/%s_covariates.txt -a %s/%s_snps.txt -lmm 2 -maf %s > %s/gn2/%s.json' % (flat_files('mapping'),
+ # this_dataset.group.name,
+ # flat_files('genotype/bimbam'),
+ # genofile_name,
+ # maf,
+ # TEMPDIR,
+ # gwa_output_filename)
+ #else:
+ # gemma_command += ' -a %s/%s_snps.txt -lmm 2 -maf %s > %s/gn2/%s.json' % (flat_files('genotype/bimbam'),
+ # genofile_name,
+ # maf,
+ # TEMPDIR,
+ # gwa_output_filename)
+
logger.debug("gemma_command:" + gemma_command)
os.system(gemma_command)
if use_loco == "True":
marker_obs = parse_loco_output(this_dataset, gwa_output_filename)
else:
+ #marker_obs = parse_loco_output(this_dataset, gwa_output_filename)
marker_obs = parse_gemma_output(genofile_name)
return marker_obs