2019/07/21-17:10:27 (Linux thalamus unknown)

Author: Philipp Benner

kmers.go kmer.gokmers.go renamed to kmer.go

@@ -21,17 +21,18 @@ package gonetics
 //import "fmt"
 import "sort"
 
-/* -------------------------------------------------------------------------- */
+/* K-mer equivalence class
+ * -------------------------------------------------------------------------- */
 
-type Kmer struct {
+type KmerClass struct {
   // K-mer ID
   K    int
   I    int
   // K-mer string representation
   Name string
 }
 
-func (obj Kmer) Equals(b Kmer) bool {
+func (obj KmerClass) Equals(b KmerClass) bool {
   if obj.K != b.K {
     return false
   }
@@ -43,7 +44,7 @@ func (obj Kmer) Equals(b Kmer) bool {
 
 /* -------------------------------------------------------------------------- */
 
-type KmerList []Kmer
+type KmerList []KmerClass
 
 func (obj KmerList) Clone() KmerList {
   r := make(KmerList, len(obj))
@@ -98,7 +99,7 @@ func (obj KmerList) Union(b ...KmerList) KmerList {
 
 /* -------------------------------------------------------------------------- */
 
-type KmerSet map[Kmer]struct{}
+type KmerSet map[KmerClass]struct{}
 
 func (obj KmerSet) AsList() KmerList {
   r := make(KmerList, len(obj))
@@ -116,7 +117,7 @@ type KmerCounts struct {
   // this is a sorted list of k-mers and might contain more entries than
   // the counts map
   Kmers  KmerList
-  Counts map[Kmer]int
+  Counts map[KmerClass]int
 }
 
 func (obj KmerCounts) Len() int {
@@ -135,83 +136,18 @@ func (obj KmerCounts) At(i int) int {
   }
 }
 
-func (obj KmerCounts) GetCount(kmer Kmer) int {
+func (obj KmerCounts) GetCount(kmer KmerClass) int {
   if c, ok := obj.Counts[kmer]; ok {
     return c
   } else {
     return 0
   }
 }
 
-func (obj KmerCounts) GetKmer(i int) Kmer {
+func (obj KmerCounts) GetKmer(i int) KmerClass {
   return obj.Kmers[i]
 }
 
 func (obj KmerCounts) Iterate() KmerCountsIterator {
   return KmerCountsIterator{obj, 0}
 }
-
-/* -------------------------------------------------------------------------- */
-
-type KmerCountsList struct {
-  Kmers    KmerList
-  Counts []map[Kmer]int
-}
-
-func NewKmerCountsList(counts ...KmerCounts) KmerCountsList {
-  r := KmerCountsList{}
-  return r.Append(counts...)
-}
-
-func (obj KmerCountsList) Append(args ...KmerCounts) KmerCountsList {
-  if len(args) == 0 {
-    return obj
-  }
-  idLists := make([]KmerList, len(args))
-  counts  := obj.Counts
-  for i, c := range args {
-    idLists[i] = c.Kmers
-    counts = append(counts, c.Counts)
-  }
-  ids := obj.Kmers.Union(idLists...)
-  return KmerCountsList{Kmers: ids, Counts: counts}
-}
-
-func (obj KmerCountsList) Len() int {
-  return len(obj.Counts)
-}
-
-func (obj KmerCountsList) At(i int) KmerCounts {
-  return KmerCounts{Kmers: obj.Kmers, Counts: obj.Counts[i]}
-}
-
-func (obj *KmerCountsList) Slice(i, j int) KmerCountsList {
-  return KmerCountsList{Kmers: obj.Kmers, Counts: obj.Counts[i:j]}
-}
-
-/* -------------------------------------------------------------------------- */
-
-type KmerCountsIterator struct {
-  KmerCounts
-  i int
-}
-
-func (obj KmerCountsIterator) Ok() bool {
-  return obj.i < obj.Len()
-}
-
-func (obj KmerCountsIterator) GetKmer() Kmer {
-  return obj.Kmers[obj.i]
-}
-
-func (obj KmerCountsIterator) GetName() string {
-  return obj.Kmers[obj.i].Name
-}
-
-func (obj KmerCountsIterator) GetCount() int {
-  return obj.At(obj.i)
-}
-
-func (obj *KmerCountsIterator) Next() {
-  obj.i++
-}

kmer_catalogue.go

@@ -0,0 +1,191 @@
+/* Copyright (C) 2019 Philipp Benner
+ *
+ * This program is free software: you can redistribute it and/or modify
+ * it under the terms of the GNU General Public License as published by
+ * the Free Software Foundation, either version 3 of the License, or
+ * (at your option) any later version.
+ *
+ * This program is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+ * GNU General Public License for more details.
+ *
+ * You should have received a copy of the GNU General Public License
+ * along with this program.  If not, see <http://www.gnu.org/licenses/>.
+ */
+
+package gonetics
+
+/* -------------------------------------------------------------------------- */
+
+//import "fmt"
+import "sort"
+import "strings"
+
+/* -------------------------------------------------------------------------- */
+
+type KmerCatalogue struct {
+  KmerEquivalenceRelation
+  names []map[int]string   // k-mer names (equivalent k-mers separated by pipe)
+  idmap []map[string]int   // unique k-mer IDs
+}
+
+/* -------------------------------------------------------------------------- */
+
+func NewKmerCatalogue(n, m int, comp, rev, rc bool, maxAmbiguous []int, al ComplementableAlphabet) (*KmerCatalogue, error) {
+  r := KmerCatalogue{}
+  if f, err := NewKmerEquivalenceRelation(n, m, comp, rev, rc, maxAmbiguous, al); err != nil {
+    return nil, err
+  } else {
+    r.KmerEquivalenceRelation = f
+  }
+  idmap := make([]map[string]int, m-n+1)
+  names := make([]map[int]string, m-n+1)
+  for k := n; k <= m; k++ {
+    idmap[k-n] = make(map[string]int)
+    names[k-n] = make(map[int]string)
+  }
+  r.names = names
+  r.idmap = idmap
+  return &r, nil
+}
+
+/* -------------------------------------------------------------------------- */
+
+func (obj *KmerCatalogue) GetId(kmer string) int {
+  k := len(kmer)
+  if k < obj.n || k > obj.m {
+    panic("k-mer has invalid length")
+  }
+  if i, ok := obj.idmap[k-obj.n][kmer]; ok {
+    return i
+  } else {
+    r := obj.EquivalenceClass(kmer)
+    for _, kmer := range strings.Split(r.Name, "|") {
+      obj.idmap[k-obj.n][kmer] = r.I
+    }
+    obj.names[k-obj.n][r.I] = r.Name
+    return r.I
+  }
+}
+
+func (obj *KmerCatalogue) GetName(kmer string) string {
+  k := len(kmer)
+  if k < obj.n || k > obj.m {
+    panic("k-mer has invalid length")
+  }
+  if i, ok := obj.idmap[k-obj.n][kmer]; ok {
+    return obj.names[k-obj.n][i]
+  } else {
+    r := obj.EquivalenceClass(kmer)
+    for _, kmer := range strings.Split(r.Name, "|") {
+      obj.idmap[k-obj.n][kmer] = r.I
+    }
+    obj.names[k-obj.n][r.I] = r.Name
+    return r.Name
+  }
+}
+
+func (obj *KmerCatalogue) IdToName(k, id int) string {
+  if name, ok := obj.names[k-obj.n][id]; ok {
+    return name
+  }
+  panic("k-mer name not found")
+}
+
+func (obj *KmerCatalogue) GetNames(kmer string) []string {
+  return strings.Split(obj.GetName(kmer), "|")
+}
+
+func (obj *KmerCatalogue) ObservedKmers() int {
+  r := 0
+  for i := 0; i < len(obj.names); i++ {
+    r += len(obj.names[i])
+  }
+  return r
+}
+
+/* -------------------------------------------------------------------------- */
+
+func (obj *KmerCatalogue) scanSubKmers_(kmer []byte, k int) []int {
+  idMap := make(map[int]struct{})
+  // loop over sequence
+  for i := 0; i < len(kmer); i++ {
+    if i+k-1 >= len(kmer) {
+      break
+    }
+    it := NewKmerInstantiationIterator(obj.al, string(kmer[i:i+k]), true)
+    for ; it.Ok(); it.Next() {
+      if id, ok := obj.idmap[k-obj.n][it.Get()]; ok {
+        idMap[id] = struct{}{}
+      }
+    }
+  }
+  ids := []int{}
+  for id, _ := range idMap {
+    ids = append(ids, id)
+  }
+  sort.Ints(ids)
+  return ids
+}
+
+func (obj *KmerCatalogue) scanSubKmers(kmer []byte) []string {
+  names := []string{}
+  for k := obj.n; k <= obj.m; k++ {
+    ids := obj.scanSubKmers_(kmer, k)
+    for _, id := range ids {
+      names = append(names , obj.IdToName(k, id))
+    }
+  }
+  return names
+}
+
+/* -------------------------------------------------------------------------- */
+
+func (obj *KmerCatalogue) relatedKmers(s []byte, m, k int) []int {
+  idMap := make(map[int]struct{})
+  // loop over positions where the k-mer can be fixed
+  for j := 0; j <= k - len(s); j++ {
+    for it := NewKmerCylinderIterator(k, obj.ma[k-obj.n] - m, obj.al, j, string(s)); it.Ok(); it.Next() {
+      if id, ok := obj.idmap[k-obj.n][it.Get()]; ok {
+        idMap[id] = struct{}{}
+      }
+    }
+  }
+  ids := []int{}
+  for id, _ := range idMap {
+    ids = append(ids, id)
+  }
+  sort.Ints(ids)
+  return ids
+}
+
+func (obj *KmerCatalogue) RelatedKmers(kmer string) []string {
+  s := []byte(obj.GetNames(kmer)[0])
+  m := obj.countAmbiguous(kmer)
+  // scan k-mer for sub-k-mers
+  names := obj.scanSubKmers(s)
+  // loop over k-mer sizes
+  for k := len(s)+1; k <= obj.m; k++ {
+    ids := obj.relatedKmers(s, m, k)
+    for _, id := range ids {
+      names = append(names , obj.IdToName(k, id))
+    }
+  }
+  return names
+}
+
+func (obj *KmerCatalogue) countAmbiguous(kmer string) int {
+  m := 0
+  s := []byte(kmer)
+  for i := 0; i < len(s); i++ {
+    if ok, err := obj.al.IsAmbiguous(s[i]); err != nil {
+      panic("internal error")
+    } else {
+      if ok {
+        m++
+      }
+    }
+  }
+  return m
+}

kmers_catalogue_iterator.go kmer_catalogue_iterator.gokmers_catalogue_iterator.go renamed to kmer_catalogue_iterator.go

@@ -22,15 +22,15 @@ package gonetics
 
 /* -------------------------------------------------------------------------- */
 
-type KmersCatalogueIterator struct {
+type KmerCatalogueIterator struct {
   names []string
   ids   []int
   i       int
 }
 
 /* -------------------------------------------------------------------------- */
 
-func NewKmersCatalogueIterator(kmersSet KmersCatalogue) KmersCatalogueIterator {
+func NewKmerCatalogueIterator(kmersSet KmerCatalogue) KmerCatalogueIterator {
   names := []string{}
   ids   := []int{}
   for k := 0; k < len(kmersSet.names); k++ {
@@ -44,23 +44,23 @@ func NewKmersCatalogueIterator(kmersSet KmersCatalogue) KmersCatalogueIterator {
     names = append(names, n...)
     ids   = append(ids,   i...)
   }
-  return KmersCatalogueIterator{names: names, ids: ids, i: 0}
+  return KmerCatalogueIterator{names: names, ids: ids, i: 0}
 }
 
 /* -------------------------------------------------------------------------- */
 
-func (obj KmersCatalogueIterator) Ok() bool {
+func (obj KmerCatalogueIterator) Ok() bool {
   return obj.i < len(obj.names)
 }
 
-func (obj KmersCatalogueIterator) Get() string {
+func (obj KmerCatalogueIterator) Get() string {
   return obj.names[obj.i]
 }
 
-func (obj KmersCatalogueIterator) GetId() int {
+func (obj KmerCatalogueIterator) GetId() int {
   return obj.ids[obj.i]
 }
 
-func (obj *KmersCatalogueIterator) Next() {
+func (obj *KmerCatalogueIterator) Next() {
   obj.i++
 }