changeset 1:509993d9fdca draft default tip

"planemo upload for repository https://github.com/galaxyproject/tools-devteam/tree/master/tools/dwt_ivc_all commit f929353ffb0623f2218d7dec459c7da62f3b0d24"
author devteam
date Mon, 06 Jul 2020 18:12:29 +0000
parents 91fad0f30fd3
children
files execute_dwt_IvC_all.R execute_dwt_IvC_all.pl execute_dwt_IvC_all.xml test-data/in1.tsv test-data/in2.tsv test-data/out2.pdf
diffstat 6 files changed, 223 insertions(+), 221 deletions(-) [+]
line wrap: on
line diff
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/execute_dwt_IvC_all.R	Mon Jul 06 18:12:29 2020 +0000
@@ -0,0 +1,163 @@
+###########################################################################################
+## code to do wavelet Indel vs. Control
+## signal is the difference I-C; function is second moment i.e. variance from zero not mean
+## to perform wavelet transf. of signal, scale-by-scale analysis of the function
+## create null bands by permuting the original data series
+## generate plots and table matrix of correlation coefficients including p-values
+############################################################################################
+library("wavethresh");
+library("waveslim");
+
+options(echo = FALSE)
+
+## normalize data
+norm <- function(data) {
+    v <- (data - mean(data)) / sd(data);
+    if (sum(is.na(v)) >= 1) {
+        v <- data;
+    }
+    return(v);
+}
+
+dwt_cor <- function(data_short, names_short, data_long, names_long, test, pdf, table, filter = 4, bc = "symmetric", wf = "haar", boundary = "reflection") {
+    print(test);
+    print(pdf);
+    print(table);
+
+    pdf(file = pdf);
+    final_pvalue <- NULL;
+    title <- NULL;
+
+    short_levels <- wavethresh::wd(data_short[, 1], filter.number = filter, bc = bc)$nlevels;
+    title <- c("motif");
+    for (i in 1:short_levels) {
+        title <- c(title, paste(i, "moment2", sep = "_"), paste(i, "pval", sep = "_"), paste(i, "test", sep = "_"));
+    }
+    print(title);
+
+    ## loop to compare a vs a
+    for (i in seq_len(length(names_short))) {
+        wave1_dwt <- NULL;
+        m2_dwt <- NULL;
+        diff <- NULL;
+        var_dwt <- NULL;
+        out <- NULL;
+        out <- vector(length = length(title));
+
+        print(names_short[i]);
+        print(names_long[i]);
+
+        ## need exit if not comparing motif(a) vs motif(a)
+        if (names_short[i] != names_long[i]) {
+            stop(paste("motif", names_short[i], "is not the same as", names_long[i], sep = " "));
+        }
+        else {
+            ## signal is the difference I-C data sets
+            diff <- data_short[, i] - data_long[, i];
+
+            ## normalize the signal
+            diff <- norm(diff);
+
+            ## function is 2nd moment
+            ## 2nd moment m_j = 1/N[sum_N(W_j + V_J)^2] = 1/N sum_N(W_j)^2 + (X_bar)^2
+            wave1_dwt <- waveslim::dwt(diff, wf = wf, short_levels, boundary = boundary);
+            var_dwt <- waveslim::wave.variance(wave1_dwt);
+            m2_dwt <- vector(length = short_levels)
+            for (level in 1:short_levels) {
+                m2_dwt[level] <- var_dwt[level, 1] + (mean(diff)^2);
+            }
+
+            ## CI bands by permutation of time series
+            feature1 <- NULL;
+            feature2 <- NULL;
+            feature1 <- data_short[, i];
+            feature2 <- data_long[, i];
+            null <- NULL;
+            results <- NULL;
+            med <- NULL;
+            m2_25 <- NULL;
+            m2_975 <- NULL;
+
+            for (k in 1:1000) {
+                nk_1 <- NULL;
+                nk_2 <- NULL;
+                m2_null <- NULL;
+                var_null <- NULL;
+                null_levels <- NULL;
+                null_wave1 <- NULL;
+                null_diff <- NULL;
+                nk_1 <- sample(feature1, length(feature1), replace = FALSE);
+                nk_2 <- sample(feature2, length(feature2), replace = FALSE);
+                null_levels <- wavethresh::wd(nk_1, filter.number = filter, bc = bc)$nlevels;
+                null_diff <- nk_1 - nk_2;
+                null_diff <- norm(null_diff);
+                null_wave1 <- waveslim::dwt(null_diff, wf = wf, short_levels, boundary = boundary);
+                var_null <- waveslim::wave.variance(null_wave1);
+                m2_null <- vector(length = null_levels);
+                for (level in 1:null_levels) {
+                    m2_null[level] <- var_null[level, 1] + (mean(null_diff)^2);
+                }
+                null <- rbind(null, m2_null);
+            }
+
+            null <- apply(null, 2, sort, na.last = TRUE);
+            m2_25 <- null[25, ];
+            m2_975 <- null[975, ];
+            med <- apply(null, 2, median, na.rm = TRUE);
+
+            ## plot
+            results <- cbind(m2_dwt, m2_25, m2_975);
+            matplot(results, type = "b", pch = "*", lty = 1, col = c(1, 2, 2), xlab = "Wavelet Scale", ylab = c("Wavelet 2nd Moment", test), main = (names_short[i]), cex.main = 0.75);
+            abline(h = 1);
+
+            ## get pvalues by comparison to null distribution
+            out <- c(names_short[i]);
+            for (m in seq_len(length(m2_dwt))) {
+                print(paste("scale", m, sep = " "));
+                print(paste("m2", m2_dwt[m], sep = " "));
+                print(paste("median", med[m], sep = " "));
+                out <- c(out, format(m2_dwt[m], digits = 4));
+                pv <- NULL;
+                if (is.na(m2_dwt[m])) {
+                    pv <- "NA";
+                }
+                else {
+                    if (m2_dwt[m] >= med[m]) {
+                        ## R tail test
+                        tail <- "R";
+                        pv <- (length(which(null[, m] >= m2_dwt[m]))) / (length(na.exclude(null[, m])));
+                    }
+                    else{
+                        if (m2_dwt[m] < med[m]) {
+                            ## L tail test
+                            tail <- "L";
+                            pv <- (length(which(null[, m] <= m2_dwt[m]))) / (length(na.exclude(null[, m])));
+                        }
+                    }
+                }
+                out <- c(out, pv);
+                print(pv);
+                out <- c(out, tail);
+            }
+            final_pvalue <- rbind(final_pvalue, out);
+            print(out);
+        }
+    }
+
+    colnames(final_pvalue) <- title;
+    write.table(final_pvalue, file = table, sep = "\t", quote = FALSE, row.names = FALSE);
+    dev.off();
+}
+## execute
+## read in data
+args <- commandArgs(trailingOnly = TRUE)
+
+input_data <- read.delim(args[1]);
+input_data_names <- colnames(input_data);
+
+control_data <- read.delim(args[2]);
+control_data_names <- colnames(control_data);
+
+## call the test function to implement IvC test
+dwt_cor(input_data, input_data_names, control_data, control_data_names, test = "IvC", pdf = args[3], table = args[4]);
+print("done with the correlation test");
--- a/execute_dwt_IvC_all.pl	Thu Jan 23 12:31:01 2014 -0500
+++ /dev/null	Thu Jan 01 00:00:00 1970 +0000
@@ -1,210 +0,0 @@
-#!/usr/bin/perl -w
-use warnings;
-use IO::Handle;
-
-$usage = "execute_dwt_IvC_all.pl [TABULAR.in] [TABULAR.in] [TABULAR.out] [PDF.out]  \n";
-die $usage unless @ARGV == 4;
-
-#get the input arguments
-my $firstInputFile = $ARGV[0];
-my $secondInputFile = $ARGV[1];
-my $firstOutputFile = $ARGV[2];
-my $secondOutputFile = $ARGV[3];
-
-open (INPUT1, "<", $firstInputFile) || die("Could not open file $firstInputFile \n");
-open (INPUT2, "<", $secondInputFile) || die("Could not open file $secondInputFile \n");
-open (OUTPUT1, ">", $firstOutputFile) || die("Could not open file $firstOutputFile \n");
-open (OUTPUT2, ">", $secondOutputFile) || die("Could not open file $secondOutputFile \n");
-open (ERROR,  ">", "error.txt")  or die ("Could not open file error.txt \n");
-
-#save all error messages into the error file $errorFile using the error file handle ERROR
-STDERR -> fdopen( \*ERROR,  "w" ) or die ("Could not direct errors to the error file error.txt \n");
-
-
-print "There are two input data files: \n";
-print "The input data file is: $firstInputFile \n";
-print "The control data file is: $secondInputFile \n";
-
-# IvC test
-$test = "IvC";
-
-# construct an R script to implement the IvC test
-print "\n";
-
-$r_script = "get_dwt_IvC_test.r"; 
-print "$r_script \n";
-
-# R script
-open(Rcmd, ">", "$r_script") or die "Cannot open $r_script \n\n";
-print Rcmd "
-        ###########################################################################################
-        # code to do wavelet Indel vs. Control
-        # signal is the difference I-C; function is second moment i.e. variance from zero not mean
-        # to perform wavelet transf. of signal, scale-by-scale analysis of the function 
-        # create null bands by permuting the original data series
-        # generate plots and table matrix of correlation coefficients including p-values
-        ############################################################################################
-        library(\"Rwave\");
-        library(\"wavethresh\");
-        library(\"waveslim\");
-        
-        options(echo = FALSE)
-        
-        # normalize data
-        norm <- function(data){
-            v <- (data - mean(data))/sd(data);
-            if(sum(is.na(v)) >= 1){
-                v <- data;
-            }
-            return(v);
-        }
-        
-        dwt_cor <- function(data.short, names.short, data.long, names.long, test, pdf, table, filter = 4, bc = \"symmetric\", wf = \"haar\", boundary = \"reflection\") {
-            print(test);
-            print(pdf);
-            print(table);
-            
-            pdf(file = pdf);
-            final_pvalue = NULL;
-            title = NULL;
-                
-            short.levels <- wd(data.short[, 1], filter.number = filter, bc = bc)\$nlevels;
-            title <- c(\"motif\");
-            for (i in 1:short.levels){
-            	title <- c(title, paste(i, \"moment2\", sep = \"_\"), paste(i, \"pval\", sep = \"_\"), paste(i, \"test\", sep = \"_\"));
-            }
-            print(title);
-        
-            # loop to compare a vs a
-            for(i in 1:length(names.short)){
-        		wave1.dwt = NULL;
-        		m2.dwt = diff = var.dwt = NULL;
-        		out = NULL;
-                out <- vector(length = length(title));
-        
-        		print(names.short[i]);
-        		print(names.long[i]);
-                        
-        		# need exit if not comparing motif(a) vs motif(a)
-        		if (names.short[i] != names.long[i]){
-                	stop(paste(\"motif\", names.short[i], \"is not the same as\", names.long[i], sep = \" \"));
-        		}
-        		else {
-                	# signal is the difference I-C data sets
-                    diff<-data.short[,i]-data.long[,i];
-        
-                    # normalize the signal
-                    diff<-norm(diff);
-        
-                    # function is 2nd moment
-                    # 2nd moment m_j = 1/N[sum_N(W_j + V_J)^2] = 1/N sum_N(W_j)^2 + (X_bar)^2 
-            		wave1.dwt <- dwt(diff, wf = wf, short.levels, boundary = boundary);
-            		var.dwt <- wave.variance(wave1.dwt);
-                	m2.dwt <- vector(length = short.levels)
-                    for(level in 1:short.levels){
-                    	m2.dwt[level] <- var.dwt[level, 1] + (mean(diff)^2);
-                    }
-                                
-            		# CI bands by permutation of time series
-            		feature1 = feature2 = NULL;
-            		feature1 = data.short[, i];
-            		feature2 = data.long[, i];
-            		null = results = med = NULL; 
-            		m2_25 = m2_975 = NULL;
-            
-            		for (k in 1:1000) {
-                		nk_1 = nk_2 = NULL;
-                		m2_null = var_null = NULL;
-                		null.levels = null_wave1 = null_diff = NULL;
-                		nk_1 <- sample(feature1, length(feature1), replace = FALSE);
-                		nk_2 <- sample(feature2, length(feature2), replace = FALSE);
-                		null.levels <- wd(nk_1, filter.number = filter, bc = bc)\$nlevels;
-                		null_diff <- nk_1-nk_2;
-                		null_diff <- norm(null_diff);
-                		null_wave1 <- dwt(null_diff, wf = wf, short.levels, boundary = boundary);
-                        var_null <- wave.variance(null_wave1);
-                		m2_null <- vector(length = null.levels);
-                		for(level in 1:null.levels){
-                        	m2_null[level] <- var_null[level, 1] + (mean(null_diff)^2);
-                		}
-                		null= rbind(null, m2_null);
-            		}
-                
-            		null <- apply(null, 2, sort, na.last = TRUE);
-            		m2_25 <- null[25,];
-            		m2_975 <- null[975,];
-            		med <- apply(null, 2, median, na.rm = TRUE);
-
-            		# plot
-            		results <- cbind(m2.dwt, m2_25, m2_975);
-            		matplot(results, type = \"b\", pch = \"*\", lty = 1, col = c(1, 2, 2), xlab = \"Wavelet Scale\", ylab = c(\"Wavelet 2nd Moment\", test), main = (names.short[i]), cex.main = 0.75);
-            		abline(h = 1);
-
-            		# get pvalues by comparison to null distribution
-            		out <- c(names.short[i]);
-            		for (m in 1:length(m2.dwt)){
-                    	print(paste(\"scale\", m, sep = \" \"));
-                        print(paste(\"m2\", m2.dwt[m], sep = \" \"));
-                        print(paste(\"median\", med[m], sep = \" \"));
-                        out <- c(out, format(m2.dwt[m], digits = 4));	
-                        pv = NULL;
-                        if(is.na(m2.dwt[m])){
-                        	pv <- \"NA\"; 
-                        } 
-                        else {
-                        	if (m2.dwt[m] >= med[m]){
-                            	# R tail test
-                                tail <- \"R\";
-                                pv <- (length(which(null[, m] >= m2.dwt[m])))/(length(na.exclude(null[, m])));
-                            }
-                            else{
-                                if (m2.dwt[m] < med[m]){
-                                	# L tail test
-                                    tail <- \"L\";
-                                    pv <- (length(which(null[, m] <= m2.dwt[m])))/(length(na.exclude(null[, m])));
-                                }
-                            }
-                        }
-                        out <- c(out, pv);
-                        print(pv);  
-                        out <- c(out, tail);
-                    }
-                    final_pvalue <-rbind(final_pvalue, out);
-                    print(out);
-                }
-            }
-            
-            colnames(final_pvalue) <- title;
-            write.table(final_pvalue, file = table, sep = \"\\t\", quote = FALSE, row.names = FALSE);
-            dev.off();
-        }\n";
-
-print Rcmd "
-        # execute
-        # read in data 
-        
-        inputData <- read.delim(\"$firstInputFile\");
-        inputDataNames <- colnames(inputData);
-        
-        controlData <- read.delim(\"$secondInputFile\");
-        controlDataNames <- colnames(controlData);
-        
-        # call the test function to implement IvC test
-        dwt_cor(inputData, inputDataNames, controlData, controlDataNames, test = \"$test\", pdf = \"$secondOutputFile\", table = \"$firstOutputFile\");
-        print (\"done with the correlation test\");
-\n";
-
-print Rcmd "#eof\n";
-
-close Rcmd;
-
-system("echo \"wavelet IvC test started on \`hostname\` at \`date\`\"\n");
-system("R --no-restore --no-save --no-readline < $r_script > $r_script.out\n");
-system("echo \"wavelet IvC test ended on \`hostname\` at \`date\`\"\n");
-
-#close the input and output and error files
-close(ERROR);
-close(OUTPUT2);
-close(OUTPUT1);
-close(INPUT2);
-close(INPUT1);
\ No newline at end of file
--- a/execute_dwt_IvC_all.xml	Thu Jan 23 12:31:01 2014 -0500
+++ b/execute_dwt_IvC_all.xml	Mon Jul 06 18:12:29 2020 +0000
@@ -1,20 +1,35 @@
-<tool id="compute_p-values_second_moments_feature_occurrences_between_two_datasets_using_discrete_wavelet_transfom" name="Compute P-values and Second Moments for Feature Occurrences" version="1.0.0">
+<tool id="compute_p-values_second_moments_feature_occurrences_between_two_datasets_using_discrete_wavelet_transfom" name="Compute P-values and Second Moments for Feature Occurrences" version="1.0.1">
   <description>between two datasets using Discrete Wavelet Transfoms</description>
-  
-  <command interpreter="perl">
-  	execute_dwt_IvC_all.pl $inputFile1 $inputFile2 $outputFile1 $outputFile2
+  <requirements>
+    <requirement type="package" version="1.7.5">r-waveslim</requirement>
+    <requirement type="package" version="4.6.8">r-wavethresh</requirement>
+  </requirements>
+  <command detect_errors="exit_code">
+    Rscript --vanilla '$__tool_directory__/execute_dwt_IvC_all.R'
+      '$inputFile1'
+      '$inputFile2'
+      '$outputFile2'
+      '$outputFile1'
   </command>
-
   <inputs>
-  	<param format="tabular" name="inputFile1" type="data" label="Select the first input file"/>	
-  	<param format="tabular" name="inputFile2" type="data" label="Select the second input file"/>
+    <param format="tabular" name="inputFile1" type="data" label="Select the first input file"/>
+    <param format="tabular" name="inputFile2" type="data" label="Select the second input file"/>
   </inputs>
-  
   <outputs>
-    <data format="tabular" name="outputFile1"/> 
-    <data format="pdf" name="outputFile2"/>
+    <data format="tabular" name="outputFile1" label="${tool.name} on ${on_string}: statistics"/> 
+    <data format="pdf" name="outputFile2" label="${tool.name} on ${on_string}: pdf"/>
   </outputs>
-  	
+  <tests>
+    <test>
+      <param ftype="tabular" name="inputFile1" value="in1.tsv"/>
+      <param ftype="tabular" name="inputFile2" value="in2.tsv"/>
+      <output name="outputFile1" ftype="tabular">
+        <assert_contents><has_line_matching expression="^motif\t1_moment2.*"/></assert_contents>
+        <assert_contents><has_line_matching expression="^translinTarget.*" /></assert_contents>
+      </output>
+      <output name="outputFile2" ftype="pdf" file="out2.pdf" compare="sim_size"/>
+    </test>
+  </tests>
   <help> 
 
 .. class:: infomark
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/test-data/in1.tsv	Mon Jul 06 18:12:29 2020 +0000
@@ -0,0 +1,17 @@
+deletionHoptspot	insertionHoptspot	dnaPolPauseFrameshift	topoisomeraseCleavageSite	translinTarget
+269	366	330	238	1129
+239	328	327	283	1188
+254	351	358	297	1151
+262	371	355	256	1107
+254	361	352	234	1192
+265	354	367	240	1182
+255	359	333	235	1217
+271	389	387	272	1241
+240	305	341	249	1159
+272	351	337	257	1169
+275	351	337	233	1158
+305	331	361	253	1172
+277	341	343	253	1113
+266	362	355	267	1162
+235	326	329	241	1230
+254	335	360	251	1172
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/test-data/in2.tsv	Mon Jul 06 18:12:29 2020 +0000
@@ -0,0 +1,17 @@
+deletionHoptspot	insertionHoptspot	dnaPolPauseFrameshift	topoisomeraseCleavageSite	translinTarget
+104	146	142	113	478
+89	146	151	94	495
+100	176	151	88	435
+96	163	128	114	468
+99	138	144	91	513
+112	126	162	106	468
+86	127	145	83	491
+104	145	171	110	496
+91	121	147	104	469
+103	141	145	98	458
+92	134	142	117	468
+97	146	145	107	471
+115	121	136	109	470
+113	135	138	101	491
+111	150	138	102	451
+94	128	151	138	481
Binary file test-data/out2.pdf has changed