bp_bin/remove_mids

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