From 56538279d3ecf2e039b995dccb610041897b2888 Mon Sep 17 00:00:00 2001 From: zsloan Date: Wed, 30 Mar 2016 18:24:07 +0000 Subject: GN1 mapping now when mapping scale is centimorgan (for R/qtl) Fixed some issues related to zooming in on chromosome X Fixed problem where it threw an error for mapping methods that didn't return additive effect Fixed problem where the endMb in the top form would be set to the endMb of the last chromosome when viewing full genome --- .../wqflask/marker_regression/marker_regression.py | 23 ++++------ .../marker_regression/marker_regression_gn1.py | 51 +++++++++++++--------- .../wqflask/templates/marker_regression_gn1.html | 6 ++- 3 files changed, 43 insertions(+), 37 deletions(-) diff --git a/wqflask/wqflask/marker_regression/marker_regression.py b/wqflask/wqflask/marker_regression/marker_regression.py index 85404f58..6ce35a6c 100644 --- a/wqflask/wqflask/marker_regression/marker_regression.py +++ b/wqflask/wqflask/marker_regression/marker_regression.py @@ -337,9 +337,6 @@ class MarkerRegression(object): count, p_values = self.parse_rqtl_output(plink_output_filename) def geno_to_rqtl_function(self): # TODO: Need to figure out why some genofiles have the wrong format and don't convert properly - print("Adding some custom helper functions to the R environment") - - ro.r(""" trim <- function( x ) { gsub("(^[[:space:]]+|[[:space:]]+$)", "", x) } @@ -371,8 +368,6 @@ class MarkerRegression(object): """ % (self.dataset.group.name + ".geno")) def run_rqtl_geno(self): - print("Calling R/qtl") - self.geno_to_rqtl_function() ## Get pointers to some common R functions @@ -397,7 +392,7 @@ class MarkerRegression(object): genofilelocation = webqtlConfig.HTMLPATH + "genotypes/" + self.dataset.group.name + ".geno" crossfilelocation = webqtlConfig.HTMLPATH + "genotypes/" + self.dataset.group.name + ".cross" - print("Conversion of geno to cross at location:", genofilelocation, " to ", crossfilelocation) + #print("Conversion of geno to cross at location:", genofilelocation, " to ", crossfilelocation) cross_object = GENOtoCSVR(genofilelocation, crossfilelocation) # TODO: Add the SEX if that is available @@ -419,7 +414,7 @@ class MarkerRegression(object): else: print("No covariates"); result_data_frame = scantwo(cross_object, pheno = "the_pheno", model=self.model, method=self.method, n_cluster = 16) - print("Pair scan results:", result_data_frame) + #print("Pair scan results:", result_data_frame) self.pair_scan_filename = webqtlUtil.genRandStr("scantwo_") + ".png" png(file=webqtlConfig.TMPDIR+self.pair_scan_filename) @@ -454,13 +449,13 @@ class MarkerRegression(object): ro.r('genotypes <- pull.geno(the_cross)') # Get the genotype matrix userinputS = self.control.replace(" ", "").split(",") # TODO: sanitize user input, Never Ever trust a user covariate_names = ', '.join('"{0}"'.format(w) for w in userinputS) - print("Marker names of selected covariates:", covariate_names) + #print("Marker names of selected covariates:", covariate_names) ro.r('covnames <- c(' + covariate_names + ')') ro.r('covInGeno <- which(covnames %in% colnames(genotypes))') ro.r('covnames <- covnames[covInGeno]') ro.r("cat('covnames (purged): ', covnames,'\n')") ro.r('covariates <- genotypes[,covnames]') # Get the covariate matrix by using the marker name as index to the genotype file - print("R/qtl matrix of covariates:", ro.r["covariates"]) + #print("R/qtl matrix of covariates:", ro.r["covariates"]) return ro.r["covariates"] def sanitize_rqtl_phenotype(self): @@ -484,7 +479,7 @@ class MarkerRegression(object): result = result[1] output = [tuple([result[j][i] for j in range(result.ncol)]) for i in range(result.nrow)] - print("R/qtl scantwo output:", output) + #print("R/qtl scantwo output:", output) for i, line in enumerate(result.iter_row()): marker = {} @@ -494,7 +489,7 @@ class MarkerRegression(object): marker['chr2'] = int(output[i][2]) pair_scan_results.append(marker) - print("pair_scan_results:", pair_scan_results) + #print("pair_scan_results:", pair_scan_results) return pair_scan_results @@ -502,7 +497,7 @@ class MarkerRegression(object): qtl_results = [] output = [tuple([result[j][i] for j in range(result.ncol)]) for i in range(result.nrow)] - print("R/qtl scanone output:", output) + #print("R/qtl scanone output:", output) for i, line in enumerate(result.iter_row()): marker = {} @@ -517,7 +512,7 @@ class MarkerRegression(object): def process_rqtl_perm_results(self, results): perm_vals = [] for line in str(results).split("\n")[1:(int(self.num_perm)+1)]: - print("R/qtl permutation line:", line.split()) + #print("R/qtl permutation line:", line.split()) perm_vals.append(float(line.split()[1])) self.suggestive = np.percentile(np.array(perm_vals), 67) @@ -532,7 +527,7 @@ class MarkerRegression(object): self.gen_pheno_txt_file_plink(pheno_filename = plink_output_filename) plink_command = PLINK_COMMAND + ' --noweb --bed %s/%s.bed --bim %s/%s.bim --fam %s/%s.fam --no-fid --no-parents --no-sex --no-pheno --pheno %s%s.txt --pheno-name %s --maf %s --missing-phenotype -9999 --out %s%s --assoc ' % (PLINK_PATH, self.dataset.group.name, PLINK_PATH, self.dataset.group.name, PLINK_PATH, self.dataset.group.name, webqtlConfig.TMPDIR, plink_output_filename, self.this_trait.name, self.maf, webqtlConfig.TMPDIR, plink_output_filename) - print("plink_command:", plink_command) + #print("plink_command:", plink_command) os.system(plink_command) diff --git a/wqflask/wqflask/marker_regression/marker_regression_gn1.py b/wqflask/wqflask/marker_regression/marker_regression_gn1.py index e3466bef..19d521a4 100644 --- a/wqflask/wqflask/marker_regression/marker_regression_gn1.py +++ b/wqflask/wqflask/marker_regression/marker_regression_gn1.py @@ -485,7 +485,7 @@ class MarkerRegression(object): ################################################################ # GeneCollection goes here ################################################################ - if self.plotScale == 'physic': + if self.plotScale == 'physic' and self.selectedChr != -1: #StartMb or EndMb if self.startMb < 0 or self.endMb < 0: self.startMb = 0 @@ -594,7 +594,7 @@ class MarkerRegression(object): if self.traitList and self.traitList[0].dataset and self.traitList[0].dataset.type == 'Geno': btminfo.append(HT.BR(), 'Mapping using genotype data as a trait will result in infinity LRS at one locus. In order to display the result properly, all LRSs higher than 100 are capped at 100.') - if self.permChecked and not self.multipleInterval and 0 1: + + #if len(self.genotype) > 1: + if self.selectedChr == -1: #ZS: If viewing full genome/all chromosomes for i, _chr in enumerate(self.genotype): thisChr = [] Locus0CM = _chr[0].cM @@ -1759,15 +1761,16 @@ class MarkerRegression(object): ChrAInfo.append(thisChr) else: for i, _chr in enumerate(self.genotype): - thisChr = [] - Locus0CM = _chr[0].cM - for _locus in _chr: - if _locus.name != ' - ': - if _locus.cM != preLpos: - distinctCount += 1 - preLpos = _locus.cM - thisChr.append([_locus.name, _locus.cM-Locus0CM]) - ChrAInfo.append(thisChr) + if _chr.name == self.ChrList[self.selectedChr][0]: + thisChr = [] + Locus0CM = _chr[0].cM + for _locus in _chr: + if _locus.name != ' - ': + if _locus.cM != preLpos: + distinctCount += 1 + preLpos = _locus.cM + thisChr.append([_locus.name, _locus.cM-Locus0CM]) + ChrAInfo.append(thisChr) stepA = (plotWidth+0.0)/distinctCount @@ -1776,8 +1779,8 @@ class MarkerRegression(object): offsetA = -stepA lineColor = pid.lightblue startPosX = xLeftOffset + for j, ChrInfo in enumerate(ChrAInfo): - if ChrInfo == self.selectedChr: preLpos = -1 for i, item in enumerate(ChrInfo): Lname,Lpos = item @@ -1809,7 +1812,7 @@ class MarkerRegression(object): xLeftOffset+offsetA,yZero+40+Zorder*(LRectWidth+3)+LRectWidth) HREF="/show_trait?trait_id=%s&dataset=%s" % (Lname, self.dataset.group.name+"Geno") #HREF="javascript:showDatabase3('%s','%s','%s','');" % (showLocusForm,fd.RISet+"Geno", Lname) - Areas=HT.Area(shape='rect',coords=COORDS,href=HREF, title="Locus : " + Lname) + Areas=HT.Area(shape='rect', coords=COORDS, href=HREF, target="_blank", title="Locus : " + Lname) gifmap.areas.append(Areas) ##piddle bug if j == 0: @@ -1923,7 +1926,8 @@ class MarkerRegression(object): if self.multipleInterval: lrsEdgeWidth = 1 else: - additiveMax = max(map(lambda X : abs(X['additive']), self.qtlresults)) + if self.additiveChecked: + additiveMax = max(map(lambda X : abs(X['additive']), self.qtlresults)) #if INTERCROSS: # dominanceMax = max(map(lambda X : abs(X.dominance), self.qtlresults[0])) #else: @@ -1957,8 +1961,13 @@ class MarkerRegression(object): #startPosX += (self.ChrLengthDistList[j]+self.GraphInterval)*plotXScale - #for j, _chr in enumerate(self.genotype): - if self.selectedChr == -1 or qtlresult['chr'] == self.selectedChr: + #for j, _chr in enumerate(self.genotype): + #ZS: This is beause the chromosome value stored in qtlresult['chr'] can be (for example) either X or 20 depending upon the mapping method/scale used + if self.plotScale == "physic": + this_chr = str(self.ChrList[self.selectedChr][0]) + else: + this_chr = str(self.ChrList[self.selectedChr][1]+1) + if self.selectedChr == -1 or str(qtlresult['chr']) == this_chr: #LRSCoordXY = [] #AdditiveCoordXY = [] #DominanceCoordXY = [] @@ -2127,7 +2136,7 @@ class MarkerRegression(object): canvas.drawRect(startPix, yTopOffset, min(startPix+spacingAmt, xLeftOffset+plotWidth), \ yBottom, edgeColor=theBackColor, fillColor=theBackColor) - drawRegionDistance = self.ChrLengthDistList[self.selectedChr] + drawRegionDistance = self.ChrLengthDistList[self.ChrList[self.selectedChr][1]] self.ChrLengthDistList = [drawRegionDistance] if self.plotScale == 'physic': plotXScale = plotWidth / (endMb-startMb) diff --git a/wqflask/wqflask/templates/marker_regression_gn1.html b/wqflask/wqflask/templates/marker_regression_gn1.html index 550173fc..b2679fbf 100644 --- a/wqflask/wqflask/templates/marker_regression_gn1.html +++ b/wqflask/wqflask/templates/marker_regression_gn1.html @@ -39,7 +39,7 @@ Location: Chr {{ this_trait.chr }} @ {{ this_trait.mb }} Mb {% endif %} -
+
@@ -56,7 +56,7 @@ @@ -68,7 +68,9 @@
View:  - to + to
+ {% if mapping_method == "reaper" %} Allele Effects
+ {% endif %} Gene Track *
Legend
Haplotype Analyst * -- cgit v1.2.3