Convert segments to families function

Author: mikessh

pom.xml

@@ -4,7 +4,7 @@
 
     <groupId>com.antigenomics</groupId>
     <artifactId>vdjtools</artifactId>
-    <version>1.1.9</version>
+    <version>1.1.10-SNAPSHOT</version>
     <packaging>jar</packaging>
 
     <name>vdjtools</name>
@@ -21,6 +21,7 @@
                 <includes>
                     <include>rscripts/*</include>
                     <include>profile/*</include>
+                    <include>vj_families.txt</include>
                 </includes>
             </resource>
         </resources>

src/main/groovy/com/antigenomics/vdjtools/VdjTools.groovy

@@ -33,6 +33,7 @@ package com.antigenomics.vdjtools
 import com.antigenomics.vdjtools.annotate.Annotate
 import com.antigenomics.vdjtools.annotate.CalcCdrAaStats
 import com.antigenomics.vdjtools.annotate.CalcDegreeStats
+import com.antigenomics.vdjtools.annotate.SegmentsToFamilies
 import com.antigenomics.vdjtools.basic.*
 import com.antigenomics.vdjtools.diversity.CalcDiversityStats
 import com.antigenomics.vdjtools.diversity.PlotQuantileStats
@@ -171,6 +172,8 @@ def getScript = { String scriptName ->
             return new CalcDegreeStats()
         case "ANNOTATE":
             return new Annotate()
+        case "SEGMENTSTOFAMILIES":
+            return new SegmentsToFamilies()
         case "SCANDATABASE":
             println "Moved to VDJdb since 1.0.5, please visit vdjdb.cdr3.net"
             System.exit(0)

src/main/groovy/com/antigenomics/vdjtools/annotate/SegmentsToFamilies.groovy

@@ -0,0 +1,98 @@
+package com.antigenomics.vdjtools.annotate
+
+import com.antigenomics.vdjtools.io.SampleWriter
+import com.antigenomics.vdjtools.misc.CommonUtil
+import com.antigenomics.vdjtools.sample.Sample
+import com.antigenomics.vdjtools.sample.SampleCollection
+import com.antigenomics.vdjtools.sample.SegmentConverter
+
+def cli = new CliBuilder(usage: "SegmentsToFamilies [options] " +
+        "[sample1 sample2 ... if not -m] output_prefix")
+cli.h("display help message")
+cli.m(longOpt: "metadata", argName: "filename", args: 1,
+        "Metadata file. First and second columns should contain file name and sample id. " +
+                "Header is mandatory and will be used to assign column names for metadata.")
+cli.s(longOpt: "species", argName: "name", args: 1,
+        "Dataset species, 'human' or 'mouse'.", required: true)
+cli.c(longOpt: "compress", "Compress output sample files.")
+
+
+def opt = cli.parse(args)
+
+if (opt == null) {
+    System.exit(2)
+}
+
+if (opt.h) {
+    cli.usage()
+    System.exit(2)
+}
+
+// Check if metadata is provided
+
+def metadataFileName = opt.m
+
+if (metadataFileName ? opt.arguments().size() != 1 : opt.arguments().size() < 2) {
+    if (metadataFileName)
+        println "Only output prefix should be provided in case of -m"
+    else
+        println "At least 1 sample files should be provided if not using -m"
+    cli.usage()
+    System.exit(2)
+}
+
+def outputFilePrefix = opt.arguments()[-1],
+    species = (String) opt.s,
+    compress = (boolean) opt.c
+
+if (!["human", "mouse"].any { it.equalsIgnoreCase(species) }) {
+    println "Should specify either human or mouse as species."
+    System.exit(2)
+}
+
+def scriptName = getClass().canonicalName.split("\\.")[-1]
+
+//
+// Batch load samples
+//
+
+println "[${new Date()} $scriptName] Reading samples"
+
+def sampleCollection = metadataFileName ?
+        new SampleCollection((String) metadataFileName) :
+        new SampleCollection(opt.arguments()[0..-2])
+
+println "[${new Date()} $scriptName] ${sampleCollection.size()} samples prepared"
+
+// Load segment conversions
+def vSegmentMap = new HashMap<String, String>(),
+    jSegmentMap = new HashMap<String, String>()
+CommonUtil.resourceStreamReader("vj_families.txt").splitEachLine("\t") {
+    if (it[0].equalsIgnoreCase(species)) {
+        if (it[2].equalsIgnoreCase("v")) {
+            vSegmentMap.put(it[3], it[4])
+        } else {
+            jSegmentMap.put(it[3], it[4])
+        }
+    }
+}
+
+def converter = new SegmentConverter(vSegmentMap, jSegmentMap)
+
+//
+// Iterate over samples and change V segments
+//
+def sw = new SampleWriter(compress)
+
+sampleCollection.eachWithIndex { sample, ind ->
+    def sampleId = sample.sampleMetadata.sampleId
+    println "[${new Date()} $scriptName] Changing segments for $sampleId.."
+
+    // print output
+    sw.writeConventional(new Sample(sample, converter), outputFilePrefix)
+}
+
+sampleCollection.metadataTable.storeWithOutput(outputFilePrefix, compress,
+        "segm2fam")
+
+println "[${new Date()} $scriptName] Finished"

src/main/java/com/antigenomics/vdjtools/sample/Clonotype.java

@@ -497,6 +497,15 @@ public int hashCode() {
         return result;
     }
 
+    Clonotype withSegments(Segment v, Segment j) {
+        return new Clonotype(parent, count, freq,
+                segmPoints, v, d, j,
+                cdr3nt, cdr3aa,
+                inFrame, noStop, isComplete,
+                annotation
+        );
+    }
+
     @Override
     public Clonotype getClonotype() {
         return this;

src/main/java/com/antigenomics/vdjtools/sample/ClonotypeConverter.java

@@ -0,0 +1,6 @@
+package com.antigenomics.vdjtools.sample;
+
+
+public interface ClonotypeConverter {
+    Clonotype convert(Clonotype clonotype);
+}

src/main/java/com/antigenomics/vdjtools/sample/Sample.java

@@ -88,6 +88,15 @@ public Sample(Sample other, Map<Clonotype, Integer> samplerMap) {
         Collections.sort(clonotypes);
     }
 
+    public Sample(Sample other, ClonotypeConverter clonotypeConverter) {
+        this.sampleMetadata = other.sampleMetadata;
+        this.annotationHeader = other.annotationHeader;
+
+        for (Clonotype clonotype : other.clonotypes) {
+            this.addClonotype(clonotypeConverter.convert(clonotype));
+        }
+    }
+
     /**
      * Creates a new sample by filtering and selecting top N clonotypes from the specified sample.
      *

src/main/java/com/antigenomics/vdjtools/sample/SegmentConverter.java

@@ -0,0 +1,31 @@
+package com.antigenomics.vdjtools.sample;
+
+import com.antigenomics.vdjtools.misc.Segment;
+import com.antigenomics.vdjtools.misc.SegmentFactory;
+
+import java.util.HashMap;
+import java.util.Map;
+
+public class SegmentConverter implements ClonotypeConverter {
+    private final Map<String, Segment> vSegmentMap = new HashMap<>(),
+            jSegmentMap = new HashMap<>();
+
+    public SegmentConverter(Map<String, String> vSegmentMap, Map<String, String> jSegmentMap) {
+        for (Map.Entry<String, String> conv : vSegmentMap.entrySet()) {
+            this.vSegmentMap.put(conv.getKey(),
+                    SegmentFactory.INSTANCE.create(conv.getValue()));
+        }
+
+        for (Map.Entry<String, String> conv : jSegmentMap.entrySet()) {
+            this.jSegmentMap.put(conv.getKey(),
+                    SegmentFactory.INSTANCE.create(conv.getValue()));
+        }
+    }
+
+    @Override
+    public Clonotype convert(Clonotype clonotype) {
+        return clonotype.withSegments(
+                vSegmentMap.getOrDefault(clonotype.getV(), clonotype.getVBinary()),
+                jSegmentMap.getOrDefault(clonotype.getJ(), clonotype.getJBinary()));
+    }
+}