binaries/src/mafft/core/mafft-homologs.tmpl

   1 #!/usr/bin/env ruby
   2
   3 localdb = "sp"
   4 # database name from which homologues are collected
   5 # by locally installed blast. Leave this if you do
   6 # not use the '-l' option.
   7
   8 mafftpath = "_BINDIR/mafft"
   9 # path of mafft. "/usr/local/bin/mafft"
  10 # if mafft is in your command path, "mafft" is ok.
  11
  12 blastpath = "blastall"
  13 # path of blastall.
  14 # if blastall is in your command path, "blastall" is ok.
  15
  16 # mafft-homologs.rb  v. 2.1 aligns sequences together with homologues
  17 # automatically collected from SwissProt via NCBI BLAST.
  18 #
  19 # mafft > 5.58 is required
  20 #
  21 # Usage:
  22 #   mafft-homologs.rb [options] input > output
  23 # Options:
  24 #   -a #      the number of collected sequences (default: 50)
  25 #   -e #      threshold value (default: 1e-10)
  26 #   -o "xxx"  options for mafft
  27 #             (default: " --op 1.53 --ep 0.123 --maxiterate 1000")
  28 #   -l        locally carries out blast searches instead of NCBI blast
  29 #             (requires locally installed blast and a database)
  30 #   -f        outputs collected homologues also (default: off)
  31 #   -w        entire sequences are subjected to BLAST search
  32 #             (default: well-aligned region only)
  33
  34 #require 'getopts'
  35 require 'optparse'
  36 require 'tempfile'
  37
  38 # mktemp
  39 GC.disable
  40 temp_vf = Tempfile.new("_vf").path
  41 temp_if = Tempfile.new("_if").path
  42 temp_pf = Tempfile.new("_pf").path
  43 temp_af = Tempfile.new("_af").path
  44 temp_qf = Tempfile.new("_qf").path
  45 temp_bf = Tempfile.new("_bf").path
  46 temp_rid = Tempfile.new("_rid").path
  47 temp_res = Tempfile.new("_res").path
  48
  49
  50 system( mafftpath + " --help > #{temp_vf} 2>&1" )
  51 pfp = File.open( "#{temp_vf}", 'r' )
  52 while pfp.gets
  53         break if $_ =~ /MAFFT v/
  54 end
  55 pfp.close
  56
  57 if( $_ ) then
  58         mafftversion = $_.sub( /^\D*/, "" ).split(" ").slice(0).strip.to_s
  59 else
  60         mafftversion = "0"
  61 end
  62 if( mafftversion < "5.58" ) then
  63         STDERR.puts ""
  64         STDERR.puts "======================================================"
  65         STDERR.puts "Install new mafft (v. >= 5.58)"
  66         STDERR.puts "======================================================"
  67         STDERR.puts ""
  68         exit
  69 end
  70
  71 srand ( 0 )
  72
  73 def readfasta( fp, name, seq )
  74         nseq = 0
  75         tmpseq = ""
  76         while fp.gets
  77                 if $_ =~ /^>/ then
  78                         name.push( $_.sub(/>/,"").strip )
  79                         seq.push( tmpseq ) if nseq > 0
  80                         nseq += 1
  81                         tmpseq = ""
  82                 else
  83                         tmpseq += $_.strip
  84                 end
  85         end
  86         seq.push( tmpseq )
  87         return nseq
  88 end
  89
  90 nadd = 50
  91 eval = 1e-10
  92 local = 0
  93 fullout = 0
  94 entiresearch = 0
  95 corewin = 50
  96 corethr = 0.3
  97 mafftopt = " --op 1.53 --ep 0.123 --localpair --maxiterate 1000 --reorder "
  98
  99
 100 #if getopts( "s", "f", "w", "l", "h", "e:", "a:", "o:", "c:", "d:" ) == nil ||  ARGV.length == 0 || $OPT_h then
 101 #       puts "Usage: #{$0} [-h -l -e# -a# -o\"[options for mafft]\"] input_file"
 102 #       exit
 103 #end
 104 params = ARGV.getopts( "sfwlhe:a:o:c:d:" )
 105
 106
 107 #if $OPT_c then
 108 if params["c"] != nil then
 109         corewin = params["c"].to_i
 110 end
 111
 112 #if $OPT_d then
 113 if params["d"] != nil then
 114         corethr = params["d"].to_f
 115 end
 116
 117 #if $OPT_w
 118 if params["w"] == true then
 119         entiresearch = 1
 120 end
 121
 122 #if $OPT_f
 123 if params["f"] == true then
 124         fullout = 1
 125 end
 126
 127 #if $OPT_s
 128 if params["s"] == true then
 129         fullout = 0
 130 end
 131
 132 #if $OPT_l
 133 if params["l"] == true then
 134         local = 1
 135 end
 136
 137 #if $OPT_e then
 138 if params["e"] != nil then
 139 #       eval = $OPT_e.to_f
 140         eval = params["e"].to_f
 141 end
 142
 143 #if $OPT_a then
 144 if params["a"] != nil then
 145         nadd = params["a"].to_i
 146 end
 147
 148 #if $OPT_o then
 149 if params["o"] != nil then
 150         mafftopt += " " + params["o"] + " "
 151 end
 152
 153 infn = ARGV[0].to_s.strip
 154
 155 system "cat " + infn + " > #{temp_if}"
 156 ar = mafftopt.split(" ")
 157 nar = ar.length
 158 for i in 0..(nar-1)
 159         if ar[i] == "--seed" then
 160                 system "cat #{ar[i+1]} >> #{temp_if}"
 161         end
 162 end
 163
 164 nseq = 0
 165 ifp = File.open( "#{temp_if}", 'r' )
 166         while ifp.gets
 167                 nseq += 1 if $_ =~ /^>/
 168         end
 169 ifp.close
 170
 171 if nseq >= 100 then
 172         STDERR.puts "The number of input sequences must be <100."
 173         exit
 174 elsif nseq == 1 then
 175         system( "cp #{temp_if}"  + " #{temp_pf}" )
 176 else
 177         STDERR.puts "Performing preliminary alignment .. "
 178         if entiresearch == 1 then
 179 #               system( mafftpath + " --maxiterate 1000 --localpair #{temp_if} > #{temp_pf}" )
 180                 system( mafftpath + " --maxiterate 0 --retree 2 #{temp_if} > #{temp_pf}" )
 181         else
 182                 system( mafftpath + " --maxiterate 1000 --localpair --core --coreext --corethr #{corethr.to_s} --corewin #{corewin.to_s} #{temp_if} > #{temp_pf}" )
 183         end
 184 end
 185
 186 pfp = File.open( "#{temp_pf}", 'r' )
 187 inname = []
 188 inseq = []
 189 slen = []
 190 act = []
 191 nin = 0
 192 nin = readfasta( pfp, inname, inseq )
 193 for i in 0..(nin-1)
 194         slen.push( inseq[i].gsub(/-/,"").length )
 195         act.push( 1 )
 196 end
 197 pfp.close
 198
 199 pfp = File.open( "#{temp_if}", 'r' )
 200 orname = []
 201 orseq = []
 202 nin = 0
 203 nin = readfasta( pfp, orname, orseq )
 204 pfp.close
 205
 206 allen = inseq[0].length
 207 for i in 0..(nin-2)
 208         for j in (i+1)..(nin-1)
 209                 next if act[i] == 0
 210                 next if act[j] == 0
 211                 pid = 0.0
 212                 total = 0
 213                 for a in 0..(allen-1)
 214                         next if inseq[i][a,1] == "-" || inseq[j][a,1] == "-"
 215                         total += 1
 216                         pid += 1.0 if inseq[i][a,1] == inseq[j][a,1]
 217                 end
 218                 pid /= total
 219 #               puts "#{i.to_s}, #{j.to_s}, #{pid.to_s}"
 220                 if pid > 0.5 then
 221                         if slen[i] < slen[j]
 222                                 act[i] = 0
 223                         else
 224                                 act[j] = 0
 225                         end
 226                 end
 227         end
 228 end
 229 #p act
 230
 231
 232 afp = File.open( "#{temp_af}", 'w' )
 233
 234 STDERR.puts "Searching .. \n"
 235 ids = []
 236 add = []
 237 sco = []
 238 for i in 0..(nin-1)
 239         inseq[i].gsub!(/-/,"")
 240         afp.puts ">" + orname[i]
 241         afp.puts orseq[i]
 242
 243 #       afp.puts ">" + inname[i]
 244 #       afp.puts inseq[i]
 245
 246         STDERR.puts "Query (#{i+1}/#{nin})\n" + inname[i]
 247         if act[i] == 0 then
 248                 STDERR.puts "Skip.\n\n"
 249                 next
 250         end
 251
 252         if local == 0 then
 253                 command = "lynx -source 'https://www.ncbi.nlm.nih.gov/blast/Blast.cgi?QUERY=" + inseq[i] + "&DATABASE=swissprot&HITLIST_SIZE=" + nadd.to_s + "&FILTER=L&EXPECT='" + eval.to_s + "'&FORMAT_TYPE=TEXT&PROGRAM=blastp&SERVICE=plain&NCBI_GI=on&PAGE=Proteins&CMD=Put' > #{temp_rid}"
 254                 system command
 255
 256                 ridp = File.open( "#{temp_rid}", 'r' )
 257                 while ridp.gets
 258                         break if $_ =~ / RID = (.*)/
 259                 end
 260                 ridp.close
 261                 rid = $1.strip
 262                 STDERR.puts "Submitted to NCBI. rid = " + rid
 263
 264                 STDERR.printf "Waiting "
 265                 while 1
 266                         STDERR.printf "."
 267                         sleep 10
 268                         command = "lynx -source 'https://www.ncbi.nlm.nih.gov/blast/Blast.cgi?RID=" + rid + "&DESCRIPTIONS=500&ALIGNMENTS=" + nadd.to_s + "&ALIGNMENT_TYPE=Pairwise&OVERVIEW=no&CMD=Get&FORMAT_TYPE=XML' > #{temp_res}"
 269                         system command
 270                         resp = File.open( "#{temp_res}", 'r' )
 271 #                       resp.gets
 272 #                       if $_ =~ /WAITING/ then
 273 #                               resp.close
 274 #                               next
 275 #                       end
 276                         while( resp.gets )
 277                                 break if $_ =~ /QBlastInfoBegin/
 278                         end
 279                         resp.gets
 280                         if $_ =~ /WAITING/ then
 281                                 resp.close
 282                                 next
 283                         else
 284                                 resp.close
 285                                 break
 286                         end
 287                 end
 288         else
 289 #               puts "Not supported"
 290 #               exit
 291                 qfp = File.open( "#{temp_qf}", 'w' )
 292                         qfp.puts "> "
 293                         qfp.puts inseq[i]
 294                 qfp.close
 295                 command = blastpath + "  -p blastp  -e #{eval} -b 1000 -m 7 -i #{temp_qf} -d #{localdb} > #{temp_res}"
 296                 system command
 297                 resp = File.open( "#{temp_res}", 'r' )
 298         end
 299         STDERR.puts " Done.\n\n"
 300
 301         resp = File.open( "#{temp_res}", 'r' )
 302         while 1
 303                 while resp.gets
 304                         break if $_ =~ /<Hit_id>(.*)<\/Hit_id>/ || $_ =~ /(<Iteration_stat>)/
 305                 end
 306                 id = $1
 307                 break if $_ =~ /<Iteration_stat>/
 308 #               p id
 309                 while resp.gets
 310                         break if $_ =~ /<Hsp_bit-score>(.*)<\/Hsp_bit-score>/
 311                 end
 312                 score = $1.to_f
 313 #               p score
 314
 315                 known = ids.index( id )
 316                 if known != nil then
 317                         if sco[known] >= score then
 318                                 next
 319                         else
 320                                 ids.delete_at( known )
 321                                 add.delete_at( known )
 322                                 sco.delete_at( known )
 323                         end
 324                 end
 325                 while resp.gets
 326                         break if $_ =~ /<Hsp_hseq>(.*)<\/Hsp_hseq>/
 327                 end
 328 #               break if $1 == nil
 329                 target = $1.sub( /-/, "" ).sub( /U/, "X" )
 330 #               p target
 331 #               STDERR.puts "adding 1 seq"
 332                 ids.push( id )
 333                 sco.push( score )
 334                 add.push( target )
 335         end
 336         resp.close
 337 end
 338
 339 n = ids.length
 340
 341 outnum = 0
 342 while n > 0 && outnum < nadd
 343         m = rand( n )
 344         afp.puts ">_addedbymaffte_" + ids[m]
 345         afp.puts add[m]
 346         ids.delete_at( m )
 347         add.delete_at( m )
 348         n -= 1
 349         outnum += 1
 350 end
 351 afp.close
 352
 353 STDERR.puts "Performing alignment .. "
 354 system( mafftpath + mafftopt + " #{temp_af} > #{temp_bf}" )
 355 STDERR.puts "done."
 356
 357 bfp = File.open( "#{temp_bf}", 'r' )
 358 outseq = []
 359 outnam = []
 360 readfasta( bfp, outnam, outseq )
 361 bfp.close
 362
 363 outseq2 = []
 364 outnam2 = []
 365
 366 len = outseq.length
 367 for i in 0..(len-1)
 368 #       p outnam[i]
 369         if fullout == 0 && outnam[i] =~ /_addedbymaffte_/ then
 370                 next
 371         end
 372         outseq2.push( outseq[i] )
 373         outnam2.push( outnam[i].sub( /_addedbymaffte_/, "_ho_" ) )
 374 end
 375
 376 nout = outseq2.length
 377 len = outseq[0].length
 378 p = len
 379 while p>0
 380         p -= 1
 381     allgap = 1
 382     for j in 0..(nout-1)
 383                 if outseq2[j][p,1] != "-" then
 384                         allgap = 0
 385                         break
 386                 end
 387     end
 388     if allgap == 1 then
 389         for j in 0..(nout-1)
 390             outseq2[j][p,1] = ""
 391         end
 392     end
 393 end
 394 for i in 0..(nout-1)
 395         puts ">" + outnam2[i]
 396         puts outseq2[i].gsub( /.{1,60}/, "\\0\n" )
 397 end
 398
 399
 400 system( "rm -rf #{temp_if} #{temp_vf} #{temp_af} #{temp_bf} #{temp_pf} #{temp_qf} #{temp_res} #{temp_rid}" )
 401 if File.exist?( "#{temp_af}.tree" ) then
 402         system( "sed 's/_addedbymaffte_/_ho_/'  #{temp_af}.tree > #{ARGV[0].to_s}.tree" )
 403         system( "rm #{temp_af}.tree" )
 404 end