diff options
6 files changed, 165 insertions, 499 deletions
diff --git a/wqflask/tests/unit/wqflask/api/test_mapping.py b/wqflask/tests/unit/wqflask/api/test_mapping.py index b094294a..159c982b 100644 --- a/wqflask/tests/unit/wqflask/api/test_mapping.py +++ b/wqflask/tests/unit/wqflask/api/test_mapping.py @@ -58,7 +58,7 @@ class TestMapping(unittest.TestCase): self.assertEqual(results_2, expected_results) - @mock.patch("wqflask.api.mapping.rqtl_mapping.run_rqtl_geno") + @mock.patch("wqflask.api.mapping.rqtl_mapping.run_rqtl") @mock.patch("wqflask.api.mapping.gemma_mapping.run_gemma") @mock.patch("wqflask.api.mapping.initialize_parameters") @mock.patch("wqflask.api.mapping.retrieve_sample_data") diff --git a/wqflask/tests/unit/wqflask/marker_regression/test_rqtl_mapping.py b/wqflask/tests/unit/wqflask/marker_regression/test_rqtl_mapping.py index 91d2c587..9d13e943 100644 --- a/wqflask/tests/unit/wqflask/marker_regression/test_rqtl_mapping.py +++ b/wqflask/tests/unit/wqflask/marker_regression/test_rqtl_mapping.py @@ -1,42 +1,43 @@ import unittest from unittest import mock -from wqflask import app -from wqflask.marker_regression.rqtl_mapping import get_trait_data_type -from wqflask.marker_regression.rqtl_mapping import sanitize_rqtl_phenotype -from wqflask.marker_regression.rqtl_mapping import sanitize_rqtl_names +from dataclasses import dataclass +from wqflask.marker_regression.rqtl_mapping import run_rqtl -class TestRqtlMapping(unittest.TestCase): +@dataclass +class MockGroup: + name: str + genofile: str + +@dataclass +class MockDataset: + group: MockGroup - def setUp(self): - self.app_context = app.app_context() - self.app_context.push() - - def tearDown(self): - self.app_context.pop() - - @mock.patch("wqflask.marker_regression.rqtl_mapping.g") - @mock.patch("wqflask.marker_regression.rqtl_mapping.logger") - def test_get_trait_data(self, mock_logger, mock_db): - """test for getting trait data_type return True""" - query_value = """SELECT value FROM TraitMetadata WHERE type='trait_data_type'""" - mock_db.db.execute.return_value.fetchone.return_value = [ - """{"type":"trait_data_type","name":"T1","traid_id":"fer434f"}"""] - results = get_trait_data_type("traid_id") - mock_db.db.execute.assert_called_with(query_value) - self.assertEqual(results, "fer434f") - - def test_sanitize_rqtl_phenotype(self): - """test for sanitizing rqtl phenotype""" - vals = ['f', "x", "r", "x", "x"] - results = sanitize_rqtl_phenotype(vals) - expected_phenotype_string = 'c(f,NA,r,NA,NA)' - - self.assertEqual(results, expected_phenotype_string) - - def test_sanitize_rqtl_names(self): - """test for sanitzing rqtl names""" - vals = ['f', "x", "r", "x", "x"] - expected_sanitized_name = "c('f',NA,'r',NA,NA)" - results = sanitize_rqtl_names(vals) - self.assertEqual(expected_sanitized_name, results) +class TestRqtlMapping(unittest.TestCase): + """Tests for functions in rqtl_mapping.py""" + @mock.patch("wqflask.marker_regression.rqtl_mapping.requests.post") + @mock.patch("wqflask.marker_regression.rqtl_mapping.locate") + @mock.patch("wqflask.marker_regression.rqtl_mapping.write_phenotype_file") + def test_run_rqtl_with_perm(self, mock_write_pheno_file, mock_locate, mock_post): + """Test for run_rqtl with permutations > 0""" + dataset_group = MockGroup("GP1", "file_geno") + dataset = MockDataset(dataset_group) + + mock_write_pheno_file.return_value = "pheno_filename" + mock_locate.return_value = "geno_filename" + mock_post.return_value = mock.Mock(ok=True) + mock_post.return_value.json.return_value = {"perm_results": [], + "suggestive": 3, + "significant": 4, + "results" : []} + + results = run_rqtl(trait_name="the_trait", vals=[], samples=[], + dataset=dataset, mapping_scale="cM", model="normal", method="hk", + num_perm=5, perm_strata_list=[], do_control="false", control_marker="", + manhattan_plot=True, cofactors="") + + mock_write_pheno_file.assert_called_once() + mock_locate.assert_called_once() + mock_post.assert_called_once() + + self.assertEqual(results, ([], 3, 4, [])) diff --git a/wqflask/tests/unit/wqflask/marker_regression/test_run_mapping.py b/wqflask/tests/unit/wqflask/marker_regression/test_run_mapping.py index 78cd3be9..c220a072 100644 --- a/wqflask/tests/unit/wqflask/marker_regression/test_run_mapping.py +++ b/wqflask/tests/unit/wqflask/marker_regression/test_run_mapping.py @@ -229,20 +229,20 @@ class TestRunMapping(unittest.TestCase): used_samples = ["S1", "S2"] sample_list = AttributeSetter({"sample_attribute_values": { "S1": { - "C1": "c1_value", - "C2": "c2_value", - "W1": "w1_value" + "c1": "c1_value", + "c2": "c2_value", + "w1": "w1_value" }, "S2": { - "W1": "w2_value", - "W2": "w2_value" + "w1": "w2_value", + "w2": "w2_value" }, "S3": { - "C1": "c1_value", - "C2": "c2_value" + "c1": "c1_value", + "c2": "c2_value" }, diff --git a/wqflask/wqflask/marker_regression/rqtl_mapping.py b/wqflask/wqflask/marker_regression/rqtl_mapping.py index 1c8477bf..09afb8d1 100644 --- a/wqflask/wqflask/marker_regression/rqtl_mapping.py +++ b/wqflask/wqflask/marker_regression/rqtl_mapping.py @@ -1,460 +1,133 @@ -import rpy2.robjects as ro -import rpy2.robjects.numpy2ri as np2r -import numpy as np -import json +import csv +import hashlib +import io +import requests +import shutil +from typing import Dict +from typing import List +from typing import Optional +from typing import TextIO -from flask import g +import numpy as np from base.webqtlConfig import TMPDIR from base.trait import create_trait -from base.data_set import create_dataset -from utility import webqtlUtil -from utility.tools import locate, TEMPDIR -from wqflask.marker_regression.gemma_mapping import generate_random_n_string -from flask import g +from utility.tools import locate import utility.logger logger = utility.logger.getLogger(__name__) -# Get a trait's type (numeric, categorical, etc) from the DB - - -def get_trait_data_type(trait_db_string): - logger.info("get_trait_data_type") - the_query = "SELECT value FROM TraitMetadata WHERE type='trait_data_type'" - logger.info("the_query done") - results_json = g.db.execute(the_query).fetchone() - logger.info("the_query executed") - results_ob = json.loads(results_json[0]) - logger.info("json results loaded") - if trait_db_string in results_ob: - logger.info("found") - return results_ob[trait_db_string] - else: - logger.info("not found") - return "numeric" - - -# Run qtl mapping using R/qtl -def run_rqtl_geno(vals, samples, dataset, mapping_scale, method, model, permCheck, num_perm, perm_strata_list, do_control, control_marker, manhattan_plot, pair_scan, cofactors): - logger.info("Start run_rqtl_geno") - # Get pointers to some common R functions - r_library = ro.r["library"] # Map the library function - r_c = ro.r["c"] # Map the c function - plot = ro.r["plot"] # Map the plot function - png = ro.r["png"] # Map the png function - dev_off = ro.r["dev.off"] # Map the device off function - - print((r_library("qtl"))) # Load R/qtl - - logger.info("QTL library loaded") +GN3_RQTL_URL = "http://localhost:8086/api/rqtl/compute" +GN3_TMP_PATH = "/export/local/home/zas1024/genenetwork3/tmp" - # Get pointers to some R/qtl functions - scanone = ro.r["scanone"] # Map the scanone function - scantwo = ro.r["scantwo"] # Map the scantwo function - # Map the calc.genoprob function - calc_genoprob = ro.r["calc.genoprob"] +def run_rqtl(trait_name, vals, samples, dataset, mapping_scale, model, method, num_perm, perm_strata_list, do_control, control_marker, manhattan_plot, cofactors): + """Run R/qtl by making a request to the GN3 endpoint and reading in the output file(s)""" - crossname = dataset.group.name - # try: - # generate_cross_from_rdata(dataset) - # read_cross_from_rdata = ro.r["generate_cross_from_rdata"] # Map the local read_cross_from_rdata function - # genofilelocation = locate(crossname + ".RData", "genotype/rdata") - # cross_object = read_cross_from_rdata(genofilelocation) # Map the local GENOtoCSVR function - # except: - - if mapping_scale == "morgan": - scale_units = "cM" - else: - scale_units = "Mb" - - generate_cross_from_geno(dataset, scale_units) - # Map the local GENOtoCSVR function - GENOtoCSVR = ro.r["GENOtoCSVR"] - crossfilelocation = TMPDIR + crossname + ".cross" + pheno_file = write_phenotype_file(trait_name, samples, vals, dataset, cofactors, perm_strata_list) if dataset.group.genofile: - genofilelocation = locate(dataset.group.genofile, "genotype") + geno_file = locate(dataset.group.genofile, "genotype") else: - genofilelocation = locate(dataset.group.name + ".geno", "genotype") - logger.info("Going to create a cross from geno") - # TODO: Add the SEX if that is available - cross_object = GENOtoCSVR(genofilelocation, crossfilelocation) - logger.info("before calc_genoprob") - if manhattan_plot: - cross_object = calc_genoprob(cross_object) + geno_file = locate(dataset.group.name + ".geno", "genotype") + + post_data = { + "pheno_file": pheno_file, + "geno_file": geno_file, + "model": model, + "method": method, + "nperm": num_perm, + "scale": mapping_scale + } + + if do_control == "true" and control_marker: + post_data["control_marker"] = control_marker + + if not manhattan_plot: + post_data["interval"] = True + if cofactors: + post_data["addcovar"] = True + + if perm_strata_list: + post_data["pstrata"] = True + + rqtl_output = requests.post(GN3_RQTL_URL, data=post_data).json() + if num_perm > 0: + return rqtl_output['perm_results'], rqtl_output['suggestive'], rqtl_output['significant'], rqtl_output['results'] else: - cross_object = calc_genoprob(cross_object, step=5, stepwidth="max") - logger.info("after calc_genoprob") + return rqtl_output['results'] - pheno_string = sanitize_rqtl_phenotype(vals) - logger.info("phenostring done") - names_string = sanitize_rqtl_names(samples) - logger.info("sanitized pheno and names") - # Add the phenotype - cross_object = add_phenotype(cross_object, pheno_string, "the_pheno") - # Add the phenotype - cross_object = add_names(cross_object, names_string, "the_names") - logger.info("Added pheno and names") - # Create the additive covariate markers - marker_covars = create_marker_covariates(control_marker, cross_object) - logger.info("Marker covars done") - if cofactors != "": - logger.info("Cofactors: " + cofactors) - # Create the covariates from selected traits - cross_object, trait_covars = add_cofactors( - cross_object, dataset, cofactors, samples) - ro.r('all_covars <- cbind(marker_covars, trait_covars)') - else: - ro.r('all_covars <- marker_covars') - covars = ro.r['all_covars'] - # DEBUG to save the session object to file - if pair_scan: - if do_control == "true": - logger.info("Using covariate") - result_data_frame = scantwo( - cross_object, pheno="the_pheno", addcovar=covars, model=model, method=method, n_cluster=16) - else: - logger.info("No covariates") - result_data_frame = scantwo( - cross_object, pheno="the_pheno", model=model, method=method, n_cluster=16) - pair_scan_filename = webqtlUtil.genRandStr("scantwo_") + ".png" - png(file=TEMPDIR + pair_scan_filename) - plot(result_data_frame) - dev_off() +def get_hash_of_textio(the_file: TextIO) -> str: + """Given a StringIO, return the hash of its contents""" - return process_pair_scan_results(result_data_frame) - else: - if do_control == "true" or cofactors != "": - logger.info("Using covariate") - ro.r(f"qtl_results = scanone(the_cross, pheno='the_pheno', addcovar=all_covars, model='{model}', method='{method}')") - result_data_frame = ro.r("qtl_results") - else: - ro.r(f"qtl_results = scanone(the_cross, pheno='the_pheno', model='{model}', method='{method}')") - result_data_frame = np.asarray(ro.r("qtl_results")).T + the_file.seek(0) + hash_of_file = hashlib.md5(the_file.read().encode()).hexdigest() - marker_names = np.asarray(ro.r("row.names(qtl_results)")) + return hash_of_file - # Do permutation (if requested by user) - if num_perm > 0 and permCheck == "ON": - # ZS: The strata list would only be populated if "Stratified" was checked on before mapping - if len(perm_strata_list) > 0: - cross_object, strata_ob = add_perm_strata( - cross_object, perm_strata_list) - if do_control == "true" or cofactors != "": - perm_data_frame = scanone(cross_object, pheno_col="the_pheno", addcovar=covars, n_perm=int( - num_perm), perm_strata=strata_ob, model=model, method=method) - else: - perm_data_frame = scanone( - cross_object, pheno_col="the_pheno", n_perm=num_perm, perm_strata=strata_ob, model=model, method=method) - else: - if do_control == "true" or cofactors != "": - perm_data_frame = scanone(cross_object, pheno_col="the_pheno", addcovar=covars, n_perm=int( - num_perm), model=model, method=method) - else: - perm_data_frame = scanone( - cross_object, pheno_col="the_pheno", n_perm=num_perm, model=model, method=method) - - # Functions that sets the thresholds for the webinterface - perm_output, suggestive, significant = process_rqtl_perm_results( - num_perm, perm_data_frame) - return perm_output, suggestive, significant, process_rqtl_results(marker_names, result_data_frame, dataset.group.species) - else: - return process_rqtl_results(marker_names, result_data_frame, dataset.group.species) +def write_phenotype_file(trait_name: str, + samples: List[str], + vals: List, + dataset_ob, + cofactors: Optional[str] = None, + perm_strata_list: Optional[List] = None) -> TextIO: + """Given trait name, sample list, value list, dataset ob, and optional string + representing cofactors, return the file's full path/name + """ + cofactor_data = cofactors_to_dict(cofactors, dataset_ob, samples) -def generate_cross_from_rdata(dataset): - rdata_location = locate(dataset.group.name + ".RData", "genotype/rdata") - ro.r(""" - generate_cross_from_rdata <- function(filename = '%s') { - load(file=filename) - cross = cunique - return(cross) - } - """ % (rdata_location)) + pheno_file = io.StringIO() + writer = csv.writer(pheno_file, delimiter="\t", quoting=csv.QUOTE_NONE) + header_row = ["Samples", trait_name] + header_row += [cofactor for cofactor in cofactor_data] + if perm_strata_list: + header_row.append("Strata") -# TODO: Need to figure out why some genofiles have the wrong format and don't convert properly -def generate_cross_from_geno(dataset, scale_units): - - cross_filename = (f"{str(dataset.group.name)}_" - f"{generate_random_n_string(6)}") - - ro.r(""" - trim <- function( x ) { gsub("(^[[:space:]]+|[[:space:]]+$)", "", x) } - getGenoCode <- function(header, name = 'unk'){ - mat = which(unlist(lapply(header,function(x){ length(grep(paste('@',name,sep=''), x)) })) == 1) - return(trim(strsplit(header[mat],':')[[1]][2])) - } - GENOtoCSVR <- function(genotypes = '%s', out = '%s.csvr', phenotype = NULL, sex = NULL, verbose = FALSE){ - header = readLines(genotypes, 40) # Assume a geno header is not longer than 40 lines - toskip = which(unlist(lapply(header, function(x){ length(grep("Chr\t", x)) })) == 1)-1 # Major hack to skip the geno headers - type <- getGenoCode(header, 'type') - if(type == '4-way'){ - genocodes <- NULL - } else { - genocodes <- c(getGenoCode(header, 'mat'), getGenoCode(header, 'het'), getGenoCode(header, 'pat')) # Get the genotype codes - } - genodata <- read.csv(genotypes, sep='\t', skip=toskip, header=TRUE, na.strings=getGenoCode(header,'unk'), colClasses='character', comment.char = '#') - cat('Genodata:', toskip, " ", dim(genodata), genocodes, '\n') - if(is.null(phenotype)) phenotype <- runif((ncol(genodata)-4)) # If there isn't a phenotype, generate a random one - if(is.null(sex)) sex <- rep('m', (ncol(genodata)-4)) # If there isn't a sex phenotype, treat all as males - outCSVR <- rbind(c('Pheno', '', '', phenotype), # Phenotype - c('sex', '', '', sex), # Sex phenotype for the mice - cbind(genodata[,c('Locus','Chr', '%s')], genodata[, 5:ncol(genodata)])) # Genotypes - write.table(outCSVR, file = out, row.names=FALSE, col.names=FALSE,quote=FALSE, sep=',') # Save it to a file - require(qtl) - if(type == '4-way'){ - cat('Loading in as 4-WAY\n') - cross = read.cross(file=out, 'csvr', genotypes=NULL, crosstype="4way") # Load the created cross file using R/qtl read.cross - }else if(type == 'f2'){ - cat('Loading in as F2\n') - cross = read.cross(file=out, 'csvr', genotypes=genocodes, crosstype="f2") # Load the created cross file using R/qtl read.cross - }else{ - cat('Loading in as normal\n') - cross = read.cross(file=out, 'csvr', genotypes=genocodes) # Load the created cross file using R/qtl read.cross - } - if(type == 'riset'){ - cat('Converting to RISELF\n') - cross <- convert2riself(cross) # If its a RIL, convert to a RIL in R/qtl - } - return(cross) - } - """ % (dataset.group.genofile, cross_filename, scale_units)) - - -def add_perm_strata(cross, perm_strata): - col_string = 'c("the_strata")' - perm_strata_string = "c(" - for item in perm_strata: - perm_strata_string += str(item) + "," - - perm_strata_string = perm_strata_string[:-1] + ")" - - cross = add_phenotype(cross, perm_strata_string, "the_strata") - - strata_ob = pull_var("perm_strata", cross, col_string) - - return cross, strata_ob - - -def sanitize_rqtl_phenotype(vals): - pheno_as_string = "c(" - for i, val in enumerate(vals): - if val == "x": - if i < (len(vals) - 1): - pheno_as_string += "NA," - else: - pheno_as_string += "NA" - else: - if i < (len(vals) - 1): - pheno_as_string += str(val) + "," - else: - pheno_as_string += str(val) - pheno_as_string += ")" - - return pheno_as_string - - -def sanitize_rqtl_names(vals): - pheno_as_string = "c(" - for i, val in enumerate(vals): - if val == "x": - if i < (len(vals) - 1): - pheno_as_string += "NA," - else: - pheno_as_string += "NA" + writer.writerow(header_row) + for i, sample in enumerate(samples): + this_row = [sample] + if vals[i] != "x": + this_row.append(vals[i]) else: - if i < (len(vals) - 1): - pheno_as_string += "'" + str(val) + "'," - else: - pheno_as_string += "'" + str(val) + "'" - pheno_as_string += ")" - - return pheno_as_string - - -def add_phenotype(cross, pheno_as_string, col_name): - ro.globalenv["the_cross"] = cross - ro.r('pheno <- data.frame(pull.pheno(the_cross))') - ro.r('the_cross$pheno <- cbind(pheno, ' + col_name + \ - ' = as.numeric(' + pheno_as_string + '))') - return ro.r["the_cross"] - - -def add_categorical_covar(cross, covar_as_string, i): - ro.globalenv["the_cross"] = cross - logger.info("cross set") - ro.r('covar <- as.factor(' + covar_as_string + ')') - logger.info("covar set") - ro.r('newcovar <- model.matrix(~covar-1)') - logger.info("model.matrix finished") - ro.r('cat("new covar columns", ncol(newcovar), "\n")') - nCol = ro.r('ncol(newcovar)') - logger.info("ncol covar done: " + str(nCol[0])) - ro.r('pheno <- data.frame(pull.pheno(the_cross))') - logger.info("pheno pulled from cross") - nCol = int(nCol[0]) - logger.info("nCol python int:" + str(nCol)) - col_names = [] - # logger.info("loop") - for x in range(1, (nCol + 1)): - #logger.info("loop" + str(x)); - col_name = "covar_" + str(i) + "_" + str(x) - #logger.info("col_name" + col_name); - ro.r('the_cross$pheno <- cbind(pheno, ' + \ - col_name + ' = newcovar[,' + str(x) + '])') - col_names.append(col_name) - #logger.info("loop" + str(x) + "done"); - - logger.info("returning from add_categorical_covar") - return ro.r["the_cross"], col_names - - -def add_names(cross, names_as_string, col_name): - ro.globalenv["the_cross"] = cross - ro.r('pheno <- data.frame(pull.pheno(the_cross))') - ro.r('the_cross$pheno <- cbind(pheno, ' + \ - col_name + ' = ' + names_as_string + ')') - return ro.r["the_cross"] - - -def pull_var(var_name, cross, var_string): - ro.globalenv["the_cross"] = cross - ro.r(var_name + ' <- pull.pheno(the_cross, ' + var_string + ')') - - return ro.r[var_name] - - -def add_cofactors(cross, this_dataset, covariates, samples): - ro.numpy2ri.activate() - - covariate_list = covariates.split(",") - covar_name_string = "c(" - for i, covariate in enumerate(covariate_list): - logger.info("Covariate: " + covariate) - this_covar_data = [] - covar_as_string = "c(" - trait_name = covariate.split(":")[0] - dataset_ob = create_dataset(covariate.split(":")[1]) - trait_ob = create_trait(dataset=dataset_ob, - name=trait_name, - cellid=None) - - this_dataset.group.get_samplelist() - trait_samples = this_dataset.group.samplelist - trait_sample_data = trait_ob.data - for index, sample in enumerate(samples): - if sample in trait_samples: - if sample in trait_sample_data: - sample_value = trait_sample_data[sample].value - this_covar_data.append(sample_value) - else: - this_covar_data.append("NA") - - for j, item in enumerate(this_covar_data): - if j < (len(this_covar_data) - 1): - covar_as_string += str(item) + "," - else: - covar_as_string += str(item) - - covar_as_string += ")" - - datatype = get_trait_data_type(covariate) - logger.info("Covariate: " + covariate + " is of type: " + datatype) - if(datatype == "categorical"): # Cat variable - logger.info("call of add_categorical_covar") - cross, col_names = add_categorical_covar( - cross, covar_as_string, i) # Expand and add it to the cross - logger.info("add_categorical_covar returned") - # Go through the additional covar names - for z, col_name in enumerate(col_names): - if i < (len(covariate_list) - 1): - covar_name_string += '"' + col_name + '", ' - else: - if(z < (len(col_names) - 1)): - covar_name_string += '"' + col_name + '", ' + this_row.append("NA") + for cofactor in cofactor_data: + this_row.append(cofactor_data[cofactor][i]) + if perm_strata_list: + this_row.append(perm_strata_list[i]) + writer.writerow(this_row) + + hash_of_file = get_hash_of_textio(pheno_file) + file_path = TMPDIR + hash_of_file + ".csv" + + with open(file_path, "w") as fd: + pheno_file.seek(0) + shutil.copyfileobj(pheno_file, fd) + + return file_path + + +def cofactors_to_dict(cofactors: str, dataset_ob, samples) -> Dict: + """Given a string of cofactors, the trait being mapped's dataset ob, + and list of samples, return cofactor data as a Dict + + """ + cofactor_dict = {} + if cofactors: + dataset_ob.group.get_samplelist() + sample_list = dataset_ob.group.samplelist + for cofactor in cofactors.split(","): + cofactor_name, cofactor_dataset = cofactor.split(":") + if cofactor_dataset == dataset_ob.name: + cofactor_dict[cofactor_name] = [] + trait_ob = create_trait(dataset=dataset_ob, + name=cofactor_name) + sample_data = trait_ob.data + for index, sample in enumerate(samples): + if sample in sample_data: + sample_value = sample_data[sample].value + cofactor_dict[cofactor_name].append(sample_value) else: - covar_name_string += '"' + col_name + '"' - else: - col_name = "covar_" + str(i) - cross = add_phenotype(cross, covar_as_string, col_name) - if i < (len(covariate_list) - 1): - covar_name_string += '"' + col_name + '", ' - else: - covar_name_string += '"' + col_name + '"' - - covar_name_string += ")" - covars_ob = pull_var("trait_covars", cross, covar_name_string) - return cross, covars_ob - - -def create_marker_covariates(control_marker, cross): - ro.globalenv["the_cross"] = cross - # Get the genotype matrix - ro.r('genotypes <- pull.geno(the_cross)') - # TODO: sanitize user input, Never Ever trust a user - userinput_sanitized = control_marker.replace(" ", "").split(",") - logger.debug(userinput_sanitized) - if len(userinput_sanitized) > 0: - covariate_names = ', '.join('"{0}"'.format(w) - for w in userinput_sanitized) - ro.r('covnames <- c(' + covariate_names + ')') - else: - ro.r('covnames <- c()') - ro.r('covInGeno <- which(covnames %in% colnames(genotypes))') - ro.r('covnames <- covnames[covInGeno]') - ro.r("cat('covnames (purged): ', covnames,'\n')") - # Get the covariate matrix by using the marker name as index to the genotype file - ro.r('marker_covars <- genotypes[,covnames]') - # TODO: Create a design matrix from the marker covars for the markers in case of an F2, 4way, etc - return ro.r["marker_covars"] - - -def process_pair_scan_results(result): - pair_scan_results = [] - - result = result[1] - output = [tuple([result[j][i] for j in range(result.ncol)]) - for i in range(result.nrow)] - - for i, line in enumerate(result.iter_row()): - marker = {} - marker['name'] = result.rownames[i] - marker['chr1'] = output[i][0] - marker['Mb'] = output[i][1] - marker['chr2'] = int(output[i][2]) - pair_scan_results.append(marker) - - return pair_scan_results - - -def process_rqtl_perm_results(num_perm, results): - perm_vals = [item[0] for item in results] - - perm_output = perm_vals - suggestive = np.percentile(np.array(perm_vals), 67) - significant = np.percentile(np.array(perm_vals), 95) - - return perm_output, suggestive, significant - - -def process_rqtl_results(marker_names, results, species_name): # TODO: how to make this a one liner and not copy the stuff in a loop - qtl_results = [] - - for i, line in enumerate(results): - marker = {} - marker['name'] = marker_names[i] - if species_name == "mouse" and line[0] == 20: - marker['chr'] = "X" - else: - try: - marker['chr'] = int(line[0]) - except: - marker['chr'] = line[0] - marker['cM'] = marker['Mb'] = line[1] - marker['lod_score'] = line[2] - qtl_results.append(marker) - - return qtl_results + cofactor_dict[cofactor_name].append("NA") + return cofactor_dict diff --git a/wqflask/wqflask/marker_regression/run_mapping.py b/wqflask/wqflask/marker_regression/run_mapping.py index a3b579ec..c5b980a7 100644 --- a/wqflask/wqflask/marker_regression/run_mapping.py +++ b/wqflask/wqflask/marker_regression/run_mapping.py @@ -242,11 +242,11 @@ class RunMapping: # if start_vars['pair_scan'] == "true": # self.pair_scan = True if self.permCheck and self.num_perm > 0: - self.perm_output, self.suggestive, self.significant, results = rqtl_mapping.run_rqtl_geno( - self.vals, self.samples, self.dataset, self.mapping_scale, self.method, self.model, self.permCheck, self.num_perm, perm_strata, self.do_control, self.control_marker, self.manhattan_plot, self.pair_scan, self.covariates) + self.perm_output, self.suggestive, self.significant, results = rqtl_mapping.run_rqtl( + self.this_trait.name, self.vals, self.samples, self.dataset, self.mapping_scale, self.model, self.method, self.num_perm, perm_strata, self.do_control, self.control_marker, self.manhattan_plot, self.covariates) else: - results = rqtl_mapping.run_rqtl_geno(self.vals, self.samples, self.dataset, self.mapping_scale, self.method, self.model, self.permCheck, - self.num_perm, perm_strata, self.do_control, self.control_marker, self.manhattan_plot, self.pair_scan, self.covariates) + results = rqtl_mapping.run_rqtl(self.this_trait.name, self.vals, self.samples, self.dataset, self.mapping_scale, self.model, self.method, + self.num_perm, perm_strata, self.do_control, self.control_marker, self.manhattan_plot, self.covariates) elif self.mapping_method == "reaper": if "startMb" in start_vars: # ZS: Check if first time page loaded, so it can default to ON if "additiveCheck" in start_vars: @@ -765,9 +765,9 @@ def get_perm_strata(this_trait, sample_list, categorical_vars, used_samples): if sample in list(sample_list.sample_attribute_values.keys()): combined_string = "" for var in categorical_vars: - if var in list(sample_list.sample_attribute_values[sample].keys()): + if var.lower() in sample_list.sample_attribute_values[sample]: combined_string += str( - sample_list.sample_attribute_values[sample][var]) + sample_list.sample_attribute_values[sample][var.lower()]) else: combined_string += "NA" else: diff --git a/wqflask/wqflask/show_trait/show_trait.py b/wqflask/wqflask/show_trait/show_trait.py index 837c7a54..a02da872 100644 --- a/wqflask/wqflask/show_trait/show_trait.py +++ b/wqflask/wqflask/show_trait/show_trait.py @@ -279,6 +279,8 @@ class ShowTrait: hddn['suggestive'] = 0 hddn['num_perm'] = 0 hddn['categorical_vars'] = "" + if categorical_var_list: + hddn['categorical_vars'] = ",".join(categorical_var_list) hddn['manhattan_plot'] = "" hddn['control_marker'] = "" if not self.temp_trait: @@ -684,23 +686,13 @@ def get_ncbi_summary(this_trait): return None -def get_categorical_variables(this_trait, sample_list): +def get_categorical_variables(this_trait, sample_list) -> list: categorical_var_list = [] if len(sample_list.attributes) > 0: for attribute in sample_list.attributes: - attribute_vals = [] - for sample_name in list(this_trait.data.keys()): - if sample_list.attributes[attribute].name in this_trait.data[sample_name].extra_attributes: - attribute_vals.append( - this_trait.data[sample_name].extra_attributes[sample_list.attributes[attribute].name]) - else: - attribute_vals.append("N/A") - num_distinct = len(set(attribute_vals)) - - if num_distinct < 10: - categorical_var_list.append( - sample_list.attributes[attribute].name) + if len(sample_list.attributes[attribute].distinct_values) < 10: + categorical_var_list.append(sample_list.attributes[attribute].name) return categorical_var_list |