Mercurial > repos > jjohnson > cummerbund

comparison cuffdiff_mds_plot.pl @ 0:da7241f92ecf

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Uploaded
author jjohnson
date Mon, 04 Feb 2013 19:50:25 -0500
parents
children
comparison
equal deleted inserted replaced
-1:000000000000 0:da7241f92ecf
1 #!/usr/bin/perl -w
2
3 ###############################################################
4 # cuffdiff_mds_plot.pl
5 # John Garbe
6 # Septmeber 2012
7 #
8 # Given a sample tracking file from cuffdiff2, convert it to
9 # the proper format for loading into R and generating an MDS plot
10 #
11 ################################################################
12
13 # check to make sure having correct files
14 my $usage = "usage: cuffdiff_mds_plot.pl [TABULAR.in] [TABULAR.out] [PLOT.png]\n";
15 die $usage unless @ARGV == 3;
16
17 #get the input arguments
18 my $inputFile = $ARGV[0];
19 my $outputFile = $ARGV[1];
20 my $plotFile = $ARGV[2];
21
22 #Open files
23 open (INPUT, "<", $inputFile) || die("Could not open file $inputFile \n");
24 open (OUTPUT, ">", $outputFile) || die("Could not open file $outputFile \n");
25 open (PLOT, ">", $plotFile) || die("Could not open file $plotFile \n");
26
27 # header looks like this:
28 # tracking_id condition replicate raw_frags internal_scaled_frags external_scaled_frags FPKM effective_length status
29 my $header = <INPUT>;
30
31 # read in the sample tracking file
32 while (<INPUT>) {
33 chomp;
34 @line = split /\t/;
35 $tracking_id{$line[0]} = 1;
36 $sample = $line[1] . "-" . $line[2];
37 $fpkm{$sample}{$line[0]} = $line[6];
38 }
39 close(INPUT);
40
41 @sorted_tracking_id = sort( keys(%tracking_id));
42
43 # print out header
44 foreach $tracking_id (@sorted_tracking_id) {
45 print OUTPUT "\t$tracking_id";
46 }
47 print OUTPUT "\n";
48
49 # print out data
50 foreach $sample (keys(%fpkm)) {
51 print OUTPUT "$sample";
52
53 foreach $tracking_id (@sorted_tracking_id) {
54 print OUTPUT "\t$fpkm{$sample}{$tracking_id}";
55 }
56
57 print OUTPUT "\n";
58 }
59 close(OUTPUT);
60
61 #variables to store the name of the R script file
62 my $r_script = "cuffinks2mdf.r";
63
64 open(Rcmd,">", $r_script) or die "Cannot open $r_script \n\n";
65 print Rcmd "
66 datat <- read.table(\"$outputFile\");
67 cmd <- cmdscale(dist(datat));
68 png(filename=\"$plotFile\");
69 plot(cmd[,1], cmd[,2], type=\"n\", ann=FALSE);
70 text(cmd[,1], cmd[,2], labels = rownames(datat));
71 title(main=\"Multidimensional Scaling Plot\");
72 title(xlab= \"Dimension 1\");
73 title(ylab= \"Dimension 2\");
74 dev.off();
75 #eof" . "\n";
76
77 close Rcmd;
78
79
80 system("R --no-restore --no-save --no-readline < $r_script > $r_script.out");
81 #close the input and output files
82 close(PLOT);
83