code_ruby/lib/maasha/sam.rb

   1 # as published by the Free Software Foundation; either version 2
   2 # of the License, or (at your option) any later version.
   3
   4 # This program is distributed in the hope that it will be useful,
   5 # but WITHOUT ANY WARRANTY; without even the implied warranty of
   6 # MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
   7 # GNU General Public License for more details.
   8
   9 # You should have received a copy of the GNU General Public License
  10 # along with this program; if not, write to the Free Software
  11 # Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA.
  12
  13 # http://www.gnu.org/copyleft/gpl.html
  14
  15 # >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>><<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
  16
  17 # This software is part of the Biopieces framework (www.biopieces.org).
  18
  19 # >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>><<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
  20
  21 # SAM format version v1.4-r962 - April 17, 2011
  22 #
  23 # http://samtools.sourceforge.net/SAM1.pdf
  24
  25 require 'pp'
  26 require 'maasha/filesys'
  27 require 'maasha/seq'
  28
  29 # Error class for all exceptions to do with Genbank.
  30 class SamError < StandardError; end
  31
  32 REGEX_HEADER  = Regexp.new(/^@[A-Za-z][A-Za-z](\t[A-Za-z][A-Za-z0-9]:[ -~]+)+$/)
  33 REGEX_COMMENT = Regexp.new(/^@CO\t.*/)
  34
  35 # Class to parse and write SAM files.
  36 class Sam < Filesys
  37   attr_accessor :io
  38
  39   # Method to initialize a Sam object.
  40   def initialize(io = nil)
  41     @io          = io
  42     @header_hash = {}
  43   end
  44
  45   # Method to parse the header of a SAM file.
  46   # Each header line should match:
  47   # /^@[A-Za-z][A-Za-z](\t[A-Za-z][A-Za-z0-9]:[ -~]+)+$/ or /^@CO\t.*/.
  48   # Tags containing lowercase letters are reserved for end users.
  49   def header
  50     @io.each_line do |line|
  51       if line =~ /^@([A-Za-z][A-Za-z])/
  52         line.chomp!
  53
  54         tag = $1
  55
  56         case tag
  57         when 'HD' then parse_header(line)
  58         when 'SQ' then parse_sequence(line)
  59         when 'RG' then parse_read_group(line)
  60         when 'PG' then parse_program(line)
  61         when 'CO' then parse_comment(line)
  62         else
  63           raise SamError, "Unknown header tag: #{tag}"
  64         end
  65       else
  66         @io.rewind
  67         break
  68       end
  69     end
  70
  71     return @header_hash.empty? ? nil : @header_hash
  72   end
  73
  74   def each
  75     @io.each_line do |line|
  76       unless line[0] == '@'
  77         entry = parse_alignment(line.chomp)
  78
  79         yield entry if block_given?
  80       end
  81     end
  82   end
  83
  84   private
  85
  86   # Method to subparse header lines.
  87   def parse_header(line)
  88     hash   = {}
  89     fields = line.split("\t")
  90
  91     if fields[1] =~ /^VN:([0-9]+\.[0-9]+)$/
  92       hash[:VN] = $1.to_f
  93     else
  94       raise SamError, "Bad version number: #{fields[1]}"
  95     end
  96
  97     if fields.size > 2
  98       if fields[2] =~ /^SO:(unknown|unsorted|queryname|coordinate)$/
  99         hash[:SO] = $1
 100       else
 101         raise SamError, "Bad sort order: #{fields[2]}"
 102       end
 103     end
 104
 105     @header_hash[:HD] = hash
 106   end
 107
 108   # Method to subparse sequence lines.
 109   def parse_sequence(line)
 110     @header_hash[:SQ] = Hash.new unless @header_hash[:SQ].is_a? Hash
 111     hash = {}
 112
 113     fields = line.split("\t")
 114
 115     if fields[1] =~ /^SN:([!-)+-<>-~][!-~]*)$/
 116       seq_name = $1.to_sym
 117     else
 118       raise SamError, "Bad sequence name: #{fields[1]}"
 119     end
 120
 121     if fields[2] =~ /^LN:(\d+)$/
 122       hash[:LN] = $1.to_i
 123     else
 124       raise SamError, "Bad sequence length: #{fields[2]}"
 125     end
 126
 127     (3 ... fields.size).each do |i|
 128       if fields[i] =~ /^(AS|M5|SP|UR):([ -~]+)$/
 129         hash[$1.to_sym] = $2
 130       else
 131         raise SamError, "Bad sequence tag: #{fields[i]}"
 132       end
 133     end
 134
 135     @header_hash[:SQ][:SN] = Hash.new unless @header_hash[:SQ][:SN].is_a? Hash
 136
 137     if @header_hash[:SQ][:SN].has_key? seq_name
 138       raise SamError, "Non-unique sequence name: #{seq_name}"
 139     else
 140       @header_hash[:SQ][:SN][seq_name] = hash
 141     end
 142   end
 143
 144   # Method to subparse read group lines.
 145   def parse_read_group(line)
 146     @header_hash[:RG] = Hash.new unless @header_hash[:RG].is_a? Hash
 147     hash = {}
 148
 149     fields = line.split("\t")
 150
 151     if fields[1] =~ /^ID:([ -~]+)$/
 152       id = $1.to_sym
 153     else
 154       raise SamError, "Bad read group identifier: #{fields[1]}"
 155     end
 156
 157     (2 ... fields.size).each do |i|
 158       if fields[i] =~ /^(CN|DS|DT|FO|KS|LB|PG|PI|PL|PU|SM):([ -~]+)$/
 159         hash[$1.to_sym] = $2
 160       else
 161         raise SamError, "Bad read group tag: #{fields[i]}"
 162       end
 163     end
 164
 165     if hash.has_key? :FO
 166       unless hash[:FO] =~ /^\*|[ACMGRSVTWYHKDBN]+$/
 167         raise SamError, "Bad flow order: #{hash[:FO]}"
 168       end
 169     end
 170
 171     if hash.has_key? :PL
 172       unless hash[:PL] =~ /^(CAPILLARY|LS454|ILLUMINA|SOLID|HELICOS|IONTORRENT|PACBIO)$/
 173         raise SamError, "Bad platform: #{hash[:PL]}"
 174       end
 175     end
 176
 177     @header_hash[:RG][:ID] = Hash.new unless @header_hash[:RG][:ID].is_a? Hash
 178
 179     if @header_hash[:RG][:ID].has_key? id
 180       raise SamError, "Non-unique read group identifier: #{id}"
 181     else
 182       @header_hash[:RG][:ID][id] = hash
 183     end
 184   end
 185
 186   # Method to subparse program lines.
 187   def parse_program(line)
 188     @header_hash[:PG] = Hash.new unless @header_hash[:PG].is_a? Hash
 189     hash = {}
 190
 191     fields = line.split("\t")
 192
 193     if fields[1] =~ /^ID:([ -~]+)$/
 194       id = $1.to_sym
 195     else
 196       raise SamError, "Bad program record identifier: #{fields[1]}"
 197     end
 198
 199     (2 ... fields.size).each do |i|
 200       if fields[i] =~ /^(PN|CL|PP|VN):([ -~]+)$/
 201         hash[$1.to_sym] = $2
 202       else
 203         raise SamError, "Bad program record tag: #{fields[i]}"
 204       end
 205     end
 206
 207     @header_hash[:PG][:ID] = Hash.new unless @header_hash[:PG][:ID].is_a? Hash
 208
 209     if @header_hash[:PG][:ID].has_key? id
 210       raise SamError, "Non-unique program record identifier: #{id}"
 211     else
 212       @header_hash[:PG][:ID][id] = hash
 213     end
 214   end
 215
 216   # Method to subparse comment lines.
 217   def parse_comment(line)
 218     @header_hash[:CO] = Array.new unless @header_hash[:CO].is_a? Array
 219
 220     if line =~ /^@CO\t(.+)/
 221       @header_hash[:CO] << $1
 222     else
 223       raise SamError, "Bad comment line: #{line}"
 224     end
 225   end
 226
 227   # Method to subparse alignment lines.
 228   def parse_alignment(line)
 229     fields = line.split("\t")
 230
 231     raise SamError, "Bad number of fields: #{fields.size} < 11" if fields.size < 11
 232
 233     qname = fields[0]
 234     flag  = fields[1].to_i
 235     rname = fields[2]
 236     pos   = fields[3].to_i
 237     mapq  = fields[4].to_i
 238     cigar = fields[5]
 239     rnext = fields[6]
 240     pnext = fields[7].to_i
 241     tlen  = fields[8].to_i
 242     seq   = fields[9]
 243     qual  = fields[10]
 244
 245     raise SamError, "Bad qname: #{qname}" unless qname =~ /^[!-?A-~]{1,255}$/
 246     raise SamError, "Bad flag: #{flag}"   unless (0 .. 2**16 - 1).include? flag
 247     raise SamError, "Bad rname: #{rname}" unless rname =~ /^(\*|[!-()+-<>-~][!-~]*)$/
 248     raise SamError, "Bad pos: #{pos}"     unless (0 .. 2**29 - 1).include? pos
 249     raise SamError, "Bad mapq: #{mapq}"   unless (0 .. 2**8 - 1).include? mapq
 250     raise SamError, "Bad cigar: #{cigar}" unless cigar =~ /^(\*|([0-9]+[MIDNSHPX=])+)$/
 251     raise SamError, "Bad rnext: #{rnext}" unless rnext =~ /^(\*|=|[!-()+-<>-~][!-~]*)$/
 252     raise SamError, "Bad pnext: #{pnext}" unless (0 .. 2**29 - 1).include? pnext
 253     raise SamError, "Bad tlen: #{tlen}"   unless (-2**29 + 1 .. 2**29 - 1).include? tlen
 254     raise SamError, "Bad seq: #{seq}"     unless seq  =~ /^(\*|[A-Za-z=.]+)$/
 255     raise SamError, "Bad qual: #{qual}"   unless qual =~ /^[!-~]+$/
 256
 257     entry = {}
 258     entry[:QNAME] = qname
 259     entry[:FLAG]  = flag
 260     entry[:RNAME] = rname
 261     entry[:POS]   = pos
 262     entry[:MAPQ]  = mapq
 263     entry[:CIGAR] = cigar
 264     entry[:RNEXT] = rnext
 265     entry[:PNEXT] = pnext
 266     entry[:TLEN]  = tlen
 267     entry[:SEQ]   = seq
 268     entry[:QUAL]  = qual
 269
 270     entry
 271   end
 272 end
 273
 274
 275 # >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>><<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
 276
 277
 278 __END__
 279