diff options
Diffstat (limited to 'gn2/wqflask/heatmap/heatmap.py')
| -rw-r--r-- | gn2/wqflask/heatmap/heatmap.py | 191 |
1 files changed, 191 insertions, 0 deletions
diff --git a/gn2/wqflask/heatmap/heatmap.py b/gn2/wqflask/heatmap/heatmap.py new file mode 100644 index 00000000..c2dd55bd --- /dev/null +++ b/gn2/wqflask/heatmap/heatmap.py @@ -0,0 +1,191 @@ +import string +import os +import random +from gn2.base import species +from gn2.base import webqtlConfig +from gn2.utility import helper_functions + +from gn2.utility.tools import flat_files, REAPER_COMMAND, TEMPDIR +from redis import Redis +from flask import Flask, g + +from gn2.wqflask.database import database_connection +from gn2.utility.tools import get_setting + +Redis = Redis() + + +class Heatmap: + + def __init__(self, db_cursor, start_vars, temp_uuid): + trait_db_list = [trait.strip() + for trait in start_vars['trait_list'].split(',')] + helper_functions.get_trait_db_obs(self, trait_db_list) + + self.temp_uuid = temp_uuid + self.num_permutations = 5000 + self.dataset = self.trait_list[0][1] + + self.json_data = {} # The dictionary that will be used to create the json object that contains all the data needed to create the figure + + self.all_sample_list = [] + self.traits = [] + + chrnames = [] + self.species = species.TheSpecies(dataset=self.trait_list[0][1]) + + with database_connection(get_setting("SQL_URI")) as conn, conn.cursor() as db_cursor: + for this_chr in self.species.chromosomes.chromosomes(db_cursor): + chrnames.append([self.species.chromosomes.chromosomes(db_cursor)[this_chr].name, + self.species.chromosomes.chromosomes(db_cursor)[this_chr].mb_length]) + + for trait_db in self.trait_list: + + this_trait = trait_db[0] + self.traits.append(this_trait.name) + this_sample_data = this_trait.data + + for sample in this_sample_data: + if sample not in self.all_sample_list: + self.all_sample_list.append(sample) + + self.sample_data = [] + for trait_db in self.trait_list: + this_trait = trait_db[0] + this_sample_data = this_trait.data + + this_trait_vals = [] + for sample in self.all_sample_list: + if sample in this_sample_data: + this_trait_vals.append(this_sample_data[sample].value) + else: + this_trait_vals.append('') + self.sample_data.append(this_trait_vals) + + self.gen_reaper_results() + + lodnames = [] + chr_pos = [] + pos = [] + markernames = [] + + for trait in list(self.trait_results.keys()): + lodnames.append(trait) + + self.dataset.group.get_markers() + for marker in self.dataset.group.markers.markers: + chr_pos.append(marker['chr']) + pos.append(marker['Mb']) + markernames.append(marker['name']) + + self.json_data['chrnames'] = chrnames + self.json_data['lodnames'] = lodnames + self.json_data['chr'] = chr_pos + self.json_data['pos'] = pos + self.json_data['markernames'] = markernames + + for trait in self.trait_results: + self.json_data[trait] = self.trait_results[trait] + + self.js_data = dict( + json_data=self.json_data + ) + + def gen_reaper_results(self): + self.trait_results = {} + for trait_db in self.trait_list: + self.dataset.group.get_markers() + this_trait = trait_db[0] + + genotype = self.dataset.group.read_genotype_file(use_reaper=False) + samples, values, variances, sample_aliases = this_trait.export_informative() + + if self.dataset.group.genofile != None: + genofile_name = self.dataset.group.genofile[:-5] + else: + genofile_name = self.dataset.group.name + + trimmed_samples = [] + trimmed_values = [] + for i in range(0, len(samples)): + if samples[i] in self.dataset.group.samplelist: + trimmed_samples.append(str(samples[i])) + trimmed_values.append(values[i]) + + trait_filename = str(this_trait.name) + "_" + \ + str(self.dataset.name) + "_pheno" + gen_pheno_txt_file(trimmed_samples, trimmed_values, trait_filename) + + output_filename = self.dataset.group.name + "_GWA_" + \ + ''.join(random.choice(string.ascii_uppercase + string.digits) + for _ in range(6)) + + reaper_command = REAPER_COMMAND + ' --geno {0}/{1}.geno --traits {2}/gn2/{3}.txt -n 1000 -o {4}{5}.txt'.format(flat_files('genotype'), + genofile_name, + TEMPDIR, + trait_filename, + webqtlConfig.GENERATED_IMAGE_DIR, + output_filename) + + os.system(reaper_command) + + reaper_results = parse_reaper_output(output_filename) + + lrs_values = [float(qtl['lrs_value']) for qtl in reaper_results] + + self.trait_results[this_trait.name] = [] + for qtl in reaper_results: + if qtl['additive'] > 0: + self.trait_results[this_trait.name].append( + -float(qtl['lrs_value'])) + else: + self.trait_results[this_trait.name].append( + float(qtl['lrs_value'])) + + +def gen_pheno_txt_file(samples, vals, filename): + """Generates phenotype file for GEMMA""" + + with open("{0}/gn2/{1}.txt".format(TEMPDIR, filename), "w") as outfile: + outfile.write("Trait\t") + + filtered_sample_list = [] + filtered_vals_list = [] + for i, sample in enumerate(samples): + if vals[i] != "x": + filtered_sample_list.append(sample) + filtered_vals_list.append(str(vals[i])) + + samples_string = "\t".join(filtered_sample_list) + outfile.write(samples_string + "\n") + outfile.write("T1\t") + values_string = "\t".join(filtered_vals_list) + outfile.write(values_string) + + +def parse_reaper_output(gwa_filename): + included_markers = [] + p_values = [] + marker_obs = [] + + with open("{}{}.txt".format(webqtlConfig.GENERATED_IMAGE_DIR, gwa_filename)) as output_file: + for line in output_file: + if line.startswith("ID\t"): + continue + else: + marker = {} + marker['name'] = line.split("\t")[1] + try: + marker['chr'] = int(line.split("\t")[2]) + except: + marker['chr'] = line.split("\t")[2] + marker['cM'] = float(line.split("\t")[3]) + marker['Mb'] = float(line.split("\t")[4]) + if float(line.split("\t")[7]) != 1: + marker['p_value'] = float(line.split("\t")[7]) + marker['lrs_value'] = float(line.split("\t")[5]) + marker['lod_score'] = marker['lrs_value'] / 4.61 + marker['additive'] = float(line.split("\t")[6]) + marker_obs.append(marker) + + return marker_obs |