diff options
| -rw-r--r-- | wqflask/wqflask/marker_regression/gemma_mapping.py | 65 |
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 |