Added substitution commands and moved gap commands

Author: fredericlemoine

align/align.go

@@ -43,6 +43,7 @@ type Alignment interface {
 	IterateAll(it func(name string, sequence []rune, comment string))
 	NbSequences() int
 	Length() int
+	Mutate(rate float64) // Adds uniform substitutions in the alignment (~sequencing errors)
 	ShuffleSequences()
 	ShuffleSites(rate float64, roguerate float64) []string
 	SimulateRogue(prop float64, proplen float64) ([]string, []string)
@@ -400,6 +401,40 @@ func (a *align) AddGaps(lenprop float64, prop float64) {
 	}
 }
 
+// Add substitutions uniformly to the alignment
+// if rate < 0 : does nothing
+// if rate > 1 : rate=1
+func (a *align) Mutate(rate float64) {
+	if rate <= 0 {
+		return
+	}
+	if rate > 1 {
+		rate = 1
+	}
+	r := 0.0
+	newchar := 0
+	leng := a.Length()
+	nb := a.NbSequences()
+	// We take a random position in the sequences between min and max
+	for i := 0; i < nb; i++ {
+		seq := a.seqs[i]
+		for j := 0; j < leng; j++ {
+			r = rand.Float64()
+			// We mutate only if rand is <= rate
+			// taking a random nucleotide or amino acid uniformly
+			if r <= rate {
+				if a.Alphabet() == AMINOACIDS {
+					newchar = rand.Intn(len(stdaminoacid))
+					seq.sequence[j] = stdaminoacid[newchar]
+				} else {
+					newchar = rand.Intn(len(stdnucleotides))
+					seq.sequence[j] = stdnucleotides[newchar]
+				}
+			}
+		}
+	}
+}
+
 // Simulate rogue taxa in the alignment:
 // take the proportion prop of sequences as rogue taxa => R
 // For each t in R

cmd/addgaps.go

@@ -1,38 +1,34 @@
 package cmd
 
 import (
+	"math/rand"
+
 	"github.com/spf13/cobra"
 )
 
