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 $msp1_lib = shift;
  11 if(!defined($msp1_lib)) { print STDERR "Error: No msp1 library provided to summarize methylseq\n"; exit(1); }
  12 my $hpa2_lib = shift;
  13 if(!defined($hpa2_lib)) { print STDERR "Error: No hpa2 library provided to summarize methylseq\n"; exit(1); }
  14 my $genome = shift;
  15
  16 my $assayed_bed = "Assayed.bed";
  17 my $methylated_bed = "Methylated.bed";
  18 my $methylseq_calls = "MethylseqCalls.tab";
  19 my $summary = "methylseq_summary";
  20 if(!defined($summary)) { return; }
  21
  22 open(FILE,$summary);
  23
  24 my $tag_count = <FILE>; chomp $tag_count;
  25 $tag_count =~ /^(\d+) MSP1 tags, (\d+) HPA2 tags$/;
  26 my($msp_count,$hpa_count) = ($1,$2);
  27
  28 my $assayed_count = <FILE>; chomp $assayed_count;
  29    $assayed_count =~ /^(\d+) of (\d+) regions assayed/;
  30    my($n_assayed,$n_regions) = ($1,$2);
  31
  32 my $methylated_count = <FILE>; chomp $methylated_count;
  33    $methylated_count =~ /^(\d+) of \d+ regions methylated/;
  34    my $n_methylated = $1;
  35
  36 my $error_count = <FILE>; chomp $error_count;
  37    $error_count =~ /^(\d+) of \d+ regions with error/;
  38    my $n_error = $1;
  39
  40 my $pcnt_assayed = $n_assayed / $n_regions;
  41 my $pcnt_methylated = $n_methylated / $n_assayed;
  42 my $pcnt_error = $n_error / $n_regions;
  43
  44 my $assayed_color;
  45 if($pcnt_assayed < 0.7) { $assayed_color = "FF6600"; }
  46 elsif($pcnt_assayed < 0.9) { $assayed_color = "4444FF"; }
  47 else {$assayed_color= "66FF66"; }
  48
  49 my $methylated_color;
  50 if($pcnt_methylated < 0.3) { $methylated_color = "FF6600"; }
  51 elsif($pcnt_methylated < 0.6) { $methylated_color = "4444FF"; }
  52 else {$methylated_color= "66FF66"; }
  53
  54 my $error_color;
  55 if($pcnt_error > 0.05) { $error_color = "FF6600"; }
  56 elsif($pcnt_error > 0.01) { $error_color = "4444FF"; }
  57 else {$error_color= "66FF66"; }
  58
  59 print "<TR>";
  60 printf "<TD>$msp1_lib</TD><TD>$hpa2_lib</TD><TD>$n_assayed (%.4f)</TD><TD>$n_methylated (%0.4f)</TD><TD>$n_error (%0.4f)</TD>",$pcnt_assayed,$pcnt_methylated,$pcnt_error;
  61 print "<TD><A HREF=../../Tasks/$task/$methylseq_calls>Calls</A><BR>\n";
  62 print "<A HREF=../../Tasks/$task/Assayed.bed>Assayed Regions</A>";
  63 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/Assayed.bed TARGET=\"_blank\">(Genome Browser)</A></BR>\n";
  64 print "<A HREF=../../Tasks/$task/Methylated.bed>Methylated Regions</A> ";
  65 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/Methylated.bed TARGET=\"_blank\">(Genome Browser)</A></TD>\n";
  66 print "</TR>";