Mercurial > repos > jjohnson > cummerbund
comparison cuffdiff_mds_plot.pl @ 3:728b0aa0df6e
Add cuffdiff_mds_plot for MultiDimensionalScaling plot for read_group_tracking
author | Jim Johnson <jj@umn.edu> |
---|---|
date | Mon, 08 Oct 2012 16:10:13 -0500 |
parents | |
children |
comparison
equal
deleted
inserted
replaced
2:de024d31e777 | 3:728b0aa0df6e |
---|---|
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 |