fixes

RevantOrg · Sep 3, 2023 · 835bc9a · 835bc9a
1 parent 53c6981
commit 835bc9a
Show file tree

Hide file tree

Showing 10 changed files with 187 additions and 52 deletions.
diff --git a/scripts/6-repeatAlphabet/1-buildAlphabet.sh b/scripts/6-repeatAlphabet/1-buildAlphabet.sh
@@ -32,7 +32,7 @@ CONCATENATE_BLOCKS=${12}  # 0=do not try to merge adjacent blocks from the same
 KEEP_PERIODIC="1"  # 1=do not remove rare characters if they are periodic. Usually good.
 # ----------------------------------------------------------------------------------------
 
-set -o pipefail; set -e; set -u;
+set -euo pipefail
 export LC_ALL=C  # To speed up the $sort$ command.
 READ_LENGTHS_FILE="${INPUT_DIR}/reads-lengths.txt"
 READ_IDS_FILE="${INPUT_DIR}/reads-ids.txt"

diff --git a/scripts/6-repeatAlphabet/2-fixEndBlocks.sh b/scripts/6-repeatAlphabet/2-fixEndBlocks.sh
@@ -23,13 +23,14 @@ N_THREADS=$6
 DELETE_TMP_FILES=$7
 GENOME_LENGTH=$8
 N_HAPLOTYPES=$9
+MAX_KMER_LENGTH_BPS=${10}
 TIGHT_MODE="0"
 SPANNING_BPS="150"  # Bps before and after a k-mer to consider it observed in a read.
 # ------------------------------------ REVANT --------------------------------------------
 REVANT_LIBRARIES="${REVANT_BINARIES}/../lib/*"
 # ----------------------------------------------------------------------------------------
 
-set -o pipefail; set -e; set -u
+set -euo pipefail
 export LC_ALL=C  # To speed up the $sort$ command.
 TMPFILE_NAME="fixEndBlocks-tmp"
 TMPFILE_PATH="${INPUT_DIR}/${TMPFILE_NAME}"
