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
