Manually create ped file for KING

Author: bbimber

SequenceAnalysis/src/org/labkey/sequenceanalysis/run/variant/KingInferenceStep.java

@@ -11,6 +11,7 @@
 import org.jetbrains.annotations.Nullable;
 import org.json.JSONObject;
 import org.labkey.api.pipeline.PipelineJobException;
+import org.labkey.api.reader.Readers;
 import org.labkey.api.sequenceanalysis.SequenceAnalysisService;
 import org.labkey.api.sequenceanalysis.pipeline.AbstractVariantProcessingStepProvider;
 import org.labkey.api.sequenceanalysis.pipeline.PedigreeToolParameterDescriptor;
@@ -24,12 +25,18 @@
 import org.labkey.api.sequenceanalysis.run.AbstractCommandPipelineStep;
 import org.labkey.api.sequenceanalysis.run.AbstractCommandWrapper;
 import org.labkey.api.util.Compress;
+import org.labkey.api.writer.PrintWriters;
 import org.labkey.sequenceanalysis.pipeline.ProcessVariantsHandler;
 
+import java.io.BufferedReader;
 import java.io.File;
 import java.io.IOException;
+import java.io.PrintWriter;
 import java.util.ArrayList;
+import java.util.Arrays;
+import java.util.HashMap;
 import java.util.List;
+import java.util.Map;
 
 public class KingInferenceStep extends AbstractCommandPipelineStep<KingInferenceStep.KingWrapper> implements VariantProcessingStep
 {
@@ -121,19 +128,6 @@ public Output processVariants(File inputVCF, File outputDirectory, ReferenceGeno
         plinkArgs.add("--max-alleles");
         plinkArgs.add("2");
 
-        String demographicsProviderName = getProvider().getParameterByName(PedigreeToolParameterDescriptor.NAME).extractValue(getPipelineCtx().getJob(), getProvider(), getStepIdx());
-        if (demographicsProviderName != null)
-        {
-            File pedFile = ProcessVariantsHandler.getPedigreeFile(getPipelineCtx().getSourceDirectory(true), demographicsProviderName);
-            if (!pedFile.exists())
-            {
-                throw new PipelineJobException("Unable to find pedigree file: " + pedFile.getPath());
-            }
-
-            plinkArgs.add("--ped");
-            plinkArgs.add(pedFile.getPath());
-        }
-
         plinkArgs.add("--out");
         plinkArgs.add(plinkOut.getPath());
 
@@ -166,6 +160,23 @@ public Output processVariants(File inputVCF, File outputDirectory, ReferenceGeno
         kingArgs.add("--prefix");
         kingArgs.add(SequenceAnalysisService.get().getUnzippedBaseName(inputVCF.getName()));
 
+        // Update the pedigree / fam file:
+        String demographicsProviderName = getProvider().getParameterByName(PedigreeToolParameterDescriptor.NAME).extractValue(getPipelineCtx().getJob(), getProvider(), getStepIdx());
+        if (demographicsProviderName != null)
+        {
+            File pedFile = ProcessVariantsHandler.getPedigreeFile(getPipelineCtx().getSourceDirectory(true), demographicsProviderName);
+            if (!pedFile.exists())
+            {
+                throw new PipelineJobException("Unable to find pedigree file: " + pedFile.getPath());
+            }
+
+            File kingFam = createFamFile(pedFile, new File(plinkOutBed.getParentFile(), "plink.fam"), kingArgs);
+            kingArgs.add("--ped");
+            kingArgs.add(kingFam.getPath());
+
+            output.addIntermediateFile(kingFam);
+        }
+
         if (threads != null)
         {
             kingArgs.add("--cpus");
@@ -202,6 +213,58 @@ public Output processVariants(File inputVCF, File outputDirectory, ReferenceGeno
         return output;
     }
 
+    private File createFamFile(File pedFile, File famFile, List<String> kingArgs) throws PipelineJobException
+    {
+        File newFamFile = new File(famFile.getParentFile(), "king.fam");
+
+        Map<String, String> pedMap = new HashMap<>();
+        try (BufferedReader reader = Readers.getReader(pedFile))
+        {
+            String line;
+            while ((line = reader.readLine()) != null)
+            {
+                String[] tokens = line.split(" ");
+                if (tokens.length != 6)
+                {
+                    throw new PipelineJobException("Improper ped line length: " + tokens.length);
+                }
+
+                pedMap.put(tokens[1], StringUtils.join(Arrays.asList("0", tokens[1], tokens[2], tokens[3], tokens[4], "-9"), "\t"));
+            }
+        }
+        catch (IOException e)
+        {
+            throw new PipelineJobException(e);
+        }
+
+        try (BufferedReader reader = Readers.getReader(famFile);PrintWriter writer = PrintWriters.getPrintWriter(newFamFile))
+        {
+            String line;
+            while ((line = reader.readLine()) != null)
+            {
+                String[] tokens = line.split("\t");
+                if (tokens.length != 6)
+                {
+                    throw new PipelineJobException("Improper ped line length: " + tokens.length);
+                }
+
+                String newRow = pedMap.get(tokens[1]);
+                if (newRow == null)
+                {
+                    throw new PipelineJobException("Unable to find pedigree entry for: " + tokens[1]);
+                }
+
+                writer.println(newRow);
+            }
+        }
+        catch (IOException e)
+        {
+            throw new PipelineJobException(e);
+        }
+
+        return newFamFile;
+    }
+
     public static class KingWrapper extends AbstractCommandWrapper
     {
         public KingWrapper(@Nullable Logger logger)