Merge branch 'release/v3.0.8'

CHANGELOG

@@ -1,4 +1,24 @@
 
+MiXCR 3.0.8 (22 Jul 2019)
+========================
+
+-- Alignments are forced to the corresponding edge if the same V/J gene is detected in the opposite 
+   PE read mate
+-- Average quality threshold (`-OaverageQualityThreshold=...`) in `assembleContigs` increased to 
+   `20`
+-- Added `exportAlignmentsForClones` action
+-- All indels in homopolimeric stretches are now shifted left in all alignment algorihtms
+-- Smarter base V/J/C hit selection in `assembleContigs`
+-- Read quality trimming (see help for `--trimming-window-size` and `--trimming-quality-threshold` 
+   options in `align`) (disabled by default until 3.1)
+-- Fix exception in `-mutationsDetailed` and similar export options
+-- Fixes in step skipping logic in `analyze` (additional automatic parameter adjustment for VDJC 
+   libraries covering only `VRegion`)
+-- Several more fixes for `assembleContigs`
+-- minor: Fixes incorrect numbers in `assemble` report
+-- minor: Fix for filterAlignmentsAction
+
+
 MiXCR 3.0.7 (20 May 2019)
 ========================
 

itests.sh

@@ -103,6 +103,11 @@ if [[ $create_standard_results == true ]]; then
     go_assemble ${s}_paired
     mixcr align -s hs -r ${s}_single.vdjca.report ${s}_R1.fastq ${s}_single.vdjca
     go_assemble ${s}_single
+
+    mixcr align -s hs -p kAligner2 -r ${s}_paired.vdjca.report ${s}_R1.fastq ${s}_R2.fastq ${s}_paired2.vdjca
+    go_assemble ${s}_paired2
+    mixcr align -s hs -p kAligner2 -r ${s}_single.vdjca.report ${s}_R1.fastq ${s}_single2.vdjca
+    go_assemble ${s}_single2
   done
 fi
 

pom.xml

@@ -33,20 +33,20 @@
 
     <groupId>com.milaboratory</groupId>
     <artifactId>mixcr</artifactId>
-    <version>3.0.7</version>
+    <version>3.0.8</version>
     <packaging>jar</packaging>
     <name>MiXCR</name>
 
     <properties>
         <project.build.sourceEncoding>UTF-8</project.build.sourceEncoding>
-        <milib.version>1.10</milib.version>
+        <milib.version>1.11</milib.version>
     </properties>
 
     <dependencies>
         <dependency>
             <groupId>io.repseq</groupId>
             <artifactId>repseqio</artifactId>
-            <version>1.3.1</version>
+            <version>1.3.2</version>
             <exclusions>
                 <exclusion>
                     <groupId>com.milaboratory</groupId>

src/main/java/com/milaboratory/mixcr/assembler/CloneAssembler.java

