Added option --log to goalign mutate gaps

Author: fredericlemoine

align/align.go

@@ -19,7 +19,8 @@ import (
 // Alignment represents a set of aligned sequences (multiple Sequence Alignment)
 type Alignment interface {
 	SeqBag
-	AddGaps(rate, lenprop float64, rand *mathrand.Rand)
+	// Add random gaps to random sequences, and output the indices of the affected sequences
+	AddGaps(rate, lenprop float64, rand *mathrand.Rand) []int
 	AddAmbiguities(rate, lenprop float64, rand *mathrand.Rand)
 	Append(Alignment) error // Appends alignment sequences to this alignment
 	AvgAllelesPerSite() float64
@@ -996,7 +997,8 @@ func (a *align) Recombine(prop float64, lenprop float64, swap bool, rand *mathra
 // Add prop*100% gaps to lenprop*100% of the sequences
 // if prop < 0 || lenprop<0 : does nothing
 // if prop > 1 || lenprop>1 : does nothing
-func (a *align) AddGaps(lenprop float64, prop float64, rand *mathrand.Rand) {
+// Returns the indices of the affected sequences
+func (a *align) AddGaps(lenprop float64, prop float64, rand *mathrand.Rand) (affected []int) {
 	if prop < 0 || prop > 1 {
 		return
 	}
@@ -1012,10 +1014,12 @@ func (a *align) AddGaps(lenprop float64, prop float64, rand *mathrand.Rand) {
 	for i := 0; i < nb; i++ {
 		permsites := rand.Perm(a.Length())
 		seq := a.seqs[permseqs[i]]
+		affected = append(affected, permseqs[i])
 		for j := 0; j < nbgaps; j++ {
 			seq.sequence[permsites[j]] = GAP
 		}
 	}
+	return
 }
 
 // Add prop*100% ambiguities to lenprop*100% of the sequences

cmd/addgaps.go

@@ -1,13 +1,16 @@
 package cmd
 
 import (
+	"fmt"
+
 	"github.com/evolbioinfo/goalign/align"
 	"github.com/evolbioinfo/goalign/io"
 	"github.com/evolbioinfo/goalign/io/utils"
 	"github.com/spf13/cobra"
 )
 
 var gapnbseqs float64
+var addgapslogfile string
 
 // rogueCmd represents the rogue command
 var addgapsCmd = &cobra.Command{
@@ -21,6 +24,7 @@ goalign mutate gaps -i align.fa -n 0.5 -r 0.5
 	RunE: func(cmd *cobra.Command, args []string) (err error) {
 		var aligns *align.AlignChannel
 		var f utils.StringWriterCloser
+		var addgapslog utils.StringWriterCloser
 
 		if aligns, err = readalign(infile); err != nil {
 			io.LogError(err)
@@ -32,9 +36,25 @@ goalign mutate gaps -i align.fa -n 0.5 -r 0.5
 		}
 		defer utils.CloseWriteFile(f, mutateOutput)
 
+		if addgapslog, err = utils.OpenWriteFile(addgapslogfile); err != nil {
+			io.LogError(err)
+			return
+		}
+		defer utils.CloseWriteFile(addgapslog, addgapslogfile)
+
 		for al := range aligns.Achan {
-			al.AddGaps(mutateRate, gapnbseqs, globalRand)
+			affected := al.AddGaps(mutateRate, gapnbseqs, globalRand)
 			writeAlign(al, f)
+			for _, s := range affected {
+				if n, ok := al.GetSequenceNameById(s); !ok {
+					err = fmt.Errorf("affected sequence not found in the alignment: %d", s)
+					io.LogError(err)
+					return
+				} else {
+					fmt.Fprintf(addgapslog, "%d\t%s\n", s, n)
+				}
+			}
+
 		}
 
 		if aligns.Err != nil {
@@ -48,4 +68,5 @@ goalign mutate gaps -i align.fa -n 0.5 -r 0.5
 func init() {
 	mutateCmd.AddCommand(addgapsCmd)
 	addgapsCmd.PersistentFlags().Float64VarP(&gapnbseqs, "prop-seq", "n", 0.5, "Proportion of the sequences in which to add gaps")
+	addgapsCmd.PersistentFlags().StringVar(&addgapslogfile, "log", "none", "Log file with the names of affected sequences")
 }