Vidjil: release 2013.04

	* First public release

Makefile

+
+all:
+	make -C src
+
+test: all
+	make -C src/tests 
+	make should
+
+should: all
+	@echo
+	@echo "*** Launching .should_get tests..."
+	src/tests/should-to-tap.sh src/tests/stanford.should_get
+	src/tests/should-to-tap.sh src/tests/clones_simul.should_get
+	src/tests/should-to-tap.sh src/tests/clones_simul_cluster.should_get
+	src/tests/should-to-tap.sh src/tests/segment_S22.should_get
+	src/tests/should-to-tap.sh src/tests/segment_lec.should_get
+	src/tests/should-to-tap.sh src/tests/segment_simul.should_get
+	@echo "*** All .should_get tests passed"
+
+data germline: %:
+	make -C $@
+
+clean:
+	make -C src clean
+
+cleanall: clean
+	make -C data $^
+	make -C germline $^
+
+.PHONY: all test should clean cleanall distrib data germline
+
+RELEASE_TAG="notag"
+RELEASE_H = src/release.h
+RELEASE_SOURCE = $(wildcard src/*.cpp) $(wildcard src/*.h)  $(wildcard src/core/*.cpp)  $(wildcard src/tests/*.cpp) $(wildcard src/core/*.h)  $(wildcard src/tests/*.h)  
+RELEASE_MAKE = ./Makefile  src/Makefile src/tests/Makefile germline/Makefile data/Makefile
+RELEASE_TESTS =  data/get-sequences $(wildcard data/*.fa) $(wildcard data/*.fq) src/tests/should-to-tap.sh $(wildcard src/tests/*.should_get) $(wildcard src/tests/bugs/*.fa)
+RELEASE_FILES = $(RELEASE_SOURCE) $(RELEASE_TESTS) $(RELEASE_MAKE)  germline/get-germline germline/split-from-imgt.py  doc/README doc/LICENSE
+RELEASE_ARCHIVE = vidjil-$(RELEASE_TAG).tgz
+
+CURRENT_DIR = vidjil
+DIST_DIR=$(CURRENT_DIR)-$(RELEASE_TAG)
+RELEASE_FILES_VID = $(RELEASE_FILES)
+
+
+# make distrib RELEASE_TAG=2013.04alpha
+distrib:	
+	$(info ==== Release $(RELEASE_TAG) ====)
+
+	# Tag the release
+	if test "$(RELEASE_TAG)" != "notag"; then \
+		git tag -f release-$(RELEASE_TAG); \
+		echo '#define RELEASE_TAG "$(RELEASE_TAG)"' > $(RELEASE_H); \
+	fi
+
+	mkdir -p release 
+	rm -rf release/$(RELEASE_ARCHIVE) release/$(DIST_DIR)
+	mkdir -p release/$(DIST_DIR)
+	for file in  $(RELEASE_FILES_VID); do\
+		dir=release/$(DIST_DIR)/`dirname "$$file"`;	\
+		mkdir -p $$dir;	\
+		cp "$$file" $$dir;	\
+	done
+	cd release && tar cvzf  $(RELEASE_ARCHIVE) $(DIST_DIR) \
+	&& rm -rf $(DIST_DIR)
+
+	# Untag the source
+	rm -f $(RELEASE_H) ; touch $(RELEASE_H)
+
+	# Check archive
+	cd release && tar xvfz $(RELEASE_ARCHIVE)
+	cd release/$(DIST_DIR) && make
+	cd release/$(DIST_DIR) && make germline
+	cd release/$(DIST_DIR) && make data
+	cd release/$(DIST_DIR) && make test
+
+
+

data/Makefile

+DATA=Stanford_S22.fasta
+
+$(DATA):
+	sh get-sequences
+
+clean:
+	$(RM) -f $(DATA)
\ No newline at end of file

data/clones_simul

+# 100 reads (200bp) from 2 first clones of S22, should be 50%/50%
+
+
+# exact
+wgsim -S42 -N100 -e0 -d0 -1 200 -2 1 -h clones_simul.fa.base clones_simul_exact.fa /dev/null
+
+# some errors
+wgsim -S42 -N100 -e0.01 -d0 -1 200 -2 1 -h clones_simul.fa.base clones_simul.fa /dev/null

data/clones_simul.fa.base

+>clone1
+ggctggagtgggtttcatacattagtagtaatagtggtgccatatactacgcagactctgtgaagggccgattcaccatctccagaaacaatgccaaggactcactgtatctgcaaatgaacagcctgagagccgaggacacggctgtgtattactgtgcgagagcgatcccccggtattactatgatactagtggcccaaacgactactggggccagggaaccctggtcaccgtctcctcag
+>clone2
+aggcttgagtggatgggatggatcaacgctggcaatggtaacacaaaatattcacagaagttccagggcagagtcaccattaccagggacacatccgcgagcacagcctacatggagctgagcagcctgagatctgaagacacggctgtgtattactgtgcgagagtgggcagcagctggtctgatgcttttgattatctggggccaagggacaatggtcaccgtctcctcag

data/get-sequences

+#!/bin/sh
+
+wget http://www.emi.unsw.edu.au/~ihmmune/Stanford_S22.zip
+unzip -d $(dirname $0) -o Stanford_S22.zip
+rm -f Stanford_S22.zip
+

data/segment_S22.fa

+>lcl|FLN1FA002RWEZA.1| 
+ggctggagtgggtttcatacattagtagtaatagtggtgccatatactacgcagactctgtgaagggccgattcaccatc
+tccagaaacaatgccaaggactcactgtatctgcaaatgaacagcctgagagccgaggacacggctgtgtattactgtgc
+gagagcgatcccccggtattactatgatactagtggcccaaacgactactggggccagggaaccctggtcaccgtctcct
+cag
+>lcl|FLN1FA001BLION.1| 
+aggcttgagtggatgggatggatcaacgctggcaatggtaacacaaaatattcacagaagttccagggcagagtcaccat
+taccagggacacatccgcgagcacagcctacatggagctgagcagcctgagatctgaagacacggctgtgtattactgtg
+cgagagtgggcagcagctggtctgatgcttttgattatctggggccaagggacaatggtcaccgtctcctcag
+>lcl|FLN1FA002SGU3N.1| 
+gggcctggagtgggtctcaggtattagttggaatagtggtagaataggctatgcggactctgtgaacggccgattcacca
+tctccagagacaacgccaagaactccctgtatctgcaaatgaacagtctgagagctgaggacacggccttgtattactgt
+gcaaaagataaaattgggtattatgattacgtttgggggagttatcgtttgggtctactggggccagggaaccctggtca
+ccgtctcctcag

data/segment_lec.fa

+>Junc#01
+TATTACTGTGCCACCTGGGCCTTATTATAAGAAACTCTTT
+>Junc#05
+TACTACTGTGCTGCGTGGAGACCCACTGGTTGGTTCAAGA 
+>Rep#01
+CTCATACACCCAGGAGGTGGAGCTGGATATTGATACTACGAAATCTAATTGAAAATGATTCTGGGGTCTATTACTGTGCCACCTGGGCCTTATTATAAGAAACTCTTTGGCAGTGGAAC
+>Rep#05
+ATCCTTACCATCAAGTCCGTAGAGAAAGAAGACATGGCCGTTTACTACTGTGCTGCGTGGAGACCCACTGGTTGGTTCAAGATATTTGCTGAAGGGACTAAGCTCATAGTAACTT
+
+>Junc_rc#05
+TCTTGAACCAACCAGTGGGTCTCCACGCAGCACAGTAGTA
+>Rep_rc#01
+GTTCCACTGCCAAAGAGTTTCTTATAATAAGGCCCAGGTGGCACAGTAATAGACCCCAGAATCATTTTCAATTAGATTTCGTAGTATCAATATCCAGCTCCACCTCCTGGGTGTATGAG
+>Rep_rc#05
+AAGTTACTATGAGCTTAGTCCCTTCAGCAAATATCTTGAACCAACCAGTGGGTCTCCACGCAGCACAGTAGTAAACGGCCATGTCTTCTTTCTCTACGGACTTGATGGTAAGGAT
\ No newline at end of file

data/segment_simul.fa

+>seq1
+aagtccatcagcactgcctacctgcagtggagcagcctgaaggcctcggacaccgccatgtattactgtgcgaCTTCataaccggaaccaGTTAtggttcgactcctggggccaaggaaccctggtcaccgtctcctcag
+
+>seq1_polyT
+TTTTTTTTTTaagtccatcagcactgcctacctgcagtggagcagcctgaaggcctcggacaccgccatgtattactgtgcgaCTTCataaccggaaccaGTTAtggttcgactcctggggccaaggaaccctggtcaccgtctcctcagTTTTTTTTTT

data/test1.fa

+>seq1
+ACAAC
+>seq2
+C
+G
+ACCCCCAA
+>seq3
+
+A
+
+>seq4
+AAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAA
+
+     
+
+>
+A
+
+A
+
+
+T
+
+N
\ No newline at end of file

data/test1.fq

+@seq1
+ACAAC
++
+11111
+@seq2
+CGACCCCCAA
++
+1111111111
+@seq3
+A
++seq3
+1
+@seq4
+AAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAA
++
+@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
+@
+AATN
++
+=5$2

doc/CHANGELOG

+2013-04-18  The Vidjil Team 
+
+	* First public release
+
+

doc/README

+# Vidjil -- V(D)J recombinations analysis <http://bioinfo.lifl.fr/vidjil>
+# Copyright (C) 2011, 2012, 2013 by Bonsai bioinformatics at LIFL (UMR CNRS 8022, Université Lille) and Inria Lille
+# Contact: mathieu.giraud@lifl.fr, mikael.salson@lifl.fr
+
+V(D)J recombinations in lymphocytes are a key for immunologic diversity
+as they have an influence on the production of antibodies and antigen
+receptors. They are also useful markers for pathologies: in many
+cases, such a clonality marker is used for patient follow-up to
+quantify the minimal residual disease in leukemias.  
+
+Vidjil process high-througput sequencing data to extract V(D)J
+junctions and gather them into clones. This analysis is based on a
+seed heuristics and is fast and scalable, as, in the first phase, no
+alignment is done with database germline sequences.  Starting from a
+set of reads, Vidjil detects the junctions in each read. This is based
+on an ultra-fast seed-based heuristic, which has be prooven to be
+reliable output the most abundant clones, based on their junctions
+
+Vidjil can also clusterize similar clones, or leave this to the user
+after a manual review. The method is described in the paper referenced
+below.
+
+Vidjil is open-source, released under GNU GPLv3 license.
+
+
+### Supported platforms
+
+Vidjil has been successfully tested on the following platforms :
+ - CentOS 6.3 amd64
+ - CentOS 6.3 i386
+ - Debian Squeeze 
+ - Fedora 17
+ - FreeBSD 9.1 amd64
+ - NetBSD 6.0.1 amd64
+ - Ubuntu 12.04 amd64
+ - Ubuntu 12.04 i386
+
+
+### Installation
+
+make data
+   # get some IGH rearrangements from a single individual, as described in:
+   # Boyd, S. D., and al. Individual variation in the germline Ig gene
+   # repertoire inferred from variable region gene rearrangements. J
+   # Immunol, 184(12), 6986–92.
+
+make germline
+   # get IMGT germline databases -- you have to agree to IMGT license: 
+   # academic research only, provided that it is referred to IMGT®,
+   # and cited as "IMGT®, the international ImMunoGeneTics information system® 
+   # http://www.imgt.org (founder and director: Marie-Paule Lefranc, Montpellier, France). 
+   # Lefranc, M.-P., IMGT®, the international ImMunoGeneTics database,
+   # Nucl. Acids Res., 29, 207-209 (2001). PMID: 11125093
+
+make                     # compile Vijil
+make test                # run self-tests
+
+./vidjil -h              # display help/usage
+
+### Optional dependencies
+
+clustalw (to compute alignments between junctions from a same clone)
+neato (to display graph of neighbors for the automatic clusterisation)
+
+### Vidjil parameters
+
+Launching vidjil with -h option provides the list of parameters that can be
+used.
+
+### List of junctions
+
+Vidjil allows to specify a list of junctions that must be followed
+(even if those junctions are 'rare', below the -r/-R/-% thresholds).
+The parameter -l is made for providing such a list in a file having
+the following format: junction label (separed by one space)
+
+The first column of the file is the junction to be followed
+while the remaining columns consist of the junction's label.
+In Vidjil output, the labels are output alongside their junctions.
+
+### Manual clustering
+
+The -e option allows to specify a file for manually clustering two junctions
+considered as similar. Such a file may be automatically produced by vidjil
+(out/edges), depending on the option provided. Only the two first columns 
+(separed by one space) are important to vidjil, they only consist of the 
+two junctions that must be clustered.
+
+
+### Examples of use
+
+All the following examples are on a IGH VDJ recombinations : they thus
+require the "-G germline/IGH" and the "-d" options.
+
+
+./vidjil -G germline/IGH -d data/Stanford_S22.fasta
+   # Extract (with an ultra-fast heuristic) all junctions
+   # Results are in out/segmented.vdj.fa, which is a FASTA file 
+   # embedding segmentation information in the headers
+   # ('.vdj' format, see below)
+
+>5--junction--1 
+TTGTAGTGGTGGTAGCTGCTACTCCTTTGACTACTGGGGC
+>5--junction--2 
+TGTAGTGGTGGTAGCTGTTACTCCCACGTCTGGGGCCAAG
+(...)
+
+   Junctions of size 40 (modifiable by -w) have been extracted.
+   These two junctions have 5 occurrences in the set of reads.
+
+./vidjil -c clones -G germline/IGH -x -r 1 -R 1 -d ./data/clones_simul.fa
+   # Extracts the junctions (-r 1, with at least 1 read each), 
+   # then gather them into clones (-R 1, with at least 1 read each:
+   # there are many 1-read clones due to sequencing errors.) 
+   # A more natural option could be -R 5.
+   # No representative selection / clustalw postprocessing (-x)
+   # Results are in out/segmented.fa, out/junctions.fa-* and out/clones*
+   # out/segmented.fa list segmented reads using the .vdj format (see below)
+
+./vidjil -c clones -G germline/IGH -x -r 1 -R 5 -n 5 -d ./data/clones_simul.fa
+   # Junction extraction + clone gathering,
+   # with automatic clusterisation, distance five (-n 5)
+
+./vidjil -c segment -G germline/IGH -d data/segment_S22.fa
+   # Segment the reads onto VDJ germline using a full comparison 
+   # (dynamic programming) with all sequences.
+   # The output is displayed in .vdj format (see below)
+
+
+### .vdj format
+
+Segmentations of V(D)J recombinations are displayed using a dedicated
+format. This format is compatible with FASTA format. A line starting
+with a > is of the following form:
+
+>name + VDJ  startV endV   startD endD   startJ  endJ   Vgene   delV/N1/delD5'   Dgene   delD3'/N2/delJ   Jgene
+
+        name          sequence name
+        +             strand on which the sequence is mapped
+        VDJ           type of segmentation (can be "VJ", "VDJ", 
+    	              or shorter tags such as "V" for incomplete sequences).	
+		      The following line are for "VDJ" recombinations :
+
+        startV endV   start and end position of the V gene in the sequence (start at 0)
+        startD endD                      ... of the D gene ...
+        startJ endJ                      ... of the J gene ...
+
+        Vgene         name of the V gene 
+
+        delV          number of deletions at the end (3') of the V
+        N1            nucleotide sequence inserted between the V and the D
+        delD5'        number of deletions at the start (5') of the D
+
+        Dgene         name of the D gene being rearranged
+
+        delD3'        number of deletions at the end (3') of the D
+        N2            nucleotide sequence inserted between the D and the J
+        delJ          number of deletions at the start (5') of the J
+
+        Jgene         name of the J gene being rearranged
+        
+
+Following such a line, the nucleotide sequence may be given, giving in
+this case a valid FASTA file.
+
+For VJ recombinations the output is similar, the fields that are not
+applicable being removed:
+>name + VJ  startV endV   startJ endJ   Vgene   delV/N1/delJ   Jgene
+

germline/Makefile

+DATA=???[VDJ].fa
+
+data:
+	sh get-germline
+
+clean:
+	rm -f $(DATA)
\ No newline at end of file

germline/get-germline

+#!/bin/sh
+
+cd $(dirname $0)
+wget -O - http://www.imgt.org/download/GENE-DB/GENEDB-ReferenceSequences.fasta-nt-WithoutGaps-F+ORF+inframeP | python split-from-imgt.py 
+

germline/split-from-imgt.py

+
+
+import sys
+
+# Parse lines such as:
+# >M12949|TRGV1*01|Homo sapiens|ORF|...
+
+open_files = {}
+current_file = None
+
+for l in sys.stdin:
+
+    if ">" in l:
+        current_file = None
+        if "Homo sapiens" in l and ("V-REGION" in l or "D-REGION" in l or "J-REGION" in l):
+            system = l.split('|')[1][:4]
+            if system.startswith('IG') or system.startswith('TR'):
+
+                if system in open_files:
+                    current_file = open_files[system]
+                else:
+                    name = '%s.fa' % system
+                    print "  ==>", name
+                    current_file = open(name, 'w')
+                    open_files[system] = current_file
+
+
+    if current_file:
+            current_file.write(l)
+
+
+

src/Makefile

+CC=g++
+OPTIM=-O2 -DNO_SPACED_SEEDS
+CFLAGS=-W -Wall $(OPTIM) $(DEBUG)
+LDFLAGS=
+EXEC=vidjil kmer kmer_count cut count detailed_affect
+SRCCORE=$(wildcard core/*.cpp) $(wildcard core/*.h) 
+MAINCORE=$(wildcard *.cpp)
+OBJCORE=$(SRCCORE:.cpp=.o)
+LIBCORE=core.a
+
+BINDIR=..
+
+v: vidjil
+
+all: $(EXEC)
+
+debug:
+	make clean
+	make DEBUG="-ggdb"
+
+$(EXEC): %: %.o $(LIBCORE) 
+	$(CC) -o $(BINDIR)/$@ $^ $(LDFLAGS) $(CFLAGS)
+
+core: $(LIBCORE)
+
+$(LIBCORE): $(OBJCORE)
+	rm -f $(LIBCORE)
+	ar rc $@ $^
+
+%.o: %.cpp %.h
+	$(CC) -o $@ -c $< $(CFLAGS)
+
+%.o: %.cpp 
+	$(CC) -o $@ -c $< $(CFLAGS)
+
+clean:
+	rm -f core/*.o *.o $(LIBCORE)
+	cd $(BINDIR) ; rm -f $(EXEC)
+	make -C tests $@
+
+forcedep:
+	g++ -M $(SRCCORE) $(MAINCORE) > dep.mk
+
+DEP=$(wildcard dep.mk)
+
+ifeq (${DEP},)
+$(shell	g++ -M $(SRCCORE) $(MAINCORE) > dep.mk)
+endif
+include dep.mk
+
+
+
+
+RELEASE_H = release.h
+RELEASE_SOURCE = $(wildcard src/*.cpp) $(wildcard src/*.h)  $(wildcard src/core/*.cpp) $(wildcard src/core/*.h)  
+RELEASE_FILES = $(RELEASE_SOURCE) src/Makefile germline/get-germline doc/README doc/LICENSE
+RELEASE_ARCHIVE = vidjil-$(RELEASE_TAG).tgz
+
+CURRENT_DIR = .
+RELEASE_FILES_VID = $(addprefix $(CURRENT_DIR)/, $(RELEASE_FILES))
+
+
+# make distrib RELEASE_TAG=2013.04alpha
+distrib:	
+	$(info ==== Release $(RELEASE_TAG) ====)
+
+	# Tag the release
+	git tag release-$(RELEASE_H)
+	echo '#define RELEASE_TAG "$(RELEASE_TAG)"' > $(RELEASE_H)
+
+	mkdir -p release 
+	rm -f release/$(RELEASE_ARCHIVE)
+	cd .. ; tar cvfz  $(CURRENT_DIR)/release/$(RELEASE_ARCHIVE) $(RELEASE_FILES_VID)
+
+	# Untag the source
+	rm $(RELEASE_H) ; touch $(RELEASE_H)
+
+	# Check archive
+	cd release ; tar xvfz $(RELEASE_ARCHIVE)
+	cd release/$(CURRENT_DIR) ; make
+
+

src/core/affectanalyser.h

+#ifndef AFFECT_ANALYSER_H
+#define AFFECT_ANALYSER_H
+
+#include "kmerstore.h"
+#include <set>
+#include <vector>
+#include <cassert>
+
+// Define two constant affectations: ambiguous and unknown.
+
+/* Declaration of types */
+
+typedef enum affect_options_e {
+  AO_NONE, AO_NO_CONSECUTIVE, AO_NO_MULTIPLICITY
+} affect_options_t;
+
+/**
+ * Class that records for every k-mer of a given sequence
+ * in which sequences this k-mer was also seen.
+ * It can either record one affectation per kmer (the only sequence where it
+ * occurs or ambiguous case if there are several possibilities, or 
+ * unknown otherwise), or all the possible affectations (ie. there is no ambiguous
+ * case, all the possibilities for one k-mer are stored).
+ *
+ * Input: Index that constitutes the k-mer sequence repertoire
+ * Input: Sequence whose k-mers must be affected
+ */
+template<class T>
+class AffectAnalyser {
+ public:
+
+  /* Queries */
+
+  /**
+   * @return the total number of affectations
+   */
+  virtual int count() const = 0;
+
+  /**
+   * @param affect: An affectation
+   * @return the number of times this affectation has been given in the read.
+   */
+  virtual int count(const T &affect)  const = 0;
+
+  /**
+   * @param i: the position to consider
+   * @pre i >= 0 && i < count()
+   * @return the affectation of the k-mer at position i.
+   */
+  virtual const T&getAffectation(int i)  const = 0;
+
+  /**
+   * @param options: options can either be AO_NONE or AO_NO_CONSECUTIVE 
+   * @return all the affectations contained in the read from left to right.
+   *         if AO_NO_CONSECUTIVE is given: two consecutive elements in the vector
+   *                                     will be different (we remove consecutive
+   *                                     duplicates)
+   */
+  virtual vector<T> getAllAffectations(affect_options_t options) const = 0;
+
+  /**
+   * @return the distinct affectations
+   */
+  virtual set<T> getDistinctAffectations() const = 0;
+
+  /**
+   * @return the sequence we are analysing
+   */
+  virtual const string &getSequence() const = 0;
+
+  /**
+   * @param affect: an affectation
+   * @return the first occurrence of this affectation in the read
+   *         or string::npos if the affectation was not found
+   * @post getAffectation(first(affect)) == affect 
+   * ==>  getAffectation(1...first(affect)-1) != affect
+   */
+  virtual int first(const T &affect) const  = 0;
+
+  /**
+   * @param affect: an affectation
+   * @return the last occurrence of this affectation in the read
+   *         or string::npos if the affectation was not found
+   * @post getAffectation(last(affect)) == affect 
+   * ==> getAffectation(last(affect)+1 ... count() -1) != affect
+   */
+  virtual int last(const T &affect) const  = 0;
+
+  /**
+   * @return a string representation of the object
+   */
+  virtual string toString() const  = 0;
+};
+
+template <class T>
+class KmerAffectAnalyser: public AffectAnalyser<T> {
+ private:
+  IKmerStore<T> &kms;
+  const string &seq;
+  vector<T> affectations;
+ public:
+  /**
+   * @param kms: the index storing the affectation for the k-mers
+   *             (parameter is not copied)
+   * @param seq: the sequence to analyse (parameter is not copied)
+   */
+  KmerAffectAnalyser(IKmerStore<T> &kms, const string &seq);
+  ~KmerAffectAnalyser();
+
+  int count() const;
+
+  int count(const T &affect) const;
+
+  const T&getAffectation(int i) const;
+
+  vector<T> getAllAffectations(affect_options_t options) const;
+
+  set<T> getDistinctAffectations() const;
+
+  const string &getSequence() const;
+
+  int first(const T &affect) const;
+
+  int last(const T &affect) const ;
+
+  string toString() const;
+};
+
+
+template <class T>
+KmerAffectAnalyser<T>::KmerAffectAnalyser(IKmerStore<T> &kms, 
+                                       const string &seq)
+  :kms(kms), seq(seq) {
+  assert(seq.length() >=  (size_t)kms.getS());
+  affectations = kms.getResults(seq, true);