code_ruby/lib/maasha/sam.rb

   1 # Copyright (C) 2007-2011 Martin A. Hansen.
   2
   3 # This program is free software; you can redistribute it and/or
   4 # modify it under the terms of the GNU General Public License
   5 # as published by the Free Software Foundation; either version 2
   6 # of the License, or (at your option) any later version.
   7
   8 # This program is distributed in the hope that it will be useful,
   9 # but WITHOUT ANY WARRANTY; without even the implied warranty of
  10 # MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
  11 # GNU General Public License for more details.
  12
  13 # You should have received a copy of the GNU General Public License
  14 # along with this program; if not, write to the Free Software
  15 # Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA.
  16
  17 # http://www.gnu.org/copyleft/gpl.html
  18
  19 # >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>><<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
  20
  21 # This software is part of the Biopieces framework (www.biopieces.org).
  22
  23 # >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>><<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
  24
  25 # SAM format version v1.4-r962 - April 17, 2011
  26 #
  27 # http://samtools.sourceforge.net/SAM1.pdf
  28
  29 require 'pp'
  30 require 'maasha/filesys'
  31 require 'maasha/seq'
  32 require 'maasha/cigar'
  33
  34 REGEX_HEADER  = Regexp.new(/^@[A-Za-z][A-Za-z](\t[A-Za-z][A-Za-z0-9]:[ -~]+)+$/)
  35 REGEX_COMMENT = Regexp.new(/^@CO\t.*/)
  36
  37 FLAG_MULTI               = 0x1   # Template having multiple fragments in sequencing
  38 FLAG_ALIGNED             = 0x2   # Each fragment properly aligned according to the aligner
  39 FLAG_UNMAPPED            = 0x4   # Fragment unmapped
  40 FLAG_NEXT_UNMAPPED       = 0x8   # Next fragment in the template unmapped
  41 FLAG_REVCOMP             = 0x10  # SEQ being reverse complemented
  42 FLAG_NEXT_REVCOMP        = 0x20  # SEQ of the next fragment in the template being reversed
  43 FLAG_FIRST               = 0x40  # The first fragment in the template
  44 FLAG_LAST                = 0x80  # The last fragment in the template
  45 FLAG_SECONDARY_ALIGNMENT = 0x100 # Secondary alignment
  46 FLAG_QUALITY_FAIL        = 0x200 # Not passing quality controls
  47 FLAG_DUPLICATES          = 0x400 # PCR or optical duplicate
  48
  49 # Error class for all exceptions to do with Genbank.
  50 class SamError < StandardError; end
  51
  52 # Class to parse and write SAM files.
  53 class Sam < Filesys
  54   attr_accessor :io, :header
  55
  56   # Class method to convert a SAM entry
  57   # to a Biopiece record.
  58   def self.to_bp(sam)
  59     bp = {}
  60
  61     bp[:REC_TYPE] = 'SAM'
  62     bp[:Q_ID]     = sam[:QNAME]
  63     bp[:STRAND]   = sam[:FLAG].revcomp? ? '-' : '+'
  64     bp[:S_ID]     = sam[:RNAME]
  65     bp[:S_BEG]    = sam[:POS]
  66     bp[:MAPQ]     = sam[:MAPQ]
  67     bp[:CIGAR]    = sam[:CIGAR].to_s
  68
  69     unless sam[:RNEXT] == '*'
  70       bp[:Q_ID2]  = sam[:RNEXT]
  71       bp[:S_BEG2] = sam[:PNEXT]
  72       bp[:TLEN]   = sam[:TLEN]
  73     end
  74
  75     bp[:SEQ]      = sam[:SEQ].seq
  76
  77     unless sam[:SEQ].qual.nil?
  78       bp[:SCORES] = sam[:SEQ].convert_phred2illumina!.qual
  79     end
  80
  81     if sam.has_key? :NM and sam[:NM].to_i > 0
  82       bp[:ALIGN] = self.align_descriptors(sam)
  83     end
  84
  85     bp
  86   end
  87
  88   # Class method to convert a Biopiece record
  89   # into a SAM entry.
  90   def self.to_sam(bp)
  91     "FISK" # FIXME
  92   end
  93
  94   # Create alignment descriptors according to the KISS
  95   # format description:
  96   # http://code.google.com/p/biopieces/wiki/KissFormat
  97   def self.align_descriptors(sam)
  98     offset = 0
  99     align  = []
 100
 101     # Insertions
 102     sam[:CIGAR].each do |len, op|
 103       if op == 'I'
 104         (0 ... len).each_with_index do |i|
 105           nt = sam[:SEQ].seq[offset + i]
 106
 107           align << [offset + i, "->#{nt}"]
 108         end
 109       end
 110
 111       offset += len
 112     end
 113
 114     offset    = 0
 115     deletions = 0
 116
 117     sam[:MD].scan(/\d+|\^[A-Z]+|[A-Z]+/).each do |m|
 118       if m =~ /\d+/      # Matches
 119         offset += m.to_i
 120       elsif m[0] == '^'  # Deletions
 121         m.each_char do |nt|
 122           unless nt == '^'
 123             align << [offset, "#{nt}>-"]
 124             deletions += 1
 125             offset    += 1
 126           end
 127         end
 128       else               # Mismatches
 129         m.each_char do |nt|
 130           nt2 = sam[:SEQ].seq[offset - deletions]
 131
 132           align << [offset, "#{nt}>#{nt2}"]
 133
 134           offset += 1
 135         end
 136       end
 137     end
 138
 139     align.sort_by { |a| a.first }.map { |k, v| "#{k}:#{v}" }.join(",")
 140   end
 141
 142   # Method to initialize a Sam object.
 143   def initialize(io = nil)
 144     @io     = io
 145     @header = {}
 146
 147     parse_header
 148   end
 149
 150   def each
 151     @io.each_line do |line|
 152       unless line[0] == '@'
 153         entry = parse_alignment(line.chomp)
 154
 155         yield entry if block_given?
 156       end
 157     end
 158   end
 159
 160   private
 161
 162   # Method to parse the header section of a SAM file.
 163   # Each header line should match:
 164   # /^@[A-Za-z][A-Za-z](\t[A-Za-z][A-Za-z0-9]:[ -~]+)+$/ or /^@CO\t.*/.
 165   # Tags containing lowercase letters are reserved for end users.
 166   def parse_header
 167     @io.each_line do |line|
 168       if line =~ /^@([A-Za-z][A-Za-z])/
 169         line.chomp!
 170
 171         tag = $1
 172
 173         case tag
 174         when 'HD' then subparse_header(line)
 175         when 'SQ' then subparse_sequence(line)
 176         when 'RG' then subparse_read_group(line)
 177         when 'PG' then subparse_program(line)
 178         when 'CO' then subparse_comment(line)
 179         else
 180           raise SamError, "Unknown header tag: #{tag}"
 181         end
 182       else
 183         @io.rewind
 184         break
 185       end
 186     end
 187
 188     return @header.empty? ? nil : @header
 189   end
 190
 191   # Method to subparse header lines.
 192   def subparse_header(line)
 193     hash   = {}
 194     fields = line.split("\t")
 195
 196     if fields[1] =~ /^VN:([0-9]+\.[0-9]+)$/
 197       hash[:VN] = $1.to_f
 198     else
 199       raise SamError, "Bad version number: #{fields[1]}"
 200     end
 201
 202     if fields.size > 2
 203       if fields[2] =~ /^SO:(unknown|unsorted|queryname|coordinate)$/
 204         hash[:SO] = $1
 205       else
 206         raise SamError, "Bad sort order: #{fields[2]}"
 207       end
 208     end
 209
 210     @header[:HD] = hash
 211   end
 212
 213   # Method to subparse sequence lines.
 214   def subparse_sequence(line)
 215     @header[:SQ] = Hash.new unless @header[:SQ].is_a? Hash
 216     hash = {}
 217
 218     fields = line.split("\t")
 219
 220     if fields[1] =~ /^SN:([!-)+-<>-~][!-~]*)$/
 221       seq_name = $1.to_sym
 222     else
 223       raise SamError, "Bad sequence name: #{fields[1]}"
 224     end
 225
 226     if fields[2] =~ /^LN:(\d+)$/
 227       hash[:LN] = $1.to_i
 228     else
 229       raise SamError, "Bad sequence length: #{fields[2]}"
 230     end
 231
 232     (3 ... fields.size).each do |i|
 233       if fields[i] =~ /^(AS|M5|SP|UR):([ -~]+)$/
 234         hash[$1.to_sym] = $2
 235       else
 236         raise SamError, "Bad sequence tag: #{fields[i]}"
 237       end
 238     end
 239
 240     @header[:SQ][:SN] = Hash.new unless @header[:SQ][:SN].is_a? Hash
 241
 242     if @header[:SQ][:SN].has_key? seq_name
 243       raise SamError, "Non-unique sequence name: #{seq_name}"
 244     else
 245       @header[:SQ][:SN][seq_name] = hash
 246     end
 247   end
 248
 249   # Method to subparse read group lines.
 250   def subparse_read_group(line)
 251     @header[:RG] = Hash.new unless @header[:RG].is_a? Hash
 252     hash = {}
 253
 254     fields = line.split("\t")
 255
 256     if fields[1] =~ /^ID:([ -~]+)$/
 257       id = $1.to_sym
 258     else
 259       raise SamError, "Bad read group identifier: #{fields[1]}"
 260     end
 261
 262     (2 ... fields.size).each do |i|
 263       if fields[i] =~ /^(CN|DS|DT|FO|KS|LB|PG|PI|PL|PU|SM):([ -~]+)$/
 264         hash[$1.to_sym] = $2
 265       else
 266         raise SamError, "Bad read group tag: #{fields[i]}"
 267       end
 268     end
 269
 270     if hash.has_key? :FO
 271       unless hash[:FO] =~ /^\*|[ACMGRSVTWYHKDBN]+$/
 272         raise SamError, "Bad flow order: #{hash[:FO]}"
 273       end
 274     end
 275
 276     if hash.has_key? :PL
 277       unless hash[:PL] =~ /^(CAPILLARY|LS454|ILLUMINA|SOLID|HELICOS|IONTORRENT|PACBIO)$/
 278         raise SamError, "Bad platform: #{hash[:PL]}"
 279       end
 280     end
 281
 282     @header[:RG][:ID] = Hash.new unless @header[:RG][:ID].is_a? Hash
 283
 284     if @header[:RG][:ID].has_key? id
 285       raise SamError, "Non-unique read group identifier: #{id}"
 286     else
 287       @header[:RG][:ID][id] = hash
 288     end
 289   end
 290
 291   # Method to subparse program lines.
 292   def subparse_program(line)
 293     @header[:PG] = Hash.new unless @header[:PG].is_a? Hash
 294     hash = {}
 295
 296     fields = line.split("\t")
 297
 298     if fields[1] =~ /^ID:([ -~]+)$/
 299       id = $1.to_sym
 300     else
 301       raise SamError, "Bad program record identifier: #{fields[1]}"
 302     end
 303
 304     (2 ... fields.size).each do |i|
 305       if fields[i] =~ /^(PN|CL|PP|VN):([ -~]+)$/
 306         hash[$1.to_sym] = $2
 307       else
 308         raise SamError, "Bad program record tag: #{fields[i]}"
 309       end
 310     end
 311
 312     @header[:PG][:ID] = Hash.new unless @header[:PG][:ID].is_a? Hash
 313
 314     if @header[:PG][:ID].has_key? id
 315       raise SamError, "Non-unique program record identifier: #{id}"
 316     else
 317       @header[:PG][:ID][id] = hash
 318     end
 319   end
 320
 321   # Method to subparse comment lines.
 322   def subparse_comment(line)
 323     @header[:CO] = Array.new unless @header[:CO].is_a? Array
 324
 325     if line =~ /^@CO\t(.+)/
 326       @header[:CO] << $1
 327     else
 328       raise SamError, "Bad comment line: #{line}"
 329     end
 330   end
 331
 332   # Method to subparse alignment lines.
 333   def parse_alignment(line)
 334     fields = line.split("\t")
 335
 336     raise SamError, "Bad number of fields: #{fields.size} < 11" if fields.size < 11
 337
 338     qname = fields[0]
 339     flag  = fields[1].to_i
 340     rname = fields[2]
 341     pos   = fields[3].to_i
 342     mapq  = fields[4].to_i
 343     cigar = fields[5]
 344     rnext = fields[6]
 345     pnext = fields[7].to_i
 346     tlen  = fields[8].to_i
 347     seq   = fields[9]
 348     qual  = fields[10]
 349
 350     check_qname(qname)
 351     check_flag(flag)
 352     check_rname(rname)
 353     check_pos(pos)
 354     check_mapq(mapq)
 355     check_rnext(rnext)
 356     check_pnext(pnext)
 357     check_tlen(tlen)
 358     check_seq(seq)
 359     check_qual(qual)
 360
 361     entry = {}
 362     entry[:QNAME] = qname
 363     entry[:FLAG]  = Flag.new(flag)
 364     entry[:RNAME] = rname
 365     entry[:POS]   = pos
 366     entry[:MAPQ]  = mapq
 367     entry[:CIGAR] = Cigar.new(cigar)
 368     entry[:RNEXT] = rnext
 369     entry[:PNEXT] = pnext
 370     entry[:TLEN]  = tlen
 371     entry[:SEQ]   = (qual == '*') ? Seq.new(qname, seq) : Seq.new(qname, seq, qual)
 372     entry[:QUAL]  = qual
 373
 374     # Optional fields - where some are really important! HATE HATE HATE SAM!!!
 375
 376     fields[11 .. -1].each do |field|
 377       tag, type, val = field.split(':')
 378
 379       raise SamError, "Non-unique optional tag: #{tag}" if entry.has_key? tag.to_sym
 380
 381       # A [!-~] Printable character
 382
 383       # i [-+]?[0-9]+ Singed 32-bit integer
 384       if type == 'i'
 385         raise SamError, "Bad tag in optional field: #{field}" unless val =~ /^[-+]?[0-9]+$/
 386         val = val.to_i
 387       end
 388
 389       # f [-+]?[0-9]*\.?[0-9]+([eE][-+]?[0-9]+)? Single-precision floating number
 390       # Z [ !-~]+ Printable string, including space
 391       # H [0-9A-F]+ Byte array in the Hex format
 392       # B [cCsSiIf](,[-+]?[0-9]*\.?[0-9]+([eE][-+]?[0-9]+)?)+ Integer or numeric array
 393
 394       entry[tag.to_sym] = val
 395     end
 396
 397     entry
 398   end
 399
 400   # Method to check qname.
 401   def check_qname(qname)
 402     raise SamError, "Bad qname: #{qname}" unless qname =~ /^[!-?A-~]{1,255}$/
 403   end
 404
 405   # Method to check flag.
 406   def check_flag(flag)
 407     raise SamError, "Bad flag: #{flag}" unless (0 .. 2**16 - 1).include? flag
 408   end
 409
 410   # Method to check if rname, when not '*' and
 411   # @SQ header lines are present, is located in
 412   # the header hash.
 413   def check_rname(rname)
 414     raise SamError, "Bad rname: #{rname}" unless rname =~ /^(\*|[!-()+-<>-~][!-~]*)$/
 415
 416     unless @header.empty? or rname == '*'
 417       unless @header[:SQ][:SN].has_key? rname.to_sym
 418         raise SamError, "rname not found in header hash: #{rname}"
 419       end
 420     end
 421   end
 422
 423   # Method to check pos.
 424   def check_pos(pos)
 425     raise SamError, "Bad pos: #{pos}" unless (0 .. 2**29 - 1).include? pos
 426   end
 427
 428   # Method to check mapq.
 429   def check_mapq(mapq)
 430     raise SamError, "Bad mapq: #{mapq}" unless (0 .. 2**8 - 1).include? mapq
 431   end
 432
 433   # Method to check if rnext, when not '*' or '='
 434   # and @SQ header lines are present, is located
 435   # in the header hash.
 436   def check_rnext(rnext)
 437     raise SamError, "Bad rnext: #{rnext}" unless rnext =~ /^(\*|=|[!-()+-<>-~][!-~]*)$/
 438
 439     unless @header.empty? or rnext == '*' or rnext == '='
 440       unless @header[:SQ][:SN].has_key? rnext.to_sym
 441         raise SamError, "rnext not found in header hash: #{rnext}"
 442       end
 443     end
 444   end
 445
 446   # Method to check pnext.
 447   def check_pnext(pnext)
 448     raise SamError, "Bad pnext: #{pnext}" unless (0 .. 2**29 - 1).include? pnext
 449   end
 450
 451   # Method to check tlen.
 452   def check_tlen(tlen)
 453     raise SamError, "Bad tlen: #{tlen}" unless (-2**29 + 1 .. 2**29 - 1).include? tlen
 454   end
 455
 456   # Method to check seq.
 457   def check_seq(seq)
 458     raise SamError, "Bad seq: #{seq}" unless seq  =~ /^(\*|[A-Za-z=.]+)$/
 459   end
 460
 461   # Method to check qual.
 462   def check_qual(qual)
 463     raise SamError, "Bad qual: #{qual}" unless qual =~ /^[!-~]+$/
 464   end
 465
 466   # Method to deconvolute the SAM flag field.
 467   class Flag
 468     attr_reader :flag
 469
 470     # Method to initialize a Flag object.
 471     def initialize(flag)
 472       @flag = flag
 473     end
 474
 475     # Method to test if template have
 476     # multiple fragments in sequencing.
 477     def multi?
 478       flag & FLAG_MULTI
 479     end
 480
 481     # Method to test if each fragment
 482     # properly aligned according to the aligner.
 483     def aligned?
 484       flag & FLAG_ALIGNED
 485     end
 486
 487     # Method to test if the fragment was unmapped.
 488     def unmapped?
 489       flag & FLAG_UNMAPPED
 490     end
 491
 492     # Method to test if the next fragment was unmapped.
 493     def next_unmapped?
 494       flag & FLAG_NEXT_UNMAPPED
 495     end
 496
 497     # Method to test if the fragment was reverse complemented.
 498     def revcomp?
 499       flag & FLAG_REVCOMP
 500     end
 501
 502     # Method to test if the next fragment was reverse complemented.
 503     def next_revcomp?
 504       flag & FLAG_NEXT_REVCOMP
 505     end
 506
 507     # Method to test if the fragment was first in the template.
 508     def first?
 509       flag & FLAG_FIRST
 510     end
 511
 512     # Method to test if the fragment was last in the template.
 513     def last?
 514       flag & FLAG_LAST
 515     end
 516
 517     # Method to test for secondary alignment.
 518     def secondary_alignment?
 519       flag & FLAG_SECONDARY_ALIGNMENT
 520     end
 521
 522     # Method to test for quality fail.
 523     def quality_fail?
 524       flag & FLAG_QUALITY_FAIL
 525     end
 526
 527     # Method to test for PCR or optical duplicates.
 528     def duplicates?
 529       flag & FLAG_DUPLICATES
 530     end
 531   end
 532 end
 533
 534
 535 # >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>><<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
 536
 537
 538 __END__
 539