Merge pull request #46 from maize-genetics/fix-maf-to-gvcf-efficiency

Fixing RAM issues with MAFToGVCF.  This should reduce RAM requirements

Author: zrm22

src/main/kotlin/biokotlin/genome/MAFToGVCF.kt

@@ -114,18 +114,19 @@ class MAFToGVCF {
             val sampleName = variantsMap.keys.first()
 
             // sort the variants by contig and position.
-            val variants = variantsMap.values.first().sortedWith(VariantContextComparator(contigList))
+            val variants = variantsMap.values.first().sortedBy { Position(it.chr,it.startPos) }
 
-            exportVariantContext(sampleName, variants, gvcfOutput, refSeqs)
+            exportVariantContext(sampleName, variants, gvcfOutput, refSeqs, outJustGT, delAsSymbolic, maxDeletionSize)
             if (compressAndIndex) {
                 // compress and index the file with bgzip and tabix.
                 compressAndIndexFile(gvcfOutput)
             }
         } else if (variantsMap.size == 2) {
             val outputNames = twoOutputFiles(gvcfOutput)
             variantsMap.entries.forEachIndexed { index, (name, variants) ->
-                val sortedVariants = variants.sortedWith(VariantContextComparator(contigList))
-                val outputFile = exportVariantContext(name, sortedVariants, outputNames[index], refSeqs)
+                val sortedVariants = variants.sortedBy { Position(it.chr,it.startPos) }
+
+                val outputFile = exportVariantContext(name, sortedVariants, outputNames[index], refSeqs, outJustGT, delAsSymbolic, maxDeletionSize)
                 if (compressAndIndex) {
                     compressAndIndexFile(outputNames[index])
                 }
@@ -149,28 +150,28 @@ class MAFToGVCF {
         delAsSymbolic: Boolean = false,
         maxDeletionSize: Int = 0,
         anchorwaveLegacy: Boolean = false
-    ): Map<String, List<VariantContext>> {
+    ): Map<String, List<AssemblyVariantInfo>> {
 
         val mafRecords = loadMAFRecords(mafFile)
         return if (twoGvcfs) {
             val splitGenomes = splitMafRecordsIntoTwoGenomes(mafRecords)
 
             if (splitGenomes.size == 1) {
                 println("getVariantContextsfromMAF:twoGvcfs is true but only 1 genome in the MAF file. Processing the single genome")
-                mapOf(sampleName to buildVariantsForAllAlignments(sampleName, mafRecords, refSeqs, fillGaps,outJustGT,outputType, delAsSymbolic, maxDeletionSize, anchorwaveLegacy))
+                mapOf(sampleName to buildVariantsForAllAlignments(mafRecords, refSeqs, fillGaps, outputType, anchorwaveLegacy))
             } else {
-                splitGenomes.mapIndexed { index, mafrecs ->
+                splitGenomes.mapIndexed { index, splitMafRecords ->
                     //append _1 and _2 to the sampleName for the split genomes
 
                     val genomeName = "${sampleName}_${index + 1}"
                     println("MAFToGVCF:getVariantContextsfromMAF: Splitting ${sampleName} into ${genomeName}")
-                    Pair(genomeName, buildVariantsForAllAlignments(genomeName, mafrecs, refSeqs, fillGaps, outJustGT, outputType, delAsSymbolic, maxDeletionSize, anchorwaveLegacy))
+                    Pair(genomeName, buildVariantsForAllAlignments(splitMafRecords, refSeqs, fillGaps, outputType, anchorwaveLegacy))
                 }.toMap()
             }
 
         } else {
             println("getVariantContextsfromMAF: Processing a single genome")
-            mapOf(sampleName to buildVariantsForAllAlignments(sampleName, mafRecords, refSeqs, fillGaps, outJustGT, outputType, delAsSymbolic, maxDeletionSize, anchorwaveLegacy))
+            mapOf(sampleName to buildVariantsForAllAlignments(mafRecords, refSeqs, fillGaps, outputType, anchorwaveLegacy))
         }
 
     }
@@ -237,16 +238,12 @@ class MAFToGVCF {
      * Function to build the variants for all the alignments.
      */
     fun buildVariantsForAllAlignments(
-        sampleName: String,
         mafRecords: List<MAFRecord>,
         refGenomeSequence: Map<String, NucSeq>,
         fillGaps: Boolean,
-        outJustGT: Boolean,
         outputType: OUTPUT_TYPE,
-        delAsSymbolic: Boolean,
-        maxDeletionSize: Int,
         anchorwaveLegacy: Boolean = false
-    ): List<VariantContext> {
+    ): List<AssemblyVariantInfo> {
         var variantInfos = mutableListOf<AssemblyVariantInfo>()
 
         for (record in mafRecords) {
@@ -263,7 +260,7 @@ class MAFToGVCF {
             variantInfos = fillInMissingVariantBlocks(variantInfos, refGenomeSequence, true)
         }
 
-        return createVariantContextsFromInfo(sampleName, variantInfos, outJustGT, delAsSymbolic, maxDeletionSize)
+        return variantInfos
     }
 
     fun removeRefBlocks(variantInfos: MutableList<AssemblyVariantInfo>) : MutableList<AssemblyVariantInfo> {
@@ -1117,9 +1114,12 @@ class MAFToGVCF {
      */
     fun exportVariantContext(
         sampleName: String,
-        variantContexts: List<VariantContext>,
+        variantContexts: List<AssemblyVariantInfo>,
         outputFileName: String,
-        refGenomeSequence: Map<String, NucSeq>
+        refGenomeSequence: Map<String, NucSeq>,
+        outputJustGT: Boolean,
+        delAsSymbolic: Boolean,
+        maxDeletionSize: Int
     ) {
         val writer = VariantContextWriterBuilder()
             .unsetOption(Options.INDEX_ON_THE_FLY)
@@ -1132,7 +1132,7 @@ class MAFToGVCF {
         addSequenceDictionary(header, refGenomeSequence)
         writer.writeHeader(header)
         for (variant in variantContexts) {
-            writer.add(variant)
+            writer.add(convertVariantInfoToContext(sampleName, variant, outputJustGT, delAsSymbolic, maxDeletionSize))
         }
 
         writer.close()

src/test/kotlin/biokotlin/genome/MAFToGVCFTest.kt

@@ -619,6 +619,24 @@ class MAFToGVCFTest : StringSpec({
 
         }
     }
+    
+    "TestSort" {
+        val chr1 = "1"
+        val chr2 = "2"
+        val chr10 = "10"
+
+        val asmInfos = listOf(AssemblyVariantInfo(chr1,10,20,"0","A","T",true, intArrayOf(10,20), "chr1", 10, 20, "+"),
+            AssemblyVariantInfo(chr2,10,20,"0","A","T",true, intArrayOf(10,20), "chr2", 10, 20, "+"),
+            AssemblyVariantInfo(chr10,10,20,"0","A","T",true, intArrayOf(10,20), "chr10", 10, 20, "+"),
+            AssemblyVariantInfo(chr2,30,40,"0","A","T",true, intArrayOf(10,20), "chr2", 10, 20, "+"),
+            AssemblyVariantInfo(chr1,1,9,"0","A","T",true, intArrayOf(10,20), "chr1", 10, 20, "+"),
+            )
+
+        val expectedInfos = listOf(Pair(chr1,1), Pair(chr1,10), Pair(chr2,10), Pair(chr2,30), Pair(chr10,10))
+
+        val sorted = asmInfos.sortedBy { Position(it.chr,it.startPos) }
+        sorted.map { Pair(it.chr,it.startPos) } shouldBe expectedInfos
+    }
 
 })