From 08185087a75e8159a6f6565d8e5ad0e3a8d4105c Mon Sep 17 00:00:00 2001 From: Guruprasad Anada Date: Thu, 3 Nov 2011 10:49:45 -0400 Subject: [PATCH] Added 3 tools: 1) to compute weighted average values based on interval intersections (under Regional variation), 2) to perform logistic regression (under Multiple regression), 3) to compute partial R-squared (under Multiple regression). --- tool_conf.xml.sample | 3 + tools/regVariation/WeightedAverage.py | 103 ++++++++++ tools/regVariation/WeightedAverage.xml | 71 +++++++ tools/regVariation/logistic_regression_vif.py | 185 ++++++++++++++++++ .../regVariation/logistic_regression_vif.xml | 74 +++++++ tools/regVariation/partialR_square.py | 146 ++++++++++++++ tools/regVariation/partialR_square.xml | 68 +++++++ 7 files changed, 650 insertions(+) create mode 100755 tools/regVariation/WeightedAverage.py create mode 100755 tools/regVariation/WeightedAverage.xml create mode 100755 tools/regVariation/logistic_regression_vif.py create mode 100755 tools/regVariation/logistic_regression_vif.xml create mode 100755 tools/regVariation/partialR_square.py create mode 100755 tools/regVariation/partialR_square.xml diff --git a/tool_conf.xml.sample b/tool_conf.xml.sample index d309da3e364..d189a3b2319 100644 --- a/tool_conf.xml.sample +++ b/tool_conf.xml.sample @@ -175,6 +175,7 @@
+ @@ -194,8 +195,10 @@
+ +
diff --git a/tools/regVariation/WeightedAverage.py b/tools/regVariation/WeightedAverage.py new file mode 100755 index 00000000000..d0299dd18ed --- /dev/null +++ b/tools/regVariation/WeightedAverage.py @@ -0,0 +1,103 @@ +#!/usr/bin/env python +""" + +usage: %prog bed_file_1 bed_file_2 out_file + -1, --cols1=N,N,N,N: Columns for chr, start, end, strand in first file + -2, --cols2=N,N,N,N,N: Columns for chr, start, end, strand, name/value in second file +""" +from galaxy import eggs + +import collections +import sys, string +#import numpy +from galaxy import eggs +import pkg_resources +pkg_resources.require( "bx-python" ) +import sys, traceback, fileinput +from warnings import warn +from galaxy.tools.util.galaxyops import * +from bx.cookbook import doc_optparse + + +#export PYTHONPATH=~/galaxy/lib/ +#running command python WeightedAverage.py interval_interpolate.bed value_interpolate.bed interpolate_result.bed + +def stop_err(msg): + sys.stderr.write(msg) + sys.exit() + + +def FindRate(chromosome,start_stop,dictType): + OverlapList=[] + for tempO in dictType[chromosome]: + DatabaseInterval=[tempO[0],tempO[1]] + Overlap=GetOverlap(start_stop,DatabaseInterval) + if Overlap>0: + OverlapList.append([Overlap,tempO[2]]) + + if len(OverlapList)>0: + + SumRecomb=0 + SumOverlap=0 + for member in OverlapList: + SumRecomb+=member[0]*member[1] + SumOverlap+=member[0] + averageRate=SumRecomb/SumOverlap + + return averageRate + + else: + return 'NA' + + + +def GetOverlap(a,b): + return min(a[1],b[1])-max(a[0],b[0]) + +options, args = doc_optparse.parse( __doc__ ) + +try: + chr_col_1, start_col_1, end_col_1, strand_col1 = parse_cols_arg( options.cols1 ) + chr_col_2, start_col_2, end_col_2, strand_col2, name_col_2 = parse_cols_arg( options.cols2 ) + input1, input2, input3 = args +except Exception, eee: + print eee + stop_err( "Data issue: click the pencil icon in the history item to correct the metadata attributes." ) + + + +fd2=open(input2) +lines2=fd2.readlines() +RecombChrDict=collections.defaultdict(list) + +skipped=0 +for line in lines2: + temp=line.strip().split() + try: + assert float(temp[int(name_col_2)]) + except: + skipped+=1 + continue + tempIndex=[int(temp[int(start_col_2)]),int(temp[int(end_col_2)]),float(temp[int(name_col_2)])] + RecombChrDict[temp[int(chr_col_2)]].append(tempIndex) + +print "Skipped %d features with invalid values" %(skipped) + +fd1=open(input1) +lines=fd1.readlines() +finalProduct='' +for line in lines: + temp=line.strip().split('\t') + chromosome=temp[int(chr_col_1)] + start=int(temp[int(start_col_1)]) + stop=int(temp[int(end_col_1)]) + start_stop=[start,stop] + RecombRate=FindRate(chromosome,start_stop,RecombChrDict) + try: + RecombRate="%.4f" %(float(RecombRate)) + except: + RecombRate=RecombRate + finalProduct+=line.strip()+'\t'+str(RecombRate)+'\n' +fdd=open(input3,'w') +fdd.writelines(finalProduct) +fdd.close() diff --git a/tools/regVariation/WeightedAverage.xml b/tools/regVariation/WeightedAverage.xml new file mode 100755 index 00000000000..9e9406985d2 --- /dev/null +++ b/tools/regVariation/WeightedAverage.xml @@ -0,0 +1,71 @@ + + of the values of features overlapping an interval + WeightedAverage.py $genomic_interval $genomic_feature $out_file1 -1 ${genomic_interval.metadata.chromCol},${genomic_interval.metadata.startCol},${genomic_interval.metadata.endCol},${genomic_interval.metadata.strandCol} -2 ${genomic_feature.metadata.chromCol},${genomic_feature.metadata.startCol},${genomic_feature.metadata.endCol},${genomic_feature.metadata.strandCol},${genomic_feature.metadata.nameCol} + + + + + + + + + + + + + + + + + + + + +.. class:: infomark + +**What it does** + +For each interval in your first dataset, this tool calculates the weighted average value of the overlapping features in your second dataset. + +- When a genomic interval partially or totally overlaps a single genomic feature, the value of that genomic feature is assigned to the genomic interval. +- When a genomic interval partially or totally overlaps with more than one genomic features, the average of the values of the overlapping genomic features weighted by the corresponding number of overlapping bases is assigned to the genomic interval. +- When a genomic interval does not overlap with any genomic feature, 'NA' will be assigned as it's value. + +----- + +.. class:: warningmark + +**Note** + +The input datasets should be in **bed** or **interval** format. Please use "edit attributes"/pencil icon to specify the column containing the values for the features in the second dataset as **name/identifier** column. + +The output will contain all the columns in the first input plus a new column containing the assigned value for each interval. + +----- + +**Example** + +- Suppose our first dataset contains the following **genomic intervals**:: + + chr start stop + chr1 1000 2000 + chr1 3000 5000 + chr1 8000 9000 + +- and our second dataset contains the following **genomic features** each having an associated value (in fourth column) :: + + chr start stop name + chr1 900 1200 0.5 + chr1 2900 3100 0.2 + chr1 4800 5100 0.8 + +- For each **genomic interval** in our first dataset, this tool calculates the weighted average value of the overlapping **genomic features** in our second dataset :: + + chr1 1000 2000 0.5 + chr1 3000 5000 0.6 + chr1 8000 9000 NA + + + + + \ No newline at end of file diff --git a/tools/regVariation/logistic_regression_vif.py b/tools/regVariation/logistic_regression_vif.py new file mode 100755 index 00000000000..b4198c6588f --- /dev/null +++ b/tools/regVariation/logistic_regression_vif.py @@ -0,0 +1,185 @@ +#!/usr/bin/env python + +from galaxy import eggs +import sys, string +from rpy import * +import numpy + +#export PYTHONPATH=~/galaxy/lib/ + +def stop_err(msg): + sys.stderr.write(msg) + sys.exit() +#infile = 'logreg_inp.tab' +#y_col=3 +#x_cols=[1,2,3] +#outfile='logreg_out.txt' +#python logistic_regression_vif.py logreg_inp.tab 4 1,2,3 logreg_out2.tabular # running test +infile = sys.argv[1] +y_col = int(sys.argv[2])-1 +x_cols = sys.argv[3].split(',') +outfile = sys.argv[4] + + +print "Predictor columns: %s; Response column: %d" %(x_cols,y_col+1) +fout = open(outfile,'w') +elems = [] +for i, line in enumerate( file ( infile )): + line = line.rstrip('\r\n') + if len( line )>0 and not line.startswith( '#' ): + elems = line.split( '\t' ) + break + if i == 30: + break # Hopefully we'll never get here... + +if len( elems )<1: + stop_err( "The data in your input dataset is either missing or not formatted properly." ) + +y_vals = [] +x_vals = [] + +for k,col in enumerate(x_cols): + x_cols[k] = int(col)-1 + x_vals.append([]) + +NA = 'NA' +for ind,line in enumerate( file( infile )): + if line and not line.startswith( '#' ): + try: + fields = line.split("\t") + try: + yval = float(fields[y_col]) + except: + yval = r('NA') + y_vals.append(yval) + for k,col in enumerate(x_cols): + try: + xval = float(fields[col]) + except: + xval = r('NA') + x_vals[k].append(xval) + except: + pass + +x_vals1 = numpy.asarray(x_vals).transpose() + +check1=0 +check0=0 +for i in y_vals: + if i == 1: + check1=1 + if i == 0: + check0=1 +if check1==0 or check0==0: + sys.exit("Warning: logistic regression must have at least two classes") + +for i in y_vals: + if i not in [1,0,r('NA')]: + print >>fout, str(i) + sys.exit("Warning: the current version of this tool can run only with two classes and need to be labeled as 0 and 1.") + + +dat= r.list(x=array(x_vals1), y=y_vals) +novif=0 +set_default_mode(NO_CONVERSION) +try: + linear_model = r.glm(r("y ~ x"), data = r.na_exclude(dat),family="binomial") + #r('library(car)') + #r.assign('dat',dat) + #r.assign('ncols',len(x_cols)) + #r.vif(r('glm(dat$y ~ ., data = na.exclude(data.frame(as.matrix(dat$x,ncol=ncols))->datx),family="binomial")')).as_py() + +except RException, rex: + stop_err("Error performing logistic regression on the input data.\nEither the response column or one of the predictor columns contain only non-numeric or invalid values.") +if len(x_cols)>1: + try: + + r('library(car)') + r.assign('dat',dat) + r.assign('ncols',len(x_cols)) + vif=r.vif(r('glm(dat$y ~ ., data = na.exclude(data.frame(as.matrix(dat$x,ncol=ncols))->datx),family="binomial")')) + except RException, rex: + print rex +else: + novif=1 + +set_default_mode(BASIC_CONVERSION) + +coeffs=linear_model.as_py()['coefficients'] +null_deviance=linear_model.as_py()['null.deviance'] +residual_deviance=linear_model.as_py()['deviance'] +yintercept= coeffs['(Intercept)'] +summary = r.summary(linear_model) +co = summary.get('coefficients', 'NA') +""" +if len(co) != len(x_vals)+1: + stop_err("Stopped performing logistic regression on the input data, since one of the predictor columns contains only non-numeric or invalid values.") +""" + +try: + yintercept = r.round(float(yintercept), digits=10) + pvaly = r.round(float(co[0][3]), digits=10) +except: + pass +print >>fout, "response column\tc%d" %(y_col+1) +tempP=[] +for i in x_cols: + tempP.append('c'+str(i+1)) +tempP=','.join(tempP) +print >>fout, "predictor column(s)\t%s" %(tempP) +print >>fout, "Y-intercept\t%s" %(yintercept) +print >>fout, "p-value (Y-intercept)\t%s" %(pvaly) + +if len(x_vals) == 1: #Simple linear regression case with 1 predictor variable + try: + slope = r.round(float(coeffs['x']), digits=10) + except: + slope = 'NA' + try: + pval = r.round(float(co[1][3]), digits=10) + except: + pval = 'NA' + print >>fout, "Slope (c%d)\t%s" %(x_cols[0]+1,slope) + print >>fout, "p-value (c%d)\t%s" %(x_cols[0]+1,pval) +else: #Multiple regression case with >1 predictors + ind=1 + while ind < len(coeffs.keys()): + try: + slope = r.round(float(coeffs['x'+str(ind)]), digits=10) + except: + slope = 'NA' + print >>fout, "Slope (c%d)\t%s" %(x_cols[ind-1]+1,slope) + try: + pval = r.round(float(co[ind][3]), digits=10) + except: + pval = 'NA' + print >>fout, "p-value (c%d)\t%s" %(x_cols[ind-1]+1,pval) + ind+=1 + +rsq = summary.get('r.squared','NA') + + +try: + rsq= r.round(float((null_deviance-residual_deviance)/null_deviance), digits=5) + null_deviance= r.round(float(null_deviance), digits=5) + residual_deviance= r.round(float(residual_deviance), digits=5) + + #rsq = r.round(float(rsq), digits=5) + +except: + pass + +print >>fout, "Null deviance\t%s" %(null_deviance) +print >>fout, "Residual deviance\t%s" %(residual_deviance) +print >>fout, "pseudo R-squared\t%s" %(rsq) +print >>fout, "\n" +print >>fout, 'vif' + +if novif==0: + py_vif=vif.as_py() + count=0 + for i in sorted(py_vif.keys()): + print >>fout,'c'+str(x_cols[count]+1) ,str(py_vif[i]) + count+=1 +elif novif==1: + print >>fout, "vif can calculate only when model have more than 1 predictor" diff --git a/tools/regVariation/logistic_regression_vif.xml b/tools/regVariation/logistic_regression_vif.xml new file mode 100755 index 00000000000..2b491bd1aac --- /dev/null +++ b/tools/regVariation/logistic_regression_vif.xml @@ -0,0 +1,74 @@ + + + + logistic_regression_vif.py + $input1 + $response_col + $predictor_cols + $out_file1 + 1>/dev/null + + + + + + + + + + + + + + rpy + + + + + + + + + + + + + +.. class:: infomark + +**TIP:** If your data is not TAB delimited, use *Edit Datasets->Convert characters* + +----- + +.. class:: infomark + +**What it does** + +This tool uses the **'glm'** function from R statistical package to perform logistic regression on the input data. It outputs one file containing the summary statistics of the performed regression. Also, it calculates VIF(Variance Inflation Factor) with **'vif'** function from library (car) in R. + + +*R Development Core Team (2010). R: A language and environment for statistical computing. R Foundation for Statistical Computing, Vienna, Austria. ISBN 3-900051-07-0, URL http://www.R-project.org.* + +----- + +.. class:: warningmark + +**Note** + +- This tool currently treats all predictor variables as continuous numeric variables and response variable as categorical variable. Currently, the response variable can have only two classes, namely 0 and 1. The program will take 0 as base class. + +- Rows containing non-numeric (or missing) data in any of the chosen columns will be skipped from the analysis. + +- The summary statistics in the output are described below: + +- Pseudo R-squared: the proportion of model improvement from null model +- p-value: p-value for the z-test of the null hypothesis that the corresponding slope is equal to zero against the two-sided alternative. +- Coefficient indicates log ratio of (probability to be class 1 / probability to be class 0) + +- This tool also provides **Variance Inflation Factor or VIF** which quantifies the level of multicollinearity. The tool will automatic generate VIF if the model has more than one predictor. The higher the VIF, the higher is the multicollinearity. Multicollinearity will inflate standard error and reduce level of significance of the predictor. In the worst case, it can reverse direction of slope for highly correlated predictors if one of them is significant. A general thumb-rule is to use those predictors having VIF lower than 10 or 5. +- **vif** is calculated by + - First, regressing each predictor over all other predictors, and recording R-squared for each regression. + - Second, computing vif as 1/(1- R_squared) + + + diff --git a/tools/regVariation/partialR_square.py b/tools/regVariation/partialR_square.py new file mode 100755 index 00000000000..326a7bef41a --- /dev/null +++ b/tools/regVariation/partialR_square.py @@ -0,0 +1,146 @@ +#!/usr/bin/env python + +from galaxy import eggs + +import sys, string +from rpy import * +import numpy + +#export PYTHONPATH=~/galaxy/lib/ +#running command python partialR_square.py reg_inp.tab 4 1,2,3 partialR_result.tabular + +def stop_err(msg): + sys.stderr.write(msg) + sys.exit() + +def sscombs(s): + if len(s) == 1: + return [s] + else: + ssc = sscombs(s[1:]) + return [s[0]] + [s[0]+comb for comb in ssc] + ssc + + +infile = sys.argv[1] +y_col = int(sys.argv[2])-1 +x_cols = sys.argv[3].split(',') +outfile = sys.argv[4] + +print "Predictor columns: %s; Response column: %d" %(x_cols,y_col+1) +fout = open(outfile,'w') + +for i, line in enumerate( file ( infile )): + line = line.rstrip('\r\n') + if len( line )>0 and not line.startswith( '#' ): + elems = line.split( '\t' ) + break + if i == 30: + break # Hopefully we'll never get here... + +if len( elems )<1: + stop_err( "The data in your input dataset is either missing or not formatted properly." ) + +y_vals = [] +x_vals = [] + +for k,col in enumerate(x_cols): + x_cols[k] = int(col)-1 + x_vals.append([]) + """ + try: + float( elems[x_cols[k]] ) + except: + try: + msg = "This operation cannot be performed on non-numeric column %d containing value '%s'." %( col, elems[x_cols[k]] ) + except: + msg = "This operation cannot be performed on non-numeric data." + stop_err( msg ) + """ +NA = 'NA' +for ind,line in enumerate( file( infile )): + if line and not line.startswith( '#' ): + try: + fields = line.split("\t") + try: + yval = float(fields[y_col]) + except Exception, ey: + yval = r('NA') + #print >>sys.stderr, "ey = %s" %ey + y_vals.append(yval) + for k,col in enumerate(x_cols): + try: + xval = float(fields[col]) + except Exception, ex: + xval = r('NA') + #print >>sys.stderr, "ex = %s" %ex + x_vals[k].append(xval) + except: + pass + +x_vals1 = numpy.asarray(x_vals).transpose() +dat= r.list(x=array(x_vals1), y=y_vals) + +set_default_mode(NO_CONVERSION) +try: + full = r.lm(r("y ~ x"), data= r.na_exclude(dat)) #full model includes all the predictor variables specified by the user +except RException, rex: + stop_err("Error performing linear regression on the input data.\nEither the response column or one of the predictor columns contain no numeric values.") +set_default_mode(BASIC_CONVERSION) + +summary = r.summary(full) +fullr2 = summary.get('r.squared','NA') + +if fullr2 == 'NA': + stop_error("Error in linear regression") + +if len(x_vals) < 10: + s = "" + for ch in range(len(x_vals)): + s += str(ch) +else: + stop_err("This tool only works with less than 10 predictors.") + +print >>fout, "#Model\tR-sq\tpartial_R_Terms\tpartial_R_Value" +all_combos = sorted(sscombs(s), key=len) +all_combos.reverse() +for j,cols in enumerate(all_combos): + #if len(cols) == len(s): #Same as the full model above + # continue + if len(cols) == 1: + x_vals1 = x_vals[int(cols)] + else: + x_v = [] + for col in cols: + x_v.append(x_vals[int(col)]) + x_vals1 = numpy.asarray(x_v).transpose() + dat= r.list(x=array(x_vals1), y=y_vals) + set_default_mode(NO_CONVERSION) + red = r.lm(r("y ~ x"), data= dat) #Reduced model + set_default_mode(BASIC_CONVERSION) + summary = r.summary(red) + redr2 = summary.get('r.squared','NA') + try: + partial_R = (float(fullr2)-float(redr2))/(1-float(redr2)) + except: + partial_R = 'NA' + col_str = "" + for col in cols: + col_str = col_str + str(int(x_cols[int(col)]) + 1) + " " + col_str.strip() + partial_R_col_str = "" + for col in s: + if col not in cols: + partial_R_col_str = partial_R_col_str + str(int(x_cols[int(col)]) + 1) + " " + partial_R_col_str.strip() + if len(cols) == len(s): #full model + partial_R_col_str = "-" + partial_R = "-" + try: + redr2 = "%.4f" %(float(redr2)) + except: + pass + try: + partial_R = "%.4f" %(float(partial_R)) + except: + pass + print >>fout, "%s\t%s\t%s\t%s" %(col_str,redr2,partial_R_col_str,partial_R) diff --git a/tools/regVariation/partialR_square.xml b/tools/regVariation/partialR_square.xml new file mode 100755 index 00000000000..4068a07a374 --- /dev/null +++ b/tools/regVariation/partialR_square.xml @@ -0,0 +1,68 @@ + + + + partialR_square.py + $input1 + $response_col + $predictor_cols + $out_file1 + 1>/dev/null + + + + + + + + + + + + + rpy + + + + + + + + + + + + + +.. class:: infomark + +**TIP:** If your data is not TAB delimited, use *Edit Datasets->Convert characters* + +----- + +.. class:: infomark + +**What it does** + +This tool computes the Partial R squared for all possible variable subsets using the following formula: + +**Partial R squared = [SSE(without i: 1,2,...,p-1) - SSE (full: 1,2,..,i..,p-1) / SSE(without i: 1,2,...,p-1)]**, which denotes the case where the 'i'th predictor is dropped. + + + +In general, **Partial R squared = [SSE(without i: 1,2,...,p-1) - SSE (full: 1,2,..,i..,p-1) / SSE(without i: 1,2,...,p-1)]**, where, + +- SSE (full: 1,2,..,i..,p-1) = Sum of Squares left out by the full set of predictors SSE(X1, X2 … Xp) +- SSE (full: 1,2,..,i..,p-1) = Sum of Squares left out by the set of predictors excluding; for example, if we omit the first predictor, it will be SSE(X2 … Xp). + + +The 4 columns in the output are described below: + +- Column 1 (Model): denotes the variables present in the model +- Column 2 (R-sq): denotes the R-squared value corresponding to the model in Column 1 +- Column 3 (Partial R squared_Terms): denotes the variable/s for which Partial R squared is computed. These are the variables that are absent in the reduced model in Column 1. A '-' in this column indicates that the model in Column 1 is the Full model. +- Column 4 (Partial R squared): denotes the Partial R squared value corresponding to the variable/s in Column 3. A '-' in this column indicates that the model in Column 1 is the Full model. + +*R Development Core Team (2010). R: A language and environment for statistical computing. R Foundation for Statistical Computing, Vienna, Austria. ISBN 3-900051-07-0, URL http://www.R-project.org.* + + +