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

Annotation of /FigKernelScripts/compute_PCA_for_metagenomes.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 cmp $b } keys(%all_subsystems);
55 :    
56 :     my $ss_count = scalar(@sorted);
57 :    
58 :     print "ss_count: $ss_count\n";
59 :    
60 :     my $header_ss_string = join("\t",@sorted);
61 :    
62 :     print OUT "metagenome\t$header_ss_string\n";
63 :    
64 :     foreach my $job_id (@jobs){
65 :     my $total_hits = $mg_ss_total_hits{$job_id};
66 :     my @cells = ();
67 :     push(@cells,$job_id);
68 :     foreach my $ss (@sorted){
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 :    
77 :     my $cell_count = scalar(@cells);
78 :     print "cell_count for $job_id: $cell_count\n";
79 :    
80 :    
81 :     my $row_string = join("\t",@cells);
82 :     print OUT "$row_string\n";
83 :     }
84 :    
85 :     close(OUT);
86 :    
87 :     open(OUT,">$output_dir/pca_script.r");
88 :     print OUT qq(m= read.table('$output_dir/data_matrix.txt', header= TRUE, row.names = 1, sep = '\\t')\n);
89 :     print OUT qq(pca_obj = prcomp(m)\n);
90 :     print OUT qq(s = summary(pca_obj)\n);
91 :     print OUT qq(write.table(s[2],'$output_dir/pca_rotational_matrix.txt', sep ='\\t')\n);
92 :     print OUT qq(write.table(s[1],'$output_dir/pca_std_deviation.txt', sep = '\\t')\n);
93 :     close(OUT);
94 :    
95 :     my $cmd = qq( R CMD BATCH $output_dir/pca_script.r );
96 :     my $errorlevel = system($cmd);
97 :    
98 :     if ($errorlevel) {
99 :     print "error with pca run\n";
100 :     }
101 :     else{
102 :     print "good pca run\n";
103 :     }
104 :    
105 :    
106 :     my @Xs;
107 :     my @Ys;
108 :    
109 :     open(IN,"$output_dir/pca_rotational_matrix.txt");
110 :     while($_ = <IN>){
111 :     chomp($_);
112 :     if($_ =~/rotation.PC1/){next;}
113 :     my @parts = split("\t",$_);
114 :     push(@Xs,$parts[1]);
115 :     push(@Ys,$parts[2]);
116 :     }
117 :     close(IN);
118 :    
119 :     my $x_string = join(",",@Xs);
120 :     my $y_string = join(",",@Ys);
121 :     open(OUT,">$output_dir/graphic_maker.r");
122 :     print OUT "pdf(file ='$output_dir/pc1_vs_pc2.pdf')\n";
123 :     print OUT "pc1<-c($x_string)\n";
124 :     print OUT "pc2<-c($y_string)\n";
125 :     print OUT "plot(pc1,pc2,xlim=c(-1,1),ylim=c(-1,1))\n";
126 :     print OUT "dev.off()\n";
127 :     close(OUT);
128 :    
129 :     my $cmd = qq( R CMD BATCH $output_dir/graphic_maker.r );
130 :     my $errorlevel = system($cmd);
131 :    
132 :     if ($errorlevel) {
133 :     print "error with making graphic\n";
134 :     }
135 :     else{
136 :     print "png made\n";
137 :     }
138 :    
139 :     #`R CMD BATCH $output_dir/pca_script.r`;
140 :    
141 :    
142 :     ############################
143 :     #
144 :     # Accessor functions
145 :     #
146 :     ############################
147 :    
148 :     sub get_mgdb {
149 :     my ($job) = @_;
150 :     unless (ref $job and $job->isa('MG_jobcache::Job')) {
151 :     warn "ERROR get_database_handle must be called with job object";
152 :     return undef;
153 :     } else {
154 :     return MGRAST::MetagenomeAnalysis->new($job);
155 :     }
156 :     }
157 :    
158 :     ############################
159 :     #
160 :     # Internal functions
161 :     #
162 :     ############################
163 :    
164 :     sub _get_job {
165 :     my ($jobnum) = @_;
166 :     my ($job) = $jobcache->Job->get_objects({id=>$jobnum})->[0];
167 :    
168 :     if (defined $job){
169 :     return $job;
170 :     } else {
171 :     print "no job returned from jobcache\n";
172 :     return undef;
173 :     }
174 :     }
175 :    
176 :     sub _jobdb {
177 :     my $handle;
178 :     print $FIG_Config::mgrast_jobcache_db."\t".$FIG_Config::mgrast_jobcache_host."\t".$FIG_Config::mgrast_jobcache_user."\n";
179 :     eval {
180 :     $handle = DBMaster->new( -database => $FIG_Config::mgrast_jobcache_db || 'JobCacheMGRast',
181 :     -host => $FIG_Config::mgrast_jobcache_host,
182 :     -user => $FIG_Config::mgrast_jobcache_user,
183 :     -password => $FIG_Config::mgrast_jobcache_password );
184 :     };
185 :     if ($@) {
186 :     warn "Unable to connect to MGRAST database: $@\n";
187 :     $handle = undef;
188 :     }
189 :     #else{ print "connected to db\n";}
190 :    
191 :     return $handle;
192 :     }
193 :    
194 :     sub get_key2taxa_mapping {
195 :    
196 :     unless (defined($_[0]->{key2taxa})) {
197 :     my $sth = $_[0]->dbh->prepare("select dbkey, str from tax_item");
198 :     $sth->execute;
199 :     $_[0]->{key2taxa} = $sth->fetchall_hashref('dbkey');
200 :     }
201 :    
202 :     return $_[0]->{key2taxa};
203 :    
204 :     }
205 :    
206 :    
207 :     sub key2taxa {
208 :     if(defined $_[1]) {
209 :     my $t = $_[0]->{key2taxa}->{$_[1]}->{str};
210 :     $t =~ s/\t+$//;
211 :     return $t;
212 :     }
213 :     return '';
214 :     }

MCS Webmaster
ViewVC Help
Powered by ViewVC 1.0.3