Mercurial > repos > devteam > dwt_var_perclass
changeset 0:cb422b6f49d2 draft
Imported from capsule None
| author | devteam |
|---|---|
| date | Mon, 27 Jan 2014 09:26:11 -0500 |
| parents | |
| children | 781e68074f84 |
| files | execute_dwt_var_perClass.pl execute_dwt_var_perClass.xml |
| diffstat | 2 files changed, 425 insertions(+), 0 deletions(-) [+] |
line wrap: on
line diff
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/execute_dwt_var_perClass.pl Mon Jan 27 09:26:11 2014 -0500 @@ -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); \ No newline at end of file
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/execute_dwt_var_perClass.xml Mon Jan 27 09:26:11 2014 -0500 @@ -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_path}/operation_icons/dwt_var_perClass.png + + </help> + +</tool>