[Bio] / FigKernelScripts / compute_PCA_with_metagenomes_as_variables.pl Repository:
ViewVC logotype

Annotation of /FigKernelScripts/compute_PCA_with_metagenomes_as_variables.pl

Parent Directory Parent Directory | Revision Log Revision Log


Revision 1.1 - (view) (download) (as text)

1 : mkubal 1.1 use FIG_Config;
2 :     use DBMaster;
3 :     use MGRAST::MetagenomeAnalysis;
4 :     use DBMaster;
5 :     use strict;
6 :     use warnings;
7 :    
8 :     if(scalar(@ARGV) == 0 ){
9 :     print "usage: perl compute_PCA_for_metagenomes.pl \$dir_for_output \$metagenome_job_id_1 \$metagenome_job_id_2 ...\n";
10 :     exit;
11 :     }
12 :    
13 :     my $output_dir = shift(@ARGV);
14 :    
15 :     my @jobs;
16 :     foreach my $arg (@ARGV){
17 :     if($arg =~/\d+/){
18 :     push(@jobs,$arg);
19 :     #print "job: $arg\n";
20 :     }
21 :     }
22 :    
23 :     my %all_subsystems;
24 :     my %mg_ss_total_hits;
25 :     my %mg_ss_hits;
26 :    
27 :     my $jobcache = &_jobdb();
28 :    
29 :     foreach my $job_id (@jobs){
30 :     my $total_hits = 0;
31 :     my $job_obj = &_get_job($job_id);
32 :     if(! $job_obj){next;}
33 :     else{print "good id: $job_id\n";}
34 :     my $mgdb = &get_mgdb($job_obj);
35 :    
36 :     my $ss = $mgdb->get_subsystem_counts("SEED:subsystem_tax");
37 :     foreach my $s (@$ss){
38 :     my $name = $mgdb->key2taxa(@$s[2]);
39 :     if($name){
40 :     $all_subsystems{$name} = 1;
41 :     my $ss_peg_count = @$s[4];
42 :     my $combo = $job_id."_".$name;
43 :     $mg_ss_hits{$combo} = $ss_peg_count;
44 :     $total_hits = $ss_peg_count + $total_hits;
45 :     #print "ss_name: $name\n";
46 :     }
47 :     }
48 :    
49 :     $mg_ss_total_hits{$job_id} = $total_hits;
50 :     #print "$job_id th: $total_hits\n";
51 :     }
52 :    
53 :     my $data_matrix = open(OUT,">$output_dir/data_matrix.txt");
54 :     my @sorted = sort { $a <=> $b } @jobs;
55 :    
56 :     my $mg_count = scalar(@sorted);
57 :    
58 :     my $header_mg_string = join("\t",@sorted);
59 :    
60 :     print OUT "metagenome\t$header_mg_string\n";
61 :    
62 :     my $count = 1;
63 :     foreach my $ss (keys(%all_subsystems)){
64 :     my @cells = ();
65 :     #push(@cells,$ss);
66 :     push(@cells,$count);
67 :     foreach my $job_id (@sorted){
68 :     my $total_hits = $mg_ss_total_hits{$job_id};
69 :     my $combo = $job_id."_".$ss;
70 :     my $ss_hits = $mg_ss_hits{$combo};
71 :     my $cell_entry;
72 :     if($ss_hits){$cell_entry = $ss_hits/$total_hits;}
73 :     else{$cell_entry = 0;}
74 :     push(@cells,$cell_entry);
75 :     }
76 :     $count++;
77 :    
78 :     #my $cell_count = scalar(@cells);
79 :     #print "cell_count for $job_id: $cell_count\n";
80 :     my $row_string = join("\t",@cells);
81 :     print OUT "$row_string\n";
82 :     }
83 :    
84 :     close(OUT);
85 :    
86 :     open(OUT,">$output_dir/pca_script.r");
87 :     print OUT qq(m= read.table('$output_dir/data_matrix.txt', header= TRUE, row.names = 1, sep = '\\t')\n);
88 :     print OUT qq(pca_obj = prcomp(m)\n);
89 :     print OUT qq(s = summary(pca_obj)\n);
90 :     print OUT qq(write.table(s[2],'$output_dir/pca_rotational_matrix.txt', sep ='\\t')\n);
91 :     print OUT qq(write.table(s[1],'$output_dir/pca_std_deviation.txt', sep = '\\t')\n);
92 :     close(OUT);
93 :    
94 :     my $cmd = qq( R CMD BATCH $output_dir/pca_script.r );
95 :     my $errorlevel = system($cmd);
96 :    
97 :     if ($errorlevel) {
98 :     print "error with pca run\n";
99 :     }
100 :     else{
101 :     print "good pca run\n";
102 :     }
103 :    
104 :    
105 :     my @Xs;
106 :     my @Ys;
107 :    
108 :     open(IN,"$output_dir/pca_rotational_matrix.txt");
109 :     while($_ = <IN>){
110 :     chomp($_);
111 :     if($_ =~/rotation.PC1/){next;}
112 :     my @parts = split("\t",$_);
113 :     push(@Xs,$parts[1]);
114 :     push(@Ys,$parts[2]);
115 :     }
116 :     close(IN);
117 :    
118 :     my $x_string = join(",",@Xs);
119 :     my $y_string = join(",",@Ys);
120 :     open(OUT,">$output_dir/graphic_maker.r");
121 :     print OUT "pdf(file ='$output_dir/pc1_vs_pc2.pdf')\n";
122 :     print OUT "pc1<-c($x_string)\n";
123 :     print OUT "pc2<-c($y_string)\n";
124 :     print OUT "plot(pc1,pc2,xlim=c(-1,1),ylim=c(-1,1))\n";
125 :     print OUT "dev.off()\n";
126 :     close(OUT);
127 :    
128 :     my $cmd = qq( R CMD BATCH $output_dir/graphic_maker.r );
129 :     my $errorlevel = system($cmd);
130 :    
131 :     if ($errorlevel) {
132 :     print "error with making graphic\n";
133 :     }
134 :     else{
135 :     print "png made\n";
136 :     }
137 :    
138 :     #`R CMD BATCH $output_dir/pca_script.r`;
139 :    
140 :    
141 :     ############################
142 :     #
143 :     # Accessor functions
144 :     #
145 :     ############################
146 :    
147 :     sub get_mgdb {
148 :     my ($job) = @_;
149 :     unless (ref $job and $job->isa('MG_jobcache::Job')) {
150 :     warn "ERROR get_database_handle must be called with job object";
151 :     return undef;
152 :     } else {
153 :     return MGRAST::MetagenomeAnalysis->new($job);
154 :     }
155 :     }
156 :    
157 :     ############################
158 :     #
159 :     # Internal functions
160 :     #
161 :     ############################
162 :    
163 :     sub _get_job {
164 :     my ($jobnum) = @_;
165 :     my ($job) = $jobcache->Job->get_objects({id=>$jobnum})->[0];
166 :    
167 :     if (defined $job){
168 :     return $job;
169 :     } else {
170 :     print "no job returned from jobcache\n";
171 :     return undef;
172 :     }
173 :     }
174 :    
175 :     sub _jobdb {
176 :     my $handle;
177 :     print $FIG_Config::mgrast_jobcache_db."\t".$FIG_Config::mgrast_jobcache_host."\t".$FIG_Config::mgrast_jobcache_user."\n";
178 :     eval {
179 :     $handle = DBMaster->new( -database => $FIG_Config::mgrast_jobcache_db || 'JobCacheMGRast',
180 :     -host => $FIG_Config::mgrast_jobcache_host,
181 :     -user => $FIG_Config::mgrast_jobcache_user,
182 :     -password => $FIG_Config::mgrast_jobcache_password );
183 :     };
184 :     if ($@) {
185 :     warn "Unable to connect to MGRAST database: $@\n";
186 :     $handle = undef;
187 :     }
188 :     #else{ print "connected to db\n";}
189 :    
190 :     return $handle;
191 :     }
192 :    
193 :     sub get_key2taxa_mapping {
194 :    
195 :     unless (defined($_[0]->{key2taxa})) {
196 :     my $sth = $_[0]->dbh->prepare("select dbkey, str from tax_item");
197 :     $sth->execute;
198 :     $_[0]->{key2taxa} = $sth->fetchall_hashref('dbkey');
199 :     }
200 :    
201 :     return $_[0]->{key2taxa};
202 :    
203 :     }
204 :    
205 :    
206 :     sub key2taxa {
207 :     if(defined $_[1]) {
208 :     my $t = $_[0]->{key2taxa}->{$_[1]}->{str};
209 :     $t =~ s/\t+$//;
210 :     return $t;
211 :     }
212 :     return '';
213 :     }

MCS Webmaster
ViewVC Help
Powered by ViewVC 1.0.3