bp_bin/find_barcodes

   1 #!/usr/bin/env ruby
   2
   3 # Copyright (C) 2007-2011 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 # Find barcodes in sequences in the stream.
  28
  29 # >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>><<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
  30
  31 require 'maasha/biopieces'
  32 require 'maasha/bits'
  33 require 'pp'
  34
  35 # Hard coded Roche Genome Sequencing MIDs
  36 GSMID_HASH = {
  37   ACGAGTGCGT: "MID1",
  38   TCTCTATGCG: "MID10",
  39   TAGACTGCAC: "MID100",
  40   TAGCGCGCGC: "MID101",
  41   TAGCTCTATC: "MID102",
  42   TATAGACATC: "MID103",
  43   TATGATACGC: "MID104",
  44   TCACTCATAC: "MID105",
  45   TCATCGAGTC: "MID106",
  46   TCGAGCTCTC: "MID107",
  47   TCGCAGACAC: "MID108",
  48   TCTGTCTCGC: "MID109",
  49   TGATACGTCT: "MID11",
  50   TGAGTGACGC: "MID110",
  51   TGATGTGTAC: "MID111",
  52   TGCTATAGAC: "MID112",
  53   TGCTCGCTAC: "MID113",
  54   ACGTGCAGCG: "MID114",
  55   ACTCACAGAG: "MID115",
  56   AGACTCAGCG: "MID116",
  57   AGAGAGTGTG: "MID117",
  58   AGCTATCGCG: "MID118",
  59   AGTCTGACTG: "MID119",
  60   TACTGAGCTA: "MID12",
  61   AGTGAGCTCG: "MID120",
  62   ATAGCTCTCG: "MID121",
  63   ATCACGTGCG: "MID122",
  64   ATCGTAGCAG: "MID123",
  65   ATCGTCTGTG: "MID124",
  66   ATGTACGATG: "MID125",
  67   ATGTGTCTAG: "MID126",
  68   CACACGATAG: "MID127",
  69   CACTCGCACG: "MID128",
  70   CAGACGTCTG: "MID129",
  71   CATAGTAGTG: "MID13",
  72   CAGTACTGCG: "MID130",
  73   CGACAGCGAG: "MID131",
  74   CGATCTGTCG: "MID132",
  75   CGCGTGCTAG: "MID133",
  76   CGCTCGAGTG: "MID134",
  77   CGTGATGACG: "MID135",
  78   CTATGTACAG: "MID136",
  79   CTCGATATAG: "MID137",
  80   CTCGCACGCG: "MID138",
  81   CTGCGTCACG: "MID139",
  82   CGAGAGATAC: "MID14",
  83   CTGTGCGTCG: "MID140",
  84   TAGCATACTG: "MID141",
  85   TATACATGTG: "MID142",
  86   TATCACTCAG: "MID143",
  87   TATCTGATAG: "MID144",
  88   TCGTGACATG: "MID145",
  89   TCTGATCGAG: "MID146",
  90   TGACATCTCG: "MID147",
  91   TGAGCTAGAG: "MID148",
  92   TGATAGAGCG: "MID149",
  93   ATACGACGTA: "MID15",
  94   TGCGTGTGCG: "MID150",
  95   TGCTAGTCAG: "MID151",
  96   TGTATCACAG: "MID152",
  97   TGTGCGCGTG: "MID153",
  98   TCACGTACTA: "MID16",
  99   CGTCTAGTAC: "MID17",
 100   TCTACGTAGC: "MID18",
 101   TGTACTACTC: "MID19",
 102   ACGCTCGACA: "MID2",
 103   ACGACTACAG: "MID20",
 104   CGTAGACTAG: "MID21",
 105   TACGAGTATG: "MID22",
 106   TACTCTCGTG: "MID23",
 107   TAGAGACGAG: "MID24",
 108   TCGTCGCTCG: "MID25",
 109   ACATACGCGT: "MID26",
 110   ACGCGAGTAT: "MID27",
 111   ACTACTATGT: "MID28",
 112   ACTGTACAGT: "MID29",
 113   AGACGCACTC: "MID3",
 114   AGACTATACT: "MID30",
 115   AGCGTCGTCT: "MID31",
 116   AGTACGCTAT: "MID32",
 117   ATAGAGTACT: "MID33",
 118   CACGCTACGT: "MID34",
 119   CAGTAGACGT: "MID35",
 120   CGACGTGACT: "MID36",
 121   TACACACACT: "MID37",
 122   TACACGTGAT: "MID38",
 123   TACAGATCGT: "MID39",
 124   AGCACTGTAG: "MID4",
 125   TACGCTGTCT: "MID40",
 126   TAGTGTAGAT: "MID41",
 127   TCGATCACGT: "MID42",
 128   TCGCACTAGT: "MID43",
 129   TCTAGCGACT: "MID44",
 130   TCTATACTAT: "MID45",
 131   TGACGTATGT: "MID46",
 132   TGTGAGTAGT: "MID47",
 133   ACAGTATATA: "MID48",
 134   ACGCGATCGA: "MID49",
 135   ATCAGACACG: "MID5",
 136   ACTAGCAGTA: "MID50",
 137   AGCTCACGTA: "MID51",
 138   AGTATACATA: "MID52",
 139   AGTCGAGAGA: "MID53",
 140   AGTGCTACGA: "MID54",
 141   CGATCGTATA: "MID55",
 142   CGCAGTACGA: "MID56",
 143   CGCGTATACA: "MID57",
 144   CGTACAGTCA: "MID58",
 145   CGTACTCAGA: "MID59",
 146   ATATCGCGAG: "MID6",
 147   CTACGCTCTA: "MID60",
 148   CTATAGCGTA: "MID61",
 149   TACGTCATCA: "MID62",
 150   TAGTCGCATA: "MID63",
 151   TATATATACA: "MID64",
 152   TATGCTAGTA: "MID65",
 153   TCACGCGAGA: "MID66",
 154   TCGATAGTGA: "MID67",
 155   TCGCTGCGTA: "MID68",
 156   TCTGACGTCA: "MID69",
 157   CGTGTCTCTA: "MID7",
 158   TGAGTCAGTA: "MID70",
 159   TGTAGTGTGA: "MID71",
 160   TGTCACACGA: "MID72",
 161   TGTCGTCGCA: "MID73",
 162   ACACATACGC: "MID74",
 163   ACAGTCGTGC: "MID75",
 164   ACATGACGAC: "MID76",
 165   ACGACAGCTC: "MID77",
 166   ACGTCTCATC: "MID78",
 167   ACTCATCTAC: "MID79",
 168   CTCGCGTGTC: "MID8",
 169   ACTCGCGCAC: "MID80",
 170   AGAGCGTCAC: "MID81",
 171   AGCGACTAGC: "MID82",
 172   AGTAGTGATC: "MID83",
 173   AGTGACACAC: "MID84",
 174   AGTGTATGTC: "MID85",
 175   ATAGATAGAC: "MID86",
 176   ATATAGTCGC: "MID87",
 177   ATCTACTGAC: "MID88",
 178   CACGTAGATC: "MID89",
 179   TAGTATCAGC: "MID9",
 180   CACGTGTCGC: "MID90",
 181   CATACTCTAC: "MID91",
 182   CGACACTATC: "MID92",
 183   CGAGACGCGC: "MID93",
 184   CGTATGCGAC: "MID94",
 185   CGTCGATCTC: "MID95",
 186   CTACGACTGC: "MID96",
 187   CTAGTCACTC: "MID97",
 188   CTCTACGCTC: "MID98",
 189   CTGTACATAC: "MID99",
 190 }
 191
 192 # Hard coded Roche Rapid Labrary MIDs
 193 RLMID_HASH = {
 194   ACACGACGACT: "RL1",
 195   ACACGTAGTAT: "RL2",
 196   ACACTACTCGT: "RL3",
 197   ACGACACGTAT: "RL4",
 198   ACGAGTAGACT: "RL5",
 199   ACGCGTCTAGT: "RL6",
 200   ACGTACACACT: "RL7",
 201   ACGTACTGTGT: "RL8",
 202   ACGTAGATCGT: "RL9",
 203   ACTACGTCTCT: "RL10",
 204   ACTATACGAGT: "RL11",
 205   ACTCGCGTCGT: "RL12"
 206 }
 207
 208 class BarCodeFinderError < StandardError; end
 209
 210 class BarCodeFinder
 211   def initialize(pos, max_mismatches)
 212     @pos            = pos
 213     @max_mismatches = max_mismatches
 214     @barcode_hash   = {}
 215     @barcode_sizes  = []
 216   end
 217
 218   def parse_barcodes(file)
 219     length = 0
 220     hash   = {}
 221
 222     File.open(file, "r") do |ios|
 223       while not ios.eof?
 224         name, barcode = ios.readline.chomp.split(/\s+/)
 225
 226         if length > 0
 227           raise BarCodeFinderError, "Uneven barcode length detected" if length != barcode.length
 228         else
 229           length = barcode.length
 230         end
 231
 232         hash[barcode.upcase.to_sym] = name
 233       end
 234     end
 235
 236     load_barcodes(hash)
 237   end
 238
 239   def load_barcodes(hash)
 240     @barcode_sizes << hash.keys.first.to_s.size
 241
 242     @barcode_hash.merge!(hash)
 243   end
 244
 245   def find_barcode(seq)
 246     raise BarCodeFinderError, "No barcodes to find" if @barcode_hash.empty?
 247
 248     @barcode_sizes.each do |size|
 249       hamming_dist = 0
 250       barcode      = seq[@pos ... @pos + size].upcase.to_sym
 251
 252       if @barcode_hash.has_key? barcode
 253         return BarCode.new(barcode, @barcode_hash[barcode], @pos, size, hamming_dist)
 254       elsif @max_mismatches > 0
 255         @barcode_hash.each_key do |key|
 256           if key.to_s.length == barcode.to_s.size
 257             hamming_dist = barcode.to_s.hamming_distance(key.to_s)
 258
 259             if hamming_dist <= @max_mismatches
 260               return BarCode.new(key, @barcode_hash[key], @pos, size, hamming_dist)
 261             end
 262           end
 263         end
 264       end
 265     end
 266
 267     nil
 268   end
 269
 270   private
 271
 272   # Class to hold a BacCode object.
 273   class BarCode
 274     attr_accessor :barcode, :name, :pos, :len, :mismatches
 275
 276     def initialize(barcode, name, pos, len, mismatches)
 277       @barcode    = barcode
 278       @name       = name
 279       @pos        = pos
 280       @len        = len
 281       @mismatches = mismatches
 282     end
 283
 284     def to_hash
 285       hash = {}
 286       hash[:BARCODE]            = @barcode
 287       hash[:BARCODE_NAME]       = @name
 288       hash[:BARCODE_POS]        = @pos
 289       hash[:BARCODE_LEN]        = @len
 290       hash[:BARCODE_MISMATCHES] = @mismatches
 291       hash
 292     end
 293
 294     def start
 295       @pos + @len
 296     end
 297   end
 298 end
 299
 300 casts = []
 301 casts << {:long=>'barcodes_in', :short=>'b', :type=>'file!', :mandatory=>false, :default=>nil, :allowed=>nil,     :disallowed=>nil}
 302 casts << {:long=>'pos',         :short=>'p', :type=>'uint',  :mandatory=>false, :default=>0,   :allowed=>nil,     :disallowed=>nil}
 303 casts << {:long=>'mismatches',  :short=>'m', :type=>'uint',  :mandatory=>false, :default=>0,   :allowed=>"0,1,2", :disallowed=>nil}
 304 casts << {:long=>'gsmids',      :short=>'g', :type=>'flag',  :mandatory=>false, :default=>nil, :allowed=>nil,     :disallowed=>nil}
 305 casts << {:long=>'rlmids',      :short=>'r', :type=>'flag',  :mandatory=>false, :default=>nil, :allowed=>nil,     :disallowed=>nil}
 306 casts << {:long=>'remove',      :short=>'R', :type=>'flag',  :mandatory=>false, :default=>nil, :allowed=>nil,     :disallowed=>nil}
 307
 308 options = Biopieces.options_parse(ARGV, casts)
 309
 310 bc_finder = BarCodeFinder.new(options[:pos], options[:mismatches])
 311 bc_finder.parse_barcodes(options[:barcodes_in]) if options[:barcodes_in]
 312 bc_finder.load_barcodes(GSMID_HASH)             if options[:gsmids]
 313 bc_finder.load_barcodes(RLMID_HASH)             if options[:rlmids]
 314
 315 Biopieces.open(options[:stream_in], options[:stream_out]) do |input, output|
 316   input.each_record do |record|
 317     if record.has_key? :SEQ
 318       if barcode = bc_finder.find_barcode(record[:SEQ])
 319         record.merge!(barcode.to_hash)
 320
 321         if options[:remove]
 322           record[:SEQ]     = record[:SEQ][barcode.start .. -1]
 323           record[:SCORES]  = record[:SCORES][barcode.start .. -1] if record[:SCORES]
 324           record[:SEQ_LEN] = record[:SEQ].length
 325         end
 326       end
 327     end
 328
 329     output.puts record
 330   end
 331 end
 332
 333
 334 # >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>><<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
 335
 336
 337 __END__