@@ -284,7 +284,7 @@ public void runClustering() {
         if (clusteredClonesAccumulators != null)
             throw new IllegalStateException("Already clustered.");
         if (!preClusteringDone)
-            throw new IllegalStateException("No preclustering done.");
+            throw new IllegalStateException("No pre-clustering is done.");
 
         @SuppressWarnings("unchecked")
         Clustering clustering = new Clustering(cloneList,
@@ -329,7 +329,7 @@ public long getAlignmentsCount() {
 
     public void buildClones() {
         if (!preClusteringDone)
-            throw new IllegalStateException("No preclustering done.");
+            throw new IllegalStateException("No pre-clustering is done.");
         ClonesBuilder builder = new ClonesBuilder();
         progressReporter = builder;
         builder.buildClones();
@@ -532,7 +532,7 @@ private final class ClonesBuilder implements CanReportProgress {
         volatile int progress;
 
         private ClonesBuilder() {
-            this.sourceSize = clusteredClonesAccumulators != null ? clusteredClonesAccumulators.size() : clones.size();
+            this.sourceSize = clusteredClonesAccumulators != null ? clusteredClonesAccumulators.size() : cloneList.size();
         }
 
         @Override
@@ -572,8 +572,14 @@ void buildClones() {
                 else {
                     for (TIntIntIterator it = idMapping.iterator(); it.hasNext(); ) {
                         it.advance();
-                        if (newIdMapping.containsKey(it.value()))
-                            it.setValue(newIdMapping.get(it.value()));
+                        int val = it.value();
+                        if (val >= 0) { // "renaming" normal clonotypes
+                            // if (newIdMapping.containsKey(val))
+                            it.setValue(newIdMapping.get(val));
+                        } else { // "renaming" clustered clonotypes
+                            // if (newIdMapping.containsKey(~val))
+                            it.setValue(newIdMapping.get(~val));
+                        }
                     }
                 }
                 source = Arrays.asList(sourceArray);
@@ -671,11 +677,11 @@ public List<CloneAccumulator> build() {
             // Score filtering step
 
             // Calculation
-            float[] maxScores = new float[3];
+            float[] maxScores = new float[2];
             for (CloneAccumulator acc : accs) {
                 if (acc == null)
                     continue;
-                for (int i = 0; i < 3; i++) {
+                for (int i = 0; i < 2; i++) {  // Only for V and J
                     GeneType gt = GeneType.VJC_REFERENCE[i];
                     maxScores[i] =
                             parameters.getSeparateBy(gt)
@@ -684,7 +690,7 @@ public List<CloneAccumulator> build() {
                 }
             }
 
-            for (int i = 0; i < 3; i++)
+            for (int i = 0; i < 2; i++) // Only for V and J
                 maxScores[i] /= parameters.preClusteringScoreFilteringRatio;
 
             // Filtering low score clonotypes
@@ -693,11 +699,12 @@ public List<CloneAccumulator> build() {
                 if (accs[i] == null)
                     continue;
 
-                for (int j = 0; j < 3; j++) {
+                for (int j = 0; j < 2; j++) { // Only for V and J
                     if (accs[i].getBestGene(GeneType.VJC_REFERENCE[j]) != null &&
                             accs[i].getBestScore(GeneType.VJC_REFERENCE[j]) < maxScores[j]) {
                         dropped.accept(accs[i]);
                         accs[i] = null;
+                        ++deleted;
                         break;
                     }
                 }
@@ -706,23 +713,23 @@ public List<CloneAccumulator> build() {
             // Filtering low quality clonotypes
             List<CloneAccumulator> result = new ArrayList<>(accs.length - deleted);
 
-            out:
             for (CloneAccumulator acc : accs) {
                 // null marks clustered clonotypes
                 if (acc == null)
                     continue;
 
                 acc.rebuildClonalSequence();
 
-                if (acc.getSequence().getConcatenated().getQuality().minValue() <
-                        parameters.minimalQuality) {
+                if (acc.getSequence().getConcatenated().getQuality().minValue() < parameters.minimalQuality) {
                     dropped.accept(acc);
-                    continue out;
+                    continue;
                 }
 
                 result.add(acc);
             }
 
+            assert result.size() == accs.length - deleted;
+
             return result;
         }
 

src/main/java/com/milaboratory/mixcr/assembler/fullseq/CoverageAccumulator.java

@@ -0,0 +1,81 @@
+/*
+ * Copyright (c) 2014-2019, Bolotin Dmitry, Chudakov Dmitry, Shugay Mikhail
+ * (here and after addressed as Inventors)
+ * All Rights Reserved
+ *
+ * Permission to use, copy, modify and distribute any part of this program for
+ * educational, research and non-profit purposes, by non-profit institutions
+ * only, without fee, and without a written agreement is hereby granted,
+ * provided that the above copyright notice, this paragraph and the following
+ * three paragraphs appear in all copies.
+ *
+ * Those desiring to incorporate this work into commercial products or use for
+ * commercial purposes should contact MiLaboratory LLC, which owns exclusive
+ * rights for distribution of this program for commercial purposes, using the
+ * following email address: licensing@milaboratory.com.
+ *
+ * IN NO EVENT SHALL THE INVENTORS BE LIABLE TO ANY PARTY FOR DIRECT, INDIRECT,
+ * SPECIAL, INCIDENTAL, OR CONSEQUENTIAL DAMAGES, INCLUDING LOST PROFITS,
+ * ARISING OUT OF THE USE OF THIS SOFTWARE, EVEN IF THE INVENTORS HAS BEEN
+ * ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
+ *
+ * THE SOFTWARE PROVIDED HEREIN IS ON AN "AS IS" BASIS, AND THE INVENTORS HAS
+ * NO OBLIGATION TO PROVIDE MAINTENANCE, SUPPORT, UPDATES, ENHANCEMENTS, OR
+ * MODIFICATIONS. THE INVENTORS MAKES NO REPRESENTATIONS AND EXTENDS NO
+ * WARRANTIES OF ANY KIND, EITHER IMPLIED OR EXPRESS, INCLUDING, BUT NOT
+ * LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY OR FITNESS FOR A
+ * PARTICULAR PURPOSE, OR THAT THE USE OF THE SOFTWARE WILL NOT INFRINGE ANY
+ * PATENT, TRADEMARK OR OTHER RIGHTS.
+ */
+package com.milaboratory.mixcr.assembler.fullseq;
+
+import com.milaboratory.core.Range;
+import com.milaboratory.core.alignment.Alignment;
+import com.milaboratory.core.sequence.NucleotideSequence;
+import com.milaboratory.mixcr.basictypes.VDJCHit;
+import gnu.trove.map.hash.TLongIntHashMap;
+
+import java.util.Arrays;
+import java.util.DoubleSummaryStatistics;
+
+public final class CoverageAccumulator {
+    public final VDJCHit hit;
+    private final long[] coverage;
+
+    public CoverageAccumulator(VDJCHit hit) {
+        this.hit = hit;
+        final int geneLength = hit.getGene().getPartitioning().getLength(hit.getAlignedFeature());
+        this.coverage = new long[geneLength];
+    }
+
+    public void accumulate(VDJCHit hit) {
+        for (int targetId = 0; targetId < hit.numberOfTargets(); targetId++) {
+            Alignment<NucleotideSequence> al = hit.getAlignment(targetId);
+            if (al == null)
+                continue;
+            Range coveredRange = al.getSequence1Range();
+            for (int i = coveredRange.getLower(); i < coveredRange.getUpper(); i++)
+                coverage[i]++;
+        }
+    }
+
+    public DoubleSummaryStatistics getStat() {
+        return Arrays.stream(coverage).mapToDouble(i -> i).summaryStatistics();
+    }
+
+    private final TLongIntHashMap nocCache = new TLongIntHashMap();
+
+    public int getNumberOfCoveredPoints(long threshold) {
+        if (nocCache.containsKey(threshold))
+            return nocCache.get(threshold);
+        else {
+            int c = (int) Arrays.stream(coverage).filter(i -> i >= threshold).count();
+            nocCache.put(threshold, c);
+            return c;
+        }
+    }
+
+    public double getFractionOfCoveredPoints(long threshold) {
+        return 1.0 * getNumberOfCoveredPoints(threshold) / coverage.length;
+    }
+}