from typing import Dict
import string
import datetime
import uuid
import requests
import json as json
from collections import OrderedDict
import numpy as np
import scipy.stats as ss
from gn2.wqflask.database import database_connection
from gn2.base import webqtlConfig
from gn2.wqflask.show_trait.SampleList import SampleList
from gn2.base.trait import create_trait
from gn2.base import data_set
from gn2.utility import helper_functions
from gn2.utility.tools import get_setting, locate_ignore_error
from gn2.utility.tools import GN_PROXY_URL
from gn2.utility.redis_tools import get_redis_conn, get_resource_id
from gn3.authentication import get_highest_user_access_role
Redis = get_redis_conn()
ONE_YEAR = 60 * 60 * 24 * 365
###############################################
#
# Todo: Put in security to ensure that user has permission to access
# confidential data sets And add i.p.limiting as necessary
#
##############################################
class ShowTrait:
def __init__(self, db_cursor, user_id, kw):
if 'trait_id' in kw and kw['dataset'] != "Temp":
self.temp_trait = False
self.trait_id = kw['trait_id']
helper_functions.get_species_dataset_trait(self, kw)
self.resource_id = get_resource_id(self.dataset,
self.trait_id)
elif 'group' in kw:
self.temp_trait = True
self.trait_id = "Temp_" + kw['species'] + "_" + kw['group'] + \
"_" + datetime.datetime.now().strftime("%m%d%H%M%S")
self.temp_species = kw['species']
self.temp_group = kw['group']
self.dataset = data_set.create_dataset(
dataset_name="Temp", dataset_type="Temp", group_name=self.temp_group)
# Put values in Redis so they can be looked up later if
# added to a collection
Redis.set(self.trait_id, kw['trait_paste'], ex=ONE_YEAR)
self.trait_vals = kw['trait_paste'].split()
self.this_trait = create_trait(dataset=self.dataset,
name=self.trait_id,
cellid=None)
else:
self.temp_trait = True
self.trait_id = kw['trait_id']
self.temp_species = self.trait_id.split("_")[1]
self.temp_group = self.trait_id.split("_")[2]
self.dataset = data_set.create_dataset(
dataset_name="Temp", dataset_type="Temp", group_name=self.temp_group)
self.this_trait = create_trait(dataset=self.dataset,
name=self.trait_id,
cellid=None)
self.trait_vals = Redis.get(self.trait_id).split()
# Get verify/rna-seq link URLs
try:
blatsequence = self.this_trait.blatseq if self.dataset.type == "ProbeSet" else self.this_trait.sequence
if not blatsequence:
# XZ, 06/03/2009: ProbeSet name is not unique among platforms. We should use ProbeSet Id instead.
seqs = ()
if self.dataset.type == "ProbeSet":
db_cursor.execute(
"SELECT Probe.Sequence, Probe.Name "
"FROM Probe, ProbeSet, ProbeSetFreeze, "
"ProbeSetXRef WHERE "
"ProbeSetXRef.ProbeSetFreezeId = ProbeSetFreeze.Id "
"AND ProbeSetXRef.ProbeSetId = ProbeSet.Id AND "
"ProbeSetFreeze.Name = %s AND "
"ProbeSet.Name = %s AND "
"Probe.ProbeSetId = ProbeSet.Id ORDER "
"BY Probe.SerialOrder",
(self.this_trait.dataset.name, self.this_trait.name,)
)
else:
db_cursor.execute(
"SELECT Geno.Sequence "
"FROM Geno, GenoXRef, GenoFreeze "
"WHERE Geno.Name = %s AND "
"Geno.Id = GenoXRef.GenoId AND "
"GenoXRef.GenoFreezeId = GenoFreeze.Id AND "
"GenoFreeze.Name = %s",
(self.this_trait.name, self.this_trait.dataset.name)
)
seqs = db_cursor.fetchall()
if not seqs:
raise ValueError
else:
blatsequence = ''
for seqt in seqs:
if int(seqt[1][-1]) % 2 == 1:
blatsequence += string.strip(seqt[0])
# --------Hongqiang add this part in order to not only blat ProbeSet, but also blat Probe
blatsequence = '%3E' + self.this_trait.name + '%0A' + blatsequence + '%0A'
# XZ, 06/03/2009: ProbeSet name is not unique among platforms. We should use ProbeSet Id instead.
seqs = ()
db_cursor.execute(
"SELECT Probe.Sequence, Probe.Name "
"FROM Probe, ProbeSet, ProbeSetFreeze, "
"ProbeSetXRef WHERE "
"ProbeSetXRef.ProbeSetFreezeId = ProbeSetFreeze.Id "
"AND ProbeSetXRef.ProbeSetId = ProbeSet.Id AND "
"ProbeSetFreeze.Name = %s AND ProbeSet.Name = %s "
"AND Probe.ProbeSetId = ProbeSet.Id "
"ORDER BY Probe.SerialOrder",
(self.this_trait.dataset.name, self.this_trait.name,)
)
seqs = db_cursor.fetchall()
for seqt in seqs:
if int(seqt[1][-1]) % 2 == 1:
blatsequence += '%3EProbe_' + \
seqt[1].strip() + '%0A' + seqt[0].strip() + '%0A'
if self.dataset.group.species == "rat":
self.UCSC_BLAT_URL = webqtlConfig.UCSC_BLAT % (
'rat', 'rn7', blatsequence)
self.UTHSC_BLAT_URL = ""
elif self.dataset.group.species == "mouse":
self.UCSC_BLAT_URL = webqtlConfig.UCSC_BLAT % (
'mouse', 'mm10', blatsequence)
self.UTHSC_BLAT_URL = webqtlConfig.UTHSC_BLAT % (
'mouse', 'mm10', blatsequence)
elif self.dataset.group.species == "human":
self.UCSC_BLAT_URL = webqtlConfig.UCSC_BLAT % (
'human', 'hg38', blatsequence)
self.UTHSC_BLAT_URL = ""
else:
self.UCSC_BLAT_URL = ""
self.UTHSC_BLAT_URL = ""
except:
self.UCSC_BLAT_URL = ""
self.UTHSC_BLAT_URL = ""
if self.dataset.type == "ProbeSet":
self.show_probes = "True"
trait_units = get_trait_units(self.this_trait)
self.get_external_links()
self.build_correlation_tools()
self.ncbi_summary = get_ncbi_summary(self.this_trait)
# Get nearest marker for composite mapping
if not self.temp_trait:
if check_if_attr_exists(self.this_trait, 'locus_chr') and self.dataset.type != "Geno" and self.dataset.type != "Publish":
self.nearest_marker = get_nearest_marker(
self.this_trait, self.dataset)
else:
self.nearest_marker = ""
self.make_sample_lists()
trait_vals_by_group = []
for sample_type in self.sample_groups:
trait_vals_by_group.append(get_trait_vals(sample_type.sample_list))
self.max_digits_by_group, self.no_decimal_place = get_max_digits(trait_vals_by_group)
self.qnorm_vals = quantile_normalize_vals(self.sample_groups, trait_vals_by_group)
self.z_scores = get_z_scores(self.sample_groups, trait_vals_by_group)
self.temp_uuid = uuid.uuid4()
self.sample_group_types = OrderedDict()
if len(self.sample_groups) > 1:
self.sample_group_types['samples_primary'] = self.dataset.group.name
self.sample_group_types['samples_other'] = "Other"
self.sample_group_types['samples_all'] = "All"
else:
self.sample_group_types['samples_primary'] = self.dataset.group.name
sample_lists = [group.sample_list for group in self.sample_groups]
self.categorical_var_list = []
self.numerical_var_list = []
if not self.temp_trait:
# ZS: Only using first samplelist, since I think mapping only uses those samples
self.categorical_var_list = get_categorical_variables(
self.this_trait, self.sample_groups[0])
self.numerical_var_list = get_numerical_variables(
self.this_trait, self.sample_groups[0])
# ZS: Get list of chromosomes to select for mapping
self.chr_list = [["All", -1]]
for i, this_chr in enumerate(self.dataset.species.chromosomes.chromosomes(db_cursor)):
self.chr_list.append(
[self.dataset.species.chromosomes.chromosomes(db_cursor)[this_chr].name, i])
self.genofiles = self.dataset.group.get_genofiles()
study_samplelist_json = self.dataset.group.get_study_samplelists()
self.study_samplelists = [study["title"] for study in study_samplelist_json]
# ZS: No need to grab scales from .geno file unless it's using
# a mapping method that reads .geno files
if "QTLReaper" or "R/qtl" in dataset.group.mapping_names:
if self.genofiles:
self.scales_in_geno = get_genotype_scales(self.genofiles)
else:
self.scales_in_geno = get_genotype_scales(
self.dataset.group.name + ".geno")
else:
self.scales_in_geno = {}
self.has_num_cases = has_num_cases(self.this_trait)
# ZS: Needed to know whether to display bar chart + get max
# sample name length in order to set table column width
self.num_values = 0
# ZS: So it knows whether to display the Binary R/qtl mapping
# method, which doesn't work unless all values are 0 or 1
self.binary = "true"
# ZS: Since we don't want to show log2 transform option for
# situations where it doesn't make sense
self.negative_vals_exist = "false"
max_samplename_width = 1
for group in self.sample_groups:
for sample in group.sample_list:
if len(sample.name) > max_samplename_width:
max_samplename_width = len(sample.name)
if sample.display_value != "x":
self.num_values += 1
if sample.display_value != 0 or sample.display_value != 1:
self.binary = "false"
if sample.value < 0:
self.negative_vals_exist = "true"
# ZS: Check whether any attributes have few enough distinct
# values to show the "Block samples by group" option
self.categorical_attr_exists = "false"
for attribute in self.sample_groups[0].attributes:
if len(self.sample_groups[0].attributes[attribute].distinct_values) <= 10:
self.categorical_attr_exists = "true"
break
sample_column_width = max_samplename_width * 8
self.stats_table_width, self.trait_table_width = get_table_widths(
self.sample_groups, sample_column_width, self.has_num_cases)
if self.num_values >= 5000:
self.maf = 0.01
else:
self.maf = 0.05
trait_symbol = None
short_description = None
if not self.temp_trait:
if self.this_trait.symbol:
trait_symbol = self.this_trait.symbol
short_description = trait_symbol
elif hasattr(self.this_trait, 'post_publication_abbreviation'):
short_description = self.this_trait.post_publication_abbreviation
elif hasattr(self.this_trait, 'pre_publication_abbreviation'):
short_description = self.this_trait.pre_publication_abbreviation
# Todo: Add back in the ones we actually need from below, as we discover we need them
hddn = OrderedDict()
if self.dataset.group.allsamples:
hddn['allsamples'] = ','.join(self.dataset.group.allsamples)
hddn['primary_samples'] = ','.join(self.primary_sample_names)
hddn['trait_id'] = self.trait_id
hddn['trait_display_name'] = self.this_trait.display_name
hddn['dataset'] = self.dataset.name
hddn['temp_trait'] = False
if self.temp_trait:
hddn['temp_trait'] = True
hddn['group'] = self.temp_group
hddn['species'] = self.temp_species
else:
hddn['group'] = self.dataset.group.name
hddn['species'] = self.dataset.group.species
hddn['use_outliers'] = False
hddn['method'] = "gemma"
hddn['mapmethod_rqtl'] = "hk"
hddn['mapmodel_rqtl'] = "normal"
hddn['pair_scan'] = ""
hddn['selected_chr'] = -1
hddn['mapping_display_all'] = True
hddn['suggestive'] = 0
hddn['study_samplelists'] = json.dumps(study_samplelist_json)
hddn['num_perm'] = 0
hddn['categorical_vars'] = ""
if self.categorical_var_list:
hddn['categorical_vars'] = ",".join(self.categorical_var_list)
hddn['manhattan_plot'] = ""
hddn['control_marker'] = ""
if not self.temp_trait:
if hasattr(self.this_trait, 'locus_chr') and self.this_trait.locus_chr != "" and self.dataset.type != "Geno" and self.dataset.type != "Publish":
hddn['control_marker'] = self.nearest_marker
hddn['do_control'] = False
hddn['maf'] = 0.05
hddn['mapping_scale'] = "physic"
hddn['compare_traits'] = []
hddn['export_data'] = ""
hddn['export_format'] = "excel"
if len(self.scales_in_geno) < 2 and bool(self.scales_in_geno):
hddn['mapping_scale'] = self.scales_in_geno[list(
self.scales_in_geno.keys())[0]][0][0]
# We'll need access to this_trait and hddn in the Jinja2
# Template, so we put it inside self
self.hddn = hddn
js_data = dict(trait_id=self.trait_id,
trait_symbol=trait_symbol,
max_digits = self.max_digits_by_group,
no_decimal_place = self.no_decimal_place,
short_description=short_description,
unit_type=trait_units,
dataset_type=self.dataset.type,
species=self.dataset.group.species,
scales_in_geno=self.scales_in_geno,
data_scale=self.dataset.data_scale,
sample_group_types=self.sample_group_types,
sample_lists=sample_lists,
se_exists=self.sample_groups[0].se_exists,
has_num_cases=self.has_num_cases,
attributes=self.sample_groups[0].attributes,
categorical_attr_exists=self.categorical_attr_exists,
categorical_vars=",".join(self.categorical_var_list),
num_values=self.num_values,
qnorm_values=self.qnorm_vals,
zscore_values=self.z_scores,
sample_column_width=sample_column_width,
temp_uuid=self.temp_uuid)
self.js_data = js_data
def get_external_links(self):
# ZS: There's some weirdness here because some fields don't
# exist while others are empty strings
self.pubmed_link = webqtlConfig.PUBMEDLINK_URL % self.this_trait.pubmed_id if check_if_attr_exists(
self.this_trait, 'pubmed_id') else None
self.ncbi_gene_link = webqtlConfig.NCBI_LOCUSID % self.this_trait.geneid if check_if_attr_exists(
self.this_trait, 'geneid') else None
self.omim_link = webqtlConfig.OMIM_ID % self.this_trait.omim if check_if_attr_exists(
self.this_trait, 'omim') else None
self.homologene_link = webqtlConfig.HOMOLOGENE_ID % self.this_trait.homologeneid if check_if_attr_exists(
self.this_trait, 'homologeneid') else None
self.genbank_link = None
if check_if_attr_exists(self.this_trait, 'genbankid'):
genbank_id = '|'.join(self.this_trait.genbankid.split('|')[0:10])
if genbank_id[-1] == '|':
genbank_id = genbank_id[0:-1]
self.genbank_link = webqtlConfig.GENBANK_ID % genbank_id
self.uniprot_link = None
if check_if_attr_exists(self.this_trait, 'uniprotid'):
self.uniprot_link = webqtlConfig.UNIPROT_URL % self.this_trait.uniprotid
self.genotation_link = self.rgd_link = self.phenogen_link = self.gtex_link = self.genebridge_link = self.ucsc_blat_link = self.biogps_link = self.protein_atlas_link = None
self.string_link = self.panther_link = self.aba_link = self.ebi_gwas_link = self.wiki_pi_link = self.genemania_link = self.ensembl_link = None
if self.this_trait.symbol:
self.genotation_link = webqtlConfig.GENOTATION_URL % self.this_trait.symbol
self.gtex_link = webqtlConfig.GTEX_URL % self.this_trait.symbol
self.string_link = webqtlConfig.STRING_URL % self.this_trait.symbol
self.panther_link = webqtlConfig.PANTHER_URL % self.this_trait.symbol
self.ebi_gwas_link = webqtlConfig.EBIGWAS_URL % self.this_trait.symbol
self.protein_atlas_link = webqtlConfig.PROTEIN_ATLAS_URL % self.this_trait.symbol
if self.dataset.group.species == "mouse" or self.dataset.group.species == "human":
self.rgd_link = webqtlConfig.RGD_URL % (
self.this_trait.symbol, self.dataset.group.species.capitalize())
if self.dataset.group.species == "mouse":
self.genemania_link = webqtlConfig.GENEMANIA_URL % (
"mus-musculus", self.this_trait.symbol)
else:
self.genemania_link = webqtlConfig.GENEMANIA_URL % (
"homo-sapiens", self.this_trait.symbol)
if self.dataset.group.species == "mouse":
self.aba_link = webqtlConfig.ABA_URL % self.this_trait.symbol
results = ()
with database_connection(get_setting("SQL_URI")) as conn, conn.cursor() as cursor:
cursor.execute(
"SELECT chromosome, txStart, txEnd FROM "
"GeneList WHERE geneSymbol = %s",
(self.this_trait.symbol,)
)
results = cursor.fetchone()
if results:
chr, transcript_start, transcript_end = results
else:
chr = transcript_start = transcript_end = None
if chr and transcript_start and transcript_end and self.this_trait.refseq_transcriptid:
transcript_start = int(transcript_start * 1000000)
transcript_end = int(transcript_end * 1000000)
self.ucsc_blat_link = webqtlConfig.UCSC_REFSEQ % (
'mm10', self.this_trait.refseq_transcriptid, chr, transcript_start, transcript_end)
if self.dataset.group.species == "rat":
self.rgd_link = webqtlConfig.RGD_URL % (
self.this_trait.symbol, self.dataset.group.species.capitalize())
self.phenogen_link = webqtlConfig.PHENOGEN_URL % (
self.this_trait.symbol)
self.genemania_link = webqtlConfig.GENEMANIA_URL % (
"rattus-norvegicus", self.this_trait.symbol)
results = ()
with database_connection(get_setting("SQL_URI")) as conn, conn.cursor() as cursor:
cursor.execute(
"SELECT kgID, chromosome, txStart, txEnd "
"FROM GeneList_rn33 WHERE geneSymbol = %s",
(self.this_trait.symbol,)
)
if results := cursor.fetchone():
kgId, chr, transcript_start, transcript_end = results
else:
kgId = chr = transcript_start = transcript_end = None
if chr and transcript_start and transcript_end and kgId:
# Convert to bases from megabases
transcript_start = int(transcript_start * 1000000)
transcript_end = int(transcript_end * 1000000)
self.ucsc_blat_link = webqtlConfig.UCSC_REFSEQ % (
'rn7', kgId, chr, transcript_start, transcript_end)
if self.this_trait.geneid and (self.dataset.group.species == "mouse" or self.dataset.group.species == "rat" or self.dataset.group.species == "human"):
self.biogps_link = webqtlConfig.BIOGPS_URL % (
self.dataset.group.species, self.this_trait.geneid)
self.gemma_link = webqtlConfig.GEMMA_URL % self.this_trait.geneid
if self.dataset.group.species == "human":
self.aba_link = webqtlConfig.ABA_URL % self.this_trait.geneid
def build_correlation_tools(self):
if self.temp_trait == True:
this_group = self.temp_group
else:
this_group = self.dataset.group.name
# We're checking a string here!
assert isinstance(this_group, str), "We need a string type thing here"
if this_group[:3] == 'BXD' and this_group not in webqtlConfig.BXD_GROUP_EXCEPTIONS:
this_group = 'BXD'
if this_group:
if self.temp_trait == True:
dataset_menu = data_set.datasets(this_group)
else:
dataset_menu = data_set.datasets(
this_group, self.dataset.group)
dataset_menu_selected = None
if len(dataset_menu):
if self.dataset:
dataset_menu_selected = self.dataset.name
return_results_menu = (100, 200, 500, 1000,
2000, 5000, 10000, 15000, 20000)
return_results_menu_selected = 500
self.corr_tools = dict(dataset_menu=dataset_menu,
dataset_menu_selected=dataset_menu_selected,
return_results_menu=return_results_menu,
return_results_menu_selected=return_results_menu_selected,)
def make_sample_lists(self):
all_samples_ordered = self.dataset.group.all_samples_ordered()
parent_f1_samples = []
if self.dataset.group.parlist and self.dataset.group.f1list:
parent_f1_samples = self.dataset.group.parlist + self.dataset.group.f1list
primary_sample_names = list(all_samples_ordered)
if not self.temp_trait:
other_sample_names = []
for sample in list(self.this_trait.data.keys()):
if (self.this_trait.data[sample].name2 != self.this_trait.data[sample].name):
if ((self.this_trait.data[sample].name2 in primary_sample_names)
and (self.this_trait.data[sample].name not in primary_sample_names)):
primary_sample_names.append(
self.this_trait.data[sample].name)
primary_sample_names.remove(
self.this_trait.data[sample].name2)
all_samples_set = set(all_samples_ordered)
if sample not in all_samples_set:
all_samples_ordered.append(sample)
other_sample_names.append(sample)
# CFW is here because the .geno file doesn't properly
# contain its full list of samples. This should probably
# be fixed.
if self.dataset.group.species == "human" or (set(primary_sample_names) == set(parent_f1_samples)) or self.dataset.group.name == "CFW":
primary_sample_names += other_sample_names
other_sample_names = []
if other_sample_names:
primary_header = "%s Only" % (self.dataset.group.name)
else:
primary_header = "Samples"
primary_samples = SampleList(dataset=self.dataset,
sample_names=primary_sample_names,
this_trait=self.this_trait,
sample_group_type='primary',
header=primary_header)
# if other_sample_names and self.dataset.group.species !=
# "human" and self.dataset.group.name != "CFW":
if len(other_sample_names) > 0:
other_sample_names.sort() # Sort other samples
other_samples = SampleList(dataset=self.dataset,
sample_names=other_sample_names,
this_trait=self.this_trait,
sample_group_type='other',
header="Other")
self.sample_groups = (primary_samples, other_samples)
else:
self.sample_groups = (primary_samples,)
else:
primary_samples = SampleList(dataset=self.dataset,
sample_names=primary_sample_names,
this_trait=self.trait_vals,
sample_group_type='primary',
header="%s Only" % (self.dataset.group.name))
self.sample_groups = (primary_samples,)
self.primary_sample_names = primary_sample_names
self.dataset.group.allsamples = all_samples_ordered
def get_trait_vals(sample_list):
trait_vals = []
for sample in sample_list:
try:
trait_vals.append(float(sample.value))
except:
continue
return trait_vals
def get_max_digits(trait_vals):
def __max_digits__(these_vals):
if not bool(these_vals):
return None
return len(str(max(these_vals))) - 1
return [__max_digits__(val) for val in trait_vals], all([val.is_integer() for val in sum(trait_vals, [])])
def normf(trait_vals):
ranked_vals = ss.rankdata(trait_vals)
p_list = []
for i, val in enumerate(trait_vals):
p_list.append(((i + 1) - 0.5) / len(trait_vals))
z = ss.norm.ppf(p_list)
normed_vals = []
for rank in ranked_vals:
normed_vals.append("%0.3f" % z[int(rank) - 1])
return normed_vals
def quantile_normalize_vals(sample_groups, trait_vals):
qnorm_by_group = []
for i, sample_type in enumerate(sample_groups):
qnorm_vals = normf(trait_vals[i])
qnorm_vals_with_x = []
counter = 0
for sample in sample_type.sample_list:
if sample.display_value == "x":
qnorm_vals_with_x.append("x")
else:
qnorm_vals_with_x.append(qnorm_vals[counter])
counter += 1
qnorm_by_group.append(qnorm_vals_with_x)
return qnorm_by_group
def get_z_scores(sample_groups, trait_vals):
zscore_by_group = []
for i, sample_type in enumerate(sample_groups):
zscores = ss.mstats.zscore(np.array(trait_vals[i])).tolist()
zscores_with_x = []
counter = 0
for sample in sample_type.sample_list:
if sample.display_value == "x":
zscores_with_x.append("x")
else:
zscores_with_x.append("%0.3f" % zscores[counter])
counter += 1
zscore_by_group.append(zscores_with_x)
return zscore_by_group
def get_nearest_marker(this_trait, this_db):
this_chr = this_trait.locus_chr
this_mb = this_trait.locus_mb
# One option is to take flanking markers, another is to take the
# two (or one) closest
with database_connection(get_setting("SQL_URI")) as conn, conn.cursor() as cursor:
cursor.execute(
"SELECT Geno.Name FROM Geno, GenoXRef, "
"GenoFreeze WHERE Geno.Chr = %s AND "
"GenoXRef.GenoId = Geno.Id AND "
"GenoFreeze.Id = GenoXRef.GenoFreezeId "
"AND GenoFreeze.Name = %s ORDER BY "
"ABS( Geno.Mb - %s) LIMIT 1",
(this_chr, f"{this_db.group.name}Geno", this_mb,)
)
if result := cursor.fetchall():
return result[0][0]
return ""
def get_table_widths(sample_groups, sample_column_width, has_num_cases=False):
stats_table_width = 250
if len(sample_groups) > 1:
stats_table_width = 450
trait_table_width = 300 + sample_column_width
if sample_groups[0].se_exists:
trait_table_width += 80
if has_num_cases:
trait_table_width += 80
trait_table_width += len(sample_groups[0].attributes) * 88
return stats_table_width, trait_table_width
def has_num_cases(this_trait):
has_n = False
if this_trait.dataset.type != "ProbeSet" and this_trait.dataset.type != "Geno":
for name, sample in list(this_trait.data.items()):
if sample.num_cases and sample.num_cases != "1":
has_n = True
break
return has_n
def get_trait_units(this_trait):
unit_type = ""
inside_brackets = False
if this_trait.description_fmt:
if ("[" in this_trait.description_fmt) and ("]" in this_trait.description_fmt):
for i in this_trait.description_fmt:
if inside_brackets:
if i != "]":
unit_type += i
else:
inside_brackets = False
if i == "[":
inside_brackets = True
if unit_type == "":
unit_type = "value"
return unit_type
def check_if_attr_exists(the_trait, id_type):
if hasattr(the_trait, id_type):
if getattr(the_trait, id_type) == None or getattr(the_trait, id_type) == "":
return False
else:
return True
else:
return False
def get_ncbi_summary(this_trait):
if check_if_attr_exists(this_trait, 'geneid'):
# ZS: Need to switch this try/except to something that checks
# the output later
try:
response = requests.get(
"https://eutils.ncbi.nlm.nih.gov/entrez/eutils/esummary.fcgi?db=gene&id=%s&retmode=json" % this_trait.geneid)
summary = json.loads(response.content)[
'result'][this_trait.geneid]['summary']
return summary
except:
return None
else:
return None
def get_categorical_variables(this_trait, sample_list) -> list:
categorical_var_list = []
if len(sample_list.attributes) > 0:
for attribute in sample_list.attributes:
if len(sample_list.attributes[attribute].distinct_values) < 10:
categorical_var_list.append(str(sample_list.attributes[attribute].id))
return categorical_var_list
def get_numerical_variables(this_trait, sample_list) -> list:
numerical_var_list = []
if len(sample_list.attributes) > 0:
for attribute in sample_list.attributes:
all_numeric = True
all_none = True
for attr_val in sample_list.attributes[attribute].distinct_values:
if not attr_val:
continue
try:
val_as_float = float(attr_val)
all_none = False
except:
all_numeric = False
break
if all_numeric and not all_none:
numerical_var_list.append(sample_list.attributes[attribute].name)
return numerical_var_list
def get_genotype_scales(genofiles):
geno_scales = {}
if isinstance(genofiles, list):
for the_file in genofiles:
file_location = the_file['location']
geno_scales[file_location] = get_scales_from_genofile(
file_location)
else:
geno_scales[genofiles] = get_scales_from_genofile(genofiles)
return geno_scales
def get_scales_from_genofile(file_location):
geno_path = locate_ignore_error(file_location, 'genotype')
# ZS: This is just to allow the code to run when
if not geno_path:
return [["physic", "Mb"]]
cm_and_mb_cols_exist = True
cm_column = None
mb_column = None
with open(geno_path, "r") as geno_fh:
for i, line in enumerate(geno_fh):
if line[0] == "#" or line[0] == "@":
# ZS: If the scale is made explicit in the metadata,
# use that
if "@scale" in line:
scale = line.split(":")[1].strip()
if scale == "morgan":
return [["morgan", "cM"]]
else:
return [["physic", "Mb"]]
else:
continue
if line[:3] == "Chr":
first_marker_line = i + 1
if line.split("\t")[2].strip() == "cM":
cm_column = 2
elif line.split("\t")[3].strip() == "cM":
cm_column = 3
if line.split("\t")[2].strip() == "Mb":
mb_column = 2
elif line.split("\t")[3].strip() == "Mb":
mb_column = 3
break
# ZS: This attempts to check whether the cM and Mb columns are
# 'real', since some .geno files have one column be a copy of
# the other column, or have one column that is all 0s
cm_all_zero = True
mb_all_zero = True
cm_mb_all_equal = True
for i, line in enumerate(geno_fh):
# ZS: I'm assuming there won't be more than 10 markers
# where the position is listed as 0
if first_marker_line <= i < first_marker_line + 10:
if cm_column:
cm_val = line.split("\t")[cm_column].strip()
if cm_val != "0":
cm_all_zero = False
if mb_column:
mb_val = line.split("\t")[mb_column].strip()
if mb_val != "0":
mb_all_zero = False
if cm_column and mb_column:
if cm_val != mb_val:
cm_mb_all_equal = False
else:
if i > first_marker_line + 10:
break
# ZS: This assumes that both won't be all zero, since if that's
# the case mapping shouldn't be an option to begin with
if mb_all_zero:
return [["morgan", "cM"]]
elif cm_mb_all_equal:
return [["physic", "Mb"]]
elif cm_and_mb_cols_exist:
return [["physic", "Mb"], ["morgan", "cM"]]
else:
return [["physic", "Mb"]]
def get_diff_of_vals(new_vals: Dict, trait_id: str) -> Dict:
""" Get the diff between current sample values and the values in the DB
Given a dict of the changed values and the trait/dataset ID, return a Dict
with keys corresponding to each sample with a changed value and a value
that is a dict with keys for the old_value and new_value
"""
trait_name = trait_id.split(":")[0]
dataset_name = trait_id.split(":")[1]
trait_ob = create_trait(name=trait_name, dataset_name=dataset_name)
old_vals = {sample : trait_ob.data[sample].value for sample in trait_ob.data}
shared_samples = set.union(set(new_vals.keys()), set(old_vals.keys()))
diff_dict = {}
for sample in shared_samples:
try:
new_val = round(float(new_vals[sample]), 3)
except:
new_val = "x"
try:
old_val = round(float(old_vals[sample]), 3)
except:
old_val = "x"
if new_val != old_val:
diff_dict[sample] = {
"new_val": new_val,
"old_val": old_val
}
return diff_dict
