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