-var gaprate float64
 var gapnbseqs float64
 
 // rogueCmd represents the rogue command
 var addgapsCmd = &cobra.Command{
 	Use:   "gaps",
-	Short: "Adds random (uniformly) gaps",
-	Long: `Adds random (uniformly) gaps on an input alignment.
+	Short: "Adds gaps uniformly in an input alignment",
+	Long: `Adds gaps uniformly in an input alignment.
 
 Example:
-goalign shuffle gaps -i align.fa -n 0.5 -r 0.5
-
+goalign mutate gaps -i align.fa -n 0.5 -r 0.5
 `,
 	Run: func(cmd *cobra.Command, args []string) {
-
-		f := openWriteFile(shuffleOutput)
-		nameFile := openWriteFile(rogueNameFile)
+		rand.Seed(mutateSeed)
+		f := openWriteFile(mutateOutput)
 		for al := range rootaligns {
-			al.AddGaps(gaprate, gapnbseqs)
+			al.AddGaps(mutateRate, gapnbseqs)
 			writeAlign(al, f)
 		}
-		nameFile.Close()
 		f.Close()
 	},
 }
 
 func init() {
-	shuffleCmd.AddCommand(addgapsCmd)
-
-	addgapsCmd.PersistentFlags().Float64VarP(&gapnbseqs, "prop-seq", "n", 0.5, "Proportion of the sequences to add gaps")
-	addgapsCmd.PersistentFlags().Float64VarP(&gaprate, "gap-rate", "r", 0.5, "Proportion of gaps to add per sequences")
+	mutateCmd.AddCommand(addgapsCmd)
+	addgapsCmd.PersistentFlags().Float64VarP(&gapnbseqs, "prop-seq", "n", 0.5, "Proportion of the sequences in which to add gaps")
 }

cmd/gaps.go

@@ -0,0 +1,31 @@
+package cmd
+
+import (
+	"math/rand"
+
+	"github.com/spf13/cobra"
+)
+
+// mutateCmd represents the mutate command
+var gapsCmd = &cobra.Command{
+	Use:   "snvs",
+	Short: "Adds substitutions uniformly in an input alignment",
+	Long: `Adds substitutions uniformly in an input alignment.
+
+      if rate <= 0 : does nothing
+      if rate > 1  : then rate = 1
+`,
+	Run: func(cmd *cobra.Command, args []string) {
+		rand.Seed(mutateSeed)
+		f := openWriteFile(mutateOutput)
+		for al := range rootaligns {
+			al.Mutate(mutateRate)
+			writeAlign(al, f)
+		}
+		f.Close()
+	},
+}
+
+func init() {
+	mutateCmd.AddCommand(gapsCmd)
+}

cmd/mutate.go

@@ -0,0 +1,26 @@
+package cmd
+
+import (
+	"time"
+
+	"github.com/spf13/cobra"
+)
+
+var mutateOutput string
+var mutateRate float64
+var mutateSeed int64
+
+// mutateCmd represents the mutate command
+var mutateCmd = &cobra.Command{
+	Use:   "mutate",
+	Short: "Adds substitutions (~sequencing errors), or gaps, uniformly in an input alignment",
+	Long: `Adds substitutions (~sequencing error), or gaps, uniformly in an input alignment.
+`,
+}
+
+func init() {
+	RootCmd.AddCommand(mutateCmd)
+	mutateCmd.PersistentFlags().Float64VarP(&mutateRate, "rate", "r", 0.1, "Mutation rate per nucleotide/amino acid")
+	mutateCmd.PersistentFlags().StringVarP(&mutateOutput, "output", "o", "stdout", "Mutated alignment output file")
+	mutateCmd.PersistentFlags().Int64VarP(&mutateSeed, "seed", "s", time.Now().UTC().UnixNano(), "Initial Random Seed")
+}

docs/commands/mutate.md

@@ -0,0 +1,120 @@
+# Goalign: toolkit and api for alignment manipulation
+
+## Commands
+
+### mutate
+This command adds different type of noises in an input alignment, with these sub-commands:
+* `goalign mutate gaps` : Adds a given proportion of gaps to a given proportion of the sequences randomly in the input alignment (uniformly).
+* `goalign mutate snvs`: Substitute nucleotides/aminoacids by random (uniform) nucleotides/aminoacids with a given rate.
+
+#### Usage
+* General command:
+```
+Usage:
+  goalign mutate [command]
+
+Available Commands:
+  gaps        Adds gaps uniformly in an input alignment
+  snvs        Adds substitutions uniformly in an input alignment
+
+Flags:
+  -h, --help            help for mutate
+  -o, --output string   Mutated alignment output file (default "stdout")
+  -r, --rate float      Mutation rate per nucleotide/amino acid (default 0.1)
+  -s, --seed int        Initial Random Seed (default 1497950022161439493)
+
+Global Flags:
+  -i, --align string   Alignment input file (default "stdin")
+  -p, --phylip         Alignment is in phylip? False=Fasta
+```
+
+* gaps command:
+```
+Usage:
+  goalign shuffle gaps [flags]
+
+Flags:
+  -r, --gap-rate float   Proportion of gaps to add per sequences (default 0.5)
+  -n, --prop-seq float   Proportion of the sequences to add gaps in (default 0.5)
+
+Global Flags:
+  -i, --align string    Alignment input file (default "stdin")
+  -o, --output string   Shuffled alignment output file (default "stdout")
+  -p, --phylip          Alignment is in phylip? False=Fasta
+  -s, --seed int        Initial Random Seed (default: number of nanoseconds elapsed since January 1, 1970 UTC)
+```
+
+* snvs command:
+```
+Usage:
+  goalign shuffle recomb [flags]
+
+Flags:
+  -l, --prop-length float   Proportion of length of sequences to recombine (default 0.5)
+  -n, --prop-seq float      Proportion of the  sequences to recombine (default 0.5)
+
+Global Flags:
+  -i, --align string    Alignment input file (default "stdin")
+  -o, --output string   Shuffled alignment output file (default "stdout")
+  -p, --phylip          Alignment is in phylip? False=Fasta
+  -s, --seed int        Initial Random Seed (default 1495021613485404707)
+```
+
+#### Examples
+* Generating a random (uniform) alignment and adding 20% gaps to 50% of the sequences:
+```
+goalign random -l 20 -s10| goalign mutate gaps -n 0.5 -r 0.2 -s10
+```
+
+Should give:
+```
+>Seq0000
+GATTAATTTGCCGTAGGCCA
+>Seq0001
+G-ATCTGAAGA-CG-A-ACT
+>Seq0002
+TTAAGTTTT-AC--CTAA-G
+>Seq0003
+GAGAGGACTAGTTCATACTT
+>Seq0004
+TT-AAACA-TTTTA-A-CGA
+>Seq0005
+TGTCGGACCTAAGTATTGAG
+>Seq0006
+TAC-A-G-TGTATT-CAGCG
+>Seq0007
+GTGGAGAGGTCTATTTTTCC
+>Seq0008
+GGTTGAAG-ACT-TA-AGC-
+>Seq0009
+GTAAAGGGTATGGCCATGTG
+```
+
+* Generating a random (uniform) alignment and adding 10% substitutions :
+```
+goalign random -l 20 -s10| goalign mutate snvs -r 0.1 -s 10
+```
+
+Should give:
+```
+>Seq0000
+GATTAATTTCCCGTAGGCCA
+>Seq0001
+GAATCTGAATATCGAACTAT
+>Seq0002
+TTAAGTTTTCACTTCTAATG
+>Seq0003
+GAGAGGACTAGTTCATAATT
+>Seq0004
+TTTTAACACTTTTACATCGA
+>Seq0005
+TGTCGGACCTAAGTTTTGTG
+>Seq0006
+TGCAACGATGTACTCCAGCG
+>Seq0007
+GTGGAGAGGTCTATTTTTGC
+>Seq0008
+GGTTAAAGGACTCTATAGCT
+>Seq0009
+GAAAAGGGTATGGCCATGTG
+```

docs/commands/shuffle.md

@@ -4,7 +4,6 @@
 
 ### shuffle
 This command adds different type of noises in an input alignment, with these sub-commands:
-* `goalign shuffle gaps` : Adds a given proportion of gaps to a given proportion of the sequences randomly in the input alignment (uniformly).
 * `goalign shuffle recomb`:Recombine a given proportion of the length of a given proportion of the sequences with other sequences (copy/paste). 
     Example:
 	```
@@ -40,7 +39,6 @@ Usage:
   goalign shuffle [command]
 
 Available Commands:
-  gaps        Adds random (uniformly) gaps
   recomb      Recombine sequences in the input alignment
   rogue       Simulate rogue taxa
   seqs        Shuffles sequence order in alignment
@@ -57,22 +55,6 @@ Global Flags:
   -p, --phylip         Alignment is in phylip? False=Fasta
 ```
 
-* gaps command:
-```
-Usage:
-  goalign shuffle gaps [flags]
-
-Flags:
-  -r, --gap-rate float   Proportion of gaps to add per sequences (default 0.5)
-  -n, --prop-seq float   Proportion of the sequences to add gaps in (default 0.5)
-
-Global Flags:
-  -i, --align string    Alignment input file (default "stdin")
-  -o, --output string   Shuffled alignment output file (default "stdout")
-  -p, --phylip          Alignment is in phylip? False=Fasta
-  -s, --seed int        Initial Random Seed (default: number of nanoseconds elapsed since January 1, 1970 UTC)
-```
-
 * recomb command:
 ```
 Usage:
@@ -152,35 +134,6 @@ Global Flags:
 ```
 
 #### Examples
-* Generating a random (uniform) alignment and adding gaps:
-```
-goalign random -l 20 -s10| goalign shuffle gaps -n 0.5 -r 0.2 -s10
-```
-
-Should give:
-```
->Seq0000
-GATTAATTTGCCGTAGGCCA
->Seq0001
-G-ATCTGAAGA-CG-A-ACT
->Seq0002
-TTAAGTTTT-AC--CTAA-G
->Seq0003
-GAGAGGACTAGTTCATACTT
->Seq0004
-TT-AAACA-TTTTA-A-CGA
->Seq0005
-TGTCGGACCTAAGTATTGAG
->Seq0006
-TAC-A-G-TGTATT-CAGCG
->Seq0007
-GTGGAGAGGTCTATTTTTCC
->Seq0008
-GGTTGAAG-ACT-TA-AGC-
->Seq0009
-GTAAAGGGTATGGCCATGTG
-```
-
 * Recombining 0.25% of the sequences on 0.5% of their length
 ```
 # Writing alignment

docs/index.md

@@ -46,6 +46,9 @@ Command                                                     | Subcommand |
 --                                                          | pssm       | Computes and prints a Position specific scoring matrix
 [concat](commands/concat.md) ([api](api/concat.md))         |            | Concatenates a set of alignment
 [divide](commands/divide.md) ([api](api/divide.md))         |            | Divide an input alignment in several output files
+[mutate](commands/mutate.md) ([api](api/mutate.md))         |            | Adds substitutions (~sequencing errors), or gaps, uniformly in an input alignment
+--                                                          | gaps       | Adds gaps uniformly in an input alignment
+--                                                          | snvs       | Adds substitutions uniformly in an input alignment
 [random](commands/random.md) ([api](api/random.md))         |            | Generate random sequences
 [reformat](commands/reformat.md) ([api](api/reformat.md))   |            | Reformats input alignment into phylip of fasta format
 --                                                          | fasta      | Reformats an input alignment into Fasta
@@ -57,7 +60,6 @@ Command                                                     | Subcommand |
 --                                                          | seqs       | Samples a subset of sequences from the input alignment
 --                                                          | sites      | Take a random subalignment
 [shuffle](commands/shuffle.md) ([api](api/shuffle.md))      |            | A set of commands to shuffle an alignment
---                                                          | gaps       | Adds random (uniformly) gaps on an input alignment
 --                                                          | recomb     | Recombines sequences in the input alignment (copy/paste)
 --                                                          | rogue      | Simulates rogue taxa
 --                                                          | seqs       | Shuffles sequence order in alignment