| 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( /-/, "" ).sub( /U/, "X" ) |
|---|
| 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}" ) |
|---|
