bp_scripts/QA_454_report.rb

   1 #!/usr/bin/env ruby
   2
   3 # Copyright (C) 2011 Martin A. Hansen (mail@maasha.dk).
   4
   5 # This program is free software; you can redistribute it and/or
   6 # modify it under the terms of the GNU General Public License
   7 # as published by the Free Software Foundation; either version 2
   8 # of the License, or (at your option) any later version.
   9
  10 # This program is distributed in the hope that it will be useful,
  11 # but WITHOUT ANY WARRANTY; without even the implied warranty of
  12 # MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
  13 # GNU General Public License for more details.
  14
  15 # You should have received a copy of the GNU General Public License
  16 # along with this program; if not, write to the Free Software
  17 # Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA.
  18
  19 # http://www.gnu.org/copyleft/gpl.html
  20
  21 require 'pp'
  22 require 'tmpdir'
  23 require 'base64'
  24 require 'erb'
  25
  26 class Report
  27   def initialize(sff_file)
  28     @sff_file = sff_file
  29     tmpdir    = Dir.mktmpdir
  30
  31     @seq_analysis   = SeqAnalyze.new(@sff_file, tmpdir)
  32     @plots          = PlotData.new(@sff_file, tmpdir)
  33     @table_mid_join = MidTable.new(@sff_file, tmpdir)
  34     @table_freq     = FreqTable.new(@sff_file, tmpdir)
  35   end
  36
  37   # Support templating of member data.
  38   def get_binding
  39     binding
  40   end
  41
  42   private
  43
  44   class SeqAnalyze
  45     attr_reader :count, :min, :max, :mean, :bases, :gc, :hard, :soft
  46
  47     def initialize(sff_file, tmpdir)
  48       @sff_file  = sff_file
  49       @out1_file = File.join(tmpdir, "out1.txt")
  50       @out2_file = File.join(tmpdir, "out2.txt")
  51       @count     = 0
  52       @min       = 0
  53       @max       = 0
  54       @mean      = 0
  55       @bases     = 0
  56       @gc        = 0
  57       @hard      = 0
  58       @soft      = 0
  59
  60       bp_seq_analyze
  61       parse_analyze_vals
  62       parse_mean_vals
  63     end
  64
  65     private
  66
  67     def bp_seq_analyze
  68       STDERR.puts "Analyzing sequences ... "
  69       system(
  70         "read_sff -i #{@sff_file} |
  71          progress_meter |
  72          analyze_vals -k SEQ -o #{@out1_file} |
  73          analyze_seq |
  74          mean_vals -k 'GC%,HARD_MASK%,SOFT_MASK%' -o #{@out2_file} -x"
  75        )
  76       STDERR.puts "done.\n"
  77     end
  78
  79     def parse_analyze_vals
  80       File.open(@out1_file, "r") do |ios|
  81         while not ios.eof?
  82           line = ios.readline.chomp
  83
  84           case line
  85           when /COUNT\s+(\d+)/; then @count = $1
  86           when /MIN\s+(\d+)/;   then @min   = $1
  87           when /MAX\s+(\d+)/;   then @max   = $1
  88           when /MEAN\s+(\d+)/;  then @mean  = $1
  89           when /SUM\s+(\d+)/;   then @bases = $1
  90           end
  91         end
  92       end
  93     end
  94
  95     def parse_mean_vals
  96       File.open(@out2_file, "r") do |ios|
  97         while not ios.eof?
  98           line = ios.readline.chomp
  99
 100           case line
 101           when /GC%_MEAN: (.+)/;        then @gc   = $1
 102           when /HARD_MASK%_MEAN: (.+)/; then @hard = $1
 103           when /SOFT_MASK%_MEAN: (.+)/; then @soft = $1
 104           end
 105         end
 106       end
 107     end
 108   end
 109
 110   class PlotData
 111     attr_reader :lendist_unclipped, :lendist_clipped, :scores_unclipped, :scores_clipped, :mean_scores
 112
 113     def initialize(sff_file, tmpdir)
 114       @sff_file = sff_file
 115       @plot1    = File.join(tmpdir, "plot1.png")
 116       @plot2    = File.join(tmpdir, "plot2.png")
 117       @plot3    = File.join(tmpdir, "plot3.png")
 118       @plot4    = File.join(tmpdir, "plot4.png")
 119       @plot5    = File.join(tmpdir, "plot5.png")
 120
 121       bp_plot
 122
 123       @lendist_unclipped = png2base64(@plot1)
 124       @lendist_clipped   = png2base64(@plot3)
 125       @scores_unclipped  = png2base64(@plot2)
 126       @scores_clipped    = png2base64(@plot4)
 127       @mean_scores       = png2base64(@plot5)
 128     end
 129
 130     def bp_plot
 131       STDERR.puts "Creating plots ... "
 132       system(
 133         "read_sff -c -i #{@sff_file} |
 134          progress_meter |
 135          plot_distribution -k SEQ_LEN -T 'Length Distribution - unclipped' -t png -o #{@plot1} |
 136          plot_scores -T 'Mean Quality Scores - unclipped' -t png -o #{@plot2} |
 137          clip_seq |
 138          plot_distribution -k SEQ_LEN -T 'Length Distribution - clipped' -t png -o #{@plot3} |
 139          plot_scores -T 'Mean Quality Scores - clipped' -t png -o #{@plot4} |
 140          mean_scores |
 141          bin_vals -k SCORES_MEAN -b 5 |
 142          plot_histogram -s num -k SCORES_MEAN_BIN -T 'Mean score bins' -X 'Bins (size 5)' -Y 'Count' -t png -o #{@plot5} -x"
 143       )
 144       STDERR.puts "done.\n"
 145     end
 146
 147     private
 148
 149     def png2base64(file)
 150       png = ""
 151
 152       File.open(file, "r") do |ios|
 153         png = ios.read
 154       end
 155
 156       "data:image/png;base64," + Base64.encode64(png)
 157     end
 158   end
 159
 160   class MidTable
 161     def initialize(sff_file, tmpdir)
 162       @sff_file  = sff_file
 163       @mid1_file = File.join(tmpdir, "mid1.tab")
 164       @mid2_file = File.join(tmpdir, "mid2.tab")
 165       @mid3_file = File.join(tmpdir, "mid3.tab")
 166       @mid4_file = File.join(tmpdir, "mid_join.tab")
 167
 168       bp_find_mids
 169       bp_merge_mid_tables
 170     end
 171
 172     def each
 173       File.open(@mid4_file, "r") do |ios|
 174         while not ios.eof? do
 175           fields = ios.readline.chomp.split("\t")
 176           yield MidRow.new(fields[0], fields[1], fields[2], fields[3], fields[4])
 177         end
 178       end
 179     end
 180
 181     private
 182
 183     def bp_find_mids
 184       STDERR.puts "Finding barcodes in raw sequences ... "
 185       system(
 186         "read_sff -i #{@sff_file} |
 187          find_barcodes -p 4 -gr |
 188          count_vals -k BARCODE_NAME |
 189          uniq_vals -k BARCODE_NAME |
 190          write_tab -c -k BARCODE_NAME,BARCODE,BARCODE_NAME_COUNT -o #{@mid1_file} -x"
 191       )
 192       STDERR.puts "done.\n"
 193       STDERR.puts "Finding barcodes in sequences >= 250 ... "
 194       system(
 195         "read_sff -i #{@sff_file} |
 196          grab -e 'SEQ_LEN >= 250' |
 197          find_barcodes -p 4 -gr |
 198          count_vals -k BARCODE_NAME |
 199          uniq_vals -k BARCODE_NAME |
 200          write_tab -c -k BARCODE_NAME,BARCODE,BARCODE_NAME_COUNT -o #{@mid2_file} -x"
 201       )
 202       STDERR.puts "done.\n"
 203       STDERR.puts "Finding barcodes in sequences >= 250 with mean score >= 20 ... "
 204       system(
 205         "read_sff -i #{@sff_file} |
 206          mean_scores |
 207          grab -e 'SEQ_LEN >= 250' |
 208          grab -e 'SCORES_MEAN >= 20' |
 209          find_barcodes -p 4 -gr |
 210          count_vals -k BARCODE_NAME |
 211          uniq_vals -k BARCODE_NAME |
 212          write_tab -c -k BARCODE_NAME,BARCODE,BARCODE_NAME_COUNT -o #{@mid3_file} -x"
 213       )
 214       STDERR.puts "done.\n"
 215     end
 216
 217     def bp_merge_mid_tables
 218       STDERR.print "Joining MID tables ... "
 219       system(
 220         "read_tab -i #{@mid1_file} |
 221          rename_keys -k BARCODE_NAME,A |
 222          rename_keys -k BARCODE_NAME_COUNT,TOTAL |
 223          read_tab -i #{@mid2_file} |
 224          rename_keys -k BARCODE_NAME,B |
 225          rename_keys -k BARCODE_NAME_COUNT,L250 |
 226          merge_records -k A,B |
 227          read_tab -i #{@mid3_file} |
 228          rename_keys -k BARCODE_NAME,C |
 229          rename_keys -k BARCODE_NAME_COUNT,L250_S20 |
 230          merge_records -k A,C |
 231          rename_keys -k A,BARCODE_NAME |
 232          sort_records -k BARCODE_NAME |
 233          write_tab -ck BARCODE_NAME,BARCODE,TOTAL,L250,L250_S20 -o #{@mid4_file} -x"
 234       )
 235       STDERR.puts "done.\n"
 236     end
 237
 238     class MidRow
 239       attr_reader :mid_num, :mid_seq, :total, :l250, :l250_s20
 240
 241       def initialize(mid_num, mid_seq, total, l250, l250_s20)
 242         @mid_num  = mid_num
 243         @mid_seq  = mid_seq
 244         @total    = total
 245         @l250     = l250
 246         @l250_s20 = l250_s20
 247       end
 248     end
 249   end
 250
 251   class FreqTable
 252     def initialize(sff_file, tmpdir)
 253       @sff_file = sff_file
 254       @tab_file = File.join(tmpdir, "freq.tab")
 255
 256       bp_frequency_table
 257     end
 258
 259     def each
 260       File.open(@tab_file, "r") do |ios|
 261         while not ios.eof? do
 262           fields = ios.readline.chomp.split("\t")
 263           yield FreqRow.new(fields[0], fields[1], fields[2], fields[3], fields[4])
 264         end
 265       end
 266     end
 267
 268     private
 269
 270     def bp_frequency_table
 271       STDERR.puts "Creating residue frequency table ... "
 272       system(
 273         "read_sff -i #{@sff_file} |
 274          progress_meter |
 275          extract_seq -l 50 |
 276          uppercase_seq |
 277          create_weight_matrix -p |
 278          flip_tab |
 279          write_tab -o #{@tab_file} -x"
 280       )
 281       STDERR.puts "done.\n"
 282     end
 283
 284     class FreqRow
 285       attr_reader :res_A, :res_C, :res_G, :res_N, :res_T
 286
 287       def initialize(res_A, res_C, res_G, res_N, res_T)
 288         @res_A = res_A
 289         @res_C = res_C
 290         @res_G = res_G
 291         @res_N = res_N
 292         @res_T = res_T
 293       end
 294     end
 295   end
 296 end
 297
 298 template = %{
 299   <html>
 300     <head>
 301       <title>QA 454 Report</title>
 302     </head>
 303     <body>
 304       <h1>QA 454 Report</h1>
 305       <p>Date: #{Time.now}</p>
 306       <p>File: <%= @sff_file %></p>
 307       <h2>Sequence analysis</h2>
 308       <ul>
 309         <li>Number of sequences in the file: <%= @seq_analysis.count %></li>
 310         <li>Minimum sequence length found: <%= @seq_analysis.min %></li>
 311         <li>Maximum sequence length found: <%= @seq_analysis.max %></li>
 312         <li>Mean sequence length found: <%= @seq_analysis.mean %></li>
 313         <li>Total number of bases in the file: <%= @seq_analysis.bases %></li>
 314         <li>Mean GC% content: <%= @seq_analysis.gc %></li>
 315         <li>Mean of hard masked sequence (i.e. % of N's): <%= @seq_analysis.hard %></li>
 316         <li>Mean of soft masked sequence (i.e. % lowercase residues = clipped sequence): <%= @seq_analysis.soft %></li>
 317       </ul>
 318       <h2>Sequence length distribution</h2>
 319       <p>The length distribution of unclipped reads:</p>
 320       <p><img alt="plot_lendist_unclipped" src="<%= @plots.lendist_unclipped %>" width="600" /></p>
 321       <p>The length distribution of clipped reads:</p>
 322       <p><img alt="plot_lendist_clipped" src="<%= @plots.lendist_clipped %>" width="600" /></p>
 323       <h2>Quality score means</h2>
 324       <p>The mean scores of the unclipped sequences:</p>
 325       <p><img alt="plot_scores_unclipped" src="<%= @plots.scores_unclipped %>" width="600" /></p>
 326       <p>The mean scores of the clipped sequences:</p>
 327       <p><img alt="plot_scores_clipped" src="<%= @plots.scores_clipped %>" width="600" /></p>
 328       <p>Histogram of bins with mean quality scores:</p>
 329       <p><img alt="plot_mean_scores" src="<%= @plots.mean_scores %>" width="600" /></p>
 330       <h2>MID tag analysis</h2>
 331       <p>The below table contains the identified MID tags and the number of times they were found:<p>
 332       <ul>
 333       <li>BARCODE_NAME is the MID tag identifier.</li>
 334       <li>BARCODE is the sequence of the MID tag.</li>
 335       <li>TOTAL is the number of times this MID tag was found.</li>
 336       <li>L250 is the a subset count of TOTAL af sequences longer than 250 bases</li>
 337       <li>L250_S20 is a subset count of L250 af sequences with a mean score above 20</li>
 338       </ul>
 339       <table>
 340       <% @table_mid_join.each do |row| %>
 341         <tr>
 342         <td><%= row.mid_num %></td>
 343         <td><%= row.mid_seq %></td>
 344         <td><%= row.total %></td>
 345         <td><%= row.l250 %></td>
 346         <td><%= row.l250_s20 %></td>
 347         </tr>
 348       <% end %>
 349       </table>
 350       <h2>Residue frequency analysis</h2>
 351       <p>The below table contains the residue frequency (in percent) of the first 50 bases:</p>
 352       <table>
 353       <% @table_freq.each do |row| %>
 354         <tr>
 355         <td><%= row.res_A %></td>
 356         <td><%= row.res_T %></td>
 357         <td><%= row.res_C %></td>
 358         <td><%= row.res_G %></td>
 359         <td><%= row.res_N %></td>
 360         </tr>
 361       <% end %>
 362       </table>
 363     </body>
 364   </html>
 365 }.gsub(/^\s+/, '')
 366
 367 html = ERB.new(template)
 368
 369 ARGV.each do |file|
 370   report = Report.new(file)
 371
 372   html.run(report.get_binding)
 373 end
 374
 375 __END__