Skip to content

Commit 669cb6e

Browse files
bbimberbbimber
committed
Initial support for velocyto
1 parent aa9547e commit 669cb6e

File tree

5 files changed

+336
-64
lines changed

5 files changed

+336
-64
lines changed

‎singlecell/src/org/labkey/singlecell/SingleCellModule.java‎

Lines changed: 4 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -44,6 +44,8 @@
4444
import org.labkey.singlecell.run.CellRangerVDJWrapper;
4545
import org.labkey.singlecell.run.NimbleAligner;
4646
import org.labkey.singlecell.run.NimbleAlignmentStep;
47+
import org.labkey.singlecell.run.VelocytoAlignmentStep;
48+
import org.labkey.singlecell.run.VelocytoPostProcessingStep;
4749

4850
import java.util.Collection;
4951
import java.util.Collections;
@@ -132,6 +134,8 @@ public static void registerPipelineSteps()
132134
SequencePipelineService.get().registerPipelineStep(new CellRangerGexCountStep.Provider());
133135
SequencePipelineService.get().registerPipelineStep(new CellRangerVDJWrapper.VDJProvider());
134136
SequencePipelineService.get().registerPipelineStep(new NimbleAligner.Provider());
137+
SequencePipelineService.get().registerPipelineStep(new VelocytoAlignmentStep.Provider());
138+
SequencePipelineService.get().registerPipelineStep(new VelocytoPostProcessingStep.Provider());
135139

136140
SequenceAnalysisService.get().registerReadsetHandler(new CellRangerFeatureBarcodeHandler());
137141

‎singlecell/src/org/labkey/singlecell/run/AbstractCellRangerDependentStep.java‎

Lines changed: 93 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,93 @@
1+
package org.labkey.singlecell.run;
2+
3+
import org.apache.commons.io.FileUtils;
4+
import org.apache.commons.lang3.StringUtils;
5+
import org.jetbrains.annotations.Nullable;
6+
import org.labkey.api.pipeline.PipelineJobException;
7+
import org.labkey.api.sequenceanalysis.model.Readset;
8+
import org.labkey.api.sequenceanalysis.pipeline.AlignmentOutputImpl;
9+
import org.labkey.api.sequenceanalysis.pipeline.AlignmentStepProvider;
10+
import org.labkey.api.sequenceanalysis.pipeline.PipelineContext;
11+
import org.labkey.api.sequenceanalysis.pipeline.ReferenceGenome;
12+
13+
import java.io.File;
14+
import java.io.IOException;
15+
import java.util.List;
16+
17+
public class AbstractCellRangerDependentStep extends CellRangerGexCountStep
18+
{
19+
public AbstractCellRangerDependentStep(AlignmentStepProvider provider, PipelineContext ctx, CellRangerWrapper wrapper)
20+
{
21+
super(provider, ctx, wrapper);
22+
}
23+
24+
@Override
25+
public IndexOutput createIndex(ReferenceGenome referenceGenome, File outputDir) throws PipelineJobException
26+
{
27+
return super.createIndex(referenceGenome, outputDir);
28+
}
29+
30+
protected File runCellRanger(AlignmentOutputImpl output, Readset rs, List<File> inputFastqs1, @Nullable List<File> inputFastqs2, File outputDirectory, ReferenceGenome referenceGenome, String basename, String readGroupId, @Nullable String platformUnit) throws PipelineJobException
31+
{
32+
File localBam = new File(outputDirectory, basename + ".cellranger.bam");
33+
File localBamIdx = new File(localBam.getPath() + ".bai");
34+
35+
36+
String idParam = StringUtils.trimToNull(getProvider().getParameterByName("id").extractValue(getPipelineCtx().getJob(), getProvider(), getStepIdx(), String.class));
37+
File cellrangerOutdir = new File(outputDirectory, CellRangerWrapper.getId(idParam, rs));
38+
39+
if (localBam.exists() && localBamIdx.exists())
40+
{
41+
getPipelineCtx().getLogger().info("Existing BAM found, re-using: " + localBam.getPath());
42+
}
43+
else
44+
{
45+
File crBam = new File(cellrangerOutdir, "outs/possorted_genome_bam.bam");
46+
if (crBam.exists())
47+
{
48+
getPipelineCtx().getLogger().info("Using previous cellranger count run");
49+
}
50+
else
51+
{
52+
getPipelineCtx().getLogger().info("Running cellranger");
53+
AlignmentOutput crOutput = super.performAlignment(rs, inputFastqs1, inputFastqs2, outputDirectory, referenceGenome, basename, readGroupId, platformUnit);
54+
crBam = crOutput.getBAM();
55+
56+
// Remove all the normal 10x outputs:
57+
output.addCommandsExecuted(crOutput.getCommandsExecuted());
58+
output.addIntermediateFiles(crOutput.getIntermediateFiles());
59+
}
60+
61+
// Remove the whole 10x folder:
62+
output.addIntermediateFile(cellrangerOutdir);
63+
64+
try
65+
{
66+
if (localBam.exists())
67+
{
68+
localBam.delete();
69+
}
70+
FileUtils.moveFile(crBam, localBam);
71+
72+
if (localBamIdx.exists())
73+
{
74+
localBamIdx.delete();
75+
}
76+
FileUtils.moveFile(new File(crBam.getPath() + ".bai"), localBamIdx);
77+
}
78+
catch (IOException e)
79+
{
80+
throw new PipelineJobException(e);
81+
}
82+
}
83+
84+
return localBam;
85+
}
86+
87+
@Override
88+
public boolean alwaysCopyIndexToWorkingDir()
89+
{
90+
return false;
91+
}
92+
93+
}

