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   AGACTCGACGT: "RL13",
 207   AGTACGAGAGT: "RL14",
 208   AGTACTACTAT: "RL15",
 209   AGTAGACGTCT: "RL16",
 210   AGTCGTACACT: "RL17",
 211   AGTGTAGTAGT: "RL18",
 212   ATAGTATACGT: "RL19",
 213   CAGTACGTACT: "RL20",
 214   CGACGACGCGT: "RL21",
 215   CGACGAGTACT: "RL22",
 216   CGATACTACGT: "RL23",
 217   CGTACGTCGAT: "RL24",
 218   CTACTCGTAGT: "RL25",
 219   GTACAGTACGT: "RL26",
 220   GTCGTACGTAT: "RL27",
 221   GTGTACGACGT: "RL28",
 222   ACACAGTGAGT: "RL29",
 223   ACACTCATACT: "RL30",
 224   ACAGACAGCGT: "RL31",
 225   ACAGACTATAT: "RL32",
 226   ACAGAGACTCT: "RL33",
 227   ACAGCTCGTGT: "RL34",
 228   ACAGTGTCGAT: "RL35",
 229   ACGAGCGCGCT: "RL36",
 230   ACGATGAGTGT: "RL37",
 231   ACGCGAGAGAT: "RL38",
 232   ACGCTCTCTCT: "RL39",
 233   ACGTCGCTGAT: "RL40",
 234   ACGTCTAGCAT: "RL41",
 235   ACTAGTGATAT: "RL42",
 236   ACTCACACTGT: "RL43",
 237   ACTCACTAGCT: "RL44",
 238   ACTCTATATAT: "RL45",
 239   ACTGATCTCGT: "RL46",
 240   ACTGCTGTACT: "RL47",
 241   ACTGTAGCGCT: "RL48",
 242   AGACACTCACT: "RL49",
 243   AGACATATAGT: "RL50"
 244 }
 245
 246 class BarCodeFinderError < StandardError; end
 247
 248 class BarCodeFinder
 249   def initialize(pos, max_mismatches)
 250     @pos            = pos
 251     @max_mismatches = max_mismatches
 252     @barcode_hash   = {}
 253     @barcode_sizes  = []
 254   end
 255
 256   def parse_barcodes(file)
 257     length = 0
 258     hash   = {}
 259
 260     File.open(file, "r") do |ios|
 261       while not ios.eof?
 262         name, barcode = ios.readline.chomp.split(/\s+/)
 263
 264         if length > 0
 265           raise BarCodeFinderError, "Uneven barcode length detected" if length != barcode.length
 266         else
 267           length = barcode.length
 268         end
 269
 270         hash[barcode.upcase.to_sym] = name
 271       end
 272     end
 273
 274     load_barcodes(hash)
 275   end
 276
 277   def load_barcodes(hash)
 278     @barcode_sizes << hash.keys.first.to_s.size
 279
 280     @barcode_hash.merge!(hash)
 281   end
 282
 283   def find_barcode(seq)
 284     raise BarCodeFinderError, "No barcodes to find" if @barcode_hash.empty?
 285
 286     @barcode_sizes.each do |size|
 287       hamming_dist = 0
 288       barcode      = seq[@pos ... @pos + size].upcase.to_sym
 289
 290       if @barcode_hash[barcode]
 291         return BarCode.new(barcode, @barcode_hash[barcode], @pos, size, hamming_dist)
 292       elsif @max_mismatches > 0
 293         @barcode_hash.each_key do |key|
 294           if key.to_s.length == barcode.to_s.size
 295             hamming_dist = barcode.to_s.hamming_distance(key.to_s)
 296
 297             if hamming_dist <= @max_mismatches
 298               return BarCode.new(key, @barcode_hash[key], @pos, size, hamming_dist)
 299             end
 300           end
 301         end
 302       end
 303     end
 304
 305     nil
 306   end
 307
 308   private
 309
 310   # Class to hold a BacCode object.
 311   class BarCode
 312     attr_accessor :barcode, :name, :pos, :len, :mismatches
 313
 314     def initialize(barcode, name, pos, len, mismatches)
 315       @barcode    = barcode
 316       @name       = name
 317       @pos        = pos
 318       @len        = len
 319       @mismatches = mismatches
 320     end
 321
 322     def to_hash
 323       hash = {}
 324       hash[:BARCODE]            = @barcode
 325       hash[:BARCODE_NAME]       = @name
 326       hash[:BARCODE_POS]        = @pos
 327       hash[:BARCODE_LEN]        = @len
 328       hash[:BARCODE_MISMATCHES] = @mismatches
 329       hash
 330     end
 331
 332     def start
 333       @pos + @len
 334     end
 335   end
 336 end
 337
 338 casts = []
 339 casts << {:long=>'barcodes_in', :short=>'b', :type=>'file!', :mandatory=>false, :default=>nil, :allowed=>nil,     :disallowed=>nil}
 340 casts << {:long=>'pos',         :short=>'p', :type=>'uint',  :mandatory=>false, :default=>0,   :allowed=>nil,     :disallowed=>nil}
 341 casts << {:long=>'mismatches',  :short=>'m', :type=>'uint',  :mandatory=>false, :default=>0,   :allowed=>"0,1,2", :disallowed=>nil}
 342 casts << {:long=>'gsmids',      :short=>'g', :type=>'flag',  :mandatory=>false, :default=>nil, :allowed=>nil,     :disallowed=>nil}
 343 casts << {:long=>'rlmids',      :short=>'r', :type=>'flag',  :mandatory=>false, :default=>nil, :allowed=>nil,     :disallowed=>nil}
 344 casts << {:long=>'remove',      :short=>'R', :type=>'flag',  :mandatory=>false, :default=>nil, :allowed=>nil,     :disallowed=>nil}
 345
 346 options = Biopieces.options_parse(ARGV, casts)
 347
 348 bc_finder = BarCodeFinder.new(options[:pos], options[:mismatches])
 349 bc_finder.parse_barcodes(options[:barcodes_in]) if options[:barcodes_in]
 350 bc_finder.load_barcodes(GSMID_HASH)             if options[:gsmids]
 351 bc_finder.load_barcodes(RLMID_HASH)             if options[:rlmids]
 352
 353 Biopieces.open(options[:stream_in], options[:stream_out]) do |input, output|
 354   input.each_record do |record|
 355     if record[:SEQ]
 356       if barcode = bc_finder.find_barcode(record[:SEQ])
 357         record.merge!(barcode.to_hash)
 358
 359         if options[:remove]
 360           record[:SEQ]     = record[:SEQ][barcode.start .. -1]
 361           record[:SCORES]  = record[:SCORES][barcode.start .. -1] if record[:SCORES]
 362           record[:SEQ_LEN] = record[:SEQ].length
 363         end
 364       end
 365     end
 366
 367     output.puts record
 368   end
 369 end
 370
 371
 372 # >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>><<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
 373
 374
 375 __END__