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 Sam < Filesys
  36   attr_accessor :io
  37
  38   def initialize(io = nil)
  39     @io          = io
  40     @header_hash = {}
  41   end
  42
  43   # Each header line should match: /^@[A-Za-z][A-Za-z](\t[A-Za-z][A-Za-z0-9]:[ -~]+)+$/ or /^@CO\t.*/.
  44   # Tags containing lowercase letters are reserved for end users.
  45   def header
  46     @io.each_line do |line|
  47       line.chomp!
  48
  49       if line =~ /^@([A-Za-z][A-Za-z])/
  50         tag = $1
  51
  52         case tag
  53         when 'HD' then parse_header(line)
  54         when 'SQ' then parse_sequence(line)
  55         when 'RG' then parse_read_group(line)
  56         when 'PG' then parse_program(line)
  57         when 'CO' then parse_comment(line)
  58         else
  59           raise SamError, "Unknown header tag: #{tag}"
  60         end
  61       else
  62         @io.rewind
  63         break
  64       end
  65     end
  66
  67     return @header_hash.empty? ? nil : @header_hash
  68   end
  69
  70   private
  71
  72   def parse_header(line)
  73     hash   = {}
  74     fields = line.split("\t")
  75
  76     if fields[1] =~ /^VN:([0-9]+\.[0-9]+)$/
  77       hash[:VN] = $1.to_f
  78     else
  79       raise SamError, "Bad version number: #{fields[1]}"
  80     end
  81
  82     if fields.size > 2
  83       if fields[2] =~ /^SO:(unknown|unsorted|queryname|coordinate)$/
  84         hash[:SO] = $1
  85       else
  86         raise SamError, "Bad sort order: #{fields[2]}"
  87       end
  88     end
  89
  90     @header_hash[:HD] = hash
  91   end
  92
  93   def parse_sequence(line)
  94     @header_hash[:SQ] = Hash.new unless @header_hash[:SQ].is_a? Hash
  95     hash = {}
  96
  97     fields = line.split("\t")
  98
  99     if fields[1] =~ /^SN:([!-)+-<>-~][!-~]*)$/
 100       seq_name = $1.to_sym
 101     else
 102       raise SamError, "Bad sequence name: #{fields[1]}"
 103     end
 104
 105     if fields[2] =~ /^LN:(\d+)$/
 106       hash[:LN] = $1.to_i
 107     else
 108       raise SamError, "Bad sequence length: #{fields[2]}"
 109     end
 110
 111     (3 ... fields.size).each do |i|
 112       if fields[i] =~ /^(AS|M5|SP|UR):([ -~]+)$/
 113         hash[$1.to_sym] = $2
 114       else
 115         raise SamError, "Bad sequence tag: #{fields[i]}"
 116       end
 117     end
 118
 119     @header_hash[:SQ][:SN] = Hash.new unless @header_hash[:SQ][:SN].is_a? Hash
 120
 121     if @header_hash[:SQ][:SN].has_key? seq_name
 122       raise SamError, "Non-unique sequence name: #{seq_name}"
 123     else
 124       @header_hash[:SQ][:SN][seq_name] = hash
 125     end
 126   end
 127
 128   def parse_read_group(line)
 129     @header_hash[:RG] = Hash.new unless @header_hash[:RG].is_a? Hash
 130     hash = {}
 131
 132     fields = line.split("\t")
 133
 134     if fields[1] =~ /^ID:([ -~]+)$/
 135       id = $1.to_sym
 136     else
 137       raise SamError, "Bad read group identifier: #{fields[1]}"
 138     end
 139
 140     (2 ... fields.size).each do |i|
 141       if fields[i] =~ /^(CN|DS|DT|FO|KS|LB|PG|PI|PL|PU|SM):([ -~]+)$/
 142         hash[$1.to_sym] = $2
 143       else
 144         raise SamError, "Bad read group tag: #{fields[i]}"
 145       end
 146     end
 147
 148     if hash.has_key? :FO
 149       unless hash[:FO] =~ /^\*|[ACMGRSVTWYHKDBN]+$/
 150         raise SamError, "Bad flow order: #{hash[:FO]}"
 151       end
 152     end
 153
 154     if hash.has_key? :PL
 155       unless hash[:PL] =~ /^(CAPILLARY|LS454|ILLUMINA|SOLID|HELICOS|IONTORRENT|PACBIO)$/
 156         raise SamError, "Bad platform: #{hash[:PL]}"
 157       end
 158     end
 159
 160     @header_hash[:RG][:ID] = Hash.new unless @header_hash[:RG][:ID].is_a? Hash
 161
 162     if @header_hash[:RG][:ID].has_key? id
 163       raise SamError, "Non-unique read group identifier: #{id}"
 164     else
 165       @header_hash[:RG][:ID][id] = hash
 166     end
 167   end
 168
 169   def parse_program(line)
 170     @header_hash[:PG] = Hash.new unless @header_hash[:PG].is_a? Hash
 171     hash = {}
 172
 173     fields = line.split("\t")
 174
 175     if fields[1] =~ /^ID:([ -~]+)$/
 176       id = $1.to_sym
 177     else
 178       raise SamError, "Bad program record identifier: #{fields[1]}"
 179     end
 180
 181     (2 ... fields.size).each do |i|
 182       if fields[i] =~ /^(PN|CL|PP|VN):([ -~]+)$/
 183         hash[$1.to_sym] = $2
 184       else
 185         raise SamError, "Bad program record tag: #{fields[i]}"
 186       end
 187     end
 188
 189     @header_hash[:PG][:ID] = Hash.new unless @header_hash[:PG][:ID].is_a? Hash
 190
 191     if @header_hash[:PG][:ID].has_key? id
 192       raise SamError, "Non-unique program record identifier: #{id}"
 193     else
 194       @header_hash[:PG][:ID][id] = hash
 195     end
 196   end
 197
 198   def parse_comment(line)
 199     @header_hash[:CO] = Array.new unless @header_hash[:CO].is_a? Array
 200
 201     if line =~ /^@CO\t(.+)/
 202       @header_hash[:CO] << $1
 203     else
 204       raise SamError, "Bad comment line: #{line}"
 205     end
 206   end
 207
 208   def get_entry
 209   end
 210 end
 211
 212
 213 # >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>><<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
 214
 215
 216 __END__
 217