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 'tempfile'
  36
  37 # mktemp
  38 GC.disable
  39 temp_vf = Tempfile.new("_vf").path
  40 temp_if = Tempfile.new("_if").path
  41 temp_pf = Tempfile.new("_pf").path
  42 temp_af = Tempfile.new("_af").path
  43 temp_qf = Tempfile.new("_qf").path
  44 temp_bf = Tempfile.new("_bf").path
  45 temp_rid = Tempfile.new("_rid").path
  46 temp_res = Tempfile.new("_res").path
  47
  48
  49 system( mafftpath + " --help > #{temp_vf} 2>&1" )
  50 pfp = File.open( "#{temp_vf}", 'r' )
  51 while pfp.gets
  52         break if $_ =~ /MAFFT v/
  53 end
  54 pfp.close
  55 if( $_ ) then
  56         mafftversion = sub( /^\D*/, "" ).split(" ").slice(0).strip.to_s
  57 else
  58         mafftversion = "0"
  59 end
  60 if( mafftversion < "5.58" ) then
  61         puts ""
  62         puts "======================================================"
  63         puts "Install new mafft (v. >= 5.58)"
  64         puts "======================================================"
  65         puts ""
  66         exit
  67 end
  68
  69 srand ( 0 )
  70
  71 def readfasta( fp, name, seq )
  72         nseq = 0
  73         tmpseq = ""
  74         while fp.gets
  75                 if $_ =~ /^>/ then
  76                         name.push( $_.sub(/>/,"").strip )
  77                         seq.push( tmpseq ) if nseq > 0
  78                         nseq += 1
  79                         tmpseq = ""
  80                 else
  81                         tmpseq += $_.strip
  82                 end
  83         end
  84         seq.push( tmpseq )
  85         return nseq
  86 end
  87
  88 nadd = 50
  89 eval = 1e-10
  90 local = 0
  91 fullout = 0
  92 entiresearch = 0
  93 corewin = 50
  94 corethr = 0.3
  95 mafftopt = " --op 1.53 --ep 0.123 --localpair --maxiterate 1000 --reorder "
  96 if getopts( "s", "f", "w", "l", "h", "e:", "a:", "o:", "c:", "d:" ) == nil ||  ARGV.length == 0 || $OPT_h then
  97         puts "Usage: #{$0} [-h -l -e# -a# -o\"[options for mafft]\"] input_file"
  98         exit
  99 end
 100
 101 if $OPT_c then
 102         corewin = $OPT_c.to_i
 103 end
 104 if $OPT_d then
 105         corethr = $OPT_d.to_f
 106 end
 107 if $OPT_w
 108         entiresearch = 1
 109 end
 110 if $OPT_f
 111         fullout = 1
 112 end
 113 if $OPT_s
 114         fullout = 0
 115 end
 116 if $OPT_l
 117         local = 1
 118 end
 119 if $OPT_e then
 120         eval = $OPT_e.to_f
 121 end
 122 if $OPT_a then
 123         nadd = $OPT_a.to_i
 124 end
 125 if $OPT_o then
 126         mafftopt += " " + $OPT_o + " "
 127 end
 128
 129 system "cat " + ARGV.to_s + " > #{temp_if}"
 130 ar = mafftopt.split(" ")
 131 nar = ar.length
 132 for i in 0..(nar-1)
 133         if ar[i] == "--seed" then
 134                 system "cat #{ar[i+1]} >> #{temp_if}"
 135         end
 136 end
 137
 138 nseq = 0
 139 ifp = File.open( "#{temp_if}", 'r' )
 140         while ifp.gets
 141                 nseq += 1 if $_ =~ /^>/
 142         end
 143 ifp.close
 144
 145 if nseq >= 100 then
 146         STDERR.puts "The number of input sequences must be <100."
 147         exit
 148 elsif nseq == 1 then
 149         system( "cp #{temp_if}"  + " #{temp_pf}" )
 150 else
 151         STDERR.puts "Performing preliminary alignment .. "
 152         if entiresearch == 1 then
 153 #               system( mafftpath + " --maxiterate 1000 --localpair #{temp_if} > #{temp_pf}" )
 154                 system( mafftpath + " --maxiterate 0 --retree 2 #{temp_if} > #{temp_pf}" )
 155         else
 156                 system( mafftpath + " --maxiterate 1000 --localpair --core --coreext --corethr #{corethr.to_s} --corewin #{corewin.to_s} #{temp_if} > #{temp_pf}" )
 157         end
 158 end
 159
 160 pfp = File.open( "#{temp_pf}", 'r' )
 161 inname = []
 162 inseq = []
 163 slen = []
 164 act = []
 165 nin = 0
 166 nin = readfasta( pfp, inname, inseq )
 167 for i in 0..(nin-1)
 168         slen.push( inseq[i].gsub(/-/,"").length )
 169         act.push( 1 )
 170 end
 171 pfp.close
 172
 173 pfp = File.open( "#{temp_if}", 'r' )
 174 orname = []
 175 orseq = []
 176 nin = 0
 177 nin = readfasta( pfp, orname, orseq )
 178 pfp.close
 179
 180 allen = inseq[0].length
 181 for i in 0..(nin-2)
 182         for j in (i+1)..(nin-1)
 183                 next if act[i] == 0
 184                 next if act[j] == 0
 185                 pid = 0.0
 186                 total = 0
 187                 for a in 0..(allen-1)
 188                         next if inseq[i][a,1] == "-" || inseq[j][a,1] == "-"
 189                         total += 1
 190                         pid += 1.0 if inseq[i][a,1] == inseq[j][a,1]
 191                 end
 192                 pid /= total
 193 #               puts "#{i.to_s}, #{j.to_s}, #{pid.to_s}"
 194                 if pid > 0.5 then
 195                         if slen[i] < slen[j]
 196                                 act[i] = 0
 197                         else
 198                                 act[j] = 0
 199                         end
 200                 end
 201         end
 202 end
 203 #p act
 204
 205
 206 afp = File.open( "#{temp_af}", 'w' )
 207
 208 STDERR.puts "Searching .. \n"
 209 ids = []
 210 add = []
 211 sco = []
 212 for i in 0..(nin-1)
 213         inseq[i].gsub!(/-/,"")
 214         afp.puts ">" + orname[i]
 215         afp.puts orseq[i]
 216
 217 #       afp.puts ">" + inname[i]
 218 #       afp.puts inseq[i]
 219
 220         STDERR.puts "Query (#{i+1}/#{nin})\n" + inname[i]
 221         if act[i] == 0 then
 222                 STDERR.puts "Skip.\n\n"
 223                 next
 224         end
 225
 226         if local == 0 then
 227                 command = "lynx -source 'http://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}"
 228                 system command
 229
 230                 ridp = File.open( "#{temp_rid}", 'r' )
 231                 while ridp.gets
 232                         break if $_ =~ / RID = (.*)/
 233                 end
 234                 ridp.close
 235                 rid = $1.strip
 236                 STDERR.puts "Submitted to NCBI. rid = " + rid
 237
 238                 STDERR.printf "Waiting "
 239                 while 1
 240                         STDERR.printf "."
 241                         sleep 10
 242                         command = "lynx -source 'http://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}"
 243                         system command
 244                         resp = File.open( "#{temp_res}", 'r' )
 245 #                       resp.gets
 246 #                       if $_ =~ /WAITING/ then
 247 #                               resp.close
 248 #                               next
 249 #                       end
 250                         while( resp.gets )
 251                                 break if $_ =~ /QBlastInfoBegin/
 252                         end
 253                         resp.gets
 254                         if $_ =~ /WAITING/ then
 255                                 resp.close
 256                                 next
 257                         else
 258                                 resp.close
 259                                 break
 260                         end
 261                 end
 262         else
 263 #               puts "Not supported"
 264 #               exit
 265                 qfp = File.open( "#{temp_qf}", 'w' )
 266                         qfp.puts "> "
 267                         qfp.puts inseq[i]
 268                 qfp.close
 269                 command = blastpath + "  -p blastp  -e #{eval} -b 1000 -m 7 -i #{temp_qf} -d #{localdb} > #{temp_res}"
 270                 system command
 271                 resp = File.open( "#{temp_res}", 'r' )
 272         end
 273         STDERR.puts " Done.\n\n"
 274
 275         resp = File.open( "#{temp_res}", 'r' )
 276         while 1
 277                 while resp.gets
 278                         break if $_ =~ /<Hit_id>(.*)<\/Hit_id>/ || $_ =~ /(<Iteration_stat>)/
 279                 end
 280                 id = $1
 281                 break if $_ =~ /<Iteration_stat>/
 282 #               p id
 283                 while resp.gets
 284                         break if $_ =~ /<Hsp_bit-score>(.*)<\/Hsp_bit-score>/
 285                 end
 286                 score = $1.to_f
 287 #               p score
 288
 289                 known = ids.index( id )
 290                 if known != nil then
 291                         if sco[known] >= score then
 292                                 next
 293                         else
 294                                 ids.delete_at( known )
 295                                 add.delete_at( known )
 296                                 sco.delete_at( known )
 297                         end
 298                 end
 299                 while resp.gets
 300                         break if $_ =~ /<Hsp_hseq>(.*)<\/Hsp_hseq>/
 301                 end
 302 #               break if $1 == nil
 303                 target = $1.sub( /-/, "" )
 304 #               p target
 305 #               STDERR.puts "adding 1 seq"
 306                 ids.push( id )
 307                 sco.push( score )
 308                 add.push( target )
 309         end
 310         resp.close
 311 end
 312
 313 n = ids.length
 314
 315 outnum = 0
 316 while n > 0 && outnum < nadd
 317         m = rand( n )
 318         afp.puts ">_addedbymaffte_" + ids[m]
 319         afp.puts add[m]
 320         ids.delete_at( m )
 321         add.delete_at( m )
 322         n -= 1
 323         outnum += 1
 324 end
 325 afp.close
 326
 327 STDERR.puts "Performing alignment .. "
 328 system( mafftpath + mafftopt + " #{temp_af} > #{temp_bf}" )
 329 STDERR.puts "done."
 330
 331 bfp = File.open( "#{temp_bf}", 'r' )
 332 outseq = []
 333 outnam = []
 334 readfasta( bfp, outnam, outseq )
 335 bfp.close
 336
 337 outseq2 = []
 338 outnam2 = []
 339
 340 len = outseq.length
 341 for i in 0..(len-1)
 342 #       p outnam[i]
 343         if fullout == 0 && outnam[i] =~ /^_addedbymaffte_/ then
 344                 next
 345         end
 346         outseq2.push( outseq[i] )
 347         outnam2.push( outnam[i].sub( /^_addedbymaffte_/, "_ho_" ) )
 348 end
 349
 350 nout = outseq2.length
 351 len = outseq[0].length
 352 p = len
 353 while p>0
 354         p -= 1
 355     allgap = 1
 356     for j in 0..(nout-1)
 357                 if outseq2[j][p,1] != "-" then
 358                         allgap = 0
 359                         break
 360                 end
 361     end
 362     if allgap == 1 then
 363         for j in 0..(nout-1)
 364             outseq2[j][p,1] = ""
 365         end
 366     end
 367 end
 368 for i in 0..(nout-1)
 369         puts ">" + outnam2[i]
 370         puts outseq2[i].gsub( /.{1,60}/, "\\0\n" )
 371 end
 372
 373
 374 system( "rm -rf #{temp_if} #{temp_vf} #{temp_af} #{temp_bf} #{temp_pf} #{temp_qf} #{temp_res} #{temp_rid}" )