mirror of
https://github.com/galaxyproject/galaxy.git
synced 2026-09-24 16:30:27 +08:00
Adding wavelet analysis tools made by Zaky from Makova lab.
This commit is contained in:
Guruprasad Anada
Binary file not shown.
|
After Width: | Height: | Size: 139 KiB |
Binary file not shown.
|
After Width: | Height: | Size: 140 KiB |
Binary file not shown.
|
After Width: | Height: | Size: 140 KiB |
Binary file not shown.
|
After Width: | Height: | Size: 142 KiB |
Binary file not shown.
|
After Width: | Height: | Size: 153 KiB |
Binary file not shown.
|
After Width: | Height: | Size: 157 KiB |
Binary file not shown.
|
After Width: | Height: | Size: 158 KiB |
Binary file not shown.
|
After Width: | Height: | Size: 156 KiB |
Binary file not shown.
|
After Width: | Height: | Size: 154 KiB |
Binary file not shown.
|
After Width: | Height: | Size: 140 KiB |
Binary file not shown.
|
After Width: | Height: | Size: 150 KiB |
Binary file not shown.
|
After Width: | Height: | Size: 165 KiB |
Binary file not shown.
|
After Width: | Height: | Size: 153 KiB |
Binary file not shown.
|
After Width: | Height: | Size: 163 KiB |
Binary file not shown.
|
After Width: | Height: | Size: 147 KiB |
Binary file not shown.
|
After Width: | Height: | Size: 169 KiB |
Binary file not shown.
|
After Width: | Height: | Size: 155 KiB |
Binary file not shown.
|
After Width: | Height: | Size: 168 KiB |
Binary file not shown.
|
After Width: | Height: | Size: 162 KiB |
Binary file not shown.
|
After Width: | Height: | Size: 152 KiB |
Binary file not shown.
|
After Width: | Height: | Size: 350 KiB |
@@ -135,6 +135,12 @@
|
||||
<tool file="regVariation/t_test_two_samples.xml" />
|
||||
<tool file="regVariation/compute_q_values.xml" />
|
||||
</section>
|
||||
<section name="Wavelet Analysis" id="dwt">
|
||||
<tool file="discreteWavelet/execute_dwt_IvC_all.xml" />
|
||||
<tool file="discreteWavelet/execute_dwt_cor_aVa_perClass.xml" />
|
||||
<tool file="discreteWavelet/execute_dwt_cor_aVb_all.xml" />
|
||||
<tool file="discreteWavelet/execute_dwt_var_perClass.xml" />
|
||||
</section>
|
||||
<section name="Graph/Display Data" id="plots">
|
||||
<tool file="plotting/histogram2.xml" />
|
||||
<tool file="plotting/scatterplot.xml" />
|
||||
|
||||
@@ -0,0 +1,210 @@
|
||||
#!/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);
|
||||
@@ -0,0 +1,112 @@
|
||||
<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">
|
||||
<description>between two datasets using Discrete Wavelet Transfoms</description>
|
||||
|
||||
<command interpreter="perl">
|
||||
execute_dwt_IvC_all.pl $inputFile1 $inputFile2 $outputFile1 $outputFile2
|
||||
</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"/>
|
||||
</inputs>
|
||||
|
||||
<outputs>
|
||||
<data format="tabular" name="outputFile1"/>
|
||||
<data format="pdf" name="outputFile2"/>
|
||||
</outputs>
|
||||
|
||||
<help>
|
||||
|
||||
.. class:: infomark
|
||||
|
||||
**What it does**
|
||||
|
||||
This program generates plots and computes table matrix of second moments, p-values, and test orientations at multiple scales for the correlation between the occurrences of features in one dataset and their occurrences in another using multiscale wavelet analysis technique.
|
||||
|
||||
The program assumes that the user has two sets of DNA sequences, S1 and S1, each of which consists of one or more sequences of equal length. Each sequence in each set is divided into the same number of multiple intervals n such that n = 2^k, where k is a positive integer and k >= 1. Thus, n could be any value of the set {2, 4, 8, 16, 32, 64, 128, ...}. k represents the number of scales.
|
||||
|
||||
The program has two input files obtained as follows:
|
||||
|
||||
For a given set of features, say motifs, the user counts the number of occurrences of each feature in each interval of each sequence in S1 and S1, and builds two tabular files representing the count results in each interval of S1 and S1. These are the input files of the program.
|
||||
|
||||
The program gives two output files:
|
||||
|
||||
- The first output file is a TABULAR format file representing the second moments, p-values, and test orientations for each feature at each scale.
|
||||
- The second output file is a PDF file consisting of as many figures as the number of features, such that each figure represents the values of the second moment for that feature at every scale.
|
||||
|
||||
-----
|
||||
|
||||
.. class:: warningmark
|
||||
|
||||
**Note**
|
||||
|
||||
In order to obtain empirical p-values, a random perumtation test is implemented by the program, which results in the fact that the program gives slightly different results each time it is run on the same input file.
|
||||
|
||||
-----
|
||||
|
||||
**Example**
|
||||
|
||||
Counting the occurrences of 5 features (motifs) in 16 intervals (one line per interval) of the DNA sequences in S1 gives the following tabular file::
|
||||
|
||||
deletionHoptspot insertionHoptspot dnaPolPauseFrameshift topoisomeraseCleavageSite translinTarget
|
||||
226 403 416 221 1165
|
||||
236 444 380 241 1223
|
||||
242 496 391 195 1116
|
||||
243 429 364 191 1118
|
||||
244 410 371 236 1063
|
||||
230 386 370 217 1087
|
||||
275 404 402 214 1044
|
||||
265 443 365 231 1086
|
||||
255 390 354 246 1114
|
||||
281 384 406 232 1102
|
||||
263 459 369 251 1135
|
||||
280 433 400 251 1159
|
||||
278 385 382 231 1147
|
||||
248 393 389 211 1162
|
||||
251 403 385 246 1114
|
||||
239 383 347 227 1172
|
||||
|
||||
And counting the occurrences of 5 features (motifs) in 16 intervals (one line per interval) of the DNA sequences in S2 gives the following tabular file::
|
||||
|
||||
deletionHoptspot insertionHoptspot dnaPolPauseFrameshift topoisomeraseCleavageSite translinTarget
|
||||
235 374 407 257 1159
|
||||
244 356 353 212 1128
|
||||
233 343 322 204 1110
|
||||
222 329 398 253 1054
|
||||
216 325 328 253 1129
|
||||
257 368 352 221 1115
|
||||
238 360 346 224 1102
|
||||
225 350 377 248 1107
|
||||
230 330 365 236 1132
|
||||
241 389 357 220 1120
|
||||
274 354 392 235 1120
|
||||
250 379 354 210 1102
|
||||
254 329 320 251 1080
|
||||
221 355 406 279 1127
|
||||
224 330 390 249 1129
|
||||
246 366 364 218 1176
|
||||
|
||||
|
||||
We notice that the number of scales here is 4 because 16 = 2^4. Runnig the program on the above input files gives the following output:
|
||||
|
||||
The first output file::
|
||||
|
||||
motif 1_moment2 1_pval 1_test 2_moment2 2_pval 2_test 3_moment2 3_pval 3_test 4_moment2 4_pval 4_test
|
||||
|
||||
deletionHoptspot 0.8751 0.376 L 1.549 0.168 R 0.6152 0.434 L 0.5735 0.488 R
|
||||
insertionHoptspot 0.902 0.396 L 1.172 0.332 R 0.6843 0.456 L 1.728 0.213 R
|
||||
dnaPolPauseFrameshift 1.65 0.013 R 0.267 0.055 L 0.1387 0.124 L 0.4516 0.498 L
|
||||
topoisomeraseCleavageSite 0.7443 0.233 L 1.023 0.432 R 1.933 0.155 R 1.09 0.3 R
|
||||
translinTarget 0.5084 0.057 L 0.8219 0.446 L 3.604 0.019 R 0.4377 0.492 L
|
||||
|
||||
The second output file:
|
||||
|
||||
.. image:: ../static/operation_icons/dwt_IvC_1.png
|
||||
.. image:: ../static/operation_icons/dwt_IvC_2.png
|
||||
.. image:: ../static/operation_icons/dwt_IvC_3.png
|
||||
.. image:: ../static/operation_icons/dwt_IvC_4.png
|
||||
.. image:: ../static/operation_icons/dwt_IvC_5.png
|
||||
|
||||
</help>
|
||||
|
||||
</tool>
|
||||
@@ -0,0 +1,221 @@
|
||||
#!/usr/bin/perl -w
|
||||
|
||||
use warnings;
|
||||
use IO::Handle;
|
||||
|
||||
$usage = "execute_dwt_cor_aVa_perClass.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 = "cor_aVa";
|
||||
|
||||
# construct an R script to implement the IvC test
|
||||
print "\n";
|
||||
|
||||
$r_script = "get_dwt_cor_aVa_test.r";
|
||||
print "$r_script \n";
|
||||
|
||||
open(Rcmd, ">", "$r_script") or die "Cannot open $r_script \n\n";
|
||||
print Rcmd "
|
||||
##################################################################################
|
||||
# code to do all correlation tests of form: motif(a) vs. motif(a)
|
||||
# add code to 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\", method = \"kendall\", 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, \"cor\", sep = \"_\"), paste(i, \"pval\", sep = \"_\"));
|
||||
}
|
||||
print(title);
|
||||
|
||||
# normalize the raw data
|
||||
data.short <- apply(data.short, 2, norm);
|
||||
data.long <- apply(data.long, 2, norm);
|
||||
|
||||
for(i in 1:length(names.short)){
|
||||
# Kendall Tau
|
||||
# DWT wavelet correlation function
|
||||
# include significance to compare
|
||||
wave1.dwt = wave2.dwt = NULL;
|
||||
tau.dwt = NULL;
|
||||
out = NULL;
|
||||
|
||||
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 {
|
||||
wave1.dwt <- dwt(data.short[, i], wf = wf, short.levels, boundary = boundary);
|
||||
wave2.dwt <- dwt(data.long[, i], wf = wf, short.levels, boundary = boundary);
|
||||
tau.dwt <- vector(length=short.levels)
|
||||
|
||||
#perform cor test on wavelet coefficients per scale
|
||||
for(level in 1:short.levels){
|
||||
w1_level = w2_level = NULL;
|
||||
w1_level <- (wave1.dwt[[level]]);
|
||||
w2_level <- (wave2.dwt[[level]]);
|
||||
tau.dwt[level] <- cor.test(w1_level, w2_level, method = method)\$estimate;
|
||||
}
|
||||
|
||||
# CI bands by permutation of time series
|
||||
feature1 = feature2 = NULL;
|
||||
feature1 = data.short[, i];
|
||||
feature2 = data.long[, i];
|
||||
null = results = med = NULL;
|
||||
cor_25 = cor_975 = NULL;
|
||||
|
||||
for (k in 1:1000) {
|
||||
nk_1 = nk_2 = NULL;
|
||||
null.levels = NULL;
|
||||
cor = NULL;
|
||||
null_wave1 = null_wave2 = 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;
|
||||
cor <- vector(length = null.levels);
|
||||
null_wave1 <- dwt(nk_1, wf = wf, short.levels, boundary = boundary);
|
||||
null_wave2 <- dwt(nk_2, wf = wf, short.levels, boundary = boundary);
|
||||
|
||||
for(level in 1:null.levels){
|
||||
null_level1 = null_level2 = NULL;
|
||||
null_level1 <- (null_wave1[[level]]);
|
||||
null_level2 <- (null_wave2[[level]]);
|
||||
cor[level] <- cor.test(null_level1, null_level2, method = method)\$estimate;
|
||||
}
|
||||
null = rbind(null, cor);
|
||||
}
|
||||
|
||||
null <- apply(null, 2, sort, na.last = TRUE);
|
||||
print(paste(\"NAs\", length(which(is.na(null))), sep = \" \"));
|
||||
cor_25 <- null[25,];
|
||||
cor_975 <- null[975,];
|
||||
med <- (apply(null, 2, median, na.rm = TRUE));
|
||||
|
||||
# plot
|
||||
results <- cbind(tau.dwt, cor_25, cor_975);
|
||||
matplot(results, type = \"b\", pch = \"*\" , lty = 1, col = c(1, 2, 2), ylim = c(-1, 1), xlab = \"Wavelet Scale\", ylab = \"Wavelet Correlation Kendall's Tau\", main = (paste(test, names.short[i], sep = \" \")), cex.main = 0.75);
|
||||
abline(h = 0);
|
||||
|
||||
# get pvalues by comparison to null distribution
|
||||
### modify pval calculation for error type II of T test ####
|
||||
out <- (names.short[i]);
|
||||
for (m in 1:length(tau.dwt)){
|
||||
print(paste(\"scale\", m, sep = \" \"));
|
||||
print(paste(\"tau\", tau.dwt[m], sep = \" \"));
|
||||
print(paste(\"med\", med[m], sep = \" \"));
|
||||
out <- c(out, format(tau.dwt[m], digits = 3));
|
||||
pv = NULL;
|
||||
if(is.na(tau.dwt[m])){
|
||||
pv <- \"NA\";
|
||||
}
|
||||
else {
|
||||
if (tau.dwt[m] >= med[m]){
|
||||
# R tail test
|
||||
print(paste(\"R\"));
|
||||
### per sv ok to use inequality not strict
|
||||
pv <- (length(which(null[, m] >= tau.dwt[m])))/(length(na.exclude(null[, m])));
|
||||
if (tau.dwt[m] == med[m]){
|
||||
print(\"tau == med\");
|
||||
print(summary(null[, m]));
|
||||
}
|
||||
}
|
||||
else if (tau.dwt[m] < med[m]){
|
||||
# L tail test
|
||||
print(paste(\"L\"));
|
||||
pv <- (length(which(null[, m] <= tau.dwt[m])))/(length(na.exclude(null[, m])));
|
||||
}
|
||||
}
|
||||
out <- c(out, pv);
|
||||
print(paste(\"pval\", pv, sep = \" \"));
|
||||
}
|
||||
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
|
||||
|
||||
inputData1 = inputData2 = NULL;
|
||||
inputData.short1 = inputData.short2 = NULL;
|
||||
inputDataNames.short1 = inputDataNames.short2 = NULL;
|
||||
|
||||
inputData1 <- read.delim(\"$firstInputFile\");
|
||||
inputData.short1 <- inputData1[, +c(1:ncol(inputData1))];
|
||||
inputDataNames.short1 <- colnames(inputData.short1);
|
||||
|
||||
inputData2 <- read.delim(\"$secondInputFile\");
|
||||
inputData.short2 <- inputData2[, +c(1:ncol(inputData2))];
|
||||
inputDataNames.short2 <- colnames(inputData.short2);
|
||||
|
||||
# cor test for motif(a) in inputData1 vs motif(a) in inputData2
|
||||
dwt_cor(inputData.short1, inputDataNames.short1, inputData.short2, inputDataNames.short2, test = \"$test\", pdf = \"$secondOutputFile\", table = \"$firstOutputFile\");
|
||||
print (\"done with the correlation test\");
|
||||
|
||||
#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);
|
||||
|
||||
@@ -0,0 +1,112 @@
|
||||
<tool id="compute_p-values_correlation_coefficients_feature_occurrences_between_two_datasets_using_discrete_wavelet_transfom" name="Compute P-values and Correlation Coefficients for Feature Occurrences" version="1.0.0">
|
||||
<description>between two datasets using Discrete Wavelet Transfoms</description>
|
||||
|
||||
<command interpreter="perl">
|
||||
execute_dwt_cor_aVa_perClass.pl $inputFile1 $inputFile2 $outputFile1 $outputFile2
|
||||
</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"/>
|
||||
</inputs>
|
||||
|
||||
<outputs>
|
||||
<data format="tabular" name="outputFile1"/>
|
||||
<data format="pdf" name="outputFile2"/>
|
||||
</outputs>
|
||||
|
||||
<help>
|
||||
|
||||
.. class:: infomark
|
||||
|
||||
**What it does**
|
||||
|
||||
This program generates plots and computes table matrix of coefficient correlations and p-values at multiple scales for the correlation between the occurrences of features in one dataset and their occurrences in another using multiscale wavelet analysis technique.
|
||||
|
||||
The program assumes that the user has two sets of DNA sequences, S1 and S1, each of which consists of one or more sequences of equal length. Each sequence in each set is divided into the same number of multiple intervals n such that n = 2^k, where k is a positive integer and k >= 1. Thus, n could be any value of the set {2, 4, 8, 16, 32, 64, 128, ...}. k represents the number of scales.
|
||||
|
||||
The program has two input files obtained as follows:
|
||||
|
||||
For a given set of features, say motifs, the user counts the number of occurrences of each feature in each interval of each sequence in S1 and S1, and builds two tabular files representing the count results in each interval of S1 and S1. These are the input files of the program.
|
||||
|
||||
The program gives two output files:
|
||||
|
||||
- The first output file is a TABULAR format file representing the coefficient correlations and p-values for each feature at each scale.
|
||||
- The second output file is a PDF file consisting of as many figures as the number of features, such that each figure represents the values of the coefficient correlation for that feature at every scale.
|
||||
|
||||
-----
|
||||
|
||||
.. class:: warningmark
|
||||
|
||||
**Note**
|
||||
|
||||
In order to obtain empirical p-values, a random perumtation test is implemented by the program, which results in the fact that the program gives slightly different results each time it is run on the same input file.
|
||||
|
||||
-----
|
||||
|
||||
**Example**
|
||||
|
||||
Counting the occurrences of 5 features (motifs) in 16 intervals (one line per interval) of the DNA sequences in S1 gives the following tabular file::
|
||||
|
||||
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
|
||||
|
||||
And counting the occurrences of 5 features (motifs) in 16 intervals (one line per interval) of the DNA sequences in S2 gives the following tabular file::
|
||||
|
||||
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
|
||||
|
||||
|
||||
We notice that the number of scales here is 4 because 16 = 2^4. Running the program on the above input files gives the following output:
|
||||
|
||||
The first output file::
|
||||
|
||||
motif 1_cor 1_pval 2_cor 2_pval 3_cor 3_pval 4_cor 4_pval
|
||||
|
||||
deletionHoptspot 0.4 0.072 0.143 0.394 -0.667 0.244 1 0.491
|
||||
insertionHoptspot 0.343 0.082 -0.0714 0.446 -1 0.12 1 0.502
|
||||
dnaPolPauseFrameshift 0.617 0.004 -0.5 0.13 0.667 0.234 1 0.506
|
||||
topoisomeraseCleavageSite -0.183 0.242 -0.286 0.256 0.333 0.353 -1 0.489
|
||||
translinTarget 0.0167 0.503 -0.0714 0.469 1 0.136 1 0.485
|
||||
|
||||
The second output file:
|
||||
|
||||
.. image:: ../static/operation_icons/dwt_cor_aVa_1.png
|
||||
.. image:: ../static/operation_icons/dwt_cor_aVa_2.png
|
||||
.. image:: ../static/operation_icons/dwt_cor_aVa_3.png
|
||||
.. image:: ../static/operation_icons/dwt_cor_aVa_4.png
|
||||
.. image:: ../static/operation_icons/dwt_cor_aVa_5.png
|
||||
|
||||
</help>
|
||||
|
||||
</tool>
|
||||
@@ -0,0 +1,223 @@
|
||||
#!/usr/bin/perl -w
|
||||
|
||||
use warnings;
|
||||
use IO::Handle;
|
||||
|
||||
$usage = "execute_dwt_cor_aVb_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 = "cor_aVb_all";
|
||||
|
||||
# construct an R script to implement the IvC test
|
||||
print "\n";
|
||||
|
||||
$r_script = "get_dwt_cor_aVa_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 all correlation tests of form: motif(a) vs. motif(b)
|
||||
# add code to 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\", method = \"kendall\", 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(\"motif1\", \"motif2\");
|
||||
for (i in 1:short.levels){
|
||||
title <- c(title, paste(i, \"cor\", sep = \"_\"), paste(i, \"pval\", sep = \"_\"));
|
||||
}
|
||||
print(title);
|
||||
|
||||
# normalize the raw data
|
||||
data.short <- apply(data.short, 2, norm);
|
||||
data.long <- apply(data.long, 2, norm);
|
||||
|
||||
# loop to compare a vs b
|
||||
for(i in 1:length(names.short)){
|
||||
for(j in 1:length(names.long)){
|
||||
if(i >= j){
|
||||
next;
|
||||
}
|
||||
else {
|
||||
# Kendall Tau
|
||||
# DWT wavelet correlation function
|
||||
# include significance to compare
|
||||
wave1.dwt = wave2.dwt = NULL;
|
||||
tau.dwt = NULL;
|
||||
out = NULL;
|
||||
|
||||
print(names.short[i]);
|
||||
print(names.long[j]);
|
||||
|
||||
# need exit if not comparing motif(a) vs motif(a)
|
||||
if (names.short[i] == names.long[j]){
|
||||
stop(paste(\"motif\", names.short[i], \"is the same as\", names.long[j], sep = \" \"));
|
||||
}
|
||||
else {
|
||||
wave1.dwt <- dwt(data.short[, i], wf = wf, short.levels, boundary = boundary);
|
||||
wave2.dwt <- dwt(data.long[, j], wf = wf, short.levels, boundary = boundary);
|
||||
tau.dwt <-vector(length = short.levels)
|
||||
|
||||
# perform cor test on wavelet coefficients per scale
|
||||
for(level in 1:short.levels){
|
||||
w1_level = w2_level = NULL;
|
||||
w1_level <- (wave1.dwt[[level]]);
|
||||
w2_level <- (wave2.dwt[[level]]);
|
||||
tau.dwt[level] <- cor.test(w1_level, w2_level, method = method)\$estimate;
|
||||
}
|
||||
|
||||
# CI bands by permutation of time series
|
||||
feature1 = feature2 = NULL;
|
||||
feature1 = data.short[, i];
|
||||
feature2 = data.long[, j];
|
||||
null = results = med = NULL;
|
||||
cor_25 = cor_975 = NULL;
|
||||
|
||||
for (k in 1:1000) {
|
||||
nk_1 = nk_2 = NULL;
|
||||
null.levels = NULL;
|
||||
cor = NULL;
|
||||
null_wave1 = null_wave2 = 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;
|
||||
cor <- vector(length = null.levels);
|
||||
null_wave1 <- dwt(nk_1, wf = wf, short.levels, boundary = boundary);
|
||||
null_wave2 <- dwt(nk_2, wf = wf, short.levels, boundary = boundary);
|
||||
|
||||
for(level in 1:null.levels){
|
||||
null_level1 = null_level2 = NULL;
|
||||
null_level1 <- (null_wave1[[level]]);
|
||||
null_level2 <- (null_wave2[[level]]);
|
||||
cor[level] <- cor.test(null_level1, null_level2, method = method)\$estimate;
|
||||
}
|
||||
null = rbind(null, cor);
|
||||
}
|
||||
|
||||
null <- apply(null, 2, sort, na.last = TRUE);
|
||||
cor_25 <- null[25, ];
|
||||
cor_975 <- null[975, ];
|
||||
med <- (apply(null, 2, median, na.rm = TRUE));
|
||||
|
||||
# plot
|
||||
results <- cbind(tau.dwt, cor_25, cor_975);
|
||||
matplot(results, type = \"b\", pch = \"*\", lty = 1, col = c(1, 2, 2), ylim = c(-1, 1), xlab = \"Wavelet Scale\", ylab = \"Wavelet Correlation Kendall's Tau\", main = (paste(test, names.short[i], \"vs.\", names.long[j], sep = \" \")), cex.main = 0.75);
|
||||
abline(h = 0);
|
||||
|
||||
# get pvalues by comparison to null distribution
|
||||
### modify pval calculation for error type II of T test ####
|
||||
out <- c(names.short[i],names.long[j]);
|
||||
for (m in 1:length(tau.dwt)){
|
||||
print(m);
|
||||
print(tau.dwt[m]);
|
||||
out <- c(out, format(tau.dwt[m], digits = 3));
|
||||
pv = NULL;
|
||||
if(is.na(tau.dwt[m])){
|
||||
pv <- \"NA\";
|
||||
}
|
||||
else{
|
||||
if (tau.dwt[m] >= med[m]){
|
||||
# R tail test
|
||||
pv <- (length(which(null[, m] >= tau.dwt[m])))/(length(na.exclude(null[, m])));
|
||||
}
|
||||
else{
|
||||
if (tau.dwt[m] < med[m]){
|
||||
# L tail test
|
||||
pv <- (length(which(null[, m] <= tau.dwt[m])))/(length(na.exclude(null[, m])));
|
||||
}
|
||||
}
|
||||
}
|
||||
out <- c(out, pv);
|
||||
print(pv);
|
||||
}
|
||||
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
|
||||
|
||||
inputData1 = inputData2 = NULL;
|
||||
inputData.short1 = inputData.short2 = NULL;
|
||||
inputDataNames.short1 = inputDataNames.short2 = NULL;
|
||||
|
||||
inputData1 <- read.delim(\"$firstInputFile\");
|
||||
inputData.short1 <- inputData1[, +c(1:ncol(inputData1))];
|
||||
inputDataNames.short1 <- colnames(inputData.short1);
|
||||
|
||||
inputData2 <- read.delim(\"$secondInputFile\");
|
||||
inputData.short2 <- inputData2[, +c(1:ncol(inputData2))];
|
||||
inputDataNames.short2 <- colnames(inputData.short2);
|
||||
|
||||
# cor test for motif(a) in inputData1 vs motif(b) in inputData2
|
||||
dwt_cor(inputData.short1, inputDataNames.short1, inputData.short2, inputDataNames.short2, test = \"$test\", pdf = \"$secondOutputFile\", table = \"$firstOutputFile\");
|
||||
print (\"done with the correlation test\");
|
||||
|
||||
#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);
|
||||
@@ -0,0 +1,123 @@
|
||||
<tool id="compute_p-values_correlation_coefficients_featureA_featureB_occurrences_between_two_datasets_using_discrete_wavelet_transfom" name="Compute P-values and Correlation Coefficients for Occurrences of Two Set of Features" version="1.0.0">
|
||||
<description>between two datasets using Discrete Wavelet Transfoms</description>
|
||||
|
||||
<command interpreter="perl">
|
||||
execute_dwt_cor_aVb_all.pl $inputFile1 $inputFile2 $outputFile1 $outputFile2
|
||||
</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"/>
|
||||
</inputs>
|
||||
|
||||
<outputs>
|
||||
<data format="tabular" name="outputFile1"/>
|
||||
<data format="pdf" name="outputFile2"/>
|
||||
</outputs>
|
||||
|
||||
<help>
|
||||
|
||||
.. class:: infomark
|
||||
|
||||
**What it does**
|
||||
|
||||
This program generates plots and computes table matrix of coefficient correlations and p-values at multiple scales for the correlation between the occurrences of features in one dataset and their occurrences in another using multiscale wavelet analysis technique.
|
||||
|
||||
The program assumes that the user has two sets of DNA sequences, S1 and S1, each of which consists of one or more sequences of equal length. Each sequence in each set is divided into the same number of multiple intervals n such that n = 2^k, where k is a positive integer and k >= 1. Thus, n could be any value of the set {2, 4, 8, 16, 32, 64, 128, ...}. k represents the number of scales.
|
||||
|
||||
The program has two input files obtained as follows:
|
||||
|
||||
For a given set of features, say motifs, the user counts the number of occurrences of each feature in each interval of each sequence in S1 and S1, and builds two tabular files representing the count results in each interval of S1 and S1. These are the input files of the program.
|
||||
|
||||
The program gives two output files:
|
||||
|
||||
- The first output file is a TABULAR format file representing the coefficient correlations and p-values for each feature at each scale.
|
||||
- The second output file is a PDF file consisting of as many figures as the number of features, such that each figure represents the values of the coefficient correlations for that feature at every scale.
|
||||
|
||||
-----
|
||||
|
||||
.. class:: warningmark
|
||||
|
||||
**Note**
|
||||
|
||||
In order to obtain empirical p-values, a random perumtation test is implemented by the program, which results in the fact that the program gives slightly different results each time it is run on the same input file.
|
||||
|
||||
-----
|
||||
|
||||
**Example**
|
||||
|
||||
Counting the occurrences of 5 features (motifs) in 16 intervals (one line per interval) of the DNA sequences in S1 gives the following tabular file::
|
||||
|
||||
deletionHoptspot insertionHoptspot dnaPolPauseFrameshift topoisomeraseCleavageSite translinTarget
|
||||
82 162 158 79 459
|
||||
111 196 154 75 459
|
||||
98 178 160 79 475
|
||||
113 201 170 113 436
|
||||
113 173 147 95 446
|
||||
107 150 155 84 436
|
||||
106 166 175 96 448
|
||||
113 176 135 106 514
|
||||
113 170 152 87 450
|
||||
95 152 167 93 467
|
||||
91 171 169 118 426
|
||||
84 139 160 100 459
|
||||
92 154 164 104 440
|
||||
100 145 154 98 472
|
||||
91 161 152 71 461
|
||||
117 164 139 97 463
|
||||
|
||||
And counting the occurrences of 5 features (motifs) in 16 intervals (one line per interval) of the DNA sequences in S2 gives the following tabular file::
|
||||
|
||||
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
|
||||
|
||||
|
||||
We notice that the number of scales here is 4 because 16 = 2^4. Running the program on the above input files gives the following output:
|
||||
|
||||
The first output file::
|
||||
|
||||
motif1 motif2 1_cor 1_pval 2_cor 2_pval 3_cor 3_pval 4_cor 4_pval
|
||||
|
||||
deletionHoptspot insertionHoptspot -0.1 0.346 -0.214 0.338 1 0.127 1 0.467
|
||||
deletionHoptspot dnaPolPauseFrameshift 0.167 0.267 -0.214 0.334 1 0.122 1 0.511
|
||||
deletionHoptspot topoisomeraseCleavageSite 0.167 0.277 0.143 0.412 -0.667 0.243 1 0.521
|
||||
deletionHoptspot translinTarget 0 0.505 0.0714 0.441 1 0.124 1 0.518
|
||||
insertionHoptspot dnaPolPauseFrameshift -0.202 0.238 0.143 0.379 -1 0.122 1 0.517
|
||||
insertionHoptspot topoisomeraseCleavageSite -0.0336 0.457 0.214 0.29 0.667 0.252 1 0.503
|
||||
insertionHoptspot translinTarget 0.0672 0.389 0.429 0.186 -1 0.119 1 0.506
|
||||
dnaPolPauseFrameshift topoisomeraseCleavageSite -0.353 0.101 0.357 0.228 0 0.612 -1 0.49
|
||||
dnaPolPauseFrameshift translinTarget -0.151 0.303 -0.571 0.09 -0.333 0.37 -1 1
|
||||
topoisomeraseCleavageSite translinTarget -0.37 0.077 -0.222 0.297 0.667 0.234 -1 0.471
|
||||
|
||||
The second output file:
|
||||
|
||||
.. image:: ../static/operation_icons/dwt_cor_aVb_all_1.png
|
||||
.. image:: ../static/operation_icons/dwt_cor_aVb_all_2.png
|
||||
.. image:: ../static/operation_icons/dwt_cor_aVb_all_3.png
|
||||
.. image:: ../static/operation_icons/dwt_cor_aVb_all_4.png
|
||||
.. image:: ../static/operation_icons/dwt_cor_aVb_all_5.png
|
||||
.. image:: ../static/operation_icons/dwt_cor_aVb_all_6.png
|
||||
.. image:: ../static/operation_icons/dwt_cor_aVb_all_7.png
|
||||
.. image:: ../static/operation_icons/dwt_cor_aVb_all_8.png
|
||||
.. image:: ../static/operation_icons/dwt_cor_aVb_all_9.png
|
||||
.. image:: ../static/operation_icons/dwt_cor_aVb_all_10.png
|
||||
|
||||
|
||||
</help>
|
||||
|
||||
</tool>
|
||||
@@ -0,0 +1,320 @@
|
||||
#!/usr/bin/perl -w
|
||||
|
||||
use warnings;
|
||||
use IO::Handle;
|
||||
use POSIX qw(floor ceil);
|
||||
|
||||
# example: perl execute_dwt_var_perClass.pl hg18_NCNR_10bp_3flanks_deletionHotspot_data_del.txt deletionHotspot 3flanks del
|
||||
|
||||
$usage = "execute_dwt_var_perClass.pl [TABULAR.in] [TABULAR.out] [TABULAR.out] [PDF.out] \n";
|
||||
die $usage unless @ARGV == 4;
|
||||
|
||||
#get the input arguments
|
||||
my $inputFile = $ARGV[0];
|
||||
my $firstOutputFile = $ARGV[1];
|
||||
my $secondOutputFile = $ARGV[2];
|
||||
my $thirdOutputFile = $ARGV[3];
|
||||
|
||||
open (INPUT, "<", $inputFile) || die("Could not open file $inputFile \n");
|
||||
open (OUTPUT1, ">", $firstOutputFile) || die("Could not open file $firstOutputFile \n");
|
||||
open (OUTPUT2, ">", $secondOutputFile) || die("Could not open file $secondOutputFile \n");
|
||||
open (OUTPUT3, ">", $thirdOutputFile) || die("Could not open file $thirdOutputFile \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");
|
||||
|
||||
# choosing meaningful names for the output files
|
||||
$max_dwt = $firstOutputFile;
|
||||
$pvalue = $secondOutputFile;
|
||||
$pdf = $thirdOutputFile;
|
||||
|
||||
# count the number of columns in the input file
|
||||
while($buffer = <INPUT>){
|
||||
#if ($buffer =~ m/interval/){
|
||||
chomp($buffer);
|
||||
$buffer =~ s/^#\s*//;
|
||||
@contrl = split(/\t/, $buffer);
|
||||
last;
|
||||
#}
|
||||
}
|
||||
print "The number of columns in the input file is: " . (@contrl) . "\n";
|
||||
print "\n";
|
||||
|
||||
# count the number of motifs in the input file
|
||||
$count = 0;
|
||||
for ($i = 0; $i < @contrl; $i++){
|
||||
$count++;
|
||||
print "# $contrl[$i]\n";
|
||||
}
|
||||
print "The number of motifs in the input file is: $count \n";
|
||||
|
||||
# check if the number of motifs is not a multiple of 12, and round up is so
|
||||
$count2 = ($count/12);
|
||||
if ($count2 =~ m/(\D)/){
|
||||
print "the number of motifs is not a multiple of 12 \n";
|
||||
$count2 = ceil($count2);
|
||||
}
|
||||
else {
|
||||
print "the number of motifs is a multiple of 12 \n";
|
||||
}
|
||||
print "There will be $count2 subfiles\n\n";
|
||||
|
||||
# split infile into subfiles only 12 motif per file for R plotting
|
||||
for ($x = 1; $x <= $count2; $x++){
|
||||
$a = (($x - 1) * 12 + 1);
|
||||
$b = $x * 12;
|
||||
|
||||
if ($x < $count2){
|
||||
print "# data.short $x <- data_test[, +c($a:$b)]; \n";
|
||||
}
|
||||
else{
|
||||
print "# data.short $x <- data_test[, +c($a:ncol(data_test)]; \n";
|
||||
}
|
||||
}
|
||||
|
||||
print "\n";
|
||||
print "There are 4 output files: \n";
|
||||
print "The first output file is a pdf file\n";
|
||||
print "The second output file is a max_dwt file\n";
|
||||
print "The third output file is a pvalues file\n";
|
||||
print "The fourth output file is a test_final_pvalues file\n";
|
||||
|
||||
# write R script
|
||||
$r_script = "get_dwt_varPermut_getMax.r";
|
||||
print "The R file name is: $r_script \n";
|
||||
|
||||
open(Rcmd, ">", "$r_script") or die "Cannot open $r_script \n\n";
|
||||
|
||||
print Rcmd "
|
||||
######################################################################
|
||||
# plot power spectra, i.e. wavelet variance by class
|
||||
# add code to create null bands by permuting the original data series
|
||||
# get class of maximum significant variance per feature
|
||||
# generate plots and table matrix of variance 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_var_permut_getMax <- function(data, names, filter = 4, bc = \"symmetric\", method = \"kendall\", wf = \"haar\", boundary = \"reflection\") {
|
||||
max_var = NULL;
|
||||
matrix = NULL;
|
||||
title = NULL;
|
||||
final_pvalue = NULL;
|
||||
short.levels = NULL;
|
||||
scale = NULL;
|
||||
|
||||
print(names);
|
||||
|
||||
par(mfcol = c(length(names), length(names)), mar = c(0, 0, 0, 0), oma = c(4, 3, 3, 2), xaxt = \"s\", cex = 1, las = 1);
|
||||
|
||||
short.levels <- wd(data[, 1], filter.number = filter, bc = bc)\$nlevels;
|
||||
|
||||
title <- c(\"motif\");
|
||||
for (i in 1:short.levels){
|
||||
title <- c(title, paste(i, \"var\", sep = \"_\"), paste(i, \"pval\", sep = \"_\"), paste(i, \"test\", sep = \"_\"));
|
||||
}
|
||||
print(title);
|
||||
|
||||
# normalize the raw data
|
||||
data<-apply(data,2,norm);
|
||||
|
||||
for(i in 1:length(names)){
|
||||
for(j in 1:length(names)){
|
||||
temp = NULL;
|
||||
results = NULL;
|
||||
wave1.dwt = NULL;
|
||||
out = NULL;
|
||||
|
||||
out <- vector(length = length(title));
|
||||
temp <- vector(length = short.levels);
|
||||
|
||||
if(i < j) {
|
||||
plot(temp, type = \"n\", axes = FALSE, xlab = NA, ylab = NA);
|
||||
box(col = \"grey\");
|
||||
grid(ny = 0, nx = NULL);
|
||||
} else {
|
||||
if (i > j){
|
||||
plot(temp, type = \"n\", axes = FALSE, xlab = NA, ylab = NA);
|
||||
box(col = \"grey\");
|
||||
grid(ny = 0, nx = NULL);
|
||||
} else {
|
||||
|
||||
wave1.dwt <- dwt(data[, i], wf = wf, short.levels, boundary = boundary);
|
||||
|
||||
temp_row = (short.levels + 1 ) * -1;
|
||||
temp_col = 1;
|
||||
temp <- wave.variance(wave1.dwt)[temp_row, temp_col];
|
||||
|
||||
#permutations code :
|
||||
feature1 = NULL;
|
||||
null = NULL;
|
||||
var_25 = NULL;
|
||||
var_975 = NULL;
|
||||
med = NULL;
|
||||
|
||||
feature1 = data[, i];
|
||||
for (k in 1:1000) {
|
||||
nk_1 = NULL;
|
||||
null.levels = NULL;
|
||||
var = NULL;
|
||||
null_wave1 = NULL;
|
||||
|
||||
nk_1 = sample(feature1, length(feature1), replace = FALSE);
|
||||
null.levels <- wd(nk_1, filter.number = filter, bc = bc)\$nlevels;
|
||||
var <- vector(length = length(null.levels));
|
||||
null_wave1 <- dwt(nk_1, wf = wf, short.levels, boundary = boundary);
|
||||
var<- wave.variance(null_wave1)[-8, 1];
|
||||
null= rbind(null, var);
|
||||
}
|
||||
null <- apply(null, 2, sort, na.last = TRUE);
|
||||
var_25 <- null[25, ];
|
||||
var_975 <- null[975, ];
|
||||
med <- (apply(null, 2, median, na.rm = TRUE));
|
||||
|
||||
# plot
|
||||
results <- cbind(temp, var_25, var_975);
|
||||
matplot(results, type = \"b\", pch = \"*\", lty = 1, col = c(1, 2, 2), axes = F);
|
||||
|
||||
# get pvalues by comparison to null distribution
|
||||
out <- (names[i]);
|
||||
for (m in 1:length(temp)){
|
||||
print(paste(\"scale\", m, sep = \" \"));
|
||||
print(paste(\"var\", temp[m], sep = \" \"));
|
||||
print(paste(\"med\", med[m], sep = \" \"));
|
||||
pv = tail = NULL;
|
||||
out <- c(out, format(temp[m], digits = 3));
|
||||
if (temp[m] >= med[m]){
|
||||
# R tail test
|
||||
print(\"R\");
|
||||
tail <- \"R\";
|
||||
pv <- (length(which(null[, m] >= temp[m])))/(length(na.exclude(null[, m])));
|
||||
|
||||
} else {
|
||||
if (temp[m] < med[m]){
|
||||
# L tail test
|
||||
print(\"L\");
|
||||
tail <- \"L\";
|
||||
pv <- (length(which(null[, m] <= temp[m])))/(length(na.exclude(null[, m])));
|
||||
}
|
||||
}
|
||||
out <- c(out, pv);
|
||||
print(pv);
|
||||
out <- c(out, tail);
|
||||
}
|
||||
final_pvalue <-rbind(final_pvalue, out);
|
||||
|
||||
|
||||
# get variances outside null bands by comparing temp to null
|
||||
## temp stores variance for each scale, and null stores permuted variances for null bands
|
||||
for (n in 1:length(temp)){
|
||||
if (temp[n] <= var_975[n]){
|
||||
temp[n] <- NA;
|
||||
} else {
|
||||
temp[n] <- temp[n];
|
||||
}
|
||||
}
|
||||
matrix <- rbind(matrix, temp)
|
||||
}
|
||||
}
|
||||
# labels
|
||||
if (i == 1){
|
||||
mtext(names[j], side = 2, line = 0.5, las = 3, cex = 0.25);
|
||||
}
|
||||
if (j == 1){
|
||||
mtext(names[i], side = 3, line = 0.5, cex = 0.25);
|
||||
}
|
||||
if (j == length(names)){
|
||||
axis(1, at = (1:short.levels), las = 3, cex.axis = 0.5);
|
||||
}
|
||||
}
|
||||
}
|
||||
colnames(final_pvalue) <- title;
|
||||
#write.table(final_pvalue, file = \"test_final_pvalue.txt\", sep = \"\\t\", quote = FALSE, row.names = FALSE, append = TRUE);
|
||||
|
||||
# get maximum variance larger than expectation by comparison to null bands
|
||||
varnames <- vector();
|
||||
for(i in 1:length(names)){
|
||||
name1 = paste(names[i], \"var\", sep = \"_\")
|
||||
varnames <- c(varnames, name1)
|
||||
}
|
||||
rownames(matrix) <- varnames;
|
||||
colnames(matrix) <- (1:short.levels);
|
||||
max_var <- names;
|
||||
scale <- vector(length = length(names));
|
||||
for (x in 1:nrow(matrix)){
|
||||
if (length(which.max(matrix[x, ])) == 0){
|
||||
scale[x] <- NA;
|
||||
}
|
||||
else{
|
||||
scale[x] <- colnames(matrix)[which.max(matrix[x, ])];
|
||||
}
|
||||
}
|
||||
max_var <- cbind(max_var, scale);
|
||||
write.table(max_var, file = \"$max_dwt\", sep = \"\\t\", quote = FALSE, row.names = FALSE, append = TRUE);
|
||||
return(final_pvalue);
|
||||
}\n";
|
||||
|
||||
print Rcmd "
|
||||
# execute
|
||||
# read in data
|
||||
|
||||
data_test = NULL;
|
||||
data_test <- read.delim(\"$inputFile\");
|
||||
|
||||
pdf(file = \"$pdf\", width = 11, height = 8);
|
||||
|
||||
# loop to read and execute on all $count2 subfiles
|
||||
final = NULL;
|
||||
for (x in 1:$count2){
|
||||
sub = NULL;
|
||||
sub_names = NULL;
|
||||
a = NULL;
|
||||
b = NULL;
|
||||
|
||||
a = ((x - 1) * 12 + 1);
|
||||
b = x * 12;
|
||||
|
||||
if (x < $count2){
|
||||
sub <- data_test[, +c(a:b)];
|
||||
sub_names <- colnames(data_test)[a:b];
|
||||
final <- rbind(final, dwt_var_permut_getMax(sub, sub_names));
|
||||
}
|
||||
else{
|
||||
sub <- data_test[, +c(a:ncol(data_test))];
|
||||
sub_names <- colnames(data_test)[a:ncol(data_test)];
|
||||
final <- rbind(final, dwt_var_permut_getMax(sub, sub_names));
|
||||
|
||||
}
|
||||
}
|
||||
|
||||
dev.off();
|
||||
|
||||
write.table(final, file = \"$pvalue\", sep = \"\\t\", quote = FALSE, row.names = FALSE);
|
||||
|
||||
#eof\n";
|
||||
|
||||
close Rcmd;
|
||||
|
||||
system("echo \"wavelet ANOVA started on \`hostname\` at \`date\`\"\n");
|
||||
system("R --no-restore --no-save --no-readline < $r_script > $r_script.out");
|
||||
system("echo \"wavelet ANOVA ended on \`hostname\` at \`date\`\"\n");
|
||||
|
||||
#close the input and output and error files
|
||||
close(ERROR);
|
||||
close(OUTPUT3);
|
||||
close(OUTPUT2);
|
||||
close(OUTPUT1);
|
||||
close(INPUT);
|
||||
@@ -0,0 +1,105 @@
|
||||
<tool id="compute_p-values_max_variances_feature_occurrences_in_one_dataset_using_discrete_wavelet_transfom" name="Compute P-values and Max Variances for Feature Occurrences" version="1.0.0">
|
||||
<description>in one dataset using Discrete Wavelet Transfoms</description>
|
||||
|
||||
<command interpreter="perl">
|
||||
execute_dwt_var_perClass.pl $inputFile $outputFile1 $outputFile2 $outputFile3
|
||||
</command>
|
||||
|
||||
<inputs>
|
||||
<param format="tabular" name="inputFile" type="data" label="Select the input file"/>
|
||||
</inputs>
|
||||
|
||||
<outputs>
|
||||
<data format="tabular" name="outputFile1"/>
|
||||
<data format="tabular" name="outputFile2"/>
|
||||
<data format="pdf" name="outputFile3"/>
|
||||
</outputs>
|
||||
|
||||
<help>
|
||||
|
||||
.. class:: infomark
|
||||
|
||||
**What it does**
|
||||
|
||||
This program generates plots and computes table matrix of maximum variances, p-values, and test orientations at multiple scales for the occurrences of a class of features in one dataset of DNA sequences using multiscale wavelet analysis technique.
|
||||
|
||||
The program assumes that the user has one set of DNA sequences, S, which consists of one or more sequences of equal length. Each sequence in S is divided into the same number of multiple intervals n such that n = 2^k, where k is a positive integer and k >= 1. Thus, n could be any value of the set {2, 4, 8, 16, 32, 64, 128, ...}. k represents the number of scales.
|
||||
|
||||
The program has one input file obtained as follows:
|
||||
|
||||
For a given set of features, say motifs, the user counts the number of occurrences of each feature in each interval of each sequence in S, and builds a tabular file representing the count results in each interval of S. This is the input file of the program.
|
||||
|
||||
The program gives three output files:
|
||||
|
||||
- The first output file is a TABULAR format file giving the scales at which each features has a maximum variances.
|
||||
- The second output file is a TABULAR format file representing the variances, p-values, and test orientation for the occurrences of features at each scale based on a random permutation test and using multiscale wavelet analysis technique.
|
||||
- The third output file is a PDF file plotting the wavelet variances of each feature at each scale.
|
||||
|
||||
-----
|
||||
|
||||
.. class:: warningmark
|
||||
|
||||
**Note**
|
||||
|
||||
- If the number of features is greater than 12, the program will divide each output file into subfiles, such that each subfile represents the results of a group of 12 features except the last subfile that will represents the results of the rest. For example, if the number of features is 17, the p-values file will consists of two subfiles, the first for the features 1-12 and the second for the features 13-17. As for the PDF file, it will consists of two pages in this case.
|
||||
- In order to obtain empirical p-values, a random perumtation test is implemented by the program, which results in the fact that the program gives slightly different results each time it is run on the same input file.
|
||||
|
||||
-----
|
||||
|
||||
|
||||
**Example**
|
||||
|
||||
Counting the occurrences of 8 features (motifs) in 16 intervals (one line per interval) of set of DNA sequences in S gives the following tabular file::
|
||||
|
||||
deletionHoptspot insertionHoptspot dnaPolPauseFrameshift indelHotspot topoisomeraseCleavageSite translinTarget vDjRecombinationSignal x-likeSite
|
||||
226 403 416 221 1165 832 749 1056
|
||||
236 444 380 241 1223 746 782 1207
|
||||
242 496 391 195 1116 643 770 1219
|
||||
243 429 364 191 1118 694 783 1223
|
||||
244 410 371 236 1063 692 805 1233
|
||||
230 386 370 217 1087 657 787 1215
|
||||
275 404 402 214 1044 697 831 1188
|
||||
265 443 365 231 1086 694 782 1184
|
||||
255 390 354 246 1114 642 773 1176
|
||||
281 384 406 232 1102 719 787 1191
|
||||
263 459 369 251 1135 643 810 1215
|
||||
280 433 400 251 1159 701 777 1151
|
||||
278 385 382 231 1147 697 707 1161
|
||||
248 393 389 211 1162 723 759 1183
|
||||
251 403 385 246 1114 752 776 1153
|
||||
239 383 347 227 1172 759 789 1141
|
||||
|
||||
We notice that the number of scales here is 4 because 16 = 2^4. Runnig the program on the above input file gives the following 3 output files:
|
||||
|
||||
The first output file::
|
||||
|
||||
motifs max_var at scale
|
||||
deletionHoptspot NA
|
||||
insertionHoptspot NA
|
||||
dnaPolPauseFrameshift NA
|
||||
indelHotspot NA
|
||||
topoisomeraseCleavageSite 3
|
||||
translinTarget NA
|
||||
vDjRecombinationSignal NA
|
||||
x.likeSite NA
|
||||
|
||||
The second output file::
|
||||
|
||||
motif 1_var 1_pval 1_test 2_var 2_pval 2_test 3_var 3_pval 3_test 4_var 4_pval 4_test
|
||||
|
||||
deletionHoptspot 0.457 0.048 L 1.18 0.334 R 1.61 0.194 R 3.41 0.055 R
|
||||
insertionHoptspot 0.556 0.109 L 1.34 0.272 R 1.59 0.223 R 2.02 0.157 R
|
||||
dnaPolPauseFrameshift 1.42 0.089 R 0.66 0.331 L 0.421 0.305 L 0.121 0.268 L
|
||||
indelHotspot 0.373 0.021 L 1.36 0.254 R 1.24 0.301 R 4.09 0.047 R
|
||||
topoisomeraseCleavageSite 0.305 0.002 L 0.936 0.489 R 3.78 0.01 R 1.25 0.272 R
|
||||
translinTarget 0.525 0.061 L 1.69 0.11 R 2.02 0.131 R 0.00891 0.069 L
|
||||
vDjRecombinationSignal 0.68 0.138 L 0.957 0.46 R 2.35 0.071 R 1.03 0.357 R
|
||||
x.likeSite 0.928 0.402 L 1.33 0.261 R 0.735 0.431 L 0.783 0.422 R
|
||||
|
||||
The third output file:
|
||||
|
||||
.. image:: ../static/operation_icons/dwt_var_perClass.png
|
||||
|
||||
</help>
|
||||
|
||||
</tool>
|
||||
Reference in New Issue
Block a user