bp_bin/order_pairs2

   1 #!/usr/bin/env ruby
   2
   3 # Copyright (C) 2007-2012 Martin A. Hansen.
   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 # >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>><<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
  22
  23 # This program is part of the Biopieces framework (www.biopieces.org).
  24
  25 # >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> DESCRIPTION <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
  26
  27 # Order records with pair end sequence data.
  28
  29 # >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>><<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
  30
  31 require 'maasha/biopieces'
  32 require 'gdbm'
  33 require 'pp'
  34
  35 ILLUMINA15 = 15
  36 ILLUMINA18 = 18
  37
  38 options = Biopieces.options_parse(ARGV)
  39
  40 tmpdir  = Biopieces.mktmpdir
  41 #tmpdir  = "Sletmig"
  42 tmpfile = File.join(tmpdir, "data.gdbm")
  43
  44 db = GDBM.new(tmpfile)
  45
  46 Biopieces.open(options[:stream_in], options[:stream_out]) do |input, output|
  47   input.each_record do |record|
  48     if record[:SEQ_NAME]
  49       db[record[:SEQ_NAME]] = Marshal.dump(record)
  50     else
  51       output.puts record
  52     end
  53   end
  54
  55   skip = {}
  56
  57   db.keys.each do |seq_name|
  58     next if skip[seq_name]
  59
  60     case seq_name
  61     when /^(.+)\/(\d)$/   # Illumina 1.5
  62       type = ILLUMINA15
  63       name = $1
  64       pair = $2.to_i
  65     when /^(.+) (\d):/    # Illumina 1.8
  66       type = ILLUMINA18
  67       name = $1.to_sym
  68       pair = $2.to_i
  69     else
  70       $stderr.puts "WARNING: Unmatched sequence name: #{record[:SEQ_NAME]}"
  71     end
  72
  73     pair2 = (pair == 1) ? 2 : 1
  74
  75     if type == ILLUMINA15
  76       seq_name2 = "#{name}/#{pair2}"
  77     else
  78       seq_name2 = seq_name.sub(/^(.+) \d:/, "#{$1} #{pair2}:")
  79     end
  80
  81     record1 = Marshal.load(db[seq_name])
  82     record2 = db[seq_name2]
  83     record2 = Marshal.load(record2) unless record2.nil?
  84
  85     if record2.nil?
  86       record1[:ORDER] = "orphan#{pair}"
  87       output.puts record1
  88     else
  89       record1[:ORDER] = "paired"
  90       record2[:ORDER] = "paired"
  91
  92       output.puts record1
  93       output.puts record2
  94
  95       skip[seq_name2] = true
  96     end
  97   end
  98 end
  99
 100 db.close
 101
 102 FileUtils.remove_entry_secure(tmpfile, true)
 103
 104 # >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>><<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
 105
 106
 107 __END__