htswanalysis/scripts/SummarizeProject2.pm

   1 #!/usr/bin/perl -w
   2 use strict;
   3 use warnings;
   4 use XML::Simple;
   5 my $root_dir = shift;
   6 my $library_info = shift;
   7 my %libs;
   8 my %quest;
   9 my $expxml = XML::Simple->new();
  10 my $xml = $expxml->XMLin("Project.xml", ForceArray => [ qw(ComparePeakCalls CompareLibraries PeakCalling ProfileReads qPCR MotifFinding Methylseq)]);
  11 my @peak_calling = ();
  12
  13 my $project_name = $xml->{Name};
  14 my $project_id= $xml->{ProjectId};
  15
  16 if(exists($xml->{PeakCalling})) {
  17 for my $i (0..scalar(@{$xml->{PeakCalling}})-1) {
  18   my $task = $xml->{PeakCalling}->[$i]->{TaskId};
  19   my $name = $xml->{PeakCalling}->[$i]->{Name};
  20   my $caller = $xml->{PeakCalling}->[$i]->{Caller};
  21   my $genome = $xml->{PeakCalling}->[$i]->{Genome};
  22   my $signal = $xml->{PeakCalling}->[$i]->{Signal}->{Library};
  23   my $bg = $xml->{PeakCalling}->[$i]->{Background}->{Library};
  24   $libs{$signal} = 0; $libs{$bg} = 0;
  25   $name =~ s/\s/_/ig;
  26   my $caller_dir = "../../Tasks/$task";
  27   my %desc = (Name => $name, Dir => $caller_dir, Caller => $caller);
  28   $desc{outfile} = "";
  29   $desc{fasta} = "";
  30   if($caller eq "QuEST") {
  31     $desc{Summary} = `$root_dir/scripts/SummarizeQuEST.pm $caller_dir`;
  32     $desc{outfile} = "$caller_dir/peak_caller.ChIP.out.bedgraph";
  33     $desc{fasta} = "$caller_dir/peak_caller.ChIP.out.fasta";
  34   } elsif($caller eq "WingPeaks") {
  35     $desc{Summary} = `$root_dir/scripts/SummarizeWingPeaks.pm $caller_dir/$name.peaks`;
  36     $desc{outfile} = "$caller_dir/$name.peaks.bed";
  37     $desc{fasta} = "$caller_dir/$name.peaks.fasta";
  38   } elsif($caller eq "MACS") {
  39     my $peakfile = $caller_dir.'/'.$name.'_peaks.xls';
  40     my $negpeakfile = $caller_dir.'/'.$name.'_negative_peaks.xls';
  41     $desc{Summary} = `$root_dir/scripts/SummarizeMACS.pm $peakfile $negpeakfile`;
  42     $desc{outfile} = "$caller_dir/".$name."_peaks.bed";
  43     $desc{fasta} = "$caller_dir/".$name."_peaks.fasta";
  44     $desc{primer_design} = "$caller_dir/ValidationPrimers.html";
  45   }
  46
  47   if($genome eq "scer") { $genome = "sacCer1"; }
  48   $desc{Genome} = $genome;
  49   $desc{Task} = $task;
  50   push @peak_calling, \%desc;
  51 }
  52 }
  53
  54 my $qPCR_Summary = "";
  55 if(exists($xml->{qPCR})) {
  56 $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";
  57 for my $i (0..scalar(@{$xml->{qPCR}})-1) {
  58   my $task = $xml->{qPCR}->[$i]->{TaskId};
  59   my $name = $xml->{qPCR}->[$i]->{Name}; $name =~ s/\s/_/g;
  60   my $genome = $xml->{qPCR}->[$i]->{Genome};
  61   my $lib = $xml->{qPCR}->[$i]->{Library};
  62   my $outfile = "../../Tasks/$task/$name.qPCR";
  63
  64   $libs{$lib} = 0;
  65
  66   if( ! -e $outfile ) {  $qPCR_Summary .= "<TR BCOLOR=#FFBBBB><TD>$name</TD><TD>Processing...</TD></TR>\n"; } else {
  67     my $summary_line = `$root_dir/scripts/Summarize_qPCR.pm $name $lib $outfile`;
  68     if($summary_line eq "") { $qPCR_Summary .= "<TR BCOLOR=#FFBBBB><TD>$name</TD><TD>Processing...</TD></TR>\n"; }
  69     else {
  70       $qPCR_Summary .= $summary_line;
  71     }
  72   }
  73 }
  74 $qPCR_Summary .= "</TABLE>\n";
  75 }
  76
  77 my $methylseq_summary = "";
  78 if(exists($xml->{Methylseq})) {
  79   $methylseq_summary .= "<H2>Methylseq Analysis</H2>\n";
  80   $methylseq_summary  .= "<TABLE BORDER=1>\n";
  81   $methylseq_summary .= "<TR>\n";
  82   $methylseq_summary .= "<TD><EM>Task</EM></TD>\n";
  83   $methylseq_summary .= "<TD><EM>msp1</EM></TD>\n";
  84   $methylseq_summary .= "<TD><EM>hpa2</EM></TD>\n";
  85   $methylseq_summary .= "<TD><EM>Assayed</EM></TD>\n";
  86   $methylseq_summary .= "<TD><EM>Methylated</EM></TD>\n";
  87   $methylseq_summary .= "<TD><EM>Errors</EM></TD>\n";
  88   $methylseq_summary .= "<TD><EM>Data</EM></TD></TR>\n";
  89
  90   for my $i (0..scalar(@{$xml->{Methylseq}})-1) {
  91     my $task = $xml->{Methylseq}->[$i]->{TaskId};
  92     my $genome = $xml->{Methylseq}->[$i]->{Genome};
  93     my $name = $xml->{Methylseq}->[$i]->{Name};
  94
  95     my $msp1   = $xml->{Methylseq}->[$i]->{Msp1}->{Library};
  96     my $hpa2   = $xml->{Methylseq}->[$i]->{Hpa2}->{Library};
  97
  98     $methylseq_summary .= `$root_dir/scripts/Summarize_Methylseq.pm $name $task ../../Tasks $msp1 $hpa2 $genome`;
  99
 100   }
 101   $methylseq_summary .= "</TABLE>\n";
 102 }
 103
 104
 105 my $profile_summary = "";
 106 if(exists($xml->{ProfileReads})) {
 107   $profile_summary = "<BR><H2>Read Profiles (<A HREF=http://genome.ucsc.edu/cgi-bin/hgGateway>Genome Browser</A>)</H2>\n";
 108   $profile_summary .= "<TABLE BORDER=1><TR><TD><EM>Library</EM></TD><TD><EM>Wiggle file</EM></TD></TR>";
 109   for my $i (0..scalar(@{$xml->{ProfileReads}})-1) {
 110     my $task = $xml->{ProfileReads}->[$i]->{TaskId};
 111     my $name = $xml->{ProfileReads}->[$i]->{Name};
 112     my $genome = $xml->{ProfileReads}->[$i]->{Genome};
 113     my $lib = $xml->{ProfileReads}->[$i]->{Library};
 114     $libs{$lib} = 0;
 115
 116     if($genome eq "scer") { $genome = "sacCer1"; }
 117     $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://illumina-mac.stanford.edu/Tasks/$task/$lib.wig.gz TARGET=\"_blank\">View in UCSC Genome Browser</A></TD>\n";
 118     $profile_summary .= "<TD><A HREF=../../Tasks/$task/$lib.profile.gif><IMG WIDTH=300 SRC=../../Tasks/$task/$lib.profile.gif></A></TD></TR>\n";
 119   }
 120   $profile_summary .= "</TABLE>\n";
 121 }
 122
 123 my $library_comparisons = "";
 124 if(exists($xml->{CompareLibraries})) {
 125 $library_comparisons .= "<H2>Library Comparisons</H2>\n";
 126 $library_comparisons .= "<TABLE BORDER=1>\n";
 127 $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";
 128 for my $i (0..scalar(@{$xml->{CompareLibraries}})-1) {
 129   my $task = $xml->{CompareLibraries}->[$i]->{TaskId};
 130   my $tf = $xml->{CompareLibraries}->[$i]->{TF};
 131   my $genome = $xml->{CompareLibraries}->[$i]->{Genome};
 132   my $features = "$root_dir/reference_data/".$genome."_uptream5k_downtream1k";
 133   my $name1 = $xml->{CompareLibraries}->[$i]->{Library}->[0]->{Library};
 134   my $name2 = $xml->{CompareLibraries}->[$i]->{Library}->[1]->{Library};
 135   my $outfile = "../../Tasks/".$task.'/'.$name1."_".$name2.".compare";
 136   $libs{$name1} = 0; $libs{$name2} = 0;
 137   my $correlation; my $color;
 138   if( !(-e $outfile) ) { $correlation = "In Progress..."; $color = "#FFBBBB"; }
 139   else {
 140     $correlation = `cat $outfile | awk '{print \$8}'`;
 141     if($correlation > 0.9) { $color = "#BBFFBB"; }
 142     elsif($correlation > 0.6) { $color = "#BBBBFF"; }
 143     else { $color = "#FFBBBB"; }
 144   }
 145   $library_comparisons .= "<TR BGCOLOR=$color><TD>$tf</TD><TD>$name1</TD><TD>$name2</TD<TD>$correlation</TD></TR>\n";
 146 }
 147 $library_comparisons .= "</TABLE>\n";
 148 }
 149
 150 my $peak_call_comparisons = "";
 151 if(exists($xml->{ComparePeakCalls})) {
 152   $peak_call_comparisons = "<H2>Peak Calling Comparisons</H2>\n";
 153   $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";
 154   for my $i (0..scalar(@{$xml->{ComparePeakCalls}})-1) {
 155     my $task = $xml->{ComparePeakCalls}->[$i]->{TaskId};
 156     my $caller1 = $xml->{ComparePeakCalls}->[$i]->{Caller1};
 157     my $caller2 = $xml->{ComparePeakCalls}->[$i]->{Caller2};
 158     my $genome = $xml->{ComparePeakCalls}->[$i]->{Genome};
 159     my $set1 = $xml->{ComparePeakCalls}->[$i]->{Set1}; $set1 =~ s/\s/_/g;
 160     my $set2 = $xml->{ComparePeakCalls}->[$i]->{Set2}; $set2 =~ s/\s/_/g;
 161     my $outfile .= "../../Tasks/".$task.'/'.$set1."-".$caller1."_v_".$set2."-".$caller2.".compare";
 162     my $correlation = -1;
 163
 164     if(! -e $outfile) {
 165       my $color = "#FFBBBB";
 166       $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";
 167     } else {
 168       my $line = `cat $outfile`; chomp $line;
 169       my($tf,$set1_id,$set2_id,$U,$I,$D,$J,$R) = split(/\t/,$line);
 170       my $color = "#FFBBBB";
 171       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 {
 172         $set1_id =~ /.+?\((\d+)\)$/; my $n1 = $1;
 173         $set2_id =~ /.+?\((\d+)\)$/; my $n2 = $1;
 174         if($R > 0.85) { $color = "#BBFFBB"; } elsif($R > 0.5) { $color="#BBBBFF"; } else { $color="#FFBBBB"; }
 175         $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";
 176       }
 177     }
 178   }
 179   $peak_call_comparisons .= "</TABLE>\n";
 180 }
 181
 182 my $motif_finding = "";
 183 if(exists($xml->{MotifFinding})) {
 184   $motif_finding = "<H2>Motif Finding</H2>\n";
 185   $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";
 186   for my $i (0..scalar(@{$xml->{MotifFinding}})-1) {
 187     my $task = $xml->{MotifFinding}->[$i]->{TaskId};
 188     my $caller = $xml->{MotifFinding}->[$i]->{Caller};
 189     my $genome = $xml->{MotifFinding}->[$i]->{Genome};
 190     my $set = $xml->{MotifFinding}->[$i]->{Set}; $set =~ s/\s/_/g;
 191     my $width = $xml->{MotifFinding}->[$i]->{Width};
 192     my $outfile .= "../../Tasks/".$task.'/'.$set."-".$caller.".w".$width.".biop";
 193
 194     if(! -e $outfile) {
 195       my $color = "#FFBBBB";
 196       $motif_finding .= "<TR BGCOLOR=$color><TD>$caller</TD><TD>$set</TD><TD COLSPAN=3><B>Processing...</B></TD></TR>\n";
 197     } else {
 198       my $color = "#BBFFBB";
 199       my $line = `cat $outfile | grep "Motif #"`; $line =~ s/\n/<BR>/g;
 200       my $score = `cat $outfile | grep "MotifScore"`; $score=~ s/\n/<BR>/g;
 201       $motif_finding .= "<TR BGCOLOR=$color><TD>$caller</TD><TD>$set</TD><TD>$width</TD><TD>$line</TD><TD>$score</TD></TR>\n";
 202     }
 203   }
 204   $motif_finding .= "</TABLE>\n";
 205 }
 206
 207 print "<HTML><HEAD><TITLE>$project_name Analysis Summary</TITLE></HEAD><BODY>\n";
 208 print "<H1><CENTER>$project_name Analysis Summary</CENTER></H1>\n";
 209 print SummarizeLibrary($library_info);
 210 print $profile_summary;
 211 print $methylseq_summary;
 212 print $library_comparisons;
 213 print $peak_call_comparisons;
 214 print $motif_finding;
 215 print $qPCR_Summary;
 216
 217 if(exists($xml->{PeakCalling})) {
 218 print "<H2>Peak Calling Summary</H2>";
 219 print "<TABLE BORDER=1><TR><TD><EM>Experiment</EM></TD><TD><EM>Peak Caller</EM></TD><TD><EM>Summary</EM></TD></TR>\n";
 220 my $row = 0;
 221 for(@peak_calling) {
 222   my %hash = %{$_};
 223   $row++;
 224   #my $color = $row % 2?"#CCCCFF":"#FFFFFF";
 225   my $color = "#EEEEFF";
 226   print "<TR BGCOLOR=$color><TD><B>$hash{Name}</B></TD>";
 227   print "<TD>$hash{Caller}</TD><TD>$hash{Summary}</TD>\n";
 228   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://illumina-mac.stanford.edu/Tasks/$hash{Task}/$hash{outfile}>View in Genome Browser</A></TD><TD><A HREF=$hash{fasta}>FASTA</A></TD>\n";
 229   if(exists($hash{primer_design})) { print "<TD><A HREF=$hash{primer_design}>Validation Primers</A></TD>\n"; }
 230   print "</TR>\n";
 231 }
 232 }
 233 print "</TABLE>\n";
 234 print "</BODY></HTML>\n";
 235 sub SummarizeLibrary {
 236   my $xml_file = shift;
 237   my $xml = XMLin($xml_file, ForceArray => 1);
 238   my $summary = "";
 239   $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";
 240   for my $i (0..scalar(@{$xml->{Library}})-1) {
 241     my $lib = $xml->{Library}->[$i]->{Name};
 242     next unless(scalar(keys %libs) == 0 || exists($libs{$lib}));
 243     my $N = scalar(@{$xml->{Library}->[$i]->{Track}});
 244     my($date,$fc,$lane,$desc);
 245     my $num_reads = 0;
 246     my $num_lanes = 0;
 247     for my $t (0..$N-1) {
 248       my $filename = $xml->{Library}->[$i]->{Track}->[$t]->{Filename};
 249       $filename = `basename $filename`;
 250       $filename =~ /^(\d+)_(.+?)_s(\d+)_(.+?)_$lib.align/;
 251       ($date,$fc,$lane,$desc) = ($1,$2,$3,$4);
 252       $num_lanes += length($lane);
 253       $num_reads += $xml->{Library}->[$i]->{Track}->[$t]->{Count};
 254     }
 255     my $bgcolor;
 256     if($num_reads < 3000000) { $bgcolor = "FF3300"; }
 257     elsif($num_reads < 5000000) { $bgcolor = "FFCC33"; }
 258     elsif($num_reads < 10000000) { $bgcolor = "00CCFF"; }
 259     else { $bgcolor = "66FF66"; }
 260
 261     $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);
 262   }
 263   $summary .= "</TABLE>\n";
 264   return $summary;
 265 }
 266