htswanalysis/scripts/Summarize_Methylseq.pm

   1 #!/usr/bin/perl -w
   2
   3 use strict;
   4 use warnings;
   5
   6 my $name = shift;
   7 if(!defined($name)) { print STDERR "Error: No name provided to summarize methylseq\n"; exit(1); }
   8 my $task = shift;
   9 if(!defined($task)) { print STDERR "Error: No task provided to summarize methylseq\n"; exit(1); }
  10 my $tasks_root = shift;
  11 if(!defined($tasks_root)) { print STDERR "Error: No root task directory provided to summarize methylseq\n"; exit(1); }
  12 my $msp1_lib = shift;
  13 if(!defined($msp1_lib)) { print STDERR "Error: No msp1 library provided to summarize methylseq\n"; exit(1); }
  14 my $hpa2_lib = shift;
  15 if(!defined($hpa2_lib)) { print STDERR "Error: No hpa2 library provided to summarize methylseq\n"; exit(1); }
  16 my $genome = shift;
  17
  18 my $methylseq_bed = "MethylseqCalls.bed";
  19 my $methylseq_calls = "MethylseqCalls.tab";
  20 my $summary = "methylseq_summary";
  21
  22 if(! -e $tasks_root."/".$task."/".$summary) {
  23   printf "<TR><TD BGCOLOR=#FFCCCC>Processing...</TD></TR>\n";
  24   exit;
  25 }
  26
  27 open(FILE,$tasks_root."/".$task."/".$summary);
  28
  29 my $tag_count = <FILE>;
  30
  31 if(!defined($tag_count)) {
  32   printf "<TR><TD BGCOLOR=#FFCCCC>Processing...</TD></TR>\n";
  33   exit;
  34 }
  35
  36 chomp $tag_count;
  37 $tag_count =~ /^(\d+) MSP1 tags, (\d+) HPA2 tags$/;
  38 my($msp_count,$hpa_count) = ($1,$2);
  39
  40 my $assayed_count = <FILE>; chomp $assayed_count;
  41    $assayed_count =~ /^(\d+) of (\d+) regions assayed/;
  42    my($n_assayed,$n_regions) = ($1,$2);
  43
  44 my $methylated_count = <FILE>; chomp $methylated_count;
  45    $methylated_count =~ /^(\d+) of \d+ regions methylated/;
  46    my $n_methylated = $1;
  47
  48 my $error_count = <FILE>; chomp $error_count;
  49    $error_count =~ /^(\d+) of \d+ regions with error/;
  50    my $n_error = $1;
  51
  52
  53 my $pcnt_assayed = $n_assayed / $n_regions;
  54 my $pcnt_methylated = $n_methylated / $n_assayed;
  55 my $pcnt_error = $n_error / $n_regions;
  56
  57 my $assayed_color;
  58 if($pcnt_assayed < 0.7) { $assayed_color = "FF6600"; }
  59 elsif($pcnt_assayed < 0.9) { $assayed_color = "4444FF"; }
  60 else {$assayed_color= "66FF66"; }
  61
  62 my $methylated_color;
  63 if($pcnt_methylated < 0.3) { $methylated_color = "FF6600"; }
  64 elsif($pcnt_methylated < 0.6) { $methylated_color = "4444FF"; }
  65 else {$methylated_color= "66FF66"; }
  66
  67 my $error_color;
  68 if($pcnt_error > 0.05) { $error_color = "FF6600"; }
  69 elsif($pcnt_error > 0.01) { $error_color = "4444FF"; }
  70 else {$error_color= "66FF66"; }
  71
  72 print "<TR>";
  73 printf "<TD>$name</TD><TD>$msp1_lib</TD><TD>$hpa2_lib</TD><TD BGCOLOR=#$assayed_color>$n_assayed (%.4f)</TD><TD BGCOLOR=#$methylated_color>$n_methylated (%0.4f)</TD><TD BGCOLOR=$error_color>$n_error (%0.4f)</TD>",$pcnt_assayed,$pcnt_methylated,$pcnt_error;
  74 print "<TD><A HREF=../../Tasks/$task/$methylseq_calls>Results</A><BR>\n";
  75 print "<A HREF=http://genome.ucsc.edu/cgi-bin/hgTracks?db=$genome&position=chr1:1000-2000&hgt.customText=http://illumina-mac.stanford.edu/Tasks/$task/$methylseq_bed TARGET=\"_blank\">(Genome Browser)</A></BR>\n";
  76 print "</TD>";
  77 print "</TR>";