@@ -50,7 +51,7 @@ function enumerateKmersThread() {
 	local LOCAL_BOUNDARIES_FILE=$3
 	local LOCAL_READ_LENGTHS_FILE=$4
 	local LOCAL_KMERS_FILE=$5
-	java ${JAVA_RUNTIME_FLAGS} -classpath "${REVANT_BINARIES}" de.mpi_cbg.revant.apps.CollectKmers 0 ${LOCAL_K} ${LOCAL_TRANSLATED_READS_FILE} ${LOCAL_BOUNDARIES_FILE} ${LOCAL_READ_LENGTHS_FILE} ${ALPHABET_FILE} null null ${LOCAL_KMERS_FILE}
+	java ${JAVA_RUNTIME_FLAGS} -classpath "${REVANT_BINARIES}" de.mpi_cbg.revant.apps.CollectKmers 0 ${LOCAL_K} ${MAX_KMER_LENGTH_BPS} 1 ${LOCAL_TRANSLATED_READS_FILE} ${LOCAL_BOUNDARIES_FILE} ${LOCAL_READ_LENGTHS_FILE} ${ALPHABET_FILE} null null ${LOCAL_KMERS_FILE}
 	if [ $? -ne 0 ]; then
 		exit
 	fi
@@ -63,7 +64,7 @@ function countKmersThread() {
 	local LOCAL_READ_LENGTHS_FILE=$4
 	local LOCAL_KMERS_FILE_INPUT=$5
     local LOCAL_KMERS_FILE_OUTPUT=$6
-	java ${JAVA_RUNTIME_FLAGS} -classpath "${REVANT_BINARIES}" de.mpi_cbg.revant.apps.CollectKmers 1 ${LOCAL_K} ${LOCAL_TRANSLATED_READS_FILE} ${LOCAL_BOUNDARIES_FILE} ${LOCAL_READ_LENGTHS_FILE} ${ALPHABET_FILE} null ${LOCAL_KMERS_FILE_INPUT} ${LOCAL_KMERS_FILE_OUTPUT}
+	java ${JAVA_RUNTIME_FLAGS} -classpath "${REVANT_BINARIES}" de.mpi_cbg.revant.apps.CollectKmers 1 ${LOCAL_K} ${MAX_KMER_LENGTH_BPS} 1 ${LOCAL_TRANSLATED_READS_FILE} ${LOCAL_BOUNDARIES_FILE} ${LOCAL_READ_LENGTHS_FILE} ${ALPHABET_FILE} null ${LOCAL_KMERS_FILE_INPUT} ${LOCAL_KMERS_FILE_OUTPUT}
 	if [ $? -ne 0 ]; then
 		exit
 	fi

diff --git a/scripts/6-repeatAlphabet/3-getUniqueSubstrings.sh b/scripts/6-repeatAlphabet/3-getUniqueSubstrings.sh
@@ -23,6 +23,7 @@ IDENTITY_THRESHOLD=$7
 DISTANCE_THRESHOLD=$8
 CHARACTER_THRESHOLD=$9
 MIN_ALIGNMENT_LENGTH=${10}  # Read-repeat
+MAX_KMER_LENGTH_BPS=${11}
 UNIQUE_MODE="1"  # Non-repetitive blocks are allowed in a k-mer, except at the first/last
 # position of the k-mer. Usually a good choice.
 SPANNING_BPS="150"  # Bps before and after a k-mer to consider it observed in a read.
@@ -31,7 +32,7 @@ REVANT_LIBRARIES="${REVANT_BINARIES}/../lib/*"
 # ----------------------------------------------------------------------------------------
 
 
-set -o pipefail; set -e; set -u
+set -euo pipefail
 export LC_ALL=C  # To speed up the $sort$ command.
 READ_IDS_FILE="${INPUT_DIR}/reads-ids.txt"
 N_READS=$(wc -l < ${READ_IDS_FILE})
@@ -58,7 +59,7 @@ function enumerateKmersThread() {
 	local LOCAL_READ_LENGTHS_FILE=$4
 	local LOCAL_K_MINUS_ONE_INTERVALS_FILE=$5
 	local LOCAL_KMERS_FILE=$6
-	java ${JAVA_RUNTIME_FLAGS} -classpath "${REVANT_BINARIES}" de.mpi_cbg.revant.apps.CollectKmers 0 ${LOCAL_K} ${LOCAL_TRANSLATED_READS_FILE} ${LOCAL_BOUNDARIES_FILE} ${LOCAL_READ_LENGTHS_FILE} ${ALPHABET_FILE} ${LOCAL_K_MINUS_ONE_INTERVALS_FILE} null ${LOCAL_KMERS_FILE}
+	java ${JAVA_RUNTIME_FLAGS} -classpath "${REVANT_BINARIES}" de.mpi_cbg.revant.apps.CollectKmers 0 ${LOCAL_K} ${MAX_KMER_LENGTH_BPS} 2 ${LOCAL_TRANSLATED_READS_FILE} ${LOCAL_BOUNDARIES_FILE} ${LOCAL_READ_LENGTHS_FILE} ${ALPHABET_FILE} ${LOCAL_K_MINUS_ONE_INTERVALS_FILE} null ${LOCAL_KMERS_FILE}
 	if [ $? -ne 0 ]; then
 		exit
 	fi
@@ -72,7 +73,7 @@ function countKmersThread() {
 	local LOCAL_K_MINUS_ONE_INTERVALS_FILE=$5
 	local LOCAL_KMERS_FILE_INPUT=$6
     local LOCAL_KMERS_FILE_OUTPUT=$7
-	java ${JAVA_RUNTIME_FLAGS} -classpath "${REVANT_BINARIES}" de.mpi_cbg.revant.apps.CollectKmers 1 ${LOCAL_K} ${LOCAL_TRANSLATED_READS_FILE} ${LOCAL_BOUNDARIES_FILE} ${LOCAL_READ_LENGTHS_FILE} ${ALPHABET_FILE} ${LOCAL_K_MINUS_ONE_INTERVALS_FILE} ${LOCAL_KMERS_FILE_INPUT} ${LOCAL_KMERS_FILE_OUTPUT}
+	java ${JAVA_RUNTIME_FLAGS} -classpath "${REVANT_BINARIES}" de.mpi_cbg.revant.apps.CollectKmers 1 ${LOCAL_K} ${MAX_KMER_LENGTH_BPS} 2 ${LOCAL_TRANSLATED_READS_FILE} ${LOCAL_BOUNDARIES_FILE} ${LOCAL_READ_LENGTHS_FILE} ${ALPHABET_FILE} ${LOCAL_K_MINUS_ONE_INTERVALS_FILE} ${LOCAL_KMERS_FILE_INPUT} ${LOCAL_KMERS_FILE_OUTPUT}
 	if [ $? -ne 0 ]; then
 		exit
 	fi
@@ -86,7 +87,7 @@ function intervalsThread() {
 	local LOCAL_UNIQUE_KMERS_FILE=$5
 	local LOCAL_K_MINUS_ONE_INTERVALS_FILE=$6
 	local LOCAL_INTERVALS_FILE=$7
-	java ${JAVA_RUNTIME_FLAGS} -classpath "${REVANT_BINARIES}:${REVANT_LIBRARIES}" de.mpi_cbg.revant.apps.GetShortestUniqueIntervals ${LOCAL_K} ${LOCAL_TRANSLATED_READS_FILE} ${LOCAL_BOUNDARIES_FILE} ${LOCAL_READ_LENGTHS_FILE} ${ALPHABET_FILE} ${LOCAL_UNIQUE_KMERS_FILE} ${N_READS} ${AVG_READ_LENGTH} ${GENOME_LENGTH} ${N_HAPLOTYPES} ${MIN_ALIGNMENT_LENGTH} ${IDENTITY_THRESHOLD} ${DISTANCE_THRESHOLD} ${CHARACTER_THRESHOLD} ${LOCAL_K_MINUS_ONE_INTERVALS_FILE} ${LOCAL_INTERVALS_FILE}
+	java ${JAVA_RUNTIME_FLAGS} -classpath "${REVANT_BINARIES}:${REVANT_LIBRARIES}" de.mpi_cbg.revant.apps.GetShortestUniqueIntervals ${LOCAL_K} ${MAX_KMER_LENGTH_BPS} ${LOCAL_TRANSLATED_READS_FILE} ${LOCAL_BOUNDARIES_FILE} ${LOCAL_READ_LENGTHS_FILE} ${ALPHABET_FILE} ${LOCAL_UNIQUE_KMERS_FILE} ${N_READS} ${AVG_READ_LENGTH} ${GENOME_LENGTH} ${N_HAPLOTYPES} ${MIN_ALIGNMENT_LENGTH} ${IDENTITY_THRESHOLD} ${DISTANCE_THRESHOLD} ${CHARACTER_THRESHOLD} ${LOCAL_K_MINUS_ONE_INTERVALS_FILE} ${LOCAL_INTERVALS_FILE}
 	if [ $? -ne 0 ]; then
 		exit
 	fi

diff --git a/scripts/6-repeatAlphabet/4-filterAlignments.sh b/scripts/6-repeatAlphabet/4-filterAlignments.sh
@@ -33,7 +33,7 @@ MIN_BLUE_INTERVAL_LENGTH="300"  # Blue intervals with fewer bps than this are no
 # considered trustworthy.
 # ----------------------------------------------------------------------------------------
 
-set -o pipefail; set -e; set -u
+set -euo pipefail
 READ_LENGTHS_FILE="${INPUT_DIR}/reads-lengths.txt"
 READ_IDS_FILE="${INPUT_DIR}/reads-ids.txt"
 N_READS=$(wc -l < ${READ_IDS_FILE})

diff --git a/scripts/6-repeatAlphabet/master.sh b/scripts/6-repeatAlphabet/master.sh
@@ -18,6 +18,8 @@ HAPLOTYPE_COVERAGE="30"  # Of one haplotype
 MAX_ALIGNMENT_ERROR="0.3"
 MIN_ALIGNMENT_LENGTH_READ_REPEAT="500"
 MIN_ALIGNMENT_LENGTH_READ_READ="500"
+AVG_READ_LENGTH="30000"
+READ_READ_ALIGNMENTS_FILE="${INPUT_DIR}/LAshow-reads-reads.txt"
 # ----------------------------- Properties of CLR reads ----------------------------------
 IS_CLR="0"
 LOW_QUALITY_LENGTH="100"  # >=100. Break reads at low-quality regions >= this length.
@@ -45,8 +47,11 @@ JAVA_RUNTIME_FLAGS="-Xms2G -Xmx10G"
 DELETE_TMP_FILES="0"
 # ----------------------------------------------------------------------------------------
 
-set -o pipefail; set -e; set -u
+set -euo pipefail
 export JAVA_RUNTIME_FLAGS
+MAX_KMER_LENGTH_BPS=$(java ${JAVA_RUNTIME_FLAGS} -classpath "${REVANT_BINARIES}" de.mpi_cbg.revant.apps.GetAlignmentLengthThreshold ${READ_READ_ALIGNMENTS_FILE} 0 ${AVG_READ_LENGTH} ${INPUT_DIR}/histogram-alignmentLength-readRead.txt)
+echo "Histogram of suffix-prefix read-read alignment lengths:"
+cat ${INPUT_DIR}/histogram-alignmentLength-readRead.txt
 if [ ${IS_CLR} -eq 1 ]; then
 	BROKEN_READS="1"
 	./0-breakReads.sh ${INPUT_DIR} ${LOW_QUALITY_LENGTH} ${N_THREADS} ${DELETE_TMP_FILES}
@@ -56,8 +61,8 @@ fi
 ./1-buildAlphabet.sh ${INPUT_DIR} ${BROKEN_READS} ${MAX_ALIGNMENT_ERROR} ${MIN_ALIGNMENT_LENGTH_READ_REPEAT} ${N_HAPLOTYPES} ${HAPLOTYPE_COVERAGE} ${N_THREADS} ${DELETE_TMP_FILES} ${MAX_SPACER_LENGTH} ${WOBBLE_LENGTH} ${FIX_TANDEM_SPACERS} ${CONCATENATE_BLOCKS}
 PERIODIC_ENDPOINTS_FIXED=$(cat ${INPUT_DIR}/buildAlphabet-tmp-return.txt)
 MIN_K_FOR_DISAMBIGUATION="2"; MAX_K_FOR_DISAMBIGUATION="4"
-./2-fixEndBlocks.sh ${INPUT_DIR} ${BROKEN_READS} ${LOW_QUALITY_TYPE} ${MIN_K_FOR_DISAMBIGUATION} ${MAX_K_FOR_DISAMBIGUATION} ${N_THREADS} ${DELETE_TMP_FILES} ${GENOME_LENGTH} ${N_HAPLOTYPES}
-./3-getUniqueSubstrings.sh ${INPUT_DIR} ${GENOME_LENGTH} ${N_HAPLOTYPES} ${MAX_K_UNIQUE} ${N_THREADS} ${DELETE_TMP_FILES} ${IDENTITY_THRESHOLD} ${DISTANCE_THRESHOLD} ${CHARACTER_THRESHOLD} ${MIN_ALIGNMENT_LENGTH_READ_REPEAT}
+./2-fixEndBlocks.sh ${INPUT_DIR} ${BROKEN_READS} ${LOW_QUALITY_TYPE} ${MIN_K_FOR_DISAMBIGUATION} ${MAX_K_FOR_DISAMBIGUATION} ${N_THREADS} ${DELETE_TMP_FILES} ${GENOME_LENGTH} ${N_HAPLOTYPES} ${MAX_KMER_LENGTH_BPS}
+./3-getUniqueSubstrings.sh ${INPUT_DIR} ${GENOME_LENGTH} ${N_HAPLOTYPES} ${MAX_K_UNIQUE} ${N_THREADS} ${DELETE_TMP_FILES} ${IDENTITY_THRESHOLD} ${DISTANCE_THRESHOLD} ${CHARACTER_THRESHOLD} ${MIN_ALIGNMENT_LENGTH_READ_REPEAT} ${MAX_KMER_LENGTH_BPS}
 ./4-filterAlignments.sh ${INPUT_DIR} ${BROKEN_READS} ${PERIODIC_ENDPOINTS_FIXED} ${MIN_ALIGNMENT_LENGTH_READ_READ} ${MIN_ALIGNMENT_LENGTH_READ_REPEAT} ${MAX_K_UNIQUE} ${ALIGNMENT_FILTERING_MODE} ${MIN_INTERSECTION_NONREPETITIVE} ${N_THREADS} ${DELETE_TMP_FILES}
 READ_LENGTHS_FILE="${INPUT_DIR}/reads-lengths.txt"
 N_READS=$(wc -l < ${READ_LENGTHS_FILE})

diff --git a/src/de/mpi_cbg/revant/apps/CollectKmers.java b/src/de/mpi_cbg/revant/apps/CollectKmers.java
@@ -16,13 +16,15 @@ public class CollectKmers {
 	public static void main(String[] args) throws IOException {
         final int MODE = Integer.parseInt(args[0]);
 		final int K = Integer.parseInt(args[1]);
-		final String TRANSLATED_FILE = args[2];
-		final String BOUNDARIES_FILE = args[3];
-		final String READ_LENGTHS_FILE = args[4];
-		final String ALPHABET_FILE = args[5];
-		final String AVOIDED_INTERVALS_FILE = args[6];  // NULL to discard it
-        final String KMERS_FILE_INPUT = args[7];  // NULL to discard it
-		final String KMERS_FILE_OUTPUT = args[8];  // Output file
+        final int MAX_KMER_LENGTH_BPS = Integer.parseInt(args[2]);  // In bps
+        final int UNIQUE_MODE = Integer.parseInt(args[3]);
+		final String TRANSLATED_FILE = args[4];
+		final String BOUNDARIES_FILE = args[5];
+		final String READ_LENGTHS_FILE = args[6];
+		final String ALPHABET_FILE = args[7];
+		final String AVOIDED_INTERVALS_FILE = args[8];  // NULL to discard it
+        final String KMERS_FILE_INPUT = args[9];  // NULL to discard it
+		final String KMERS_FILE_OUTPUT = args[10];  // Output file
 
 		boolean INTERVALS_FILE_EXISTS = !AVOIDED_INTERVALS_FILE.equalsIgnoreCase("null");
 
@@ -69,7 +71,7 @@ public static void main(String[] args) throws IOException {
 			else lastAvoidedInterval=-1;
 			RepeatAlphabet.loadBoundaries(str3);
 			readLength=Integer.parseInt(str4);
-			RepeatAlphabet.getKmers(MODE,str1,K,kmers,null,avoidedIntervals,lastAvoidedInterval,readLength,-1/*argument not used*/,-1/*argument not used*/,-1/*argument not used*/,-1/*argument not used*/,-1/*argument not used*/,-1/*argument not used*/,RepeatAlphabet.boundaries,-1/*argument not used*/,-1/*argument not used*/,-1.0/*argument not used*/,tmpKmer,tmpArray2,tmpArray3,tmpMap,tmpChar);
+			RepeatAlphabet.getKmers(MODE,str1,K,kmers,null,avoidedIntervals,lastAvoidedInterval,UNIQUE_MODE,MAX_KMER_LENGTH_BPS,readLength,-1/*argument not used*/,-1/*argument not used*/,-1/*argument not used*/,-1/*argument not used*/,-1/*argument not used*/,-1/*argument not used*/,RepeatAlphabet.boundaries,-1/*argument not used*/,-1/*argument not used*/,-1.0/*argument not used*/,tmpKmer,tmpArray2,tmpArray3,tmpMap,tmpChar);
 			str1=br1.readLine(); str2=INTERVALS_FILE_EXISTS?br2.readLine():null; 
 			str3=br3.readLine(); str4=br4.readLine(); row++;
 		}

diff --git a/src/de/mpi_cbg/revant/apps/GetAlignmentLengthThreshold.java b/src/de/mpi_cbg/revant/apps/GetAlignmentLengthThreshold.java
@@ -0,0 +1,66 @@
+package de.mpi_cbg.revant.apps;
+
+import java.io.*;
+import de.mpi_cbg.revant.factorize.Alignments;
+import de.mpi_cbg.revant.util.Math;
+import de.mpi_cbg.revant.util.IO;
+
+/**
+ * Prints to STDOUT the max length of an alignment of a given type. This is used for later
+ * enumerating only k-mers that are short enough to be contained in a suffix-prefix
+ * alignment.
+ *
+ * Remark: one could return instead the min L such that, say, 90% of all alignments have
+ * length <=L. However, using such a threshold to avoid considering unique k-mers does
+ * make the assembly graph less connected in practice.
+ */
+public class GetAlignmentLengthThreshold {
+	/**
+	 * @param args 
+	 * 1: alignment type: 0=suffix-prefix overlap; 1=local substring; 2=full containment
+     * or full identity;
+     * 4: histogram of all alignment lengths.
+	 */
+	public static void main(String[] args) throws IOException {
+		final String ALIGNMENTS_FILE = args[0];
+        final int TYPE = Integer.parseInt(args[1]);
+        final int AVG_READ_LENGTH = Integer.parseInt(args[2]);
+        final String OUTPUT_HISTOGRAM = args[3];
+
+        final int QUANTUM = IO.quantum;
+        final int N_CELLS = AVG_READ_LENGTH/QUANTUM;
+        final int IDENTITY_THRESHOLD = QUANTUM;
+
+		int i;
+        int length, lengthA, lengthB, max;
+        long nAlignments;
+		String str;
+        BufferedReader br;
+        BufferedWriter bw;
+        long[] histogram;
+
+        histogram = new long[N_CELLS];
+		br = new BufferedReader(new FileReader(ALIGNMENTS_FILE));
+		str=br.readLine(); str=br.readLine();  // Skipping header
+        str=br.readLine();
+        max=0;
+		while (str!=null)  {
+			Alignments.readAlignmentFile(str);
+            if (Alignments.readAlignmentFile_getType(IDENTITY_THRESHOLD,str)==TYPE) {
+                lengthA=Alignments.endA-Alignments.startA+1;
+                lengthB=Alignments.endB-Alignments.startB+1;
+                length=Math.max(lengthA,lengthB);
+                max=Math.max(max,length);
+                histogram[Math.min(length/QUANTUM,N_CELLS-1)]++;
+            }
+            str=br.readLine();
+        }
+        br.close();
+        bw = new BufferedWriter(new FileWriter(OUTPUT_HISTOGRAM));
+        for (i=0; i<N_CELLS; i++) bw.write((i*QUANTUM)+","+histogram[i]+"\n");
+        bw.write(max+",-1\n");
+        bw.close();
+        System.out.println(max+"");
+	}
+
+}
diff --git a/src/de/mpi_cbg/revant/apps/GetShortestUniqueIntervals.java b/src/de/mpi_cbg/revant/apps/GetShortestUniqueIntervals.java
@@ -21,24 +21,26 @@ public class GetShortestUniqueIntervals {
 
 	public static void main(String[] args) throws IOException {
 		final int K = Integer.parseInt(args[0]);
-		final String TRANSLATED_FILE = args[1];
-		final String BOUNDARIES_FILE = args[2];
-		final String READ_LENGTHS_FILE = args[3];
-		final String ALPHABET_FILE = args[4];
-		final String UNIQUE_KMERS_FILE = args[5];
-        final int N_READS = Integer.parseInt(args[6]);
-        final int AVG_READ_LENGTH = Integer.parseInt(args[7]);
-        final long GENOME_LENGTH = Long.parseLong(args[8]);  // One haplotype
-        final int N_HAPLOTYPES = Integer.parseInt(args[9]);
-        final int MIN_ALIGNMENT_LENGTH = Integer.parseInt(args[10]);  // Read-repeat   
-		final int IDENTITY_THRESHOLD = Integer.parseInt(args[11]);
-		final int DISTANCE_THRESHOLD = Integer.parseInt(args[12]);
-		final double CHARACTER_FRACTION = Double.parseDouble(args[13]);
-		final String OLD_INTERVALS_FILE = args[14];  // NULL to discard it
-		final String NEW_INTERVALS_FILE = args[15];  // Output
+        final int MAX_KMER_LENGTH_BPS = Integer.parseInt(args[1]);
+		final String TRANSLATED_FILE = args[2];
+		final String BOUNDARIES_FILE = args[3];
+		final String READ_LENGTHS_FILE = args[4];
+		final String ALPHABET_FILE = args[5];
+		final String UNIQUE_KMERS_FILE = args[6];
+        final int N_READS = Integer.parseInt(args[7]);
+        final int AVG_READ_LENGTH = Integer.parseInt(args[8]);
+        final long GENOME_LENGTH = Long.parseLong(args[9]);  // One haplotype
+        final int N_HAPLOTYPES = Integer.parseInt(args[10]);
+        final int MIN_ALIGNMENT_LENGTH = Integer.parseInt(args[11]);  // Read-repeat   
+		final int IDENTITY_THRESHOLD = Integer.parseInt(args[12]);
+		final int DISTANCE_THRESHOLD = Integer.parseInt(args[13]);
+		final double CHARACTER_FRACTION = Double.parseDouble(args[14]);
+		final String OLD_INTERVALS_FILE = args[15];  // NULL to discard it
+		final String NEW_INTERVALS_FILE = args[16];  // Output
 
 		boolean OLD_INTERVALS_FILE_EXISTS = !OLD_INTERVALS_FILE.equalsIgnoreCase("null");
         final int MIN_MISSING_LENGTH = IO.quantum;  // Arbitrary
+        final int UNIQUE_MODE = 2;  // Non-repetitive characters not allowed in k-mers
 
 		int i, p;
 		int row, nBlocks, nPairs, lastUniqueInterval, readLength;
@@ -98,7 +100,7 @@ public static void main(String[] args) throws IOException {
 			else lastUniqueInterval=-1;
 			RepeatAlphabet.loadBoundaries(str3);
 			readLength=Integer.parseInt(str4);
-			lastUniqueInterval=RepeatAlphabet.getKmers(2,str1,K,null,kmers,uniqueIntervals,lastUniqueInterval,readLength,N_READS,AVG_READ_LENGTH,GENOME_LENGTH,N_HAPLOTYPES,MIN_ALIGNMENT_LENGTH,MIN_MISSING_LENGTH,RepeatAlphabet.boundaries,IDENTITY_THRESHOLD,DISTANCE_THRESHOLD,CHARACTER_FRACTION,tmpKmer,tmpArray2,tmpArray3,null,tmpChar);            
+			lastUniqueInterval=RepeatAlphabet.getKmers(2,str1,K,null,kmers,uniqueIntervals,lastUniqueInterval,UNIQUE_MODE,MAX_KMER_LENGTH_BPS,readLength,N_READS,AVG_READ_LENGTH,GENOME_LENGTH,N_HAPLOTYPES,MIN_ALIGNMENT_LENGTH,MIN_MISSING_LENGTH,RepeatAlphabet.boundaries,IDENTITY_THRESHOLD,DISTANCE_THRESHOLD,CHARACTER_FRACTION,tmpKmer,tmpArray2,tmpArray3,null,tmpChar);            
 			if (lastUniqueInterval>0) {
 				nPairs=(lastUniqueInterval+1)/3;
 				if (pairs.length<nPairs) {