5 my $library_info = shift;
8 my $expxml = XML::Simple->new();
9 my $xml = $expxml->XMLin("Project.xml", ForceArray => [ qw(ComparePeakCalls CompareLibraries PeakCalling ProfileReads qPCR MotifFinding)]);
10 my @peak_calling = ();
12 my $project_name = $xml->{Name};
13 my $project_id= $xml->{ProjectId};
15 if(exists($xml->{PeakCalling})) {
16 for my $i (0..scalar(@{$xml->{PeakCalling}})-1) {
17 my $task = $xml->{PeakCalling}->[$i]->{TaskId};
18 my $name = $xml->{PeakCalling}->[$i]->{Name};
19 my $caller = $xml->{PeakCalling}->[$i]->{Caller};
20 my $genome = $xml->{PeakCalling}->[$i]->{Genome};
21 my $signal = $xml->{PeakCalling}->[$i]->{Signal}->{Library};
22 my $bg = $xml->{PeakCalling}->[$i]->{Background}->{Library};
23 $libs{$signal} = 0; $libs{$bg} = 0;
25 my $caller_dir = "../../Tasks/$task";
26 my %desc = (Name => $name, Dir => $caller_dir, Caller => $caller);
29 if($caller eq "QuEST") {
30 $desc{Summary} = `~/EXPTRACK/ProjectScripts/SummarizeQuEST.pm $caller_dir`;
31 $desc{outfile} = "$caller_dir/peak_caller.ChIP.out.bedgraph";
32 $desc{fasta} = "$caller_dir/peak_caller.ChIP.out.fasta";
33 } elsif($caller eq "WingPeaks") {
34 $desc{Summary} = `~/EXPTRACK/ProjectScripts/SummarizeWingPeaks.pm $caller_dir/$name.peaks`;
35 $desc{outfile} = "$caller_dir/$name.peaks.bed";
36 $desc{fasta} = "$caller_dir/$name.peaks.fasta";
37 } elsif($caller eq "MACS") {
38 my $peakfile = $caller_dir.'/'.$name.'_peaks.xls';
39 my $negpeakfile = $caller_dir.'/'.$name.'_negative_peaks.xls';
40 $desc{Summary} = `~/EXPTRACK/ProjectScripts/SummarizeMACS.pm $peakfile $negpeakfile`;
41 $desc{outfile} = "$caller_dir/".$name."_peaks.bed";
42 $desc{fasta} = "$caller_dir/".$name."_peaks.fasta";
45 if($genome eq "scer") { $genome = "sacCer1"; }
46 $desc{Genome} = $genome;
48 push @peak_calling, \%desc;
52 my $qPCR_Summary = "";
53 if(exists($xml->{qPCR})) {
54 $qPCR_Summary .= "<H2><EM>In-silico</EM> qPCR Summary</H2><BR><TABLE BORDER=1><TR><TD><EM>Library</EM></TD><TD><EM>Factor</EM></TD><TD><EM>Fold Enrichment</EM></TD></TR>\n";
55 for my $i (0..scalar(@{$xml->{qPCR}})-1) {
56 my $task = $xml->{qPCR}->[$i]->{TaskId};
57 my $name = $xml->{qPCR}->[$i]->{Name}; $name =~ s/\s/_/g;
58 my $genome = $xml->{qPCR}->[$i]->{Genome};
59 my $lib = $xml->{qPCR}->[$i]->{Library};
60 my $outfile = "../../Tasks/$task/$name.qPCR";
64 if( ! -e $outfile ) { $qPCR_Summary .= "<TR BCOLOR=#FFBBBB><TD>$name</TD><TD>Processing...</TD></TR>\n"; } else {
65 my $summary_line = `~/EXPTRACK/ProjectScripts/Summarize_qPCR.pm $name $lib $outfile`;
66 if($summary_line eq "") { $qPCR_Summary .= "<TR BCOLOR=#FFBBBB><TD>$name</TD><TD>Processing...</TD></TR>\n"; }
68 $qPCR_Summary .= $summary_line;
72 $qPCR_Summary .= "</TABLE>\n";
75 my $profile_summary = "";
76 if(exists($xml->{ProfileReads})) {
77 $profile_summary = "<BR><H2>Read Profiles (<A HREF=http://genome.ucsc.edu/cgi-bin/hgGateway>Genome Browser</A>)</H2>\n";
78 $profile_summary .= "<TABLE BORDER=1><TR><TD><EM>Library</EM></TD><TD><EM>Wiggle file</EM></TD></TR>";
79 for my $i (0..scalar(@{$xml->{ProfileReads}})-1) {
80 my $task = $xml->{ProfileReads}->[$i]->{TaskId};
81 my $name = $xml->{ProfileReads}->[$i]->{Name};
82 my $genome = $xml->{ProfileReads}->[$i]->{Genome};
83 my $lib = $xml->{ProfileReads}->[$i]->{Library};
86 if($genome eq "scer") { $genome = "sacCer1"; }
87 $profile_summary .= "<TR><TD>$name ($lib)</TD><TD><A HREF=../../Tasks/$task/$lib.wig.gz>$lib.wig.gz</A><BR><A HREF=http://genome.ucsc.edu/cgi-bin/hgTracks?db=$genome&position=chr1:1000-2000&hgt.customText=http://171.65.76.194/Tasks/$task/$lib.wig.gz TARGET=\"_blank\">View in UCSC Genome Browser</A></TD>\n";
88 $profile_summary .= "<TD><A HREF=../../Tasks/$task/$lib.profile.gif><IMG WIDTH=300 SRC=../../Tasks/$task/$lib.profile.gif></A></TD></TR>\n";
90 $profile_summary .= "</TABLE>\n";
93 my $library_comparisons = "";
94 if(exists($xml->{CompareLibraries})) {
95 $library_comparisons .= "<H2>Library Comparisons</H2>\n";
96 $library_comparisons .= "<TABLE BORDER=1>\n";
97 $library_comparisons .= "<TR><TD><EM>Factor</EM></TD><TD><EM>Library 1</EM></TD><TD><EM>Library 2</EM></TD><TD><EM>Correlation</EM></TD></TR>\n";
98 for my $i (0..scalar(@{$xml->{CompareLibraries}})-1) {
99 my $task = $xml->{CompareLibraries}->[$i]->{TaskId};
100 my $tf = $xml->{CompareLibraries}->[$i]->{TF};
101 my $genome = $xml->{CompareLibraries}->[$i]->{Genome};
102 my $features = "~/EXPTRACK/compare_peak_calls/".$genome."_uptream5k_downtream1k";
103 my $name1 = $xml->{CompareLibraries}->[$i]->{Library}->[0]->{Library};
104 my $name2 = $xml->{CompareLibraries}->[$i]->{Library}->[1]->{Library};
105 my $outfile = "../../Tasks/".$task.'/'.$name1."_".$name2.".compare";
106 $libs{$name1} = 0; $libs{$name2} = 0;
108 if(-e $outfile) { $correlation = `cat $outfile | awk '{print \$8}'`; }
110 if( !(-e $outfile) ) { $correlation = "In Progress..."; $color = "#FFBBBB"; }
111 elsif($correlation > 0.9) { $color = "#BBFFBB"; }
112 elsif($correlation > 0.6) { $color = "#BBBBFF"; }
113 else { $color = "#FFBBBB"; }
114 $library_comparisons .= "<TR BGCOLOR=$color><TD>$tf</TD><TD>$name1</TD><TD>$name2</TD<TD>$correlation</TD></TR>\n";
116 $library_comparisons .= "</TABLE>\n";
119 my $peak_call_comparisons = "";
120 if(exists($xml->{ComparePeakCalls})) {
121 $peak_call_comparisons = "<H2>Peak Calling Comparisons</H2>\n";
122 $peak_call_comparisons .= "<TABLE BORDER=1><TR><TD><EM>Caller1</EM></TD><TD><EM>Set1</EM></TD><TD><EM>Caller2</EM></TD><TD><EM>Set2</EM></TD><TD>|Union|</TD><TD>|Intersect|</TD><TD>|Difference|</TD><TD>Jaccard</TD><TD>Correlation</TD></TR>\n";
123 for my $i (0..scalar(@{$xml->{ComparePeakCalls}})-1) {
124 my $task = $xml->{ComparePeakCalls}->[$i]->{TaskId};
125 my $caller1 = $xml->{ComparePeakCalls}->[$i]->{Caller1};
126 my $caller2 = $xml->{ComparePeakCalls}->[$i]->{Caller2};
127 my $genome = $xml->{ComparePeakCalls}->[$i]->{Genome};
128 my $set1 = $xml->{ComparePeakCalls}->[$i]->{Set1}; $set1 =~ s/\s/_/g;
129 my $set2 = $xml->{ComparePeakCalls}->[$i]->{Set2}; $set2 =~ s/\s/_/g;
130 my $outfile .= "../../Tasks/".$task.'/'.$set1."-".$caller1."_v_".$set2."-".$caller2.".compare";
131 my $correlation = -1;
134 my $color = "#FFBBBB";
135 $peak_call_comparisons .= "<TR BGCOLOR=$color><TD>$caller1</TD><TD>$set1</TD><TD>$caller2</TD><TD>$set2</TD><TD COLSPAN=5><B>Processing...</B></TD></TR>\n";
137 my $line = `cat $outfile`; chomp $line;
138 my($tf,$set1_id,$set2_id,$U,$I,$D,$J,$R) = split(/\t/,$line);
139 my $color = "#FFBBBB";
140 if(!defined($R)) { $peak_call_comparisons .= "<TR BGCOLOR=$color><TD>$caller1</TD><TD>$set1</TD><TD>$caller2</TD><TD>$set2</TD><TD COLSPAN=5><B>Processing...</B></TD></TR>\n"; } else {
141 $set1_id =~ /.+?\((\d+)\)$/; my $n1 = $1;
142 $set2_id =~ /.+?\((\d+)\)$/; my $n2 = $1;
143 if($R > 0.85) { $color = "#BBFFBB"; } elsif($R > 0.5) { $color="#BBBBFF"; } else { $color="#FFBBBB"; }
144 $peak_call_comparisons .= "<TR BGCOLOR=$color><TD>$caller1</TD><TD>$set1 ($n1)</TD><TD>$caller2</TD><TD>$set2 ($n2)</TD><TD>$U</TD><TD>$I</TD><TD>$D</TD><TD>$J</TD><TD>$R</TD></TR>\n";
148 $peak_call_comparisons .= "</TABLE>\n";
151 my $motif_finding = "";
152 if(exists($xml->{MotifFinding})) {
153 $motif_finding = "<H2>Motif Finding</H2>\n";
154 $motif_finding .= "<TABLE BORDER=1><TR><TD><EM>Caller</EM></TD><TD><EM>Set</EM></TD><TD><EM>Width</EM></TD><TD>Motifs</TD><TD><EM>Motif Information</EM></TD></TR>\n";
155 for my $i (0..scalar(@{$xml->{MotifFinding}})-1) {
156 my $task = $xml->{MotifFinding}->[$i]->{TaskId};
157 my $caller = $xml->{MotifFinding}->[$i]->{Caller};
158 my $genome = $xml->{MotifFinding}->[$i]->{Genome};
159 my $set = $xml->{MotifFinding}->[$i]->{Set}; $set =~ s/\s/_/g;
160 my $width = $xml->{MotifFinding}->[$i]->{Width};
161 my $outfile .= "../../Tasks/".$task.'/'.$set."-".$caller.".w".$width.".biop";
164 my $color = "#FFBBBB";
165 $motif_finding .= "<TR BGCOLOR=$color><TD>$caller</TD><TD>$set</TD><TD COLSPAN=3><B>Processing...</B></TD></TR>\n";
167 my $color = "#BBFFBB";
168 my $line = `cat $outfile | grep "Motif #"`; $line =~ s/\n/<BR>/g;
169 my $score = `cat $outfile | grep "MotifScore"`; $score=~ s/\n/<BR>/g;
170 $motif_finding .= "<TR BGCOLOR=$color><TD>$caller</TD><TD>$set</TD><TD>$width</TD><TD>$line</TD><TD>$score</TD></TR>\n";
173 $motif_finding .= "</TABLE>\n";
176 print "<HTML><HEAD><TITLE>$project_name Analysis Summary</TITLE></HEAD><BODY>\n";
177 print "<H1><CENTER>$project_name Analysis Summary</CENTER></H1>\n";
178 print SummarizeLibrary($library_info);
179 print $profile_summary;
180 print $library_comparisons;
181 print $peak_call_comparisons;
182 print $motif_finding;
185 if(exists($xml->{PeakCalling})) {
186 print "<H2>Peak Calling Summary</H2>";
187 print "<TABLE BORDER=1><TR><TD><EM>Experiment</EM></TD><TD><EM>Peak Caller</EM></TD><TD><EM>Summary</EM></TD></TR>\n";
192 #my $color = $row % 2?"#CCCCFF":"#FFFFFF";
193 my $color = "#EEEEFF";
194 print "<TR BGCOLOR=$color><TD><B>$hash{Name}</B></TD>";
195 print "<TD>$hash{Caller}</TD><TD>$hash{Summary}</TD>\n";
196 print "<TD><A HREF=$hash{outfile}>BED file</A><BR><A HREF=http://genome.ucsc.edu/cgi-bin/hgTracks?db=$hash{Genome}&hgt.customText=http://171.65.76.194/Tasks/$hash{Task}/$hash{outfile}>View in Genome Browser</A></TD><TD><A HREF=$hash{fasta}>FASTA</A></TD>\n";
201 print "</BODY></HTML>\n";
202 sub SummarizeLibrary {
203 my $xml_file = shift;
204 my $xml = XMLin($xml_file, ForceArray => 1);
206 $summary .= "<H2>Libraries Used</H2><TABLE BORDER=1><TR><TD><EM>Library</EM></TD><TD><EM>Name</EM></TD><TD><EM>Num Lanes</EM></TD><TD><EM>Num Reads</EM></TD></TR>\n";
207 for my $i (0..scalar(@{$xml->{Library}})-1) {
208 my $lib = $xml->{Library}->[$i]->{Name};
209 next unless(scalar(keys %libs) == 0 || exists($libs{$lib}));
210 my $N = scalar(@{$xml->{Library}->[$i]->{Track}});
211 my($date,$fc,$lane,$desc);
214 for my $t (0..$N-1) {
215 my $filename = $xml->{Library}->[$i]->{Track}->[$t]->{Filename};
216 $filename = `basename $filename`;
217 $filename =~ /^(\d+)_(.+?)_s(\d+)_(.+?)_$lib.align/;
218 ($date,$fc,$lane,$desc) = ($1,$2,$3,$4);
219 $num_lanes += length($lane);
220 $num_reads += $xml->{Library}->[$i]->{Track}->[$t]->{Count};
223 if($num_reads < 3000000) { $bgcolor = "FF3300"; }
224 elsif($num_reads < 5000000) { $bgcolor = "FFCC33"; }
225 elsif($num_reads < 10000000) { $bgcolor = "00CCFF"; }
226 else { $bgcolor = "66FF66"; }
228 $summary .= sprintf("<TR BGCOLOR=#%s><TD>%s</TD><TD>%s</TD><TD>%d</TD><TD>%0.2fM</TD></TR>\n",$bgcolor,$lib,$desc,$num_lanes,$num_reads/1000000.0);
230 $summary .= "</TABLE>\n";