annotate cuffdiff_mds_plot.pl @ 0:da7241f92ecf

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