‎singlecell/src/org/labkey/singlecell/run/NimbleAligner.java‎

Lines changed: 3 additions & 64 deletions
Original file line numberDiff line numberDiff line change
@@ -45,7 +45,7 @@
4545
import java.util.Map;
4646
import java.util.stream.Collectors;
4747

48-
public class NimbleAligner extends CellRangerGexCountStep
48+
public class NimbleAligner extends AbstractCellRangerDependentStep
4949
{
5050
public static final String REF_GENOMES = "refGenomes";
5151

@@ -73,66 +73,11 @@ public AlignmentStep create(PipelineContext context)
7373
}
7474
}
7575

76-
@Override
77-
public IndexOutput createIndex(ReferenceGenome referenceGenome, File outputDir) throws PipelineJobException
78-
{
79-
return super.createIndex(referenceGenome, outputDir);
80-
}
81-
8276
@Override
8377
public AlignmentOutput performAlignment(Readset rs, List<File> inputFastqs1, @Nullable List<File> inputFastqs2, File outputDirectory, ReferenceGenome referenceGenome, String basename, String readGroupId, @Nullable String platformUnit) throws PipelineJobException
8478
{
85-
File localBam = new File(outputDirectory, basename + ".cellranger.bam");
86-
File localBamIdx = new File(localBam.getPath() + ".bai");
87-
AlignmentOutputImpl output = new AlignmentOutputImpl();;
88-
89-
String idParam = StringUtils.trimToNull(getProvider().getParameterByName("id").extractValue(getPipelineCtx().getJob(), getProvider(), getStepIdx(), String.class));
90-
File cellrangerOutdir = new File(outputDirectory, CellRangerWrapper.getId(idParam, rs));
91-
92-
if (localBam.exists() && localBamIdx.exists())
93-
{
94-
getPipelineCtx().getLogger().info("Existing BAM found, re-using: " + localBam.getPath());
95-
}
96-
else
97-
{
98-
File crBam = new File(cellrangerOutdir, "outs/possorted_genome_bam.bam");
99-
if (crBam.exists())
100-
{
101-
getPipelineCtx().getLogger().info("Using previous cellranger count run");
102-
}
103-
else
104-
{
105-
getPipelineCtx().getLogger().info("Running cellranger");
106-
AlignmentOutput crOutput = super.performAlignment(rs, inputFastqs1, inputFastqs2, outputDirectory, referenceGenome, basename, readGroupId, platformUnit);
107-
crBam = crOutput.getBAM();
108-
109-
// Remove all the normal 10x outputs:
110-
output.addCommandsExecuted(crOutput.getCommandsExecuted());
111-
output.addIntermediateFiles(crOutput.getIntermediateFiles());
112-
}
113-
114-
// Remove the whole 10x folder:
115-
output.addIntermediateFile(cellrangerOutdir);
116-
117-
try
118-
{
119-
if (localBam.exists())
120-
{
121-
localBam.delete();
122-
}
123-
FileUtils.moveFile(crBam, localBam);
124-
125-
if (localBamIdx.exists())
126-
{
127-
localBamIdx.delete();
128-
}
129-
FileUtils.moveFile(new File(crBam.getPath() + ".bai"), localBamIdx);
130-
}
131-
catch (IOException e)
132-
{
133-
throw new PipelineJobException(e);
134-
}
135-
}
79+
AlignmentOutputImpl output = new AlignmentOutputImpl();
80+
File localBam = runCellRanger(output, rs, inputFastqs1, inputFastqs2, outputDirectory, referenceGenome, basename, readGroupId, platformUnit);
13681

13782
// Now run nimble itself:
13883
doNimbleAlign(localBam, output, rs, basename);
@@ -141,12 +86,6 @@ public AlignmentOutput performAlignment(Readset rs, List<File> inputFastqs1, @Nu
14186
return output;
14287
}
14388

144-
@Override
145-
public boolean alwaysCopyIndexToWorkingDir()
146-
{
147-
return false;
148-
}
149-
15089
@Override
15190
public void init(SequenceAnalysisJobSupport support) throws PipelineJobException
15291
{

‎singlecell/src/org/labkey/singlecell/run/VelocytoAlignmentStep.java‎

Lines changed: 147 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,147 @@
1+
package org.labkey.singlecell.run;
2+
3+
import org.apache.logging.log4j.Logger;
4+
import org.jetbrains.annotations.Nullable;
5+
import org.json.JSONObject;
6+
import org.labkey.api.pipeline.PipelineJobException;
7+
import org.labkey.api.sequenceanalysis.model.Readset;
8+
import org.labkey.api.sequenceanalysis.pipeline.AbstractAlignmentStepProvider;
9+
import org.labkey.api.sequenceanalysis.pipeline.AlignmentOutputImpl;
10+
import org.labkey.api.sequenceanalysis.pipeline.AlignmentStep;
11+
import org.labkey.api.sequenceanalysis.pipeline.AlignmentStepProvider;
12+
import org.labkey.api.sequenceanalysis.pipeline.PipelineContext;
13+
import org.labkey.api.sequenceanalysis.pipeline.ReferenceGenome;
14+
import org.labkey.api.sequenceanalysis.pipeline.SequencePipelineService;
15+
import org.labkey.api.sequenceanalysis.pipeline.ToolParameterDescriptor;
16+
import org.labkey.api.sequenceanalysis.run.AbstractCommandWrapper;
17+
import org.labkey.api.sequenceanalysis.run.SimpleScriptWrapper;
18+
import org.labkey.api.util.PageFlowUtil;
19+
20+
import java.io.File;
21+
import java.util.ArrayList;
22+
import java.util.Arrays;
23+
import java.util.LinkedHashSet;
24+
import java.util.List;
25+
26+
public class VelocytoAlignmentStep extends AbstractCellRangerDependentStep
27+
{
28+
public VelocytoAlignmentStep(AlignmentStepProvider provider, PipelineContext ctx, CellRangerWrapper wrapper)
29+
{
30+
super(provider, ctx, wrapper);
31+
}
32+
33+
public static class Provider extends AbstractAlignmentStepProvider<AlignmentStep>
34+
{
35+
public Provider()
36+
{
37+
super("velocyto", "This will run velocyto to generate a supplemental feature count matrix", getCellRangerGexParams(Arrays.asList(
38+
ToolParameterDescriptor.createExpDataParam("gtf", "Gene File", "This is the ID of a GTF file containing genes from this genome.", "sequenceanalysis-genomefileselectorfield", new JSONObject()
39+
{{
40+
put("extensions", Arrays.asList("gtf"));
41+
put("width", 400);
42+
put("allowBlank", false);
43+
}}, null),
44+
ToolParameterDescriptor.createExpDataParam("mask", "Mask File", "This is the ID of an optional GTF file containing repetitive regions to mask.", "sequenceanalysis-genomefileselectorfield", new JSONObject()
45+
{{
46+
put("extensions", Arrays.asList("gtf"));
47+
put("width", 400);
48+
put("allowBlank", false);
49+
}}, null)
50+
)), new LinkedHashSet<>(PageFlowUtil.set("sequenceanalysis/field/GenomeFileSelectorField.js")), null, true, false, ALIGNMENT_MODE.MERGE_THEN_ALIGN);
51+
}
52+
53+
@Override
54+
public AlignmentStep create(PipelineContext context)
55+
{
56+
return new VelocytoAlignmentStep(this, context, new CellRangerWrapper(context.getLogger()));
57+
}
58+
}
59+
60+
@Override
61+
public AlignmentOutput performAlignment(Readset rs, List<File> inputFastqs1, @Nullable List<File> inputFastqs2, File outputDirectory, ReferenceGenome referenceGenome, String basename, String readGroupId, @Nullable String platformUnit) throws PipelineJobException
62+
{
63+
AlignmentOutputImpl output = new AlignmentOutputImpl();
64+
File localBam = runCellRanger(output, rs, inputFastqs1, inputFastqs2, outputDirectory, referenceGenome, basename, readGroupId, platformUnit);
65+
66+
File gtf = getPipelineCtx().getSequenceSupport().getCachedData(getProvider().getParameterByName("gtf").extractValue(getPipelineCtx().getJob(), getProvider(), getStepIdx(), Integer.class));
67+
if (gtf == null)
68+
{
69+
throw new PipelineJobException("Missing GTF file param");
70+
}
71+
else if (!gtf.exists())
72+
{
73+
throw new PipelineJobException("File not found: " + gtf.getPath());
74+
}
75+
76+
File mask = null;
77+
if (getProvider().getParameterByName("mask").extractValue(getPipelineCtx().getJob(), getProvider(), getStepIdx(), Integer.class) != null)
78+
{
79+
mask = getPipelineCtx().getSequenceSupport().getCachedData(getProvider().getParameterByName("mask").extractValue(getPipelineCtx().getJob(), getProvider(), getStepIdx(), Integer.class));
80+
if (!mask.exists())
81+
{
82+
throw new PipelineJobException("Missing file: " + mask.getPath());
83+
}
84+
}
85+
86+
File loom = new VelocytoWrapper(getPipelineCtx().getLogger()).runVelocytoFor10x(localBam, gtf, outputDirectory, mask);
87+
output.addSequenceOutput(loom, rs.getName() + ": velocyto", "Velocyto Counts", rs.getReadsetId(), null, referenceGenome.getGenomeId(), null);
88+
89+
return output;
90+
}
91+
92+
@Override
93+
public IndexOutput createIndex(ReferenceGenome referenceGenome, File outputDir) throws PipelineJobException
94+
{
95+
return super.createIndex(referenceGenome, outputDir);
96+
}
97+
98+
public static class VelocytoWrapper extends AbstractCommandWrapper
99+
{
100+
public VelocytoWrapper(Logger log)
101+
{
102+
super(log);
103+
}
104+
105+
public File runVelocytoFor10x(File localBam, File gtf, File outputFolder, @Nullable File mask) throws PipelineJobException
106+
{
107+
// https://velocyto.org/velocyto.py/tutorial/cli.html#run10x-run-on-10x-chromium-samples
108+
// velocyto run10x -m repeat_msk.gtf mypath/sample01 somepath/refdata-cellranger-mm10-1.2.0/genes/genes.gtf
109+
110+
SimpleScriptWrapper wrapper = new SimpleScriptWrapper(getLogger());
111+
List<String> args = new ArrayList<>();
112+
args.add(SequencePipelineService.get().getExeForPackage("VELOCYTOPATH", "velocyto").getPath());
113+
args.add("run10x");
114+
115+
args.add("-o");
116+
args.add(outputFolder.getPath());
117+
118+
Integer threads = SequencePipelineService.get().getMaxThreads(getLogger());
119+
if (threads != null && threads > 1)
120+
{
121+
args.add("--samtools-threads");
122+
args.add(String.valueOf(threads - 1));
123+
}
124+
125+
if (mask != null)
126+
{
127+
args.add("-m");
128+
args.add(mask.getPath());
129+
}
130+
131+
// Input 10x. This is the top-level project output
132+
args.add(localBam.getParentFile().getParentFile().getPath());
133+
134+
args.add(gtf.getPath());
135+
136+
wrapper.execute(args);
137+
138+
File loom = new File(outputFolder, "sample.loom");
139+
if (!loom.exists())
140+
{
141+
throw new PipelineJobException("Missing expected file: " + loom.getPath());
142+
}
143+
144+
return loom;
145+
}
146+
}
147+
}

0 commit comments

Comments
 (0)