Word size for megablast

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 {

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 <<Options for step (\ref{eq:fur3}), Pr. \ref{pr:fur}>>=
   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,15 +147,16 @@ 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 <<Parse options, Pr. \ref{pr:fur}>>=
   flag.Parse()
   //<<Respond to \ty{-v}, Pr. \ref{pr:fur}>>
   //<<Respond to \ty{-d}, Pr. \ref{pr:fur}>>
   //<<Respond to \ty{-t}, Pr. \ref{pr:fur}>>
   //<<Respond to \ty{-f}, Pr. \ref{pr:fur}>>
+  //<<Respond to \ty{-W}, Pr. \ref{pr:fur}>>
 #+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 <<Respond to \ty{-f}, Pr. \ref{pr:fur}>>=
   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 <<Respond to \ty{-W}, Pr. \ref{pr:fur}>>=
+  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 <<Construct Blast commands, Pr. \ref{pr:fur}>>=
   //<<Construct Blast options, Pr. \ref{pr:fur}>>
-  //<<Construct Blast template, Pr. \ref{pr:fur}>>
   //<<Construct megablast command, Pr. \ref{pr:fur}>>
   if !*optM {
 	  //<<Construct blastn command, Pr. \ref{pr:fur}>>
   }
 #+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 <<Construct Blast options, Pr. \ref{pr:fur}>>=
   da := *optD + "/n"
-  th := *optT
-  ev := *optE
+  th := strconv.Itoa(*optT)
+  ev := fmt.Sprintf("%g", *optE)
   ta := "megablast"
   ma := ""
+  ws := ""
   if *optMM {
 	  //<<Look up masking algorithm, Pr. \ref{pr:fur}>>
   }
+  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 <<Construct Blast template, Pr. \ref{pr:fur}>>=
-  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 <<Construct megablast command, Pr. \ref{pr:fur}>>=
-  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 <<Construct blastn command, Pr. \ref{pr:fur}>>=
   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 <<Construct tests, Pr. \ref{pr:fur}>>=
   //<<Tests without repeat masking, Pr. \ref{pr:fur}>>
   //<<Tests with repeat masking, Pr. \ref{pr:fur}>>
   //<<Test with partial intersection, Pr. \ref{pr:fur}>>
+  //<<Tests with word size, Pr. \ref{pr:fur}>>
 #+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 <<Tests with word size, Pr. \ref{pr:fur}>>=
+  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

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 {

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

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