|
| 1 | +package org.labkey.sequenceanalysis.run.variant; |
| 2 | + |
| 3 | +import htsjdk.samtools.SAMSequenceDictionary; |
| 4 | +import htsjdk.samtools.SAMSequenceRecord; |
| 5 | +import htsjdk.samtools.util.Interval; |
| 6 | +import htsjdk.variant.utils.SAMSequenceDictionaryExtractor; |
| 7 | +import org.apache.commons.io.FileUtils; |
| 8 | +import org.apache.commons.lang3.StringUtils; |
| 9 | +import org.apache.commons.lang3.math.NumberUtils; |
| 10 | +import org.apache.logging.log4j.Logger; |
| 11 | +import org.json.old.JSONObject; |
| 12 | +import org.labkey.api.pipeline.PipelineJobException; |
| 13 | +import org.labkey.api.sequenceanalysis.SequenceAnalysisService; |
| 14 | +import org.labkey.api.sequenceanalysis.pipeline.AbstractVariantProcessingStepProvider; |
| 15 | +import org.labkey.api.sequenceanalysis.pipeline.PipelineContext; |
| 16 | +import org.labkey.api.sequenceanalysis.pipeline.PipelineStepProvider; |
| 17 | +import org.labkey.api.sequenceanalysis.pipeline.ReferenceGenome; |
| 18 | +import org.labkey.api.sequenceanalysis.pipeline.SequencePipelineService; |
| 19 | +import org.labkey.api.sequenceanalysis.pipeline.ToolParameterDescriptor; |
| 20 | +import org.labkey.api.sequenceanalysis.pipeline.VariantProcessingStep; |
| 21 | +import org.labkey.api.sequenceanalysis.pipeline.VariantProcessingStepOutputImpl; |
| 22 | +import org.labkey.api.sequenceanalysis.run.AbstractCommandPipelineStep; |
| 23 | +import org.labkey.api.sequenceanalysis.run.AbstractCommandWrapper; |
| 24 | + |
| 25 | +import javax.annotation.Nullable; |
| 26 | +import java.io.File; |
| 27 | +import java.io.IOException; |
| 28 | +import java.util.ArrayList; |
| 29 | +import java.util.Arrays; |
| 30 | +import java.util.List; |
| 31 | + |
| 32 | +public class KingInferenceStep extends AbstractCommandPipelineStep<KingInferenceStep.KingWrapper> implements VariantProcessingStep |
| 33 | +{ |
| 34 | + public KingInferenceStep(PipelineStepProvider<?> provider, PipelineContext ctx) |
| 35 | + { |
| 36 | + super(provider, ctx, new KingInferenceStep.KingWrapper(ctx.getLogger())); |
| 37 | + } |
| 38 | + |
| 39 | + public static class Provider extends AbstractVariantProcessingStepProvider<KingInferenceStep> |
| 40 | + { |
| 41 | + public Provider() |
| 42 | + { |
| 43 | + super("KingInferenceStep", "KING/Relatedness", "", "This will run KING to infer kinship from a VCF", Arrays.asList( |
| 44 | + ToolParameterDescriptor.create("limitToChromosomes", "Limit to Chromosomes", "If checked, the analysis will include only the primary chromosomes", "checkbox", new JSONObject(){{ |
| 45 | + put("checked", true); |
| 46 | + }}, true) |
| 47 | + ), null, "https://www.kingrelatedness.com/manual.shtml"); |
| 48 | + } |
| 49 | + |
| 50 | + @Override |
| 51 | + public KingInferenceStep create(PipelineContext ctx) |
| 52 | + { |
| 53 | + return new KingInferenceStep(this, ctx); |
| 54 | + } |
| 55 | + } |
| 56 | + |
| 57 | + @Override |
| 58 | + public Output processVariants(File inputVCF, File outputDirectory, ReferenceGenome genome, @Nullable List<Interval> intervals) throws PipelineJobException |
| 59 | + { |
| 60 | + VariantProcessingStepOutputImpl output = new VariantProcessingStepOutputImpl(); |
| 61 | + |
| 62 | + output.addInput(inputVCF, "Input VCF"); |
| 63 | + output.addInput(genome.getWorkingFastaFile(), "Reference Genome"); |
| 64 | + |
| 65 | + File plinkOut = new File(outputDirectory, "plink"); |
| 66 | + output.addIntermediateFile(new File(plinkOut.getPath() + ".bed")); |
| 67 | + output.addIntermediateFile(new File(plinkOut.getPath() + ".fam")); |
| 68 | + output.addIntermediateFile(new File(plinkOut.getPath() + ".bim")); |
| 69 | + output.addIntermediateFile(new File(plinkOut.getPath() + ".log")); |
| 70 | + output.addIntermediateFile(new File(plinkOut.getPath() + "-temporary.psam")); |
| 71 | + |
| 72 | + PlinkPcaStep.PlinkWrapper plink = new PlinkPcaStep.PlinkWrapper(getPipelineCtx().getLogger()); |
| 73 | + List<String> plinkArgs = new ArrayList<>(); |
| 74 | + plinkArgs.add(plink.getExe().getPath()); |
| 75 | + plinkArgs.add("--vcf"); |
| 76 | + plinkArgs.add(inputVCF.getPath()); |
| 77 | + |
| 78 | + plinkArgs.add("--make-bed"); |
| 79 | + |
| 80 | + boolean limitToChromosomes = getProvider().getParameterByName("limitToChromosomes").extractValue(getPipelineCtx().getJob(), getProvider(), getStepIdx(), Boolean.class, true); |
| 81 | + if (limitToChromosomes) |
| 82 | + { |
| 83 | + SAMSequenceDictionary dict = SAMSequenceDictionaryExtractor.extractDictionary(genome.getSequenceDictionary().toPath()); |
| 84 | + List<String> toKeep = dict.getSequences().stream().filter(s -> { |
| 85 | + String name = StringUtils.replaceIgnoreCase(s.getSequenceName(), "^chr", ""); |
| 86 | + |
| 87 | + return NumberUtils.isCreatable(name) || "X".equalsIgnoreCase(name) || "Y".equalsIgnoreCase(name); |
| 88 | + }).map(SAMSequenceRecord::getSequenceName).toList(); |
| 89 | + |
| 90 | + plinkArgs.add("--chr"); |
| 91 | + plinkArgs.add(StringUtils.join(toKeep, ",")); |
| 92 | + } |
| 93 | + else |
| 94 | + { |
| 95 | + plinkArgs.add("--allow-extra-chr"); |
| 96 | + } |
| 97 | + |
| 98 | + plinkArgs.add("--silent"); |
| 99 | + |
| 100 | + plinkArgs.add("--max-alleles"); |
| 101 | + plinkArgs.add("2"); |
| 102 | + |
| 103 | + plinkArgs.add("--out"); |
| 104 | + plinkArgs.add(plinkOut.getPath()); |
| 105 | + |
| 106 | + Integer threads = SequencePipelineService.get().getMaxThreads(getPipelineCtx().getLogger()); |
| 107 | + if (threads != null) |
| 108 | + { |
| 109 | + plinkArgs.add("--threads"); |
| 110 | + plinkArgs.add(threads.toString()); |
| 111 | + } |
| 112 | + |
| 113 | + //TODO: consider --memory (in MB) |
| 114 | + |
| 115 | + plink.execute(plinkArgs); |
| 116 | + |
| 117 | + File plinkOutBed = new File(plinkOut.getPath() + ".bed"); |
| 118 | + if (!plinkOutBed.exists()) |
| 119 | + { |
| 120 | + throw new PipelineJobException("Unable to find file: " + plinkOutBed.getPath()); |
| 121 | + } |
| 122 | + |
| 123 | + KingWrapper wrapper = new KingWrapper(getPipelineCtx().getLogger()); |
| 124 | + wrapper.setWorkingDir(outputDirectory); |
| 125 | + |
| 126 | + List<String> kingArgs = new ArrayList<>(); |
| 127 | + kingArgs.add(wrapper.getExe().getPath()); |
| 128 | + |
| 129 | + kingArgs.add("-b"); |
| 130 | + kingArgs.add(plinkOutBed.getPath()); |
| 131 | + |
| 132 | + kingArgs.add("--prefix"); |
| 133 | + kingArgs.add(SequenceAnalysisService.get().getUnzippedBaseName(inputVCF.getName())); |
| 134 | + |
| 135 | + if (threads != null) |
| 136 | + { |
| 137 | + kingArgs.add("--cpus"); |
| 138 | + kingArgs.add(threads.toString()); |
| 139 | + } |
| 140 | + |
| 141 | + kingArgs.add("--kinship"); |
| 142 | + |
| 143 | + File kinshipOutput = new File(outputDirectory, SequenceAnalysisService.get().getUnzippedBaseName(inputVCF.getName()) + ".kin"); |
| 144 | + wrapper.execute(kingArgs); |
| 145 | + if (!kinshipOutput.exists()) |
| 146 | + { |
| 147 | + throw new PipelineJobException("Unable to find file: " + kinshipOutput.getPath()); |
| 148 | + } |
| 149 | + |
| 150 | + File kinshipOutputTxt = new File(kinshipOutput.getPath() + ".txt"); |
| 151 | + if (kinshipOutputTxt.exists()) |
| 152 | + { |
| 153 | + kinshipOutputTxt.delete(); |
| 154 | + } |
| 155 | + |
| 156 | + try |
| 157 | + { |
| 158 | + FileUtils.moveFile(kinshipOutput, kinshipOutputTxt); |
| 159 | + } |
| 160 | + catch (IOException e) |
| 161 | + { |
| 162 | + throw new PipelineJobException(e); |
| 163 | + } |
| 164 | + |
| 165 | + output.addSequenceOutput(kinshipOutputTxt, "King Relatedness: " + inputVCF.getName(), "KING Relatedness", null, null, genome.getGenomeId(), null); |
| 166 | + |
| 167 | + return output; |
| 168 | + } |
| 169 | + |
| 170 | + public static class KingWrapper extends AbstractCommandWrapper |
| 171 | + { |
| 172 | + public KingWrapper(@Nullable Logger logger) |
| 173 | + { |
| 174 | + super(logger); |
| 175 | + } |
| 176 | + |
| 177 | + public File getExe() |
| 178 | + { |
| 179 | + return SequencePipelineService.get().getExeForPackage("KINGPATH", "king"); |
| 180 | + } |
| 181 | + } |
| 182 | +} |
0 commit comments