diff --git a/fur/fur.go b/fur/fur.go index d473db5..e418817 100644 --- a/fur/fur.go +++ b/fur/fur.go @@ -61,6 +61,8 @@ func main() { optMM := flag.Bool("M", false, "activate masking (recommended for mammalian genomes)") optN := flag.Int("n", 100, "number of nucleotides in region") + optWW := flag.Int("W", 0, "word size for megablast mode " + + "(implies -m)") flag.Parse() if *optV { util.PrintInfo("fur") @@ -101,8 +103,18 @@ func main() { (*optT) = ncpu } if *optF > 1 || *optF <= 0 { - log.Fatalf("can't use %f as a sensitivity threshold\n" + - "please use a positive value not exceeding 1", *optF) + m := "can't use %f as a sensitivity threshold\n" + + "please use a positive value not exceeding 1" + log.Fatalf(m, *optF) + } + if *optWW > 0 { + if *optWW < 4 { + m := "couldn't set the Blast word size " + + "to %d; " + + "please use a word size of >= 4" + log.Fatalf(m, *optWW) + } + (*optM) = true } regions := make([]*fasta.Sequence, 0) rw := tabwriter.NewWriter(os.Stderr, 0, 0, 2, ' ', @@ -387,10 +399,11 @@ func main() { if len(regions) > 0 { cmds := make([]*exec.Cmd, 0) da := *optD + "/n" - th := *optT - ev := *optE + th := strconv.Itoa(*optT) + ev := fmt.Sprintf("%g", *optE) ta := "megablast" ma := "" + ws := "" if *optMM { cmd := exec.Command("blastdbcmd", "-info", "-db", *optD + "/n") out, err := cmd.CombinedOutput() @@ -409,32 +422,38 @@ func main() { fmt.Fprintf(os.Stderr, m) } } + if *optWW > 0 { + ws = strconv.Itoa(*optWW) + } of := "6 qaccver qstart qend" - tm := "blastn -db %s -num_threads %d " - tm += "-evalue %g -task %s " - if *optMM && ma != "" { - tm += "-db_soft_mask %s " + args := []string{ + "-db", da, + "-num_threads", th, + "-evalue", ev, + "-task", ta, } - tm += "-outfmt " - as := fmt.Sprintf(tm, da, th, ev, ta) if *optMM && ma != "" { - as = fmt.Sprintf(tm, da, th, ev, ta, ma) + args = append(args, "-db_soft_mask", ma) } - args := strings.Fields(as) - args = append(args, of) - cmd := exec.Command("blastn") - cmd.Args = args + if *optWW > 0 { + args = append(args, "-word_size", ws) + } + args = append(args, "-outfmt", of) + cmd := exec.Command("blastn", args...) cmds = append(cmds, cmd) if !*optM { ta = "blastn" - as := fmt.Sprintf(tm, da, th, ev, ta) + args := []string{ + "-db", da, + "-num_threads", th, + "-evalue", ev, + "-task", ta, + } if *optMM && ma != "" { - as = fmt.Sprintf(tm, da, th, ev, ta, ma) + args = append(args, "-db_soft_mask", ma) } - args = strings.Fields(as) - args = append(args, of) - cmd = exec.Command("blastn") - cmd.Args = args + args = append(args, "-outfmt", of) + cmd := exec.Command("blastn", args...) cmds = append(cmds, cmd) } for _, cmd := range cmds { diff --git a/fur/fur.org b/fur/fur.org index 944568b..21c3d60 100644 --- a/fur/fur.org +++ b/fur/fur.org @@ -120,9 +120,10 @@ The user can also print the result of Step~(\ref{eq:fur2}) and exit. #+begin_export latex Step~(\ref{eq:fur3}) is implemented in \ty{blastn}, where we expose the E-value, the number of threads, a switch to megablast mode instead -of the default blastn mode, and a switch for running with masking. We -also set the minimum length of the final regions at this point. The -number of threads is initialized to the number of CPUs. +of the default blastn mode, a switch for running with masking, and the +word size for megablast. We also set the minimum length of the final +regions at this point. The number of threads is initialized to the +number of CPUs. #+end_export #+begin_src go <>= optE := flag.Float64("e", 1e-5, "E-value for Blast") @@ -134,6 +135,8 @@ number of threads is initialized to the number of CPUs. optMM := flag.Bool("M", false, "activate masking (recommended for mammalian genomes)") optN := flag.Int("n", 100, "number of nucleotides in region") + optWW := flag.Int("W", 0, "word size for megablast mode " + + "(implies -m)") #+end_src #+begin_export latex We import \ty{runtime}. @@ -144,8 +147,8 @@ We import \ty{runtime}. #+begin_export latex We parse the options and first respond to the version request, \ty{-v}, as this might stop the program. We also respond to the -database name, \ty{-d}, the number of threads, \ty{-t}, and the -sensitivity threshold, \ty{-f}. +database name, \ty{-d}, the number of threads, \ty{-t}, the +sensitivity threshold, \ty{-f}, and the Blast word length, \ty{-W}. #+end_export #+begin_src go <>= flag.Parse() @@ -153,6 +156,7 @@ sensitivity threshold, \ty{-f}. //<> //<> //<> + //<> #+end_src #+begin_export latex We import \ty{fmt}. @@ -250,8 +254,28 @@ sets $f$ outside the range, we notify the user and exit. #+end_export #+begin_src go <>= if *optF > 1 || *optF <= 0 { - log.Fatalf("can't use %f as a sensitivity threshold\n" + - "please use a positive value not exceeding 1", *optF) + m := "can't use %f as a sensitivity threshold\n" + + "please use a positive value not exceeding 1" + log.Fatalf(m, *optF) + } +#+end_src +#+begin_export latex +The user can override megablast-specific Blast word size by explicitly +setting it with the \ty{-W} flag. Nucleotide Blast allows the word +size of 4 and greater, so we check the provided value. If it is +invalid, we politely ask the user to set a valid one. Setting the word +size only makes sense in megablast mode, so if \ty{-W} is correctly +set, we also set \ty{-m}. +#+end_export +#+begin_src go <>= + if *optWW > 0 { + if *optWW < 4 { + m := "couldn't set the Blast word size " + + "to %d; " + + "please use a word size of >= 4" + log.Fatalf(m, *optWW) + } + (*optM) = true } #+end_src #+begin_export latex @@ -1129,36 +1153,41 @@ We import \ty{exec}. "os/exec" #+end_src #+begin_export latex -We construct the Blast options, construct a template for the Blast -commands, and construct the megablast command. Then we construct the -blastn command, unless the user opted for megablast only. +We construct the Blast options and construct the megablast +command. Then we construct the blastn command, unless the user opted +for megablast only. #+end_export #+begin_src go <>= //<> - //<> //<> if !*optM { //<> } #+end_src #+begin_export latex -We set six options in Blast, the values of which we first need to +We set seven options in Blast, the values of which we first need to construct. The options are the database path, the number of threads, -the E-value, the task, the masking algorithm, and the output -format. The masking algorithm is initialized as an empty string which -we fill if the user requested masking. As output format we set the -query accession, start, and end---the coordinates for homology masking -later on. +the E-value, the task, the masking algorithm, the word size, and the +output format. The masking algorithm is initialized as an empty +string, which we fill if the user requested masking. The word size is +initialized as an empty string, which we set if the user has provided +the word length. As output format we set the query accession, start, +and end---the coordinates for homology masking later on. We convert +all numerical options values to strings. #+end_export #+begin_src go <>= da := *optD + "/n" - th := *optT - ev := *optE + th := strconv.Itoa(*optT) + ev := fmt.Sprintf("%g", *optE) ta := "megablast" ma := "" + ws := "" if *optMM { //<> } + if *optWW > 0 { + ws = strconv.Itoa(*optWW) + } of := "6 qaccver qstart qend" #+end_src #+begin_export latex @@ -1210,46 +1239,43 @@ database. In that case we should warn the user. } #+end_src #+begin_export latex -The Blast template has space for five Blast options, six with -masking. The output format is a composite string, so we append it -later to the arguments slice. -#+end_export -#+begin_src go <>= - tm := "blastn -db %s -num_threads %d " - tm += "-evalue %g -task %s " - if *optMM && ma != "" { - tm += "-db_soft_mask %s " - } - tm += "-outfmt " -#+end_src -#+begin_export latex We generate the arguments for the Blast command, append the output format, set the arguments, and store the command. #+end_export #+begin_src go <>= - as := fmt.Sprintf(tm, da, th, ev, ta) + args := []string{ + "-db", da, + "-num_threads", th, + "-evalue", ev, + "-task", ta, + } if *optMM && ma != "" { - as = fmt.Sprintf(tm, da, th, ev, ta, ma) + args = append(args, "-db_soft_mask", ma) } - args := strings.Fields(as) - args = append(args, of) - cmd := exec.Command("blastn") - cmd.Args = args + if *optWW > 0 { + args = append(args, "-word_size", ws) + } + args = append(args, "-outfmt", of) + cmd := exec.Command("blastn", args...) cmds = append(cmds, cmd) #+end_src #+begin_export latex -We repeat this for the blastn command, only with a different task. +We repeat this for the blastn command, only with a different task and +without word size. #+end_export #+begin_src go <>= ta = "blastn" - as := fmt.Sprintf(tm, da, th, ev, ta) + args := []string{ + "-db", da, + "-num_threads", th, + "-evalue", ev, + "-task", ta, + } if *optMM && ma != "" { - as = fmt.Sprintf(tm, da, th, ev, ta, ma) + args = append(args, "-db_soft_mask", ma) } - args = strings.Fields(as) - args = append(args, of) - cmd = exec.Command("blastn") - cmd.Args = args + args = append(args, "-outfmt", of) + cmd := exec.Command("blastn", args...) cmds = append(cmds, cmd) #+end_src #+begin_export latex @@ -1531,12 +1557,13 @@ We import \ty{exec}. #+end_src #+begin_export latex We construct a set of tests without repeat masking, one with repeat -masking, and one with partial intersection. +masking, one with partial intersection, and two tests with word size. #+end_export #+begin_src go <>= //<> //<> //<> + //<> #+end_src #+begin_export latex We construct four tests without repeat masking. The first runs with @@ -1593,6 +1620,16 @@ five targets lacks the marker. We set the sensitivity threshold to tests = append(tests, test) #+end_src #+begin_export latex +We run \ty{fur} in implied megablast mode by setting the word size to +4. In this way, we adjust megablast to be actually more stringent than +blastn. +#+end_export +#+begin_src go <>= + d = "test.db" + test = exec.Command("./fur", "-d", d, "-W", "4") + tests = append(tests, test) +#+end_src +#+begin_export latex For each test we compare the result we get with the result we want, which is contained in files \ty{r1.txt}, \ty{r2.txt}, and so on. #+end_export diff --git a/fur/fur_test.go b/fur/fur_test.go index ab1b4ac..c854fbe 100644 --- a/fur/fur_test.go +++ b/fur/fur_test.go @@ -32,6 +32,9 @@ func TestFur(t *testing.T) { d = "testPartial.db" test = exec.Command("./fur", "-d", d, "-f", "0.8") tests = append(tests, test) + d = "test.db" + test = exec.Command("./fur", "-d", d, "-W", "4") + tests = append(tests, test) for i, test := range tests { get, err := test.CombinedOutput() if err != nil { diff --git a/fur/makeRes.sh b/fur/makeRes.sh index 839a1e5..6576383 100644 --- a/fur/makeRes.sh +++ b/fur/makeRes.sh @@ -5,3 +5,5 @@ ./fur -d test.db -M &> r5.txt ./fur -d masked.db &> r6.txt ./fur -d masked.db -M &> r7.txt +./fur -d testPartial.db -f 0.8 &> r8.txt +./fur -d test.db -m -W 4 &> r9.txt diff --git a/fur/r9.txt b/fur/r9.txt new file mode 100644 index 0000000..e7db604 --- /dev/null +++ b/fur/r9.txt @@ -0,0 +1,21 @@ + Step Sequences Length Ns + ------------- --------- ------ -- + Subtraction_1 1 1052 0 + Intersection 1 1051 5 + Subtraction_2 1 1000 5 +>t1_1 +TGTGGCTCTGGGAGGTCCGACCAAGTGCTACTGCGGCATTAGGGCGCAGCGATCCCCTCATGCGGCAATG +AGGTCATATAATACTACGTTATATTTTGAAGATGAATAANGACTTCCGTCAAATCAGCAGTCAGTAAATA +TTGAACCANTTTAAACCGGCGCCGCGCGGTAACCTAGGACTCTCCATGCTAAACCCCAGCGCCGTAGAGT +GATCCCTAAGCTAATCAAACCGTTCGTTTGTGCCCTATGTCAGGAAGCACGACGTCGGTAATCGAGTCAC +CTGTGCCCCGTTTATCCCAGATGTTAGGATAACATTCGTACGACGTAAAGTTAGACATGGCGCAGGNCAC +GCCTGCAAACGTTGTTCCCCTTCATGTTTGCACTACTCTTAGTCGCCCTCCCTTTTCCGCGTCTCCAGGG +CGCGCCCCTCAGGACTGTCTTTGCGGTGGCTCGCTCCAGCTGGCTGAGCATTATAAGAGGATACTATTAA +TTAGTCTTGTTAACTCGCGGTTTGCGAAAGCGTATGTAGTCTGTTGTTTGCCGGGGTATTTGAGATAATC +ACACCTGGGCTCTAGTGGCTCTCAGCAAAATCGGGATTGGCCGNGAAAAACCTAACTAACTTCTGTTTGC +CGCACGTCTATTAATCAACCGCATATCGGGGACCCAGGATTCTCCTCAGACAGCAATAACGCGGACGGAG +AGATCTAATGGCAGTGCCTCGGGCTCTTGGAAAGTGAGTCAACTTNTTCCCGACAAATTAAACAGGTGGC +CACCGTAACTGCGGCATCTACGGAACATGACTCGACACTAGTGACAGATGGACTGCACGACTTCGCGGCG +GTATCGGGCTATTTTTTACGTACGTTAAACGGAAAGGGACTCTAGAGCCGTGTCAGGCTGACTGACTTCC +CGAGCCAACGCGCTTAACCGAACGCCACAGACAAACTGTTAAGGCTTGCATGGAGTATACATGGACAAGG +ATCCCAACTGAATCTGAAAG