binaries/src/mafft/scripts/mafft-sparsecore.rb

   1 #! /usr/bin/env ruby
   2 require 'optparse'
   3
   4 mafftpath = "/usr/local/bin/mafft"
   5
   6 def cleartempfiles( filenames )
   7         for f in filenames
   8                 system( "rm -rf #{f}" )
   9         end
  10 end
  11
  12
  13 seed = 0
  14 scand = "50%"
  15 npick = 500
  16 infn = ""
  17 reorderoption = "--reorder"
  18 pickoptions = " --retree 1 "
  19 coreoptions = " --globalpair --maxiterate 100 "
  20 corelastarg = " "
  21 addoptions = " "
  22 directionoptions = " --retree 0 --pileup "
  23 markcore = ""
  24 randompickup = true
  25 outnum = false
  26
  27 begin
  28         params = ARGV.getopts('m:s:n:p:i:C:L:A:o:MhuD:')
  29 rescue => e
  30     STDERR.puts e
  31         STDERR.puts "See #{$0} -h"
  32     exit 1
  33 end
  34
  35 #p params
  36
  37 mafftpath = params["m"]         if params["m"]
  38 seed = params["s"].to_i         if params["s"]
  39 scand = params["n"].to_s        if params["n"]
  40 npick = params["p"].to_i        if params["p"]
  41 infn = params["i"]              if params["i"]
  42 #pickoptions += params["P"]     if params["P"]
  43 coreoptions += params["C"]      if params["C"] # tsuikagaki!
  44 corelastarg += params["L"]      if params["L"] # tsuikagaki!
  45 addoptions  += params["A"]      if params["A"]
  46 directionoptions += params["D"] if params["D"] # tsuikagaki
  47 markcore = "*"                  if params["M"]
  48 #randompickup = false           if params["S"]
  49 reorderoption = ""              if params["o"] =~ /^i/
  50 outnum = true                   if params["u"]
  51
  52 if params["h"] then
  53         STDERR.puts "Usage: #{$0} -i inputfile [options]"
  54         STDERR.puts "Options:"
  55         STDERR.puts "   -i string     Input file."
  56         STDERR.puts "   -m string     Mafft command.  Default: mafft"
  57         STDERR.puts "   -s int        Seed.  Default:0"
  58         STDERR.puts "   -n int        Number of candidates for core sequences.  Default: upper 50% in length"
  59         STDERR.puts "   -p int        Number of core sequences.  Default: 500"
  60 #       STDERR.puts "   -P \"string\"   Mafft options for the PICKUP stage."
  61 #       STDERR.puts "                 Default: \"--retree 1\""
  62 #       STDERR.puts "   -S            Tree-based pickup.  Default: off"
  63         STDERR.puts "   -C \"string\"   Mafft options for the CORE stage."
  64         STDERR.puts "                 Default: \"--globalpair --maxiterate 100\""
  65         STDERR.puts "   -A \"string\"   Mafft options for the ADD stage."
  66         STDERR.puts "                 Default: \"\""
  67         STDERR.puts "   -D \"string\"   Mafft options for inferring the direction of nucleotide sequences."
  68         STDERR.puts "                 Default: \"\""
  69         STDERR.puts "   -o r or i     r: Reorder the sequences based on similarity.  Default"
  70         STDERR.puts "                 i: Same as input."
  71         exit 1
  72 end
  73
  74 if infn == "" then
  75         STDERR.puts "Give input file with -i."
  76         exit 1
  77 end
  78
  79
  80
  81 pid = $$.to_s
  82 tmpdir = ENV["TMPDIR"]
  83 tmpdir = "/tmp" if tmpdir == nil
  84 tempfiles = []
  85 tempfiles.push( temp_pf = tmpdir + "/_pf" + pid )
  86 tempfiles.push( temp_nf = tmpdir + "/_nf" + pid )
  87 tempfiles.push( temp_cf = tmpdir + "/_cf" + pid )
  88 tempfiles.push( temp_of = tmpdir + "/_of" + pid )
  89
  90 Signal.trap(:INT){cleartempfiles( tempfiles ); exit 1}
  91 at_exit{ cleartempfiles( tempfiles )}
  92
  93 system "#{mafftpath} --version > #{temp_of} 2>&1"
  94
  95 fp = File.open( temp_of, "r" )
  96         line = fp.gets
  97 fp.close
  98
  99
 100 versionnum = line.split(' ')[0].sub(/v/,"").to_f
 101
 102 if versionnum < 7.210 then
 103         STDERR.puts "\n"
 104         STDERR.puts "Please use mafft version >= 7.210\n"
 105         STDERR.puts "\n"
 106         exit
 107 end
 108
 109 srand( seed )
 110
 111 def readfasta( fp, name, seq )
 112         nseq = 0
 113         tmpseq = ""
 114         while fp.gets
 115                 if $_ =~ /^>/ then
 116                         name.push( $_.sub(/>/,"").chop )
 117                         seq.push( tmpseq ) if nseq > 0
 118                         nseq += 1
 119                         tmpseq = ""
 120                 else
 121                         tmpseq += $_.strip
 122                 end
 123         end
 124         seq.push( tmpseq )
 125         return nseq
 126 end
 127
 128
 129
 130 begin
 131         infp = File.open( infn, "r" )
 132 rescue => e
 133     STDERR.puts e
 134     exit 1
 135 end
 136 infp.close
 137
 138 if directionoptions =~ /--adjustdirection/ then
 139         system( mafftpath + "#{directionoptions} #{infn} > #{temp_of}" )
 140 else
 141         system( "cp #{infn} #{temp_of}" )
 142 end
 143
 144 tname = []
 145 tseq = []
 146 infp = File.open( temp_of, "r" )
 147 tin = readfasta( infp, tname, tseq )
 148 infp.close
 149 lenhash = {}
 150
 151 if outnum then
 152         for i in 0..(tin-1)
 153                 tname[i] = "_numo_s_#{i}_numo_e_" + tname[i]
 154         end
 155 end
 156
 157 npick = 0 if npick == 1
 158 npick = tin if npick > tin
 159
 160
 161 if scand =~ /%$/ then
 162         ncand = (tin * scand.to_f * 0.01 ).to_i
 163 else
 164         ncand = scand.to_i
 165 end
 166
 167 if ncand < 0 || ncand > tin then
 168         STDERR.puts "Error.  -n #{scand}?"
 169         exit 1
 170 end
 171
 172 ncand = npick if ncand < npick
 173 ncand = tin if ncand > tin
 174
 175 STDERR.puts "ncand = #{ncand}, npick = #{npick}"
 176
 177
 178 sai = []
 179 for i in 0..(tin-1)
 180         lenhash[i] = tseq[i].gsub(/-/,"").length
 181 end
 182
 183 i = 0
 184 sorted = lenhash.sort_by{|key, value| [-value, i+=1]}
 185 #for i in 0..(ncand-1)
 186 #       sai[sorted[i][0]] = 1
 187 #end
 188 #for i in ncand..(tin-1)
 189 #       sai[sorted[i][0]] = 0
 190 #end
 191
 192 ncandres = 0
 193 ntsukau = 0
 194 for i in 0..(tin-1)
 195         cand = sorted[i][0]
 196         if tname[cand] =~ /^_focus_/ then
 197                 sai[cand] = 0
 198                 ntsukau += 1
 199         elsif ncandres < ncand  then
 200                 unless  tname[cand] =~ /^_tsukawanai_/ then
 201                         sai[cand] = 1
 202                         ncandres += 1
 203                 else
 204                         sai[cand] = 0
 205                 end
 206         else
 207                 sai[cand] = 0
 208         end
 209 end
 210
 211 if ncandres+ntsukau < npick
 212         STDERR.puts "ncandres = #{ncandres}"
 213         STDERR.puts "ncand = #{ncand}"
 214         STDERR.puts "ntsukau = #{ntsukau}"
 215         STDERR.puts "npick = #{npick}"
 216         STDERR.puts "Too many _tsukawanai_ sequences."
 217         exit 1
 218 end
 219
 220 if ntsukau > npick
 221         STDERR.puts "ntsukau = #{ntsukau}"
 222         STDERR.puts "npick = #{npick}"
 223         STDERR.puts "Too many _focus_ sequences."
 224         exit 1
 225 end
 226
 227 #p sai
 228 #for i in 0..(tin-1)
 229 #       puts sai[i].to_s + " " + tname[i]
 230 #end
 231
 232 npickrand = npick - ntsukau
 233
 234 if randompickup
 235         pick = []
 236         for i in 0..(npickrand-1)
 237                 pick[i] = 1
 238         end
 239         for i in npickrand..(ncandres-1)
 240                 pick[i] = 0
 241         end
 242         pick2 = pick.sort_by{rand}
 243         pick = pick2
 244 #       p pick
 245 #       p sai
 246
 247         ipick = 0
 248         for i in 0..(tin-1)
 249                 if sai[i] == 1 then
 250                         if pick[ipick] == 0 then
 251                                 sai[i] = 0
 252                         end
 253                         ipick += 1
 254                 end
 255         end
 256 #       p sai
 257
 258         for i in 0..(tin-1)
 259                 if tname[i] =~ /^_focus_/ then
 260                         sai[i] = 1
 261                 end
 262         end
 263 #       p sai
 264
 265         pfp = File.open( temp_pf, 'w' )
 266         nfp = File.open( temp_nf, 'w' )
 267
 268         i = 0
 269         while i < tin
 270                 if sai[i] == 1 then
 271                         pfp.puts ">" + i.to_s + " " + ">" + markcore + tname[i]
 272                         pfp.puts tseq[i]
 273                 else
 274                         nfp.puts ">" + i.to_s + " " + ">" + tname[i]
 275                         nfp.puts tseq[i]
 276                 end
 277                 i += 1
 278         end
 279
 280         nfp.close
 281         pfp.close
 282
 283 else   # yamerukamo
 284         STDERR.puts "Not supported in this version"
 285         exit 1
 286 end
 287
 288 if npick > 1 then
 289         if npick < tin then
 290                 system( mafftpath + " #{coreoptions} #{temp_pf} #{corelastarg} > #{temp_cf}" ) # add de sort
 291         else
 292                 system( mafftpath + " #{coreoptions} #{reorderoption} #{temp_pf} #{corelastarg} > #{temp_cf}" ) # ima sort
 293         end
 294         res = ( File::stat(temp_cf).size == 0 )
 295 else
 296         system( "cat /dev/null > #{temp_cf}" )
 297         res = false
 298 end
 299
 300 if res == true then
 301         STDERR.puts "\n\nError in the core alignment stage.\n\n"
 302         exit 1
 303 end
 304
 305
 306 if npick < tin
 307         system( mafftpath + " #{addoptions} #{reorderoption} --add #{temp_nf} #{temp_cf} > #{temp_of}" )
 308         res = ( File::stat(temp_of).size == 0 )
 309 else
 310         system( "cp #{temp_cf} #{temp_of}" )
 311         res = false
 312 end
 313
 314 if res == true then
 315         STDERR.puts "\n\nError in the add stage.\n\n"
 316         exit 1
 317 end
 318
 319 resname = []
 320 resseq = []
 321 resfp = File.open( temp_of, "r" )
 322 nres = readfasta( resfp, resname, resseq )
 323 resfp.close
 324
 325 if reorderoption =~ /--reorder/ then
 326         for i in 0..(nres-1)
 327                 puts ">" + resname[i].sub(/^[0-9]* >/,"")
 328                 puts resseq[i]
 329         end
 330 else
 331         seqhash = {}
 332         namehash = {}
 333         seqlast = []
 334         namelast = []
 335         nlast = 0
 336         for i in 0..(nres-1)
 337                 if resname[i] =~ /^[0-9]* >/
 338                         key = resname[i].split(' ')[0]
 339                         seqhash[key] = resseq[i]
 340                         namehash[key] = resname[i]
 341                 else
 342                         seqlast.push( resseq[i] )
 343                         namelast.push( resname[i] )
 344                         nlast += 1
 345                 end
 346         end
 347         for i in 0..(nlast-1)
 348                 puts ">" + namelast[i]
 349                 puts seqlast[i]
 350         end
 351         for i in 0..(nres-nlast-1)
 352                 key = i.to_s
 353                 puts ">" + namehash[key].sub(/^[0-9]* >/,"")
 354                 puts seqhash[key]
 355         end
 356 end
 357
 358