# HG changeset patch # User Jim Johnson # Date 1349730613 18000 # Node ID 728b0aa0df6e26f7558a4ced0d87f07ece1b9c69 # Parent de024d31e777176388a1d0f36c31e5218fea81aa Add cuffdiff_mds_plot for MultiDimensionalScaling plot for read_group_tracking diff -r de024d31e777 -r 728b0aa0df6e cuffdiff_mds_plot.pl --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/cuffdiff_mds_plot.pl Mon Oct 08 16:10:13 2012 -0500 @@ -0,0 +1,83 @@ +#!/usr/bin/perl -w + +############################################################### +# cuffdiff_mds_plot.pl +# John Garbe +# Septmeber 2012 +# +# Given a sample tracking file from cuffdiff2, convert it to +# the proper format for loading into R and generating an MDS plot +# +################################################################ + +# check to make sure having correct files +my $usage = "usage: cuffdiff_mds_plot.pl [TABULAR.in] [TABULAR.out] [PLOT.png]\n"; +die $usage unless @ARGV == 3; + +#get the input arguments +my $inputFile = $ARGV[0]; +my $outputFile = $ARGV[1]; +my $plotFile = $ARGV[2]; + +#Open files +open (INPUT, "<", $inputFile) || die("Could not open file $inputFile \n"); +open (OUTPUT, ">", $outputFile) || die("Could not open file $outputFile \n"); +open (PLOT, ">", $plotFile) || die("Could not open file $plotFile \n"); + +# header looks like this: +# tracking_id condition replicate raw_frags internal_scaled_frags external_scaled_frags FPKM effective_length status +my $header = ; + +# read in the sample tracking file +while () { + chomp; + @line = split /\t/; + $tracking_id{$line[0]} = 1; + $sample = $line[1] . "-" . $line[2]; + $fpkm{$sample}{$line[0]} = $line[6]; +} +close(INPUT); + +@sorted_tracking_id = sort( keys(%tracking_id)); + +# print out header +foreach $tracking_id (@sorted_tracking_id) { + print OUTPUT "\t$tracking_id"; +} +print OUTPUT "\n"; + +# print out data +foreach $sample (keys(%fpkm)) { + print OUTPUT "$sample"; + + foreach $tracking_id (@sorted_tracking_id) { + print OUTPUT "\t$fpkm{$sample}{$tracking_id}"; + } + + print OUTPUT "\n"; +} +close(OUTPUT); + +#variables to store the name of the R script file +my $r_script = "cuffinks2mdf.r"; + +open(Rcmd,">", $r_script) or die "Cannot open $r_script \n\n"; +print Rcmd " + datat <- read.table(\"$outputFile\"); + cmd <- cmdscale(dist(datat)); + png(filename=\"$plotFile\"); + plot(cmd[,1], cmd[,2], type=\"n\", ann=FALSE); + text(cmd[,1], cmd[,2], labels = rownames(datat)); + title(main=\"Multidimensional Scaling Plot\"); + title(xlab= \"Dimension 1\"); + title(ylab= \"Dimension 2\"); + dev.off(); + #eof" . "\n"; + +close Rcmd; + + +system("R --no-restore --no-save --no-readline < $r_script > $r_script.out"); +#close the input and output files +close(PLOT); + diff -r de024d31e777 -r 728b0aa0df6e cuffdiff_mds_plot.xml --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/cuffdiff_mds_plot.xml Mon Oct 08 16:10:13 2012 -0500 @@ -0,0 +1,20 @@ + + multidimensional scaling (MDS) plot from cuffdiff read_group_tracking + + cuffdiff_mds_plot.pl $input data.tmp $plot + + + + + + + + + + + + + Generates a multidimensional scaling (MDS) plot of a read_group_tracking output file from cuffdiff2. + +