From 62274aa2cdac180c67a1b2ed3d908697d53df303 Mon Sep 17 00:00:00 2001 From: bbimber Date: Mon, 27 Jan 2020 22:41:13 -0800 Subject: [PATCH 01/25] Improve logging --- .../tcrdb/pipeline/CellRangerVDJUtils.java | 18 ++++++++++++------ 1 file changed, 12 insertions(+), 6 deletions(-) diff --git a/tcrdb/src/org/labkey/tcrdb/pipeline/CellRangerVDJUtils.java b/tcrdb/src/org/labkey/tcrdb/pipeline/CellRangerVDJUtils.java index dfc1a8d1a..b19822da7 100644 --- a/tcrdb/src/org/labkey/tcrdb/pipeline/CellRangerVDJUtils.java +++ b/tcrdb/src/org/labkey/tcrdb/pipeline/CellRangerVDJUtils.java @@ -541,7 +541,7 @@ else if ("Negative".equals(hto)) } else { - _log.info("skipping cell barcode without HTO call: " + barcode); + //_log.info("skipping cell barcode without HTO call: " + barcode); totalSkipped++; } continue; @@ -591,7 +591,7 @@ else if ("Negative".equals(hto)) _log.info("total rows not cells: " + nonCell); _log.info("total rows marked as cells: " + totalCells); _log.info("total clonotype rows without CDR3: " + noCDR3); - _log.info("total clonotype rows skipped for unknown barcodes: " + totalSkipped + " (" + (NumberFormat.getPercentInstance().format(totalSkipped / totalCells)) + ")"); + _log.info("total clonotype rows skipped for unknown barcodes: " + totalSkipped + " (" + (NumberFormat.getPercentInstance().format(totalSkipped / (double)totalCells)) + "%)"); _log.info("total clonotype rows skipped because they are doublets: " + doubletSkipped); _log.info("unique known cell barcodes: " + knownBarcodes.size()); _log.info("total clonotypes: " + countMapBySample.size()); @@ -611,8 +611,9 @@ else if ("Negative".equals(hto)) try (CSVReader reader = new CSVReader(Readers.getReader(consensusCsv), ',')) { String[] line; + Set uniqueClones = new HashSet<>(); Set clonesInspected = new HashSet<>(); - int clonesWithoutCounts = 0; + Set clonesWithoutCounts = new HashSet<>(); int idx = 0; while ((line = reader.readNext()) != null) { @@ -624,11 +625,16 @@ else if ("Negative".equals(hto)) } String cloneId = line[0]; + uniqueClones.add(cloneId); + Map countData = countMapBySample.get(cloneId); if (countData == null) { - _log.warn("No count data for clone: " + cloneId); - clonesWithoutCounts++; + if (!clonesWithoutCounts.contains(cloneId)) + { + _log.warn("No count data for clone: " + cloneId); + clonesWithoutCounts.add(cloneId); + } continue; } @@ -647,7 +653,7 @@ else if ("Negative".equals(hto)) totalCells += processRow(countData, cDNAMap, model, runId, am, totalCellsMapBySample, sequenceMap, rows); } - _log.info("total clones without count data: " + clonesWithoutCounts); + _log.info("total clones without count data: " + clonesWithoutCounts.size() + " (" + (NumberFormat.getPercentInstance().format(clonesWithoutCounts.size() / (double)uniqueClones.size())) + "%)"); } catch (Exception e) { From 328583871044fc66bdc8240e2d98a86beada5a33 Mon Sep 17 00:00:00 2001 From: bbimber Date: Tue, 28 Jan 2020 22:05:02 -0800 Subject: [PATCH 02/25] Preserve HTO sort order --- tcrdb/src/org/labkey/tcrdb/pipeline/CellRangerVDJUtils.java | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/tcrdb/src/org/labkey/tcrdb/pipeline/CellRangerVDJUtils.java b/tcrdb/src/org/labkey/tcrdb/pipeline/CellRangerVDJUtils.java index b19822da7..a762c381b 100644 --- a/tcrdb/src/org/labkey/tcrdb/pipeline/CellRangerVDJUtils.java +++ b/tcrdb/src/org/labkey/tcrdb/pipeline/CellRangerVDJUtils.java @@ -83,7 +83,7 @@ public void prepareVDJHashingFilesIfNeeded(PipelineJob job, SequenceAnalysisJobS _log.debug("preparing cDNA and cell hashing files"); - SequenceAnalysisService.get().writeAllCellHashingBarcodes(_sourceDir); + SequenceAnalysisService.get().writeAllCellHashingBarcodes(_sourceDir, job.getUser(), job.getContainer()); Map colMap = QueryService.get().getColumns(cDNAs, PageFlowUtil.set( FieldKey.fromString("rowid"), @@ -900,7 +900,7 @@ public static void prepareCellHashingFiles(PipelineJob job, SequenceAnalysisJobS FieldKey.fromString("hashingReadsetId")) ); - SequenceAnalysisService.get().writeAllCellHashingBarcodes(outputDir); + SequenceAnalysisService.get().writeAllCellHashingBarcodes(outputDir, job.getUser(), job.getContainer()); CellRangerVDJUtils utils = new CellRangerVDJUtils(job.getLogger(), outputDir); File barcodeOutput = utils.getValidHashingBarcodeFile(); From 613c5529394c7f6e7df619f18181b54a831a9432 Mon Sep 17 00:00:00 2001 From: bbimber Date: Wed, 29 Jan 2020 12:40:59 -0800 Subject: [PATCH 03/25] Improve reporting for discordant HTO calls --- .../CellRangerVDJCellHashingHandler.java | 2 +- .../tcrdb/pipeline/CellRangerVDJUtils.java | 17 ++++++++++++++++- 2 files changed, 17 insertions(+), 2 deletions(-) diff --git a/tcrdb/src/org/labkey/tcrdb/pipeline/CellRangerVDJCellHashingHandler.java b/tcrdb/src/org/labkey/tcrdb/pipeline/CellRangerVDJCellHashingHandler.java index ca5250b53..0e6c08b86 100644 --- a/tcrdb/src/org/labkey/tcrdb/pipeline/CellRangerVDJCellHashingHandler.java +++ b/tcrdb/src/org/labkey/tcrdb/pipeline/CellRangerVDJCellHashingHandler.java @@ -278,7 +278,7 @@ public static void processMetrics(SequenceOutputFile so, PipelineJob job, boolea String delim = description.length() > 0 ? "\n" : ""; DecimalFormat fmt = new DecimalFormat("##.##%"); - for (String metricName : Arrays.asList("InputBarcodes", "TotalCalled", "TotalCounts", "TotalSinglet", "FractionOfInputCalled", "FractionOfInputSinglet", "FractionOfInputDoublet", "FractionCalledNotInInput", "SeuratNonNegative", "MultiSeqNonNegative", "UniqueHtos", "UnknownHtoMatchingKnown")) + for (String metricName : Arrays.asList("InputBarcodes", "TotalCalled", "TotalCounts", "TotalSinglet", "FractionOfInputCalled", "FractionOfInputSinglet", "FractionOfInputDoublet", "FractionOfInputDiscordant", "FractionCalledNotInInput", "SeuratNonNegative", "MultiSeqNonNegative", "UniqueHtos", "UnknownHtoMatchingKnown")) { if (valueMap.get(metricName) != null) { diff --git a/tcrdb/src/org/labkey/tcrdb/pipeline/CellRangerVDJUtils.java b/tcrdb/src/org/labkey/tcrdb/pipeline/CellRangerVDJUtils.java index a762c381b..522a1fda2 100644 --- a/tcrdb/src/org/labkey/tcrdb/pipeline/CellRangerVDJUtils.java +++ b/tcrdb/src/org/labkey/tcrdb/pipeline/CellRangerVDJUtils.java @@ -433,6 +433,7 @@ public void importAssayData(PipelineJob job, AnalysisModel model, File outDir, I Map cellBarcodeToCDNAMap = new HashMap<>(); Set doubletBarcodes = new HashSet<>(); + Set discordantBarcodes = new HashSet<>(); if (useCellHashing) { File cellbarcodeToHtoFile = getCellToHtoFile(run); @@ -446,6 +447,7 @@ public void importAssayData(PipelineJob job, AnalysisModel model, File outDir, I //cellbarcode -> HTO name String[] line; int doublet = 0; + int discordant = 0; int negative = 0; while ((line = reader.readNext()) != null) { @@ -462,6 +464,12 @@ public void importAssayData(PipelineJob job, AnalysisModel model, File outDir, I doubletBarcodes.add(line[0]); continue; } + else if ("Discordant".equals(hto)) + { + discordant++; + discordantBarcodes.add(line[0]); + continue; + } else if ("Negative".equals(hto)) { negative++; @@ -479,6 +487,7 @@ else if ("Negative".equals(hto)) } _log.info("total doublets: " + doublet); + _log.info("total discordant: " + discordant); _log.info("total negatives: " + negative); } catch (IOException e) @@ -507,6 +516,7 @@ else if ("Negative".equals(hto)) int nonCell = 0; int totalSkipped = 0; int doubletSkipped = 0; + int discordantSkipped= 0; int hasCDR3NoClonotype = 0; Set knownBarcodes = new HashSet<>(); while ((line = reader.readNext()) != null) @@ -539,6 +549,10 @@ else if ("Negative".equals(hto)) { doubletSkipped++; } + else if (discordantBarcodes.contains(barcode)) + { + discordantSkipped++; + } else { //_log.info("skipping cell barcode without HTO call: " + barcode); @@ -592,7 +606,8 @@ else if ("Negative".equals(hto)) _log.info("total rows marked as cells: " + totalCells); _log.info("total clonotype rows without CDR3: " + noCDR3); _log.info("total clonotype rows skipped for unknown barcodes: " + totalSkipped + " (" + (NumberFormat.getPercentInstance().format(totalSkipped / (double)totalCells)) + "%)"); - _log.info("total clonotype rows skipped because they are doublets: " + doubletSkipped); + _log.info("total clonotype rows skipped because they are doublets: " + doubletSkipped + " (" + (NumberFormat.getPercentInstance().format(doubletSkipped / (double)totalCells)) + "%)"); + _log.info("total clonotype rows skipped because they are discordant calls: " + discordantSkipped + " (" + (NumberFormat.getPercentInstance().format(discordantSkipped / (double)totalCells)) + "%)"); _log.info("unique known cell barcodes: " + knownBarcodes.size()); _log.info("total clonotypes: " + countMapBySample.size()); _log.info("total cells with CDR3, lacking clonotype: " + hasCDR3NoClonotype); From 43c1ad14a5094c0c49ef6bf072cf207ef8e26258 Mon Sep 17 00:00:00 2001 From: bbimber Date: Fri, 31 Jan 2020 13:57:16 -0800 Subject: [PATCH 04/25] Support multiple types of scatter/gather (#3) * Support multiple types of scatter/gather * Reconcile AnnotationStep --- .../labkey/mgap/pipeline/AnnotationStep.java | 41 +++++++++++-------- .../RemoveAnnotationsForMgapStep.java | 14 ++++--- .../pipeline/RenameSamplesForMgapStep.java | 29 ++++++++++--- 3 files changed, 56 insertions(+), 28 deletions(-) diff --git a/mGAP/src/org/labkey/mgap/pipeline/AnnotationStep.java b/mGAP/src/org/labkey/mgap/pipeline/AnnotationStep.java index e5116a564..c8d393ff0 100644 --- a/mGAP/src/org/labkey/mgap/pipeline/AnnotationStep.java +++ b/mGAP/src/org/labkey/mgap/pipeline/AnnotationStep.java @@ -121,7 +121,7 @@ public void init(PipelineJob job, SequenceAnalysisJobSupport support, List intervals) throws PipelineJobException { VariantProcessingStepOutputImpl output = new VariantProcessingStepOutputImpl(); @@ -144,7 +144,7 @@ public Output processVariants(File inputVCF, File outputDirectory, ReferenceGeno totalSubjects = reader.getFileHeader().getSampleNamesInOrder().size(); } - boolean needToSubsetToInterval = interval != null; + boolean needToSubsetToInterval = intervals != null && !intervals.isEmpty(); boolean dropGenotypes = totalSubjects > 10; boolean dropFiltered = getProvider().getParameterByName("dropFiltered").extractValue(getPipelineCtx().getJob(), getProvider(), getStepIdx(), Boolean.class); @@ -170,8 +170,11 @@ public Output processVariants(File inputVCF, File outputDirectory, ReferenceGeno if (needToSubsetToInterval) { - selectArgs.add("-L"); - selectArgs.add(interval.getContig() + ":" + interval.getStart() + "-" + interval.getEnd()); + for (Interval interval : intervals) + { + selectArgs.add("-L"); + selectArgs.add(interval.getContig() + ":" + interval.getStart() + "-" + interval.getEnd()); + } needToSubsetToInterval = false; } @@ -202,26 +205,29 @@ public Output processVariants(File inputVCF, File outputDirectory, ReferenceGeno { List selectArgs = new ArrayList<>(); getPipelineCtx().getLogger().info("subsetting VCF by interval"); - selectArgs.add("-L"); - selectArgs.add(interval.getContig() + ":" + interval.getStart() + "-" + interval.getEnd()); + for (Interval interval : intervals) + { + selectArgs.add("-L"); + selectArgs.add(interval.getContig() + ":" + interval.getStart() + "-" + interval.getEnd()); + } needToSubsetToInterval = false; - File subset = new File(outputDirectory, SequenceAnalysisService.get().getUnzippedBaseName(inputVCF.getName()) + "." + interval.getContig() + ".subset.vcf.gz"); - if (!indexExists(subset)) + File intervalSubset = new File(outputDirectory, SequenceAnalysisService.get().getUnzippedBaseName(inputVCF.getName()) + ".intervalSubset.vcf.gz"); + if (!indexExists(intervalSubset)) { SelectVariantsWrapper wrapper = new SelectVariantsWrapper(getPipelineCtx().getLogger()); - wrapper.execute(originalGenome.getWorkingFastaFile(), inputVCF, subset, selectArgs); + wrapper.execute(originalGenome.getWorkingFastaFile(), inputVCF, intervalSubset, selectArgs); } else { - getPipelineCtx().getLogger().info("resuming with existing file: " + subset.getPath()); + getPipelineCtx().getLogger().info("resuming with existing file: " + intervalSubset.getPath()); } - output.addOutput(subset, "VCF Subset"); - output.addIntermediateFile(subset); - output.addIntermediateFile(new File(subset.getPath() + ".tbi")); + output.addOutput(intervalSubset, "VCF Subset"); + output.addIntermediateFile(intervalSubset); + output.addIntermediateFile(new File(intervalSubset.getPath() + ".tbi")); - currentVcf = subset; + currentVcf = intervalSubset; getPipelineCtx().getJob().getLogger().info("total variants: " + SequenceAnalysisService.get().getVCFLineCount(currentVcf, getPipelineCtx().getJob().getLogger(), false)); getPipelineCtx().getJob().getLogger().info("passing variants: " + SequenceAnalysisService.get().getVCFLineCount(currentVcf, getPipelineCtx().getJob().getLogger(), true)); @@ -333,8 +339,11 @@ public Output processVariants(File inputVCF, File outputDirectory, ReferenceGeno List options = new ArrayList<>(); if (needToSubsetToInterval) { - options.add("-L"); - options.add(interval.getContig() + ":" + interval.getStart() + "-" + interval.getEnd()); + for (Interval interval : intervals) + { + options.add("-L"); + options.add(interval.getContig() + ":" + interval.getStart() + "-" + interval.getEnd()); + } needToSubsetToInterval = false; } diff --git a/mGAP/src/org/labkey/mgap/pipeline/RemoveAnnotationsForMgapStep.java b/mGAP/src/org/labkey/mgap/pipeline/RemoveAnnotationsForMgapStep.java index 9c14f88c8..50dbe8666 100644 --- a/mGAP/src/org/labkey/mgap/pipeline/RemoveAnnotationsForMgapStep.java +++ b/mGAP/src/org/labkey/mgap/pipeline/RemoveAnnotationsForMgapStep.java @@ -47,7 +47,7 @@ public PipelineStep create(PipelineContext context) } @Override - public Output processVariants(File inputVCF, File outputDirectory, ReferenceGenome genome, @Nullable Interval interval) throws PipelineJobException + public Output processVariants(File inputVCF, File outputDirectory, ReferenceGenome genome, @Nullable List intervals) throws PipelineJobException { VariantProcessingStepOutputImpl output = new VariantProcessingStepOutputImpl(); @@ -58,7 +58,7 @@ public Output processVariants(File inputVCF, File outputDirectory, ReferenceGeno } else { - getWrapper().execute(inputVCF, outputFile, genome.getWorkingFastaFile(), interval); + getWrapper().execute(inputVCF, outputFile, genome.getWorkingFastaFile(), intervals); } output.setVcf(outputFile); @@ -80,7 +80,7 @@ public RemoveAnnotationsWrapper(Logger log) super(log); } - public void execute(File input, File outputFile, File referenceFasta, @Nullable Interval interval) throws PipelineJobException + public void execute(File input, File outputFile, File referenceFasta, @Nullable List intervals) throws PipelineJobException { List args = new ArrayList<>(getBaseArgs()); args.add("RemoveAnnotations"); @@ -91,10 +91,12 @@ public void execute(File input, File outputFile, File referenceFasta, @Nullable args.add("-O"); args.add(outputFile.getPath()); - if (interval != null) + if (intervals != null) { - args.add("-L"); - args.add(interval.getContig() + ":" + interval.getStart() + "-" + interval.getEnd()); + intervals.forEach(interval -> { + args.add("-L"); + args.add(interval.getContig() + ":" + interval.getStart() + "-" + interval.getEnd()); + }); } for (String key : ALLOWABLE_ANNOTATIONS) diff --git a/mGAP/src/org/labkey/mgap/pipeline/RenameSamplesForMgapStep.java b/mGAP/src/org/labkey/mgap/pipeline/RenameSamplesForMgapStep.java index 53acf0d4d..d333bcb15 100644 --- a/mGAP/src/org/labkey/mgap/pipeline/RenameSamplesForMgapStep.java +++ b/mGAP/src/org/labkey/mgap/pipeline/RenameSamplesForMgapStep.java @@ -83,11 +83,11 @@ public void init(PipelineJob job, SequenceAnalysisJobSupport support, List intervals) throws PipelineJobException { VariantProcessingStepOutputImpl output = new VariantProcessingStepOutputImpl(); - File outputFile = renameSamples(inputVCF, genome, interval); + File outputFile = renameSamples(inputVCF, genome, intervals); output.setVcf(outputFile); output.addIntermediateFile(outputFile); @@ -106,7 +106,7 @@ private File getSampleNameFile(File outputDir) return new File(outputDir, "sampleMapping.txt"); } - private File renameSamples(File currentVCF, ReferenceGenome genome, @Nullable Interval interval) throws PipelineJobException + private File renameSamples(File currentVCF, ReferenceGenome genome, @Nullable List intervals) throws PipelineJobException { getPipelineCtx().getLogger().info("renaming samples in VCF"); @@ -154,11 +154,28 @@ else if (!allSamples.contains(sample)) } writer.writeHeader(new VCFHeader(header.getMetaDataInInputOrder(), remappedSamples)); - try (CloseableIterator it = (interval == null ? reader.iterator() : reader.query(interval.getContig(), interval.getStart(), interval.getEnd()))) + if (intervals == null) { - while (it.hasNext()) + try (CloseableIterator it = reader.iterator()) { - writer.add(it.next()); + while (it.hasNext()) + { + writer.add(it.next()); + } + } + + } + else + { + for (Interval interval : intervals) + { + try (CloseableIterator it = reader.query(interval.getContig(), interval.getStart(), interval.getEnd())) + { + while (it.hasNext()) + { + writer.add(it.next()); + } + } } } } From ab6c2f39a790a08e83889ec3ca273a702e7ad799 Mon Sep 17 00:00:00 2001 From: bbimber Date: Fri, 31 Jan 2020 14:02:49 -0800 Subject: [PATCH 05/25] Prepare to refactor Seurat/Hashing --- .../pipeline/CellRangerVDJCellHashingHandler.java | 14 ++++++-------- .../labkey/tcrdb/pipeline/CellRangerVDJUtils.java | 7 ++++++- .../tcrdb/pipeline/SeuratCellHashingHandler.java | 2 +- 3 files changed, 13 insertions(+), 10 deletions(-) diff --git a/tcrdb/src/org/labkey/tcrdb/pipeline/CellRangerVDJCellHashingHandler.java b/tcrdb/src/org/labkey/tcrdb/pipeline/CellRangerVDJCellHashingHandler.java index 0e6c08b86..c6ad1a7f6 100644 --- a/tcrdb/src/org/labkey/tcrdb/pipeline/CellRangerVDJCellHashingHandler.java +++ b/tcrdb/src/org/labkey/tcrdb/pipeline/CellRangerVDJCellHashingHandler.java @@ -198,16 +198,14 @@ private void processVloupeFile(JobContext ctx, File perCellTsv, Readset rs, Reco int editDistance = ctx.getParams().optInt("editDistance", 2); File cellToHto = utils.runRemoteCellHashingTasks(output, CATEGORY, perCellTsv, rs, ctx.getSequenceSupport(), extraParams, ctx.getWorkingDirectory(), ctx.getSourceDirectory(), editDistance, scanEditDistances, genomeId); - ctx.getFileManager().addStepOutputs(action, output); - - boolean useCellHashing = utils.useCellHashing(ctx.getSequenceSupport()); - if (useCellHashing) + if (utils.useCellHashing(ctx.getSequenceSupport()) && cellToHto == null) { - if (cellToHto == null) - { - throw new PipelineJobException("Missing cell to HTO file"); - } + throw new PipelineJobException("Missing cell to HTO file"); + } + + ctx.getFileManager().addStepOutputs(action, output); + } } diff --git a/tcrdb/src/org/labkey/tcrdb/pipeline/CellRangerVDJUtils.java b/tcrdb/src/org/labkey/tcrdb/pipeline/CellRangerVDJUtils.java index 522a1fda2..50e13edeb 100644 --- a/tcrdb/src/org/labkey/tcrdb/pipeline/CellRangerVDJUtils.java +++ b/tcrdb/src/org/labkey/tcrdb/pipeline/CellRangerVDJUtils.java @@ -216,6 +216,11 @@ public File runRemoteCellHashingTasks(PipelineStepOutput output, String outputCa _log.info("No cached hashing readsets, skipping"); return null; } + else if (readsetToHashing.size() == 1) + { + _log.info("Only a single hashing readset exists, will not use hashing"); + return null; + } _log.debug("total cached readset/HTO pairs: " + readsetToHashing.size()); @@ -906,7 +911,7 @@ public static void prepareCellHashingFiles(PipelineJob job, SequenceAnalysisJobS job.getLogger().debug("preparing cell hashing files"); Container target = job.getContainer().isWorkbook() ? job.getContainer().getParent() : job.getContainer(); UserSchema tcr = QueryService.get().getUserSchema(job.getUser(), target, TCRdbSchema.NAME); - TableInfo cDNAs = tcr.getTable(TCRdbSchema.TABLE_CDNAS); + TableInfo cDNAs = tcr.getTable(TCRdbSchema.TABLE_CDNAS, null); Map colMap = QueryService.get().getColumns(cDNAs, PageFlowUtil.set( FieldKey.fromString("rowid"), diff --git a/tcrdb/src/org/labkey/tcrdb/pipeline/SeuratCellHashingHandler.java b/tcrdb/src/org/labkey/tcrdb/pipeline/SeuratCellHashingHandler.java index 09dc0b10c..d0ce232bf 100644 --- a/tcrdb/src/org/labkey/tcrdb/pipeline/SeuratCellHashingHandler.java +++ b/tcrdb/src/org/labkey/tcrdb/pipeline/SeuratCellHashingHandler.java @@ -22,7 +22,7 @@ public class SeuratCellHashingHandler extends AbstractParameterizedOutputHandler { private FileType _fileType = new FileType(".seurat.rds", false); - private static final String CATEGORY = "Seurat Cell Hashing Calls"; + public static final String CATEGORY = "Seurat Cell Hashing Calls"; public SeuratCellHashingHandler() { From 32d6857d66c72e7a6200f2fa9ed1cf5288adc11b Mon Sep 17 00:00:00 2001 From: bbimber Date: Sun, 2 Feb 2020 16:02:19 -0800 Subject: [PATCH 06/25] Allow seurat pipeline to automatically call and store HTO calls --- tcrdb/src/org/labkey/tcrdb/TCRdbModule.java | 2 + .../CellRangerCellHashingHandler.java | 19 +- .../pipeline/CellRangerSeuratHandler.java | 595 ++++++++++++++++++ .../CellRangerVDJCellHashingHandler.java | 2 +- .../tcrdb/pipeline/CellRangerVDJUtils.java | 108 +--- .../tcrdb/pipeline/CellRangerVDJWrapper.java | 2 +- .../pipeline/SeuratCellHashingHandler.java | 2 +- 7 files changed, 625 insertions(+), 105 deletions(-) create mode 100644 tcrdb/src/org/labkey/tcrdb/pipeline/CellRangerSeuratHandler.java diff --git a/tcrdb/src/org/labkey/tcrdb/TCRdbModule.java b/tcrdb/src/org/labkey/tcrdb/TCRdbModule.java index c88761935..2062f7bd0 100644 --- a/tcrdb/src/org/labkey/tcrdb/TCRdbModule.java +++ b/tcrdb/src/org/labkey/tcrdb/TCRdbModule.java @@ -27,6 +27,7 @@ import org.labkey.api.sequenceanalysis.SequenceAnalysisService; import org.labkey.api.sequenceanalysis.pipeline.SequencePipelineService; import org.labkey.tcrdb.pipeline.CellRangerCellHashingHandler; +import org.labkey.tcrdb.pipeline.CellRangerSeuratHandler; import org.labkey.tcrdb.pipeline.CellRangerVDJCellHashingHandler; import org.labkey.tcrdb.pipeline.CellRangerVDJWrapper; import org.labkey.tcrdb.pipeline.MiXCRAnalysis; @@ -120,6 +121,7 @@ public PipelineStartup() SequenceAnalysisService.get().registerFileHandler(new CellRangerCellHashingHandler()); SequenceAnalysisService.get().registerFileHandler(new CellRangerVDJCellHashingHandler()); SequenceAnalysisService.get().registerFileHandler(new SeuratCellHashingHandler()); + SequenceAnalysisService.get().registerFileHandler(new CellRangerSeuratHandler()); _hasRegistered = true; } diff --git a/tcrdb/src/org/labkey/tcrdb/pipeline/CellRangerCellHashingHandler.java b/tcrdb/src/org/labkey/tcrdb/pipeline/CellRangerCellHashingHandler.java index c07f6584a..c63db9210 100644 --- a/tcrdb/src/org/labkey/tcrdb/pipeline/CellRangerCellHashingHandler.java +++ b/tcrdb/src/org/labkey/tcrdb/pipeline/CellRangerCellHashingHandler.java @@ -94,7 +94,7 @@ public class Processor implements SequenceOutputHandler.SequenceOutputProcessor @Override public void init(PipelineJob job, SequenceAnalysisJobSupport support, List inputFiles, JSONObject params, File outputDir, List actions, List outputsToCreate) throws UnsupportedOperationException, PipelineJobException { - CellRangerVDJUtils.prepareCellHashingFiles(job, support, outputDir, "readsetId", true); + new CellRangerVDJUtils(job.getLogger(), outputDir).prepareHashingFilesIfNeeded(job, support, "readsetId"); } @Override @@ -197,6 +197,12 @@ public void complete(PipelineJob job, List inputs, List commandArgs, boolean writeLoupe, String category) throws PipelineJobException + { + CellRangerVDJUtils utils = new CellRangerVDJUtils(ctx.getLogger(), ctx.getSourceDirectory()); + return processBarcodeFile(ctx, perCellTsv, rs, genomeId, action, commandArgs, writeLoupe, category, true, utils.getValidHashingBarcodeFile()); + } + + public static File processBarcodeFile(SequenceOutputHandler.JobContext ctx, File perCellTsv, Readset rs, int genomeId, RecordedAction action, List commandArgs, boolean writeLoupe, String category, boolean createOutputFiles, File htoBarcodeWhitelist) throws PipelineJobException { ctx.getLogger().debug("inspecting file: " + perCellTsv.getPath()); @@ -251,7 +257,6 @@ public static File processBarcodeFile(SequenceOutputHandler.JobContext ctx, File } //prepare whitelist of barcodes, based on cDNA records - File htoBarcodeWhitelist = utils.getValidHashingBarcodeFile(); if (!htoBarcodeWhitelist.exists()) { throw new PipelineJobException("Unable to find file: " + htoBarcodeWhitelist.getPath()); @@ -304,7 +309,15 @@ public static File processBarcodeFile(SequenceOutputHandler.JobContext ctx, File { throw new PipelineJobException(e); } - ctx.getFileManager().addSequenceOutput(forLoupe, rs.getName() + ": Cell Hashing Calls", "10x GEX Cell Hashing Calls (Loupe)", rs.getReadsetId(), null, genomeId, null); + + if (createOutputFiles) + { + ctx.getFileManager().addSequenceOutput(forLoupe, rs.getName() + ": Cell Hashing Calls", "10x GEX Cell Hashing Calls (Loupe)", rs.getReadsetId(), null, genomeId, null); + } + else + { + ctx.getLogger().debug("Output file creation will be skipped"); + } } return cellToHto; diff --git a/tcrdb/src/org/labkey/tcrdb/pipeline/CellRangerSeuratHandler.java b/tcrdb/src/org/labkey/tcrdb/pipeline/CellRangerSeuratHandler.java new file mode 100644 index 000000000..602180f77 --- /dev/null +++ b/tcrdb/src/org/labkey/tcrdb/pipeline/CellRangerSeuratHandler.java @@ -0,0 +1,595 @@ +package org.labkey.tcrdb.pipeline; + +import au.com.bytecode.opencsv.CSVReader; +import au.com.bytecode.opencsv.CSVWriter; +import htsjdk.samtools.util.IOUtil; +import org.apache.commons.io.FileUtils; +import org.apache.commons.lang3.StringUtils; +import org.json.JSONObject; +import org.labkey.api.data.DbSchema; +import org.labkey.api.data.DbSchemaType; +import org.labkey.api.data.Table; +import org.labkey.api.data.TableInfo; +import org.labkey.api.module.ModuleLoader; +import org.labkey.api.pipeline.PipelineJob; +import org.labkey.api.pipeline.PipelineJobException; +import org.labkey.api.pipeline.RecordedAction; +import org.labkey.api.reader.Readers; +import org.labkey.api.sequenceanalysis.SequenceAnalysisService; +import org.labkey.api.sequenceanalysis.SequenceOutputFile; +import org.labkey.api.sequenceanalysis.model.Readset; +import org.labkey.api.sequenceanalysis.pipeline.AbstractParameterizedOutputHandler; +import org.labkey.api.sequenceanalysis.pipeline.SequenceAnalysisJobSupport; +import org.labkey.api.sequenceanalysis.pipeline.SequenceOutputHandler; +import org.labkey.api.sequenceanalysis.pipeline.SequencePipelineService; +import org.labkey.api.sequenceanalysis.pipeline.ToolParameterDescriptor; +import org.labkey.api.sequenceanalysis.run.SimpleScriptWrapper; +import org.labkey.api.util.FileType; +import org.labkey.api.util.FileUtil; +import org.labkey.api.writer.PrintWriters; +import org.labkey.tcrdb.TCRdbModule; + +import java.io.File; +import java.io.IOException; +import java.io.PrintWriter; +import java.util.Arrays; +import java.util.Collections; +import java.util.HashMap; +import java.util.HashSet; +import java.util.LinkedHashSet; +import java.util.List; +import java.util.Map; +import java.util.Set; + +public class CellRangerSeuratHandler extends AbstractParameterizedOutputHandler +{ + private FileType _fileType = new FileType("cloupe", false); + + public CellRangerSeuratHandler() + { + super(ModuleLoader.getInstance().getModule(TCRdbModule.class), "Run Seurat", "This will run a standard seurat-based pipeline on the selected 10x/cellranger data and save the resulting Seurat object as an rds file for external use.", new LinkedHashSet<>(), Arrays.asList( + ToolParameterDescriptor.create("projectName", "Output Name", "This will be used as the final sample/file name. If blank, the readset name will be used. The latter cannot be used when merging multiple inputs.", "textfield", new JSONObject(){{ + + }}, null), + ToolParameterDescriptor.create("doSplitJobs", "Run Separately", "If checked, each input dataset will be run separately. Otherwise they will be merged", "checkbox", new JSONObject(){{ + put("checked", true); + }}, false), + ToolParameterDescriptor.create("useOutputFileContainer", "Submit to Source File Workbook", "If checked, each job will be submitted to the same workbook as the input file, as opposed to submitting all jobs to the same workbook. This is primarily useful if submitting a large batch of files to process separately. This only applies if 'Run Separately' is selected.", "checkbox", new JSONObject(){{ + put("checked", false); + }}, false), + ToolParameterDescriptor.create("dimsToUse", "PCs To Use", "If non-blank, this is the number of PCs that seurat will use for dim reduction steps.", "ldk-integerfield", new JSONObject(){{ + + }}, null), + ToolParameterDescriptor.create("minDimsToUse", "Minimum PCs To Use", "If non-blank, the pipeline will attempt to infer the number of PCs to use for dim reduction, but will not use fewer than this value.", "ldk-integerfield", new JSONObject(){{ + + }}, 10), + ToolParameterDescriptor.create("doCellFilter", "Perform Cell Filtering", "If selected, cells will be filtered on pct.mito and number of unique genes.", "checkbox", new JSONObject(){{ + put("checked", true); + }}, true), + ToolParameterDescriptor.create("doCellCycle", "Perform Cell Cycle Correction", "If selected, the pipeline will attempt to correct for cell cycle.", "checkbox", new JSONObject(){{ + put("checked", true); + }}, true), + ToolParameterDescriptor.create("useSCTransform", "Use SCTransform", "If selected, the pipeline will use the newer SCtransform method instead of the standard Seurat pipeline.", "checkbox", new JSONObject(){{ + put("checked", false); + }}, false), + ToolParameterDescriptor.create("runSingleR", "Run SingleR", "If selected, SingleR will be run after Seurat processing.", "checkbox", new JSONObject(){{ + put("checked", true); + }}, true), + ToolParameterDescriptor.create("mergeMethod", "Merge Method", "This determines whether any batch correction will be applied when merging datasets.", "ldk-simplecombo", new JSONObject(){{ + put("storeValues", "simple;cca"); + }}, "simple"), + ToolParameterDescriptor.create("scanEditDistances", "Scan Edit Distances (Hashing)", "If checked, CITE-seq-count will be run using edit distances from 0-3 and the iteration with the highest singlets will be used.", "checkbox", new JSONObject(){{ + put("checked", true); + }}, true), + ToolParameterDescriptor.create("editDistance", "Edit Distance (Hashing)", null, "ldk-integerfield", null, 1) + )); + } + + @Override + public boolean canProcess(SequenceOutputFile o) + { + return o.getFile() != null && _fileType.isType(o.getFile()); + } + + @Override + public List validateParameters(List outputFiles, JSONObject params) + { + if (!params.optBoolean("doSplitJobs", false) && StringUtils.trimToNull(params.optString("projectName")) == null && outputFiles.size() > 1) + { + return Collections.singletonList("Must provide the output name when merging multiple inputs"); + } + + return null; + } + + @Override + public boolean doRunRemote() + { + return true; + } + + @Override + public boolean doRunLocal() + { + return false; + } + + @Override + public boolean requiresSingleGenome() + { + return false; + } + + @Override + public SequenceOutputProcessor getProcessor() + { + return new CellRangerSeuratHandler.Processor(); + } + + @Override + public boolean doSplitJobs() + { + return false; + } + + public class Processor implements SequenceOutputProcessor + { + @Override + public void init(PipelineJob job, SequenceAnalysisJobSupport support, List inputFiles, JSONObject params, File outputDir, List actions, List outputsToCreate) throws UnsupportedOperationException, PipelineJobException + { + for (SequenceOutputFile so : inputFiles) + { + if (so.getReadset() != null) + { + support.cacheReadset(so.getReadset(), job.getUser()); + } + else + { + job.getLogger().error("Output file lacks a readset and will be skipped: " + so.getRowid()); + } + } + + new CellRangerVDJUtils(job.getLogger(), outputDir).prepareHashingFilesIfNeeded(job, support,"readsetId"); + } + + @Override + public void processFilesOnWebserver(PipelineJob job, SequenceAnalysisJobSupport support, List inputFiles, JSONObject params, File outputDir, List actions, List outputsToCreate) throws UnsupportedOperationException, PipelineJobException + { + + } + + @Override + public void processFilesRemote(List inputFiles, JobContext ctx) throws UnsupportedOperationException, PipelineJobException + { + RecordedAction action = new RecordedAction(getName()); + ctx.addActions(action); + + Set rsNames = new HashSet<>(); + for (SequenceOutputFile so : inputFiles) + { + ctx.getFileManager().addInput(action, "CellRanger Loupe", so.getFile()); + if (so.getReadset() != null) + { + rsNames.add(ctx.getSequenceSupport().getCachedReadset(so.getReadset()).getName()); + } + } + + String outPrefix = StringUtils.trimToNull(ctx.getParams().getString("projectName")); + if (outPrefix == null) + { + if (rsNames.size() == 1) + { + outPrefix = rsNames.iterator().next(); + } + else + { + throw new PipelineJobException("Must provide the output prefix when merging more than one output file"); + } + } + outPrefix = FileUtil.makeLegalName(outPrefix); + + File seuratObj = new File(ctx.getWorkingDirectory(), outPrefix + ".seurat.rds"); + File doneFile = new File(seuratObj.getPath() + ".done"); + boolean seuratHasRun = doneFile.exists(); + if (seuratHasRun) + { + ctx.getLogger().info("Seurat has already run, will not repeat"); + } + + Map dataMap = new HashMap<>(); + + File pr = ctx.getFolderPipeRoot().getRootPath().getParentFile(); //drop the @files or @pipeline + for (SequenceOutputFile so : inputFiles) + { + //start with seurat 3 + File subDir = new File(so.getFile().getParentFile(), "raw_feature_bc_matrix"); + if (!subDir.exists()) + { + //try 2 + subDir = new File(so.getFile().getParentFile(), "raw_gene_bc_matrices"); + if (subDir.exists()) + { + //now infer subdir: + for (File f : subDir.listFiles()) + { + if (f.isDirectory()) + { + subDir = f; + break; + } + } + } + } + + if (!subDir.exists()) + { + throw new PipelineJobException("Unable to find raw data for input: " + so.getFile().getPath()); + } + + try + { + String subDirRel = FileUtil.relativize(pr, subDir, true); + ctx.getLogger().debug("pipe root: " + pr.getPath()); + ctx.getLogger().debug("file path: " + subDir.getPath()); + ctx.getLogger().debug("relative path: " + subDirRel); + + //Copy raw data directory locally to avoid docker permission issues + String dirName = so.getRowid() + "_RawData"; + File copyDir = new File(ctx.getWorkingDirectory(), dirName); + if (!seuratHasRun) + { + if (copyDir.exists()) + { + ctx.getLogger().debug("Deleting directory: " + copyDir.getPath()); + FileUtils.deleteDirectory(copyDir); + } + + ctx.getLogger().debug("Copying raw data directory: " + subDir.getPath()); + ctx.getLogger().debug("To: " + copyDir.getPath()); + FileUtils.copyDirectory(subDir, copyDir); + } + ctx.getFileManager().addIntermediateFile(copyDir); + + dataMap.put(so, dirName); + } + catch (IOException e) + { + throw new PipelineJobException(e); + } + } + + File rmdScript = new File(SequenceAnalysisService.get().getScriptPath(TCRdbModule.NAME, "external/scRNAseq/Seurat3.rmd")); + if (!rmdScript.exists()) + { + throw new PipelineJobException("Unable to find script: " + rmdScript.getPath()); + } + + File wrapperScript = new File(SequenceAnalysisService.get().getScriptPath(TCRdbModule.NAME, "external/scRNAseq/seuratWrapper.sh")); + if (!wrapperScript.exists()) + { + throw new PipelineJobException("Unable to find script: " + wrapperScript.getPath()); + } + + File tmpScript = new File(ctx.getWorkingDirectory(), "script.R"); + File outHtml = new File(ctx.getWorkingDirectory(), outPrefix + ".html"); + + try (PrintWriter writer = PrintWriters.getPrintWriter(tmpScript)) + { + File scriptCopy = new File(ctx.getWorkingDirectory(), rmdScript.getName()); + if (scriptCopy.exists()) + { + scriptCopy.delete(); + } + + IOUtil.copyFile(rmdScript, scriptCopy); + rmdScript = scriptCopy; + ctx.getFileManager().addIntermediateFile(rmdScript); + + scriptCopy = new File(ctx.getWorkingDirectory(), wrapperScript.getName()); + if (scriptCopy.exists()) + { + scriptCopy.delete(); + } + + IOUtil.copyFile(wrapperScript, scriptCopy); + ctx.getFileManager().addIntermediateFile(scriptCopy); + + writer.println("outPrefix <- '" + outPrefix + "'"); + writer.println("resolutionToUse <- 0.6"); + for (String v : new String[]{"dimsToUse", "minDimsToUse"}) + { + String val = StringUtils.trimToNull(ctx.getParams().optString(v)); + val = val == null ? "NULL" : val; + + writer.println(v + " <- " + val); + } + + String mergeMethod = StringUtils.trimToNull(ctx.getParams().optString("mergeMethod")); + mergeMethod = mergeMethod == null ? "NULL" : "'" + mergeMethod + "'"; + writer.println("mergeMethod <- " + mergeMethod); + + boolean doCellFilter = ctx.getParams().optBoolean("doCellFilter", true); + writer.println("doCellFilter <- " + String.valueOf(doCellFilter).toUpperCase()); + + boolean runSingleR = ctx.getParams().optBoolean("runSingleR", true); + writer.println("runSingleR <- " + String.valueOf(runSingleR).toUpperCase()); + + boolean doCellCycle = ctx.getParams().optBoolean("doCellCycle", true); + writer.println("doCellCycle <- " + String.valueOf(doCellCycle).toUpperCase()); + + boolean useSCTransform = ctx.getParams().optBoolean("useSCTransform", false); + writer.println("useSCTransform <- " + String.valueOf(useSCTransform).toUpperCase()); + + writer.println("data <- list("); + String delim = ""; + for (SequenceOutputFile so : dataMap.keySet()) + { + writer.println("\t" + delim + "'" + so.getRowid() + "'='" + dataMap.get(so) + "'"); + delim = ","; + } + writer.println(")"); + writer.println(); + writer.println(); + writer.println("setwd('/work')"); + + writer.println("rmarkdown::render('" + rmdScript.getName() + "', clean=TRUE, output_file='" + outHtml.getName() + "')"); + } + catch (IOException e) + { + throw new PipelineJobException(e); + } + + if (!seuratHasRun) + { + SimpleScriptWrapper wrapper = new SimpleScriptWrapper(ctx.getLogger()); + wrapper.setWorkingDir(ctx.getWorkingDirectory()); + wrapper.execute(Arrays.asList("/bin/bash", wrapperScript.getName(), pr.getPath())); + + try + { + FileUtils.touch(doneFile); + ctx.getFileManager().addIntermediateFile(doneFile); + } + catch (IOException e) + { + throw new PipelineJobException(e); + } + } + + if (!seuratObj.exists()) + { + throw new PipelineJobException("Unable to find expected file: " + seuratObj.getPath()); + } + + String dimsToUse = StringUtils.trimToNull(ctx.getParams().optString("dimsToUse")); + String minDimsToUse = StringUtils.trimToNull(ctx.getParams().optString("minDimsToUse")); + String mergeMethod = StringUtils.trimToNull(ctx.getParams().optString("mergeMethod")); + + String description = StringUtils.join(new String[]{ + "Correct Cell Cycle: " + ctx.getParams().optBoolean("doCellCycle", true), + "Perform Cell Filtering: " + ctx.getParams().optBoolean("doCellFilter", true), + "Min. Dims To Use: " + (minDimsToUse == null ? "NA" : minDimsToUse), + "Dims To Use: " + (dimsToUse == null ? "automatic" : dimsToUse), + "Use SCTransform: " + ctx.getParams().optBoolean("useSCTransform", false), + "Merge method: " + mergeMethod + }, "\n"); + + ctx.getFileManager().addSequenceOutput(seuratObj, "Seurat Object: " + outPrefix, "Seurat Data", (inputFiles.size() == 1 ? inputFiles.iterator().next().getReadset() : null), null, getGenomeId(inputFiles), description); + ctx.getFileManager().addOutput(action, "Seurat Object", seuratObj); + + if (!outHtml.exists()) + { + throw new PipelineJobException("Unable to find summary report"); + } + ctx.getFileManager().addOutput(action, "Seurat Report", outHtml); + ctx.getFileManager().addSequenceOutput(outHtml, "Seurat Report: " + outPrefix, "Seurat Report", (inputFiles.size() == 1 ? inputFiles.iterator().next().getReadset() : null), null, getGenomeId(inputFiles), description); + + File seuratObjRaw = new File(ctx.getWorkingDirectory(), outPrefix + ".rawData.rds"); + if (seuratObjRaw.exists()) + { + ctx.getFileManager().addIntermediateFile(seuratObjRaw); + } + + CellRangerVDJUtils utils = new CellRangerVDJUtils(ctx.getLogger(), ctx.getSourceDirectory()); + if (utils.useCellHashing(ctx.getSequenceSupport())) + { + File allCellBarcodes = new File(seuratObj.getParentFile(), seuratObj.getName().replaceAll("seurat.rds", "cellBarcodes.csv")); + Map finalCalls = new HashMap<>(); + if (!allCellBarcodes.exists()) + { + throw new PipelineJobException("Unable to find expected cell barcodes file. This might indicate the seurat object was created with an older version of the pipeline. Expected: " + allCellBarcodes.getPath()); + } + + for (SequenceOutputFile so : inputFiles) + { + String barcodePrefix = so.getRowid().toString(); + Readset rs = ctx.getSequenceSupport().getCachedReadset(so.getReadset()); + if (rs == null) + { + throw new PipelineJobException("Unable to find readset for outputfile: " + so.getRowid()); + } + else if (rs.getReadsetId() == null) + { + throw new PipelineJobException("Readset lacks a rowId for outputfile: " + so.getRowid()); + } + + //Subset barcodes by dataset: + File barcodes = new File(allCellBarcodes.getParentFile(), "cellBarcodeWhitelist." + barcodePrefix + ".txt"); + try (CSVReader reader = new CSVReader(Readers.getReader(allCellBarcodes), '\t'); CSVWriter writer = new CSVWriter(PrintWriters.getPrintWriter(barcodes), '\t', CSVWriter.NO_QUOTE_CHARACTER)) + { + String[] line; + while ((line = reader.readNext()) != null) + { + String barcode = line[0]; + if (barcode.startsWith(barcodePrefix + "_")) + { + barcode = barcode.split("_")[1]; + writer.writeNext(new String[]{barcode}); + } + } + } + catch (IOException e) + { + throw new PipelineJobException(e); + } + + // write readset-specific HTO list + Integer hashingReadsetId = CellRangerVDJUtils.getCachedReadsetMap(ctx.getSequenceSupport()).get(rs.getReadsetId()); + if (hashingReadsetId == null) + { + throw new PipelineJobException("Unable to find hashing readset Id for: " + rs.getReadsetId()); + } + + File perReadsetHtos = new File(allCellBarcodes.getParentFile(), "allowableHtos." + barcodePrefix + ".txt"); + int htosForReadset = 0; + try (CSVReader reader = new CSVReader(Readers.getReader(utils.getCDNAInfoFile()), '\t'); CSVWriter bcWriter = new CSVWriter(PrintWriters.getPrintWriter(perReadsetHtos), ',', CSVWriter.NO_QUOTE_CHARACTER)) + { + String[] line; + while ((line = reader.readNext()) != null) + { + if (hashingReadsetId.toString().equals(line[6])) + { + htosForReadset++; + bcWriter.writeNext(new String[]{line[5], line[4]}); + } + } + } + catch (IOException e) + { + throw new PipelineJobException(e); + } + + if (htosForReadset > 0) + { + ctx.getLogger().info("Total HTOs for readset: " + htosForReadset); + finalCalls.put(barcodePrefix, CellRangerCellHashingHandler.processBarcodeFile(ctx, barcodes, rs, so.getLibrary_id(), action, getClientCommandArgs(ctx.getParams()), false, SeuratCellHashingHandler.CATEGORY, false, perReadsetHtos)); + } + else + { + ctx.getLogger().info("No HTOs found for readset"); + } + } + + if (!finalCalls.isEmpty()) + { + ctx.getLogger().info("Storing cell hashing calls in seurat object"); + appendCallsToSeurat(ctx, seuratObj, finalCalls); + } + else + { + ctx.getLogger().info("Cell hashing was not used. will not append to seurat"); + } + } + } + + private void appendCallsToSeurat(JobContext ctx, File seuratObj, Map finalCalls) throws PipelineJobException + { + File rScript = new File(seuratObj.getParentFile(), "appendHashing.R"); + File bashScript = new File(seuratObj.getParentFile(), "runDocker.R"); + + try (PrintWriter rWriter = PrintWriters.getPrintWriter(rScript); PrintWriter bashWriter = PrintWriters.getPrintWriter(bashScript)) + { + rWriter.println("library(OOSAP)"); + rWriter.println("seuratObj <- readRDS('" + seuratObj.getName() + "')"); + rWriter.println("callsFiles <- list("); + finalCalls.forEach((x, y) -> { + rWriter.println("'" + x + "' = '" + y.getName() + "'"); + }); + + rWriter.println(")"); + rWriter.println(""); + rWriter.println("for (barcodePrefix in names(callsFiles)) {"); + rWriter.println(" seuratObj <- OOSAP:::AppendCellHashing(seuratObj = seuratObj, barcodeCallFile = callsFiles[[barcodePrefix]], barcodePrefix = barcodePrefix)"); + rWriter.println("}"); + rWriter.println("saveRDS(seuratObj, file = '" + seuratObj.getName() + "')"); + + bashWriter.println("#!/bin/bash"); + bashWriter.println("set -e"); + bashWriter.println("set -x"); + bashWriter.println("DOCKER=/opt/acc/sbin/exadocker"); + bashWriter.println("WD=`pwd`"); + bashWriter.println("HOME=`echo ~/`"); + + Integer maxRam = SequencePipelineService.get().getMaxRam(); + String ramOpts = ""; + if (maxRam != null) + { + ramOpts = " --memory=" +maxRam +"g"; + } + + bashWriter.println("sudo $DOCKER pull bimberlab/oosap"); + bashWriter.println("sudo $DOCKER run --rm=true " + ramOpts + "-v \"${WD}:/work\" -v \"${HOME}:/homeDir\" -u $UID -e USERID=$UID -w /work -e HOME=/homeDir bimberlab/oosap Rscript --vanilla " + rScript.getName()); + + SimpleScriptWrapper wrapper = new SimpleScriptWrapper(ctx.getLogger()); + wrapper.setWorkingDir(seuratObj.getParentFile()); + wrapper.execute(Arrays.asList("/bin/bash", bashScript.getName())); + } + catch (IOException e) + { + throw new PipelineJobException(e); + } + } + + @Override + public void complete(PipelineJob job, List inputs, List outputsCreated, SequenceAnalysisJobSupport support) throws PipelineJobException + { + for (SequenceOutputFile so : outputsCreated) + { + if (so.getFile() != null && so.getFile().getPath().endsWith(".seurat.rds")) + { + File metrics = new File(so.getFile().getPath().replaceAll(".seurat.rds", ".summary.txt")); + if (metrics.exists()) + { + job.getLogger().info("Loading metrics"); + TableInfo ti = DbSchema.get("sequenceanalysis", DbSchemaType.Module).getTable("quality_metrics"); + int total = 0; + try (CSVReader reader = new CSVReader(Readers.getReader(metrics), '\t')) + { + String[] line; + while ((line = reader.readNext()) != null) + { + if ("Category".equals(line[0])) + { + continue; + } + + Map r = new HashMap<>(); + r.put("category", "Seurat"); + r.put("metricname", line[1]); + r.put("metricvalue", line[2]); + r.put("analysis_id", so.getAnalysis_id()); + r.put("dataid", so.getDataId()); + r.put("readset", so.getReadset()); + r.put("container", job.getContainer()); + r.put("createdby", job.getUser().getUserId()); + + Table.insert(job.getUser(), ti, r); + total++; + } + } + catch (IOException e) + { + throw new PipelineJobException(e); + } + + job.getLogger().info("total metrics: " + total); + } + else + { + job.getLogger().warn("Unable to find metrics file: " + metrics.getPath()); + } + } + } + } + } + + private Integer getGenomeId(List inputFiles) + { + Set genomeIds = new HashSet<>(); + inputFiles.forEach(x -> { + genomeIds.add(x.getLibrary_id()); + }); + + return genomeIds.size() == 1 ? genomeIds.iterator().next() : null; + } +} diff --git a/tcrdb/src/org/labkey/tcrdb/pipeline/CellRangerVDJCellHashingHandler.java b/tcrdb/src/org/labkey/tcrdb/pipeline/CellRangerVDJCellHashingHandler.java index c6ad1a7f6..2152e29ba 100644 --- a/tcrdb/src/org/labkey/tcrdb/pipeline/CellRangerVDJCellHashingHandler.java +++ b/tcrdb/src/org/labkey/tcrdb/pipeline/CellRangerVDJCellHashingHandler.java @@ -105,7 +105,7 @@ public class Processor implements SequenceOutputHandler.SequenceOutputProcessor public void init(PipelineJob job, SequenceAnalysisJobSupport support, List inputFiles, JSONObject params, File outputDir, List actions, List outputsToCreate) throws UnsupportedOperationException, PipelineJobException { CellRangerVDJUtils utils = new CellRangerVDJUtils(job.getLogger(), outputDir); - utils.prepareVDJHashingFilesIfNeeded(job, support); + utils.prepareHashingFilesIfNeeded(job, support, "enrichedReadsetId"); } @Override diff --git a/tcrdb/src/org/labkey/tcrdb/pipeline/CellRangerVDJUtils.java b/tcrdb/src/org/labkey/tcrdb/pipeline/CellRangerVDJUtils.java index 50e13edeb..d031ccf6b 100644 --- a/tcrdb/src/org/labkey/tcrdb/pipeline/CellRangerVDJUtils.java +++ b/tcrdb/src/org/labkey/tcrdb/pipeline/CellRangerVDJUtils.java @@ -75,11 +75,11 @@ public CellRangerVDJUtils(Logger log, File sourceDir) _sourceDir = sourceDir; } - public void prepareVDJHashingFilesIfNeeded(PipelineJob job, SequenceAnalysisJobSupport support) throws PipelineJobException + public void prepareHashingFilesIfNeeded(PipelineJob job, SequenceAnalysisJobSupport support, String filterField) throws PipelineJobException { Container target = job.getContainer().isWorkbook() ? job.getContainer().getParent() : job.getContainer(); UserSchema tcr = QueryService.get().getUserSchema(job.getUser(), target, TCRdbSchema.NAME); - TableInfo cDNAs = tcr.getTable(TCRdbSchema.TABLE_CDNAS); + TableInfo cDNAs = tcr.getTable(TCRdbSchema.TABLE_CDNAS, null); _log.debug("preparing cDNA and cell hashing files"); @@ -101,15 +101,15 @@ public void prepareVDJHashingFilesIfNeeded(PipelineJob job, SequenceAnalysisJobS HashMap readsetToHashingMap = new HashMap<>(); try (CSVWriter writer = new CSVWriter(PrintWriters.getPrintWriter(output), '\t', CSVWriter.NO_QUOTE_CHARACTER); CSVWriter bcWriter = new CSVWriter(PrintWriters.getPrintWriter(barcodeOutput), ',', CSVWriter.NO_QUOTE_CHARACTER)) { - writer.writeNext(new String[]{"TCR_ReadsetId", "CDNA_ID", "AnimalId", "Stim", "Population", "HTO_Name", "HTO_Seq", "HashingReadsetId"}); + writer.writeNext(new String[]{"ReadsetId", "CDNA_ID", "AnimalId", "Stim", "Population", "HTO_Name", "HTO_Seq", "HashingReadsetId"}); List cachedReadsets = support.getCachedReadsets(); Set distinctHTOs = new HashSet<>(); Set hashingStatus = new HashSet<>(); for (Readset rs : cachedReadsets) { AtomicBoolean hasError = new AtomicBoolean(false); - //find cDNA records using this as enrichedReadset - new TableSelector(cDNAs, colMap.values(), new SimpleFilter(FieldKey.fromString("enrichedReadsetId"), rs.getRowId()), null).forEachResults(results -> { + //find cDNA records using this readset + new TableSelector(cDNAs, colMap.values(), new SimpleFilter(FieldKey.fromString(filterField), rs.getRowId()), null).forEachResults(results -> { if (results.getObject(FieldKey.fromString("status")) != null) { _log.info("skipping cDNA with non-null status: " + results.getString(FieldKey.fromString("rowid"))); @@ -160,6 +160,10 @@ else if (useCellHashing) { throw new PipelineJobException("The selected readsets/cDNA records use a mixture of cell hashing and non-hashing."); } + else if (hashingStatus.isEmpty()) + { + throw new PipelineJobException("There were no readsets found."); + } } // if distinct HTOs is 1, no point in running hashing. note: presence of hashing readsets is a trigger downstream @@ -906,100 +910,6 @@ public boolean useCellHashing(SequenceAnalysisJobSupport support) throws Pipelin return getCachedReadsetMap(support).size() > 1; } - public static void prepareCellHashingFiles(PipelineJob job, SequenceAnalysisJobSupport support, File outputDir, String filterFieldName, boolean throwOnZeroHto) throws PipelineJobException - { - job.getLogger().debug("preparing cell hashing files"); - Container target = job.getContainer().isWorkbook() ? job.getContainer().getParent() : job.getContainer(); - UserSchema tcr = QueryService.get().getUserSchema(job.getUser(), target, TCRdbSchema.NAME); - TableInfo cDNAs = tcr.getTable(TCRdbSchema.TABLE_CDNAS, null); - - Map colMap = QueryService.get().getColumns(cDNAs, PageFlowUtil.set( - FieldKey.fromString("rowid"), - FieldKey.fromString("sortId/hto"), - FieldKey.fromString("sortId/hto/sequence"), - FieldKey.fromString("hashingReadsetId")) - ); - - SequenceAnalysisService.get().writeAllCellHashingBarcodes(outputDir, job.getUser(), job.getContainer()); - - CellRangerVDJUtils utils = new CellRangerVDJUtils(job.getLogger(), outputDir); - File barcodeOutput = utils.getValidHashingBarcodeFile(); - HashMap readsetToHashingMap = new HashMap<>(); - try (CSVWriter bcWriter = new CSVWriter(PrintWriters.getPrintWriter(barcodeOutput), ',', CSVWriter.NO_QUOTE_CHARACTER)) - { - List cachedReadsets = support.getCachedReadsets(); - job.getLogger().debug("total cached readsets: " + cachedReadsets.size()); - Set distinctHTOs = new HashSet<>(); - Set hashingStatus = new HashSet<>(); - for (Readset rs : cachedReadsets) - { - AtomicBoolean hasError = new AtomicBoolean(false); - new TableSelector(cDNAs, colMap.values(), new SimpleFilter(FieldKey.fromString(filterFieldName), rs.getRowId()), null).forEachResults(results -> { - boolean useCellHashing = results.getObject(FieldKey.fromString("sortId/hto")) != null; - hashingStatus.add(useCellHashing); - - if (!useCellHashing) - { - return; - } - - if (results.getObject(FieldKey.fromString("hashingReadsetId")) == null || results.getInt(FieldKey.fromString("hashingReadsetId")) == 0) - { - hasError.set(true); - } - - if (results.getObject(FieldKey.fromString("sortId/hto/sequence")) == null) - { - hasError.set(true); - } - - support.cacheReadset(results.getInt(FieldKey.fromString("hashingReadsetId")), job.getUser()); - readsetToHashingMap.put(rs.getReadsetId(), results.getInt(FieldKey.fromString("hashingReadsetId"))); - - String hto = results.getString(FieldKey.fromString("sortId/hto")) + "<>" + results.getString(FieldKey.fromString("sortId/hto/sequence")); - if (!distinctHTOs.contains(hto) && !StringUtils.isEmpty(results.getString(FieldKey.fromString("sortId/hto/sequence")))) - { - distinctHTOs.add(hto); - bcWriter.writeNext(new String[]{results.getString(FieldKey.fromString("sortId/hto/sequence")), results.getString(FieldKey.fromString("sortId/hto"))}); - } - }); - - if (hasError.get()) - { - throw new PipelineJobException("No cell hashing readset or HTO found for one or more cDNAs. see the file"); - } - } - - if (hashingStatus.size() > 1) - { - throw new PipelineJobException("The selected readsets/cDNA records use a mixture of cell hashing and non-hashing."); - } - else if (hashingStatus.isEmpty()) - { - throw new PipelineJobException("There were no readsets found."); - } - - boolean useCellHashing = hashingStatus.iterator().next(); - if (useCellHashing && distinctHTOs.isEmpty()) - { - throw new PipelineJobException("Cell hashing was selected, but no HTOs were found"); - } - - job.getLogger().info("distinct HTOs: " + distinctHTOs.size()); - - support.cacheObject(READSET_TO_HASHING_MAP, readsetToHashingMap); - - if (throwOnZeroHto && distinctHTOs.isEmpty()) - { - throw new PipelineJobException("None of the provided samples use cell hashing"); - } - } - catch (IOException e) - { - throw new PipelineJobException(e); - } - } - public static class CDNA { private int _rowId; diff --git a/tcrdb/src/org/labkey/tcrdb/pipeline/CellRangerVDJWrapper.java b/tcrdb/src/org/labkey/tcrdb/pipeline/CellRangerVDJWrapper.java index 7a215554d..3275a4d9b 100644 --- a/tcrdb/src/org/labkey/tcrdb/pipeline/CellRangerVDJWrapper.java +++ b/tcrdb/src/org/labkey/tcrdb/pipeline/CellRangerVDJWrapper.java @@ -216,7 +216,7 @@ else if (nt.getLineage().contains("D")) } } - getUtils().prepareVDJHashingFilesIfNeeded(getPipelineCtx().getJob(), getPipelineCtx().getSequenceSupport()); + getUtils().prepareHashingFilesIfNeeded(getPipelineCtx().getJob(), getPipelineCtx().getSequenceSupport(), "enrichedReadsetId"); } private File getGenomeFasta() diff --git a/tcrdb/src/org/labkey/tcrdb/pipeline/SeuratCellHashingHandler.java b/tcrdb/src/org/labkey/tcrdb/pipeline/SeuratCellHashingHandler.java index d0ce232bf..59eaceee2 100644 --- a/tcrdb/src/org/labkey/tcrdb/pipeline/SeuratCellHashingHandler.java +++ b/tcrdb/src/org/labkey/tcrdb/pipeline/SeuratCellHashingHandler.java @@ -79,7 +79,7 @@ public class Processor implements SequenceOutputHandler.SequenceOutputProcessor @Override public void init(PipelineJob job, SequenceAnalysisJobSupport support, List inputFiles, JSONObject params, File outputDir, List actions, List outputsToCreate) throws UnsupportedOperationException, PipelineJobException { - CellRangerVDJUtils.prepareCellHashingFiles(job, support, outputDir, "readsetId", true); + new CellRangerVDJUtils(job.getLogger(), outputDir).prepareHashingFilesIfNeeded(job, support, "readsetId"); } @Override From 75c2d7997d558e09f2aae189e1da5e1a1f7f408d Mon Sep 17 00:00:00 2001 From: bbimber Date: Sun, 2 Feb 2020 20:14:35 -0800 Subject: [PATCH 07/25] Move seurat scripts to tcrdb module --- tcrdb/resources/external/scRNAseq/Seurat3.rmd | 168 ++++++++++++++++++ .../external/scRNAseq/seuratWrapper.sh | 27 +++ 2 files changed, 195 insertions(+) create mode 100644 tcrdb/resources/external/scRNAseq/Seurat3.rmd create mode 100644 tcrdb/resources/external/scRNAseq/seuratWrapper.sh diff --git a/tcrdb/resources/external/scRNAseq/Seurat3.rmd b/tcrdb/resources/external/scRNAseq/Seurat3.rmd new file mode 100644 index 000000000..0112ab541 --- /dev/null +++ b/tcrdb/resources/external/scRNAseq/Seurat3.rmd @@ -0,0 +1,168 @@ +```{r Setup} + +knitr::opts_chunk$set(message=FALSE, warning=FALSE,echo=TRUE,error = FALSE) +library(knitr) +library(OOSAP) + +cores <- Sys.getenv('SEQUENCEANALYSIS_MAX_THREADS') +if (cores != ''){ + print(paste0('Setting future::plan to ', cores, ' cores')) + future::plan("multiprocess", workers = as.integer(cores)) +} else { + print('SEQUENCEANALYSIS_MAX_THREADS not set, will not set cores') +} + +print('Updating future.globals.maxSize') +options(future.globals.maxSize = Inf) + +print('Global variables: ') +for (v in c('outPrefix', 'resolutionToUse', 'dimsToUse', 'minDimsToUse', 'doCellFilter', 'doCellCycle', 'useSCTransform', 'runSingleR', 'mergeMethod')){ + if (exists(v)){ + print(paste0(v, ': ', get(v))) + } else { + print(paste0(v, ': not defined')) + } +} + +``` + +## Prepare data + +```{r PreparingData} + +rawDataSaveFile <- paste0(outPrefix, '.rawData.rds') +seuratObjs <- list() +if (file.exists(rawDataSaveFile)) { + print('resuming from file') + seuratObjs <- readRDS(rawDataSaveFile) +} else { + for (datasetName in names(data)) { + print(paste0('Loading dataset: ', datasetName)) + seuratObjs[[datasetName]] <- ReadAndFilter10xData(dataDir = data[[datasetName]], datasetName = datasetName) + + print(seuratObjs[[datasetName]]) + } + + saveRDS(seuratObjs, file = rawDataSaveFile) +} + +``` + +## Merge data + +```{r MergeDatasets} + +seuratObj <- NULL +saveFile <- paste0(outPrefix, '.seurat.rds') +if (file.exists(saveFile)) { + print('resuming from file') + seuratObj <- readRDS(saveFile) +} else { + seuratObj <- MergeSeuratObjs(seuratObjs, metadata = data, method = mergeMethod) + saveRDS(seuratObj, file = saveFile) + rm(seuratObjs) +} + +print(seuratObj) + +``` + +## Initial Processing + +```{r InitialProcessing} + +seuratObj <- ProcessSeurat1(seuratObj, variableGeneTable = paste0(outPrefix, '.variableGenes.txt'), doCellFilter = doCellFilter, doCellCycle = doCellCycle, useSCTransform = useSCTransform, saveFile = saveFile) + +print(seuratObj) + +``` + +## DimRedux + +```{r DimRedux} + +seuratObj <- FindClustersAndDimRedux(seuratObj, dimsToUse = dimsToUse, minDimsToUse = minDimsToUse, saveFile = saveFile) + +Find_Markers(seuratObj, resolutionToUse = resolutionToUse, outFile = paste0(outPrefix, '.markers.txt'), saveFileMarkers = paste0(outPrefix, '.markers.rds')) + +print(seuratObj) + +``` + +## SingleR + +```{r SingleR} + +if (runSingleR) { + tryCatch({ + seuratObj <- RunSingleR(seuratObj = seuratObj, resultTableFile = paste0(outPrefix, '.singleR.txt')) + saveRDS(seuratObj, file = saveFile) + + DimPlot_SingleRClassLabs(seuratObj, plotIndividually = T) + + Tabulate_SingleRClassLabs(seuratObj, plotIndividually = T) + }, error = function(e){ + print('There was an error in SingleR') + + saveRDS(e, file = 'error.rds') + }) + + print(seuratObj) +} else { + print('SingleR will not be run') +} + +``` + +## Phenotypes + +```{r Phenotypes} + +PlotImmuneMarkers(seuratObj, reduction = 'tsne') + +PlotImmuneMarkers(seuratObj, reduction = 'umap') + +if (length(unique(seuratObj$BarcodePrefix)) > 1) { + print(Seurat::DimPlot(seuratObj, reduction = 'pca', group.by = 'BarcodePrefix', label = T)) + print(Seurat::DimPlot(seuratObj, reduction = 'tsne', group.by = 'BarcodePrefix', label = T)) + print(Seurat::DimPlot(seuratObj, reduction = 'umap', group.by = 'BarcodePrefix', label = T)) + + t <- table(Cluster = Seurat::Idents(seuratObj), Dataset = seuratObj$BarcodePrefix) + t <- round(t / colSums(t), 2) + knitr::kable(t) +} + +``` + +## Activation + +```{r ActivationScore} + +seuratObj <- ClassifySGSAndApply(seuratObj = seuratObj, geneSetName = 'HighlyActivated', geneList = OOSAP::Phenotyping_GeneList()$HighlyActivated, positivityThreshold = 0.5, saveFilePath = paste0(outPrefix, '.ha.txt')) +saveRDS(seuratObj, file = saveFile) + +``` + +## Write Summary + +```{r Summary} + +saveRDS(seuratObj, file = saveFile) +unlink(rawDataSaveFile) + +WriteSummaryMetrics(seuratObj, file = paste0(outPrefix, '.summary.txt')) + +SaveDimRedux(seuratObj, file = paste0(outPrefix, '.DimReduxComps.csv')) + +WriteCellBarcodes(seuratObj, file = paste0(outPrefix, '.cellBarcodes.csv')) + +``` + +## Print Session Info + +```{r SessionInfo} + +sessionInfo() + +``` + diff --git a/tcrdb/resources/external/scRNAseq/seuratWrapper.sh b/tcrdb/resources/external/scRNAseq/seuratWrapper.sh new file mode 100644 index 000000000..5a6f2d81b --- /dev/null +++ b/tcrdb/resources/external/scRNAseq/seuratWrapper.sh @@ -0,0 +1,27 @@ +#!/bin/bash + +set -e +set -u +set -x + +WD=`pwd` +HOME=`echo ~/` + +DOCKER=/opt/acc/sbin/exadocker +LK_ROOT=$1 + +RAM_OPTS="" +ENV_OPTS="" +if [ ! -z $SEQUENCEANALYSIS_MAX_RAM ];then + RAM_OPTS=" --memory=${SEQUENCEANALYSIS_MAX_RAM}g" + + ENV_OPTS=" -e SEQUENCEANALYSIS_MAX_RAM" +fi + +if [ ! -z SEQUENCEANALYSIS_MAX_THREADS ];then + ENV_OPTS=${ENV_OPTS}" -e SEQUENCEANALYSIS_MAX_THREADS="${SEQUENCEANALYSIS_MAX_THREADS} +fi + +sudo $DOCKER pull bimberlab/oosap + +sudo $DOCKER run --rm=true $RAM_OPTS $ENV_OPTS -v "${WD}:/work" -v "${HOME}:/homeDir" -u $UID -e USERID=$UID -w /work -e HOME=/homeDir bimberlab/oosap Rscript --vanilla script.R \ No newline at end of file From 03ec6b9932a615f96403712c8b4285fc7289a68e Mon Sep 17 00:00:00 2001 From: bbimber Date: Sun, 2 Feb 2020 20:15:24 -0800 Subject: [PATCH 08/25] Update cDNA TSV parsing --- tcrdb/src/org/labkey/tcrdb/pipeline/CellRangerVDJUtils.java | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/tcrdb/src/org/labkey/tcrdb/pipeline/CellRangerVDJUtils.java b/tcrdb/src/org/labkey/tcrdb/pipeline/CellRangerVDJUtils.java index d031ccf6b..547594e7b 100644 --- a/tcrdb/src/org/labkey/tcrdb/pipeline/CellRangerVDJUtils.java +++ b/tcrdb/src/org/labkey/tcrdb/pipeline/CellRangerVDJUtils.java @@ -390,7 +390,7 @@ public void importAssayData(PipelineJob job, AnalysisModel model, File outDir, I while ((line = reader.readNext()) != null) { //header - if (line[0].startsWith("TCR_ReadsetId")) + if (line[0].startsWith("ReadsetId")) { continue; } From 6d7cf2e1a104c002d76abe62c8fd8002a3adf514 Mon Sep 17 00:00:00 2001 From: bbimber Date: Mon, 3 Feb 2020 06:14:37 -0800 Subject: [PATCH 09/25] Abort cell hashing calls if only single HTO used --- .../tcrdb/pipeline/CellRangerVDJUtils.java | 23 +++++++++++-------- 1 file changed, 13 insertions(+), 10 deletions(-) diff --git a/tcrdb/src/org/labkey/tcrdb/pipeline/CellRangerVDJUtils.java b/tcrdb/src/org/labkey/tcrdb/pipeline/CellRangerVDJUtils.java index 547594e7b..52c77ecb0 100644 --- a/tcrdb/src/org/labkey/tcrdb/pipeline/CellRangerVDJUtils.java +++ b/tcrdb/src/org/labkey/tcrdb/pipeline/CellRangerVDJUtils.java @@ -220,13 +220,23 @@ public File runRemoteCellHashingTasks(PipelineStepOutput output, String outputCa _log.info("No cached hashing readsets, skipping"); return null; } - else if (readsetToHashing.size() == 1) + + //prepare whitelist of barcodes, based on cDNA records + File htoBarcodeWhitelist = getValidHashingBarcodeFile(); + if (!htoBarcodeWhitelist.exists()) + { + throw new PipelineJobException("Unable to find file: " + htoBarcodeWhitelist.getPath()); + } + + long lineCount = SequencePipelineService.get().getLineCount(htoBarcodeWhitelist); + if (lineCount == 1) { - _log.info("Only a single hashing readset exists, will not use hashing"); + _log.info("Only a HTO is used, will not use hashing"); return null; } - _log.debug("total cached readset/HTO pairs: " + readsetToHashing.size()); + _log.debug("total cached readset/hashing readset pairs: " + readsetToHashing.size()); + _log.debug("unique HTOs: " + lineCount); //prepare whitelist of cell indexes File cellBarcodeWhitelist = getValidCellIndexFile(); @@ -277,13 +287,6 @@ else if (readsetToHashing.size() == 1) throw new PipelineJobException(e); } - //prepare whitelist of barcodes, based on cDNA records - File htoBarcodeWhitelist = getValidHashingBarcodeFile(); - if (!htoBarcodeWhitelist.exists()) - { - throw new PipelineJobException("Unable to find file: " + htoBarcodeWhitelist.getPath()); - } - Readset htoReadset = support.getCachedReadset(readsetToHashing.get(rs.getReadsetId())); if (htoReadset == null) { From 2acb30dd3e2632b39d32af0fb9d1c46cbd65074a Mon Sep 17 00:00:00 2001 From: bbimber Date: Mon, 3 Feb 2020 10:29:30 -0800 Subject: [PATCH 10/25] Bugfix cite-seq --- .../org/labkey/tcrdb/pipeline/CellRangerSeuratHandler.java | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/tcrdb/src/org/labkey/tcrdb/pipeline/CellRangerSeuratHandler.java b/tcrdb/src/org/labkey/tcrdb/pipeline/CellRangerSeuratHandler.java index 602180f77..8a1083f1e 100644 --- a/tcrdb/src/org/labkey/tcrdb/pipeline/CellRangerSeuratHandler.java +++ b/tcrdb/src/org/labkey/tcrdb/pipeline/CellRangerSeuratHandler.java @@ -447,10 +447,10 @@ else if (rs.getReadsetId() == null) String[] line; while ((line = reader.readNext()) != null) { - if (hashingReadsetId.toString().equals(line[6])) + if (hashingReadsetId.toString().equals(line[7])) { htosForReadset++; - bcWriter.writeNext(new String[]{line[5], line[4]}); + bcWriter.writeNext(new String[]{line[6], line[5]}); } } } From 79433de850fcac7ad7b14aabc1e9a9550953f190 Mon Sep 17 00:00:00 2001 From: bbimber Date: Mon, 3 Feb 2020 10:50:37 -0800 Subject: [PATCH 11/25] When total cell barcodes with CDR3 is low, allow all valid cells to be included for hashing --- .../tcrdb/pipeline/CellRangerVDJUtils.java | 35 +++++++++++++++---- 1 file changed, 28 insertions(+), 7 deletions(-) diff --git a/tcrdb/src/org/labkey/tcrdb/pipeline/CellRangerVDJUtils.java b/tcrdb/src/org/labkey/tcrdb/pipeline/CellRangerVDJUtils.java index 52c77ecb0..8ca670b72 100644 --- a/tcrdb/src/org/labkey/tcrdb/pipeline/CellRangerVDJUtils.java +++ b/tcrdb/src/org/labkey/tcrdb/pipeline/CellRangerVDJUtils.java @@ -241,6 +241,7 @@ public File runRemoteCellHashingTasks(PipelineStepOutput output, String outputCa //prepare whitelist of cell indexes File cellBarcodeWhitelist = getValidCellIndexFile(); Set uniqueBarcodes = new HashSet<>(); + Set uniqueBarcodesIncludingNoCDR3 = new HashSet<>(); _log.debug("writing cell barcodes"); try (CSVWriter writer = new CSVWriter(PrintWriters.getPrintWriter(cellBarcodeWhitelist), ',', CSVWriter.NO_QUOTE_CHARACTER); CSVReader reader = new CSVReader(Readers.getReader(perCellTsv), ',')) { @@ -254,32 +255,36 @@ public File runRemoteCellHashingTasks(PipelineStepOutput output, String outputCa rowIdx++; if (rowIdx > 1) { - if (row.length >= 13 && "None".equals(row[12])) + if ("False".equalsIgnoreCase(row[1])) { - noCallRows++; + nonCell++; continue; } - if ("False".equalsIgnoreCase(row[1])) + //NOTE: allow these to pass for cell-hashing under some conditions + boolean hasCDR3 = !"None".equals(row[12]); + if (!hasCDR3) { - nonCell++; - continue; + noCallRows++; } //NOTE: 10x appends "-1" to barcodes String barcode = row[0].split("-")[0]; - if (!uniqueBarcodes.contains(barcode)) + if (hasCDR3 && !uniqueBarcodes.contains(barcode)) { writer.writeNext(new String[]{barcode}); uniqueBarcodes.add(barcode); } + + uniqueBarcodesIncludingNoCDR3.add(barcode); } } _log.debug("rows inspected: " + (rowIdx - 1)); _log.debug("rows without CDR3: " + noCallRows); _log.debug("rows not called as cells: " + nonCell); - _log.debug("unique cell barcodes: " + uniqueBarcodes.size()); + _log.debug("unique cell barcodes (with CDR3): " + uniqueBarcodes.size()); + _log.debug("unique cell barcodes (including no CDR3): " + uniqueBarcodesIncludingNoCDR3.size()); output.addIntermediateFile(cellBarcodeWhitelist); } catch (IOException e) @@ -287,6 +292,22 @@ public File runRemoteCellHashingTasks(PipelineStepOutput output, String outputCa throw new PipelineJobException(e); } + if (uniqueBarcodes.size() < 500 && uniqueBarcodesIncludingNoCDR3.size() > uniqueBarcodes.size()) + { + _log.info("Total cell barcodes with CDR3s is low, so cell hashing will be performing using an input that includes valid cells that lacked CDR3 data."); + try (CSVWriter writer = new CSVWriter(PrintWriters.getPrintWriter(cellBarcodeWhitelist), ',', CSVWriter.NO_QUOTE_CHARACTER)) + { + for (String barcode : uniqueBarcodesIncludingNoCDR3) + { + writer.writeNext(new String[]{barcode}); + } + } + catch (IOException e) + { + throw new PipelineJobException(e); + } + } + Readset htoReadset = support.getCachedReadset(readsetToHashing.get(rs.getReadsetId())); if (htoReadset == null) { From b5aa826a4646fb7f277a3fe0a903e5fe3ef9bce7 Mon Sep 17 00:00:00 2001 From: bbimber Date: Mon, 3 Feb 2020 11:54:00 -0800 Subject: [PATCH 12/25] remove duplicate percent signs --- .../src/org/labkey/tcrdb/pipeline/CellRangerVDJUtils.java | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/tcrdb/src/org/labkey/tcrdb/pipeline/CellRangerVDJUtils.java b/tcrdb/src/org/labkey/tcrdb/pipeline/CellRangerVDJUtils.java index 8ca670b72..2f4f7ec5d 100644 --- a/tcrdb/src/org/labkey/tcrdb/pipeline/CellRangerVDJUtils.java +++ b/tcrdb/src/org/labkey/tcrdb/pipeline/CellRangerVDJUtils.java @@ -638,9 +638,9 @@ else if (discordantBarcodes.contains(barcode)) _log.info("total rows not cells: " + nonCell); _log.info("total rows marked as cells: " + totalCells); _log.info("total clonotype rows without CDR3: " + noCDR3); - _log.info("total clonotype rows skipped for unknown barcodes: " + totalSkipped + " (" + (NumberFormat.getPercentInstance().format(totalSkipped / (double)totalCells)) + "%)"); - _log.info("total clonotype rows skipped because they are doublets: " + doubletSkipped + " (" + (NumberFormat.getPercentInstance().format(doubletSkipped / (double)totalCells)) + "%)"); - _log.info("total clonotype rows skipped because they are discordant calls: " + discordantSkipped + " (" + (NumberFormat.getPercentInstance().format(discordantSkipped / (double)totalCells)) + "%)"); + _log.info("total clonotype rows skipped for unknown barcodes: " + totalSkipped + " (" + (NumberFormat.getPercentInstance().format(totalSkipped / (double)totalCells)) + ")"); + _log.info("total clonotype rows skipped because they are doublets: " + doubletSkipped + " (" + (NumberFormat.getPercentInstance().format(doubletSkipped / (double)totalCells)) + ")"); + _log.info("total clonotype rows skipped because they are discordant calls: " + discordantSkipped + " (" + (NumberFormat.getPercentInstance().format(discordantSkipped / (double)totalCells)) + ")"); _log.info("unique known cell barcodes: " + knownBarcodes.size()); _log.info("total clonotypes: " + countMapBySample.size()); _log.info("total cells with CDR3, lacking clonotype: " + hasCDR3NoClonotype); @@ -701,7 +701,7 @@ else if (discordantBarcodes.contains(barcode)) totalCells += processRow(countData, cDNAMap, model, runId, am, totalCellsMapBySample, sequenceMap, rows); } - _log.info("total clones without count data: " + clonesWithoutCounts.size() + " (" + (NumberFormat.getPercentInstance().format(clonesWithoutCounts.size() / (double)uniqueClones.size())) + "%)"); + _log.info("total clones without count data: " + clonesWithoutCounts.size() + " (" + (NumberFormat.getPercentInstance().format(clonesWithoutCounts.size() / (double)uniqueClones.size())) + ")"); } catch (Exception e) { From 59a4b7ce2d71af72843e26c2158f4293a2f7a7e5 Mon Sep 17 00:00:00 2001 From: bbimber Date: Mon, 3 Feb 2020 15:19:18 -0800 Subject: [PATCH 13/25] Update logic to determine if cell hashing is used --- .../pipeline/CellRangerSeuratHandler.java | 4 ++++ .../tcrdb/pipeline/CellRangerVDJUtils.java | 22 ++++++++++++++----- 2 files changed, 21 insertions(+), 5 deletions(-) diff --git a/tcrdb/src/org/labkey/tcrdb/pipeline/CellRangerSeuratHandler.java b/tcrdb/src/org/labkey/tcrdb/pipeline/CellRangerSeuratHandler.java index 8a1083f1e..7ed6e818f 100644 --- a/tcrdb/src/org/labkey/tcrdb/pipeline/CellRangerSeuratHandler.java +++ b/tcrdb/src/org/labkey/tcrdb/pipeline/CellRangerSeuratHandler.java @@ -480,6 +480,10 @@ else if (rs.getReadsetId() == null) ctx.getLogger().info("Cell hashing was not used. will not append to seurat"); } } + else + { + ctx.getLogger().info("Cell hashing was not used"); + } } private void appendCallsToSeurat(JobContext ctx, File seuratObj, Map finalCalls) throws PipelineJobException diff --git a/tcrdb/src/org/labkey/tcrdb/pipeline/CellRangerVDJUtils.java b/tcrdb/src/org/labkey/tcrdb/pipeline/CellRangerVDJUtils.java index 2f4f7ec5d..3d02ec0b1 100644 --- a/tcrdb/src/org/labkey/tcrdb/pipeline/CellRangerVDJUtils.java +++ b/tcrdb/src/org/labkey/tcrdb/pipeline/CellRangerVDJUtils.java @@ -93,6 +93,7 @@ public void prepareHashingFilesIfNeeded(PipelineJob job, SequenceAnalysisJobSupp FieldKey.fromString("sortId/hto"), FieldKey.fromString("sortId/hto/sequence"), FieldKey.fromString("hashingReadsetId"), + FieldKey.fromString("hashingReadsetId/totalFiles"), FieldKey.fromString("status")) ); @@ -101,7 +102,7 @@ public void prepareHashingFilesIfNeeded(PipelineJob job, SequenceAnalysisJobSupp HashMap readsetToHashingMap = new HashMap<>(); try (CSVWriter writer = new CSVWriter(PrintWriters.getPrintWriter(output), '\t', CSVWriter.NO_QUOTE_CHARACTER); CSVWriter bcWriter = new CSVWriter(PrintWriters.getPrintWriter(barcodeOutput), ',', CSVWriter.NO_QUOTE_CHARACTER)) { - writer.writeNext(new String[]{"ReadsetId", "CDNA_ID", "AnimalId", "Stim", "Population", "HTO_Name", "HTO_Seq", "HashingReadsetId"}); + writer.writeNext(new String[]{"ReadsetId", "CDNA_ID", "AnimalId", "Stim", "Population", "HTO_Name", "HTO_Seq", "HashingReadsetId", "HasHashingReads"}); List cachedReadsets = support.getCachedReadsets(); Set distinctHTOs = new HashSet<>(); Set hashingStatus = new HashSet<>(); @@ -124,7 +125,8 @@ public void prepareHashingFilesIfNeeded(PipelineJob job, SequenceAnalysisJobSupp results.getString(FieldKey.fromString("sortId/population")), results.getString(FieldKey.fromString("sortId/hto")), results.getString(FieldKey.fromString("sortId/hto/sequence")), - String.valueOf(results.getObject(FieldKey.fromString("hashingReadsetId")) == null ? "" : results.getInt(FieldKey.fromString("hashingReadsetId"))) + String.valueOf(results.getObject(FieldKey.fromString("hashingReadsetId")) == null ? "" : results.getInt(FieldKey.fromString("hashingReadsetId"))), + String.valueOf(results.getObject(FieldKey.fromString("hashingReadsetId/totalFiles")) != null && results.getInt(FieldKey.fromString("hashingReadsetId/totalFiles")) > 0) }); boolean useCellHashing = results.getObject(FieldKey.fromString("sortId/hto")) != null; @@ -158,7 +160,7 @@ else if (useCellHashing) if (hashingStatus.size() > 1) { - throw new PipelineJobException("The selected readsets/cDNA records use a mixture of cell hashing and non-hashing."); + _log.info("The selected readsets/cDNA records use a mixture of cell hashing and non-hashing."); } else if (hashingStatus.isEmpty()) { @@ -176,7 +178,7 @@ else if (distinctHTOs.size() == 1) job.getLogger().info("There is only a single HTO in this pool, will not use hashing"); } - boolean useCellHashing = hashingStatus.iterator().next(); + boolean useCellHashing = hashingStatus.size() > 1 ? true : hashingStatus.iterator().next(); if (useCellHashing && distinctHTOs.isEmpty()) { throw new PipelineJobException("Cell hashing was selected, but no HTOs were found"); @@ -929,9 +931,19 @@ public static Map getCachedReadsetMap(SequenceAnalysisJobSuppo return support.getCachedObject(CellRangerVDJUtils.READSET_TO_HASHING_MAP, PipelineJob.createObjectMapper().getTypeFactory().constructParametricType(Map.class, Integer.class, Integer.class)); } + //NOTE: if readset ID is null, this will be interpreted as any readset using hashing public boolean useCellHashing(SequenceAnalysisJobSupport support) throws PipelineJobException { - return getCachedReadsetMap(support).size() > 1; + if (getCachedReadsetMap(support).isEmpty()) + return false; + + File htoBarcodeWhitelist = getValidHashingBarcodeFile(); + if (!htoBarcodeWhitelist.exists()) + { + throw new PipelineJobException("Unable to find file: " + htoBarcodeWhitelist.getPath()); + } + + return SequencePipelineService.get().getLineCount(htoBarcodeWhitelist) > 1; } public static class CDNA From e930d50ec050d43593d81d2fd79a23ae68cfff19 Mon Sep 17 00:00:00 2001 From: bbimber Date: Mon, 3 Feb 2020 15:46:18 -0800 Subject: [PATCH 14/25] Dont perform column/cell filtering on HTOs when a whitelist of HTOs is used --- .../org/labkey/tcrdb/pipeline/CellRangerCellHashingHandler.java | 2 +- tcrdb/src/org/labkey/tcrdb/pipeline/CellRangerVDJUtils.java | 2 +- 2 files changed, 2 insertions(+), 2 deletions(-) diff --git a/tcrdb/src/org/labkey/tcrdb/pipeline/CellRangerCellHashingHandler.java b/tcrdb/src/org/labkey/tcrdb/pipeline/CellRangerCellHashingHandler.java index c63db9210..ceb767a8a 100644 --- a/tcrdb/src/org/labkey/tcrdb/pipeline/CellRangerCellHashingHandler.java +++ b/tcrdb/src/org/labkey/tcrdb/pipeline/CellRangerCellHashingHandler.java @@ -278,7 +278,7 @@ public static File processBarcodeFile(SequenceOutputHandler.JobContext ctx, File PipelineStepOutput output = new DefaultPipelineStepOutput(); String basename = FileUtil.makeLegalName(rs.getName()); - File cellToHto = SequencePipelineService.get().runCiteSeqCount(output, category, htoReadset, htoBarcodeWhitelist, cellBarcodeWhitelist, ctx.getWorkingDirectory(), basename, ctx.getLogger(), extraParams, false, ctx.getSourceDirectory(), editDistance, scanEditDistances, rs, genomeId); + File cellToHto = SequencePipelineService.get().runCiteSeqCount(output, category, htoReadset, htoBarcodeWhitelist, cellBarcodeWhitelist, ctx.getWorkingDirectory(), basename, ctx.getLogger(), extraParams, false, false, ctx.getSourceDirectory(), editDistance, scanEditDistances, rs, genomeId); ctx.getFileManager().addStepOutputs(action, output); ctx.getFileManager().addOutput(action, category, cellToHto); diff --git a/tcrdb/src/org/labkey/tcrdb/pipeline/CellRangerVDJUtils.java b/tcrdb/src/org/labkey/tcrdb/pipeline/CellRangerVDJUtils.java index 3d02ec0b1..769d07b03 100644 --- a/tcrdb/src/org/labkey/tcrdb/pipeline/CellRangerVDJUtils.java +++ b/tcrdb/src/org/labkey/tcrdb/pipeline/CellRangerVDJUtils.java @@ -318,7 +318,7 @@ public File runRemoteCellHashingTasks(PipelineStepOutput output, String outputCa //run CiteSeqCount. this will use Multiseq to make calls per cell String basename = FileUtil.makeLegalName(rs.getName()); - File hashtagCalls = SequencePipelineService.get().runCiteSeqCount(output, outputCategory, htoReadset, htoBarcodeWhitelist, cellBarcodeWhitelist, workingDir, basename, _log, extraParams, false, sourceDir, editDistance, scanEditDistances, rs, genomeId); + File hashtagCalls = SequencePipelineService.get().runCiteSeqCount(output, outputCategory, htoReadset, htoBarcodeWhitelist, cellBarcodeWhitelist, workingDir, basename, _log, extraParams, false, false, sourceDir, editDistance, scanEditDistances, rs, genomeId); if (!hashtagCalls.exists()) { throw new PipelineJobException("Unable to find expected file: " + hashtagCalls.getPath()); From 55a64e80a6c66f738730c432793f531d0f544e76 Mon Sep 17 00:00:00 2001 From: bbimber Date: Mon, 3 Feb 2020 16:27:33 -0800 Subject: [PATCH 15/25] More directly set min reads/cell for cell hashing --- .../labkey/tcrdb/pipeline/CellRangerCellHashingHandler.java | 4 +++- .../org/labkey/tcrdb/pipeline/CellRangerSeuratHandler.java | 3 ++- .../tcrdb/pipeline/CellRangerVDJCellHashingHandler.java | 4 +++- tcrdb/src/org/labkey/tcrdb/pipeline/CellRangerVDJUtils.java | 4 ++-- tcrdb/src/org/labkey/tcrdb/pipeline/CellRangerVDJWrapper.java | 3 ++- .../org/labkey/tcrdb/pipeline/SeuratCellHashingHandler.java | 1 + 6 files changed, 13 insertions(+), 6 deletions(-) diff --git a/tcrdb/src/org/labkey/tcrdb/pipeline/CellRangerCellHashingHandler.java b/tcrdb/src/org/labkey/tcrdb/pipeline/CellRangerCellHashingHandler.java index ceb767a8a..5274dd0ef 100644 --- a/tcrdb/src/org/labkey/tcrdb/pipeline/CellRangerCellHashingHandler.java +++ b/tcrdb/src/org/labkey/tcrdb/pipeline/CellRangerCellHashingHandler.java @@ -47,6 +47,7 @@ public CellRangerCellHashingHandler() put("checked", true); }}, true), ToolParameterDescriptor.create("editDistance", "Edit Distance", null, "ldk-integerfield", null, 1), + ToolParameterDescriptor.create("minCountPerCell", "Min Reads/Cell", null, "ldk-integerfield", null, 3), ToolParameterDescriptor.create("useOutputFileContainer", "Submit to Source File Workbook", "If checked, each job will be submitted to the same workbook as the input file, as opposed to submitting all jobs to the same workbook. This is primarily useful if submitting a large batch of files to process separately. This only applies if 'Run Separately' is selected.", "checkbox", new JSONObject(){{ put("checked", true); }}, false) @@ -275,10 +276,11 @@ public static File processBarcodeFile(SequenceOutputHandler.JobContext ctx, File boolean scanEditDistances = ctx.getParams().optBoolean("scanEditDistances", false); int editDistance = ctx.getParams().optInt("editDistance", 2); + int minCountPerCell = ctx.getParams().optInt("minCountPerCell", 3); PipelineStepOutput output = new DefaultPipelineStepOutput(); String basename = FileUtil.makeLegalName(rs.getName()); - File cellToHto = SequencePipelineService.get().runCiteSeqCount(output, category, htoReadset, htoBarcodeWhitelist, cellBarcodeWhitelist, ctx.getWorkingDirectory(), basename, ctx.getLogger(), extraParams, false, false, ctx.getSourceDirectory(), editDistance, scanEditDistances, rs, genomeId); + File cellToHto = SequencePipelineService.get().runCiteSeqCount(output, category, htoReadset, htoBarcodeWhitelist, cellBarcodeWhitelist, ctx.getWorkingDirectory(), basename, ctx.getLogger(), extraParams, false, minCountPerCell, ctx.getSourceDirectory(), editDistance, scanEditDistances, rs, genomeId); ctx.getFileManager().addStepOutputs(action, output); ctx.getFileManager().addOutput(action, category, cellToHto); diff --git a/tcrdb/src/org/labkey/tcrdb/pipeline/CellRangerSeuratHandler.java b/tcrdb/src/org/labkey/tcrdb/pipeline/CellRangerSeuratHandler.java index 7ed6e818f..bf23386e2 100644 --- a/tcrdb/src/org/labkey/tcrdb/pipeline/CellRangerSeuratHandler.java +++ b/tcrdb/src/org/labkey/tcrdb/pipeline/CellRangerSeuratHandler.java @@ -81,7 +81,8 @@ public CellRangerSeuratHandler() ToolParameterDescriptor.create("scanEditDistances", "Scan Edit Distances (Hashing)", "If checked, CITE-seq-count will be run using edit distances from 0-3 and the iteration with the highest singlets will be used.", "checkbox", new JSONObject(){{ put("checked", true); }}, true), - ToolParameterDescriptor.create("editDistance", "Edit Distance (Hashing)", null, "ldk-integerfield", null, 1) + ToolParameterDescriptor.create("editDistance", "Edit Distance (Hashing)", null, "ldk-integerfield", null, 1), + ToolParameterDescriptor.create("minCountPerCell", "Min Reads/Cell", null, "ldk-integerfield", null, 3) )); } diff --git a/tcrdb/src/org/labkey/tcrdb/pipeline/CellRangerVDJCellHashingHandler.java b/tcrdb/src/org/labkey/tcrdb/pipeline/CellRangerVDJCellHashingHandler.java index 2152e29ba..1b60e2319 100644 --- a/tcrdb/src/org/labkey/tcrdb/pipeline/CellRangerVDJCellHashingHandler.java +++ b/tcrdb/src/org/labkey/tcrdb/pipeline/CellRangerVDJCellHashingHandler.java @@ -57,6 +57,7 @@ public CellRangerVDJCellHashingHandler() put("checked", true); }}, true), ToolParameterDescriptor.create("editDistance", "Edit Distance", null, "ldk-integerfield", null, 1), + ToolParameterDescriptor.create("minCountPerCell", "Min Reads/Cell", null, "ldk-integerfield", null, 3), ToolParameterDescriptor.create("useOutputFileContainer", "Submit to Source File Workbook", "If checked, each job will be submitted to the same workbook as the input file, as opposed to submitting all jobs to the same workbook. This is primarily useful if submitting a large batch of files to process separately. This only applies if 'Run Separately' is selected.", "checkbox", new JSONObject(){{ put("checked", true); }}, false) @@ -195,9 +196,10 @@ private void processVloupeFile(JobContext ctx, File perCellTsv, Readset rs, Reco //prepare whitelist of cell indexes AlignmentOutputImpl output = new AlignmentOutputImpl(); boolean scanEditDistances = ctx.getParams().optBoolean("scanEditDistances", false); + int minCountPerCell = ctx.getParams().optInt("minCountPerCell", 3); int editDistance = ctx.getParams().optInt("editDistance", 2); - File cellToHto = utils.runRemoteCellHashingTasks(output, CATEGORY, perCellTsv, rs, ctx.getSequenceSupport(), extraParams, ctx.getWorkingDirectory(), ctx.getSourceDirectory(), editDistance, scanEditDistances, genomeId); + File cellToHto = utils.runRemoteCellHashingTasks(output, CATEGORY, perCellTsv, rs, ctx.getSequenceSupport(), extraParams, ctx.getWorkingDirectory(), ctx.getSourceDirectory(), editDistance, scanEditDistances, genomeId, minCountPerCell); if (utils.useCellHashing(ctx.getSequenceSupport()) && cellToHto == null) { throw new PipelineJobException("Missing cell to HTO file"); diff --git a/tcrdb/src/org/labkey/tcrdb/pipeline/CellRangerVDJUtils.java b/tcrdb/src/org/labkey/tcrdb/pipeline/CellRangerVDJUtils.java index 769d07b03..a7402e002 100644 --- a/tcrdb/src/org/labkey/tcrdb/pipeline/CellRangerVDJUtils.java +++ b/tcrdb/src/org/labkey/tcrdb/pipeline/CellRangerVDJUtils.java @@ -214,7 +214,7 @@ public File getPerCellCsv(File outDir) return new File(outDir, "all_contig_annotations.csv"); } - public File runRemoteCellHashingTasks(PipelineStepOutput output, String outputCategory, File perCellTsv, Readset rs, SequenceAnalysisJobSupport support, List extraParams, File workingDir, File sourceDir, Integer editDistance, boolean scanEditDistances, Integer genomeId) throws PipelineJobException + public File runRemoteCellHashingTasks(PipelineStepOutput output, String outputCategory, File perCellTsv, Readset rs, SequenceAnalysisJobSupport support, List extraParams, File workingDir, File sourceDir, Integer editDistance, boolean scanEditDistances, Integer genomeId, Integer minCountPerCell) throws PipelineJobException { Map readsetToHashing = getCachedReadsetMap(support); if (readsetToHashing.isEmpty()) @@ -318,7 +318,7 @@ public File runRemoteCellHashingTasks(PipelineStepOutput output, String outputCa //run CiteSeqCount. this will use Multiseq to make calls per cell String basename = FileUtil.makeLegalName(rs.getName()); - File hashtagCalls = SequencePipelineService.get().runCiteSeqCount(output, outputCategory, htoReadset, htoBarcodeWhitelist, cellBarcodeWhitelist, workingDir, basename, _log, extraParams, false, false, sourceDir, editDistance, scanEditDistances, rs, genomeId); + File hashtagCalls = SequencePipelineService.get().runCiteSeqCount(output, outputCategory, htoReadset, htoBarcodeWhitelist, cellBarcodeWhitelist, workingDir, basename, _log, extraParams, false, minCountPerCell, sourceDir, editDistance, scanEditDistances, rs, genomeId); if (!hashtagCalls.exists()) { throw new PipelineJobException("Unable to find expected file: " + hashtagCalls.getPath()); diff --git a/tcrdb/src/org/labkey/tcrdb/pipeline/CellRangerVDJWrapper.java b/tcrdb/src/org/labkey/tcrdb/pipeline/CellRangerVDJWrapper.java index 3275a4d9b..9bf7d9bcb 100644 --- a/tcrdb/src/org/labkey/tcrdb/pipeline/CellRangerVDJWrapper.java +++ b/tcrdb/src/org/labkey/tcrdb/pipeline/CellRangerVDJWrapper.java @@ -347,8 +347,9 @@ public AlignmentStep.AlignmentOutput performAlignment(Readset rs, File inputFast { boolean scanEditDistances = getProvider().getParameterByName("scanEditDistances").extractValue(getPipelineCtx().getJob(), getProvider(), getStepIdx(), Boolean.class, false); int editDistance = getProvider().getParameterByName("editDistance").extractValue(getPipelineCtx().getJob(), getProvider(), getStepIdx(), Integer.class, 2); + int minCountPerCell = getProvider().getParameterByName("minCountPerCell").extractValue(getPipelineCtx().getJob(), getProvider(), getStepIdx(), Integer.class, 3); - getUtils().runRemoteCellHashingTasks(output, CellRangerVDJCellHashingHandler.CATEGORY, getUtils().getPerCellCsv(output.getBAM().getParentFile()), rs, getPipelineCtx().getSequenceSupport(), null, getPipelineCtx().getWorkingDirectory(), getPipelineCtx().getSourceDirectory(), editDistance, scanEditDistances, referenceGenome.getGenomeId()); + getUtils().runRemoteCellHashingTasks(output, CellRangerVDJCellHashingHandler.CATEGORY, getUtils().getPerCellCsv(output.getBAM().getParentFile()), rs, getPipelineCtx().getSequenceSupport(), null, getPipelineCtx().getWorkingDirectory(), getPipelineCtx().getSourceDirectory(), editDistance, scanEditDistances, referenceGenome.getGenomeId(), minCountPerCell); } else { diff --git a/tcrdb/src/org/labkey/tcrdb/pipeline/SeuratCellHashingHandler.java b/tcrdb/src/org/labkey/tcrdb/pipeline/SeuratCellHashingHandler.java index 59eaceee2..6449299ba 100644 --- a/tcrdb/src/org/labkey/tcrdb/pipeline/SeuratCellHashingHandler.java +++ b/tcrdb/src/org/labkey/tcrdb/pipeline/SeuratCellHashingHandler.java @@ -31,6 +31,7 @@ public SeuratCellHashingHandler() put("checked", true); }}, true), ToolParameterDescriptor.create("editDistance", "Edit Distance", null, "ldk-integerfield", null, 1), + ToolParameterDescriptor.create("minCountPerCell", "Min Reads/Cell", null, "ldk-integerfield", null, 3), ToolParameterDescriptor.create("useOutputFileContainer", "Submit to Source File Workbook", "If checked, each job will be submitted to the same workbook as the input file, as opposed to submitting all jobs to the same workbook. This is primarily useful if submitting a large batch of files to process separately..", "checkbox", new JSONObject() {{ put("checked", true); From c72c5a9431e018bf951db079e326b2f6d2496b3d Mon Sep 17 00:00:00 2001 From: bbimber Date: Mon, 3 Feb 2020 19:48:38 -0800 Subject: [PATCH 16/25] Make TSV reading more tolerant to different headers --- .../org/labkey/tcrdb/pipeline/CellRangerSeuratHandler.java | 2 ++ tcrdb/src/org/labkey/tcrdb/pipeline/MiXCRAnalysis.java | 4 +++- 2 files changed, 5 insertions(+), 1 deletion(-) diff --git a/tcrdb/src/org/labkey/tcrdb/pipeline/CellRangerSeuratHandler.java b/tcrdb/src/org/labkey/tcrdb/pipeline/CellRangerSeuratHandler.java index bf23386e2..6d7f503b7 100644 --- a/tcrdb/src/org/labkey/tcrdb/pipeline/CellRangerSeuratHandler.java +++ b/tcrdb/src/org/labkey/tcrdb/pipeline/CellRangerSeuratHandler.java @@ -496,6 +496,7 @@ private void appendCallsToSeurat(JobContext ctx, File seuratObj, Map { rWriter.println("'" + x + "' = '" + y.getName() + "'"); @@ -506,6 +507,7 @@ private void appendCallsToSeurat(JobContext ctx, File seuratObj, Map Date: Tue, 4 Feb 2020 06:54:11 -0800 Subject: [PATCH 17/25] Allow mixed cell hashing / non for combo seurat objects --- tcrdb/resources/external/scRNAseq/Seurat3.rmd | 4 ++++ .../org/labkey/tcrdb/pipeline/CellRangerSeuratHandler.java | 3 ++- 2 files changed, 6 insertions(+), 1 deletion(-) diff --git a/tcrdb/resources/external/scRNAseq/Seurat3.rmd b/tcrdb/resources/external/scRNAseq/Seurat3.rmd index 0112ab541..48eb50371 100644 --- a/tcrdb/resources/external/scRNAseq/Seurat3.rmd +++ b/tcrdb/resources/external/scRNAseq/Seurat3.rmd @@ -1,3 +1,7 @@ +--- +title: 'Seurat scRNA-seq Analysis' +--- + ```{r Setup} knitr::opts_chunk$set(message=FALSE, warning=FALSE,echo=TRUE,error = FALSE) diff --git a/tcrdb/src/org/labkey/tcrdb/pipeline/CellRangerSeuratHandler.java b/tcrdb/src/org/labkey/tcrdb/pipeline/CellRangerSeuratHandler.java index 6d7f503b7..a6331bd9c 100644 --- a/tcrdb/src/org/labkey/tcrdb/pipeline/CellRangerSeuratHandler.java +++ b/tcrdb/src/org/labkey/tcrdb/pipeline/CellRangerSeuratHandler.java @@ -438,7 +438,8 @@ else if (rs.getReadsetId() == null) Integer hashingReadsetId = CellRangerVDJUtils.getCachedReadsetMap(ctx.getSequenceSupport()).get(rs.getReadsetId()); if (hashingReadsetId == null) { - throw new PipelineJobException("Unable to find hashing readset Id for: " + rs.getReadsetId()); + ctx.getLogger().info("No hashing readset for: " + rs.getReadsetId() + ", this probably indicates either hashing is not used or the hashing data is not available."); + return; } File perReadsetHtos = new File(allCellBarcodes.getParentFile(), "allowableHtos." + barcodePrefix + ".txt"); From d6ed86753c4298c5966eb2bd96c682e867767c83 Mon Sep 17 00:00:00 2001 From: bbimber Date: Tue, 4 Feb 2020 11:37:49 -0800 Subject: [PATCH 18/25] continue, not return --- .../src/org/labkey/tcrdb/pipeline/CellRangerSeuratHandler.java | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/tcrdb/src/org/labkey/tcrdb/pipeline/CellRangerSeuratHandler.java b/tcrdb/src/org/labkey/tcrdb/pipeline/CellRangerSeuratHandler.java index a6331bd9c..52dc9016d 100644 --- a/tcrdb/src/org/labkey/tcrdb/pipeline/CellRangerSeuratHandler.java +++ b/tcrdb/src/org/labkey/tcrdb/pipeline/CellRangerSeuratHandler.java @@ -439,7 +439,7 @@ else if (rs.getReadsetId() == null) if (hashingReadsetId == null) { ctx.getLogger().info("No hashing readset for: " + rs.getReadsetId() + ", this probably indicates either hashing is not used or the hashing data is not available."); - return; + continue; } File perReadsetHtos = new File(allCellBarcodes.getParentFile(), "allowableHtos." + barcodePrefix + ".txt"); From 1893a067de7ab13b931c08701f57c8c2b7148102 Mon Sep 17 00:00:00 2001 From: bbimber Date: Tue, 4 Feb 2020 11:43:48 -0800 Subject: [PATCH 19/25] Clarify parameter --- .../src/org/labkey/tcrdb/pipeline/SeuratCellHashingHandler.java | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/tcrdb/src/org/labkey/tcrdb/pipeline/SeuratCellHashingHandler.java b/tcrdb/src/org/labkey/tcrdb/pipeline/SeuratCellHashingHandler.java index 6449299ba..a09d7f91a 100644 --- a/tcrdb/src/org/labkey/tcrdb/pipeline/SeuratCellHashingHandler.java +++ b/tcrdb/src/org/labkey/tcrdb/pipeline/SeuratCellHashingHandler.java @@ -31,7 +31,7 @@ public SeuratCellHashingHandler() put("checked", true); }}, true), ToolParameterDescriptor.create("editDistance", "Edit Distance", null, "ldk-integerfield", null, 1), - ToolParameterDescriptor.create("minCountPerCell", "Min Reads/Cell", null, "ldk-integerfield", null, 3), + ToolParameterDescriptor.create("minCountPerCell", "Min Reads/Cell (Cell Hashing)", null, "ldk-integerfield", null, 3), ToolParameterDescriptor.create("useOutputFileContainer", "Submit to Source File Workbook", "If checked, each job will be submitted to the same workbook as the input file, as opposed to submitting all jobs to the same workbook. This is primarily useful if submitting a large batch of files to process separately..", "checkbox", new JSONObject() {{ put("checked", true); From a593d8ed70b3e693672f7a5f805dc1b4493994d3 Mon Sep 17 00:00:00 2001 From: bbimber Date: Thu, 6 Feb 2020 12:21:21 -0800 Subject: [PATCH 20/25] capture changes to prime-seq LK install script --- primeseq/tools/installLabkey.sh | 37 ++++++++++++++++++++++++++------- 1 file changed, 29 insertions(+), 8 deletions(-) diff --git a/primeseq/tools/installLabkey.sh b/primeseq/tools/installLabkey.sh index 321da2d9c..f7d44ad87 100644 --- a/primeseq/tools/installLabkey.sh +++ b/primeseq/tools/installLabkey.sh @@ -10,23 +10,39 @@ labkey_home=/usr/local/labkey cd /usr/local/src #NOTE: corresponding changes must be made in javaWrapper.sh -MAJOR=18 -MINOR=2 -BRANCH=Discvr${MAJOR}${MINOR}_Installers -ARTIFACT=LabKey${MAJOR}.${MINOR} -MODULE_DIST_NAME=prime-seq-modules -PREMIUM=premium-${MAJOR}.${MINOR}.module +MAJOR=19 +MINOR_FULL="3.4" +MINOR_SHORT=3 +BRANCH=LabKey_Discvr_Discvr${MAJOR}${MINOR_SHORT}_Premuim_Installers TOMCAT_HOME=/usr/share/tomcat +TEAMCITY_USERNAME=username +MODULE_DIST_NAME=prime-seq-modules + +ARTIFACT=LabKey${MAJOR}.${MINOR_FULL} +PREMIUM=premium-${MAJOR}.${MINOR_SHORT}.module +DATAINTEGRATION=dataintegration-${MAJOR}.${MINOR_SHORT}.module + +isGzZip() { + RET=`file $1 | grep -E 'gzip compressed|Zip archive data' | wc -l` + if [ $RET == 0 ];then + echo "Not GZIP!" + exit 1 + else + echo "Is GZIP!" + fi +} #first download DATE=$(date +"%Y%m%d%H%M") MODULE_ZIP=${ARTIFACT}-ExtraModules-${DATE}.zip rm -Rf $MODULE_ZIP -wget --trust-server-names --no-check-certificate -O $MODULE_ZIP http://teamcity.labkey.org/guestAuth/repository/download/LabKey_${BRANCH}/.lastSuccessful/${MODULE_DIST_NAME}/${ARTIFACT}-{build.number}-ExtraModules.zip +wget -O $MODULE_ZIP https://${TEAMCITY_USERNAME}@teamcity.labkey.org/repository/download/${BRANCH}/.lastSuccessful/${MODULE_DIST_NAME}/${ARTIFACT}-{build.number}-ExtraModules.zip +isGzZip $MODULE_ZIP GZ=${ARTIFACT}-${DATE}-discvr-bin.tar.gz rm -Rf $GZ -wget --trust-server-names --no-check-certificate -O $GZ http://teamcity.labkey.org/guestAuth/repository/download/Labkey_${BRANCH}/.lastSuccessful/discvr/${ARTIFACT}-{build.number}-discvr-bin.tar.gz +wget -O $GZ https://${TEAMCITY_USERNAME}@teamcity.labkey.org/repository/download/${BRANCH}/.lastSuccessful/discvr/${ARTIFACT}-{build.number}-discvr-bin.tar.gz +isGzZip $GZ #extract, find name tar -xf $GZ @@ -55,6 +71,11 @@ if [ -e $PREMIUM ];then cp $PREMIUM ${labkey_home}/externalModules fi +#DataIntegration +if [ -e $DATAINTEGRATION ];then + cp $DATAINTEGRATION ${labkey_home}/externalModules +fi + #main server echo "Installing LabKey using: $GZ" cd $DIR From ea8a75ea48fb962c316d9dc95113f3ad92a6b3dc Mon Sep 17 00:00:00 2001 From: bbimber Date: Thu, 6 Feb 2020 21:47:49 -0800 Subject: [PATCH 21/25] Allow import of sorts/cDNA across workbooks --- .../web/tcrdb/panel/PoolImportPanel.js | 52 ++++++++++++------- .../src/org/labkey/tcrdb/TCRdbController.java | 32 ++++++++---- tcrdb/src/org/labkey/tcrdb/TCRdbProvider.java | 2 +- 3 files changed, 58 insertions(+), 28 deletions(-) diff --git a/tcrdb/resources/web/tcrdb/panel/PoolImportPanel.js b/tcrdb/resources/web/tcrdb/panel/PoolImportPanel.js index da4a68acd..cdaec03fa 100644 --- a/tcrdb/resources/web/tcrdb/panel/PoolImportPanel.js +++ b/tcrdb/resources/web/tcrdb/panel/PoolImportPanel.js @@ -1,10 +1,9 @@ Ext4.define('TCRdb.panel.PoolImportPanel', { extend: 'Ext.panel.Panel', - //TODO: replicate, buffer? COLUMNS: [{ - name: 'expt', - labels: ['Expt', 'Expt #', 'Experiment', 'Exp#', 'Exp #', 'Workbook'], + name: 'workbook', + labels: ['Experiment/Workbook', 'Expt', 'Expt #', 'Experiment', 'Exp#', 'Exp #', 'Workbook', 'Workbook #'], allowRowSpan: true, alwaysShow: true, transform: 'expt', @@ -27,7 +26,7 @@ Ext4.define('TCRdb.panel.PoolImportPanel', { labels: ['Animal Id', 'SubjectId', 'Subject Id'], allowRowSpan: true, allowBlank: false, - transform: 'animal', + transform: 'animal' },{ name: 'sampleDate', labels: ['Sample Date', 'Date'], @@ -149,13 +148,15 @@ Ext4.define('TCRdb.panel.PoolImportPanel', { return; } + var barcodeSeries = panel.down('#barcodeSeries').getValue(); val = val.toUpperCase(); - if (!val.match(/^SI-GA-/)) { + var re = new RegExp('^' + barcodeSeries + '-', 'i'); + if (!val.match(re)) { if (val.length > 3) { //errorMsgs.push('Every row must have name, application and proper barcodes'); } else { - val = 'SI-GA-' + val; + val = barcodeSeries + '-' + val; } } @@ -178,9 +179,10 @@ Ext4.define('TCRdb.panel.PoolImportPanel', { return val ? Ext4.data.Types.INTEGER.convert(val) : val; }, - pool: function(val, panel){ - if (panel.EXPERIMENT && Ext4.isNumeric(val) && panel.EXPERIMENT !== val){ - return panel.EXPERIMENT + '-' + val; + pool: function(val, panel, row){ + var workbook = row.workbook || panel.EXPERIMENT; + if (workbook && Ext4.isNumeric(val) && workbook !== val){ + return workbook + '-' + val; } return val; @@ -355,6 +357,13 @@ Ext4.define('TCRdb.panel.PoolImportPanel', { field.up('panel').down('#requireHTO').setValue(!val); } } + },{ + xtype: 'ldk-simplecombo', + fieldLabel: '10x Barcode Series', + itemId: 'barcodeSeries', + forceSelection: true, + storeValues: ['SI-GA'], + value: 'SI-GA' },{ xtype: 'textarea', fieldLabel: 'Paste Data Below', @@ -401,13 +410,15 @@ Ext4.define('TCRdb.panel.PoolImportPanel', { var parsedRows = this.parseRows(colArray, rows); var stimRows = []; Ext4.Array.forEach(parsedRows, function(r){ + LDK.Assert.assertNotEmpty('Expected non-null workbook', r.workbook); stimRows.push({ animalId: r.animalId, date: r.sampleDate, stim: r.stim, treatment: r.treatment || 'None', objectId: r.objectId, - population: r.population + population: r.population, + workbook: r.workbook }); }, this); @@ -484,7 +495,7 @@ Ext4.define('TCRdb.panel.PoolImportPanel', { Ext4.Array.forEach(rows, function(row, rowIdx){ var data = { objectId: LABKEY.Utils.generateUUID() - } + }; Ext4.Array.forEach(colArray, function(col, colIdx){ var cell = Ext4.isDefined(col.dataIdx) ? row[col.dataIdx] : ''; @@ -541,6 +552,7 @@ Ext4.define('TCRdb.panel.PoolImportPanel', { stimMap[key] = guid; stimIdxs[key] = stimIdx; + LDK.Assert.assertNotEmpty('Expected non-null workbook', row.workbook); ret.stimRows.push({ rowId: row.stimId || null, animalId: row.animalId, @@ -549,7 +561,7 @@ Ext4.define('TCRdb.panel.PoolImportPanel', { effector: row.effector, treatment: row.treatment || 'None', objectId: guid, - container: LABKEY.Security.currentContainer.id + workbook: row.workbook }); } @@ -559,6 +571,7 @@ Ext4.define('TCRdb.panel.PoolImportPanel', { //sorts: var sortMap = {}; Ext4.Array.forEach(parsedRows, function(row){ + LDK.Assert.assertNotEmpty('Expected non-null workbook', row.workbook); var stimGUID = stimMap[this.getStimKey(row)]; var key = this.getSortKey(row); if (!sortMap[key]){ @@ -575,7 +588,7 @@ Ext4.define('TCRdb.panel.PoolImportPanel', { hto: row.hto, buffer: row.buffer, objectId: guid, - container: LABKEY.Security.currentContainer.id + workbook: row.workbook }); } }, this); @@ -626,12 +639,13 @@ Ext4.define('TCRdb.panel.PoolImportPanel', { Ext4.Array.forEach(rowArr, function(row) { var sortKey = this.getSortKey(row); + LDK.Assert.assertNotEmpty('Expected non-null workbook', row.workbook); var cDNA = Ext4.apply({ sortGUID: sortMap[sortKey], chemistry: null, plateId: row.plateId, well: row.well || 'Pool', - container: LABKEY.Security.currentContainer.id + workbook: row.workbook }, readsetGUIDs); ret.cDNARows.push(cDNA); @@ -651,6 +665,7 @@ Ext4.define('TCRdb.panel.PoolImportPanel', { var idxValues = this.getUniqueValues(rowArr, prefix + '_library_index'); var conc = this.getUniqueValues(rowArr, prefix + '_library_conc'); var fragment = this.getUniqueValues(rowArr, prefix + '_library_fragment'); + var workbook = this.getUniqueValues(rowArr, 'workbook'); var subjectid = this.getUniqueValues(rowArr, 'animalId'); subjectid = subjectid.length === 1 ? subjectid[0] : null; @@ -662,6 +677,7 @@ Ext4.define('TCRdb.panel.PoolImportPanel', { } var guid = LABKEY.Utils.generateUUID(); + LDK.Assert.assertNotEmpty('Expected non-null workbook', workbook); readsetRows.push({ name: poolName + '-' + type, barcode5: idxValues[0], @@ -673,7 +689,7 @@ Ext4.define('TCRdb.panel.PoolImportPanel', { subjectid: subjectid, sampleType: 'mRNA', objectId: guid, - container: LABKEY.Security.currentContainer.id + workbook: workbook }); return guid; @@ -714,7 +730,7 @@ Ext4.define('TCRdb.panel.PoolImportPanel', { var data = []; var missingValues = false; - var requireHTO = this.down('#requireHTO').getValue() || this.down('#requireHashTag').getValue(); + var requireHTO = this.down('#requireHTO').getValue() || (this.down('#requireHashTag') && this.down('#requireHashTag').getValue()); Ext4.Array.forEach(parsedRows, function(row, rowIdx){ var toAdd = [rowIdx + 1]; Ext4.Array.forEach(colIdxs, function(colIdx){ @@ -771,14 +787,14 @@ Ext4.define('TCRdb.panel.PoolImportPanel', { onSubmit: function(e, dt, node, config){ Ext4.Msg.wait('Saving...'); LABKEY.Ajax.request({ - url: LABKEY.ActionURL.buildURL('tcrdb', 'importTenx'), + url: LABKEY.ActionURL.buildURL('tcrdb', 'importTenx', Laboratory.Utils.getQueryContainerPath()), method: 'POST', jsonData: config.rowData.groupedRows, scope: this, success: function(){ Ext4.Msg.hide(); Ext4.Msg.alert('Success', 'Data Saved', function(){ - window.location = LABKEY.ActionURL.buildURL('query', 'executeQuery.view', null, {'query.queryName': 'cdnas', schemaName: 'tcrdb'}) + window.location = LABKEY.ActionURL.buildURL('query', 'executeQuery.view', Laboratory.Utils.getQueryContainerPath(), {'query.queryName': 'cdnas', schemaName: 'tcrdb'}) }, this); }, failure: LDK.Utils.getErrorCallback() diff --git a/tcrdb/src/org/labkey/tcrdb/TCRdbController.java b/tcrdb/src/org/labkey/tcrdb/TCRdbController.java index 5e7179db1..c679a61c7 100644 --- a/tcrdb/src/org/labkey/tcrdb/TCRdbController.java +++ b/tcrdb/src/org/labkey/tcrdb/TCRdbController.java @@ -35,6 +35,8 @@ import org.labkey.api.data.ColumnInfo; import org.labkey.api.data.CompareType; import org.labkey.api.data.Container; +import org.labkey.api.data.ContainerManager; +import org.labkey.api.data.ContainerType; import org.labkey.api.data.DbScope; import org.labkey.api.data.Selector; import org.labkey.api.data.SimpleFilter; @@ -977,10 +979,10 @@ public static class ImportTenXAction extends MutatingApiAction> stimRows = parseRows(form, "stimRows"); - List> sortRows = parseRows(form, "sortRows"); - List> readsetRows = parseRows(form, "readsetRows"); - List> cDNARows = parseRows(form, "cDNARows"); + List> stimRows = parseRows(form, "stimRows", getContainer()); + List> sortRows = parseRows(form, "sortRows", getContainer()); + List> readsetRows = parseRows(form, "readsetRows", getContainer()); + List> cDNARows = parseRows(form, "cDNARows", getContainer()); UserSchema tcrdb = QueryService.get().getUserSchema(getUser(), getContainer(), TCRdbSchema.NAME); UserSchema sequenceAnalysis = QueryService.get().getUserSchema(getUser(), getContainer(), "sequenceanalysis"); @@ -1008,7 +1010,7 @@ public Object execute(SimpleApiJsonForm form, BindException errors) throws Excep }); - List> insertedStimRows = tcrdb.getTable(TCRdbSchema.TABLE_STIMS).getUpdateService().insertRows(getUser(), getContainer(), stimRowsToInsert, bve, null, new HashMap<>()); + List> insertedStimRows = tcrdb.getTable(TCRdbSchema.TABLE_STIMS, null).getUpdateService().insertRows(getUser(), getContainer(), stimRowsToInsert, bve, null, new HashMap<>()); if (bve.hasErrors()) { throw bve; @@ -1042,7 +1044,7 @@ public Object execute(SimpleApiJsonForm form, BindException errors) throws Excep } }); - sortRows = tcrdb.getTable(TCRdbSchema.TABLE_SORTS).getUpdateService().insertRows(getUser(), getContainer(), sortRowsToInsert, bve, null, new HashMap<>()); + sortRows = tcrdb.getTable(TCRdbSchema.TABLE_SORTS, null).getUpdateService().insertRows(getUser(), getContainer(), sortRowsToInsert, bve, null, new HashMap<>()); if (bve.hasErrors()) { throw bve; @@ -1057,7 +1059,7 @@ public Object execute(SimpleApiJsonForm form, BindException errors) throws Excep sortMap.put((String) r.get("objectId"), (Integer) r.get("rowId")); }); - readsetRows = sequenceAnalysis.getTable("sequence_readsets").getUpdateService().insertRows(getUser(), getContainer(), readsetRows, bve, null, new HashMap<>()); + readsetRows = sequenceAnalysis.getTable("sequence_readsets", null).getUpdateService().insertRows(getUser(), getContainer(), readsetRows, bve, null, new HashMap<>()); if (bve.hasErrors()) { throw bve; @@ -1096,7 +1098,7 @@ public Object execute(SimpleApiJsonForm form, BindException errors) throws Excep } } - private static List> parseRows(SimpleApiJsonForm form, String propName) throws ApiUsageException + private static List> parseRows(SimpleApiJsonForm form, String propName, Container container) throws ApiUsageException { if (!form.getJsonObject().containsKey(propName)) { @@ -1109,6 +1111,18 @@ private static List> parseRows(SimpleApiJsonForm form, Strin Arrays.stream(arr.toJSONObjectArray()).forEach(m -> { Map map = new CaseInsensitiveHashMap<>(); map.putAll(m); + + if (map.containsKey("workbook")) + { + Container parent = container.getContainerFor(ContainerType.DataType.folderManagement); + Container workbook = ContainerManager.getForPath(parent.getPath() + "/" + map.get("workbook").toString()); + if (workbook == null) + { + throw new IllegalArgumentException("Unable to identify matching workbook for: " + map.get("workbook")); + } + + map.put("container", workbook == null ? null : workbook.getId()); + } ret.add(map); }); @@ -1125,7 +1139,7 @@ public Object execute(SimpleApiJsonForm form, BindException errors) throws Excep { ApiSimpleResponse resp = new ApiSimpleResponse(); - List> stimRows = parseRows(form, "stimRows"); + List> stimRows = parseRows(form, "stimRows", getContainer()); UserSchema us = QueryService.get().getUserSchema(getUser(), getContainer(), TCRdbSchema.NAME); if (us == null) diff --git a/tcrdb/src/org/labkey/tcrdb/TCRdbProvider.java b/tcrdb/src/org/labkey/tcrdb/TCRdbProvider.java index 274b75751..286add55a 100644 --- a/tcrdb/src/org/labkey/tcrdb/TCRdbProvider.java +++ b/tcrdb/src/org/labkey/tcrdb/TCRdbProvider.java @@ -81,7 +81,7 @@ public List getDataNavItems(Container c, User u) return Collections.emptyList(); } - TCRdbImportNavItem item = new TCRdbImportNavItem(this, "TCR Stims/Sorts", LaboratoryService.NavItemCategory.data, NAME); + TCRdbImportNavItem item = new TCRdbImportNavItem(this, "TCR Stims/Sorts (SMART-seq)", LaboratoryService.NavItemCategory.data, NAME); item.setQueryCache(cache); items.add(item); From a5b7d14e958ce070a89b1b2717062fec076ea653 Mon Sep 17 00:00:00 2001 From: bbimber Date: Fri, 7 Feb 2020 12:10:38 -0800 Subject: [PATCH 22/25] Drop ##META lines from cassandra header for HTSJDK compatibility --- mGAP/src/org/labkey/mgap/pipeline/CassandraRunner.java | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/mGAP/src/org/labkey/mgap/pipeline/CassandraRunner.java b/mGAP/src/org/labkey/mgap/pipeline/CassandraRunner.java index 8c2f801ae..db04e37e7 100644 --- a/mGAP/src/org/labkey/mgap/pipeline/CassandraRunner.java +++ b/mGAP/src/org/labkey/mgap/pipeline/CassandraRunner.java @@ -104,7 +104,7 @@ private void correctHeaderAndBGzip(File inputUnzip, File outputGzip) throws Pipe writer.write("set -x\n"); writer.write("set -e\n"); writer.write("{\n"); - writer.write("cat " + inputUnzip.getPath() + " | head -n 50000 | grep -e '^#' | sed 's/Number=0,Type=String/Number=1,Type=String/';\n"); + writer.write("cat " + inputUnzip.getPath() + " | head -n 50000 | grep -e '^#' | grep -v '^##META' | sed 's/Number=0,Type=String/Number=1,Type=String/';\n"); writer.write("cat " + inputUnzip.getPath() + " | grep -v '^#';\n"); writer.write("} | bgzip > " + outputGzip + "\n"); } From 732e903d48fe281f5f7717a45554c0284617581d Mon Sep 17 00:00:00 2001 From: bbimber Date: Fri, 7 Feb 2020 12:59:24 -0800 Subject: [PATCH 23/25] Add action to auto-create branches to match LabKey release branches --- .github/workflows/sync-release-branches.yml | 16 ++++++++++++++++ 1 file changed, 16 insertions(+) create mode 100644 .github/workflows/sync-release-branches.yml diff --git a/.github/workflows/sync-release-branches.yml b/.github/workflows/sync-release-branches.yml new file mode 100644 index 000000000..4ccbde0c3 --- /dev/null +++ b/.github/workflows/sync-release-branches.yml @@ -0,0 +1,16 @@ +# Designed to keep develop branch as a perfect copy of LabKey fork +on: + schedule: + - cron: "*/15 * * * *" +jobs: + sync-develop: + runs-on: ubuntu-latest + steps: + - name: "Sync Release Branches" + uses: bimberlabinternal/DevOps/githubActions/branch-create@master + with: + source_repo: "labkey/DiscvrLabkeyModules" + source_branch_prefix: "release" + destination_repo: "BimberLab/DiscvrLabkeyModules" + destination_branch_prefix: "discvr-" + github_token: ${{ secrets.GITHUB_TOKEN }} From 6b8329d14dc204d356540ef82fbb5dd6d011e43a Mon Sep 17 00:00:00 2001 From: bbimber Date: Fri, 7 Feb 2020 13:05:23 -0800 Subject: [PATCH 24/25] Add more github workflows --- .github/workflows/sync-develop.yml | 15 +++++++++++++++ .github/workflows/sync-release-branches.yml | 4 ++-- 2 files changed, 17 insertions(+), 2 deletions(-) create mode 100644 .github/workflows/sync-develop.yml diff --git a/.github/workflows/sync-develop.yml b/.github/workflows/sync-develop.yml new file mode 100644 index 000000000..d690d89b4 --- /dev/null +++ b/.github/workflows/sync-develop.yml @@ -0,0 +1,15 @@ +# Designed to keep develop branch as a perfect copy of LabKey fork +on: + schedule: + - cron: "*/15 * * * *" +jobs: + sync-develop: + runs-on: ubuntu-latest + steps: + - name: "Sync Develop Branch" + uses: bimberlabinternal/DevOps/githubActions/git-sync@master + with: + source_repo: "labkey/BimberLabKeyModules" + source_branch: "develop" + destination_branch: "develop" + github_token: ${{ secrets.GITHUB_TOKEN }} diff --git a/.github/workflows/sync-release-branches.yml b/.github/workflows/sync-release-branches.yml index 4ccbde0c3..2c414779e 100644 --- a/.github/workflows/sync-release-branches.yml +++ b/.github/workflows/sync-release-branches.yml @@ -9,8 +9,8 @@ jobs: - name: "Sync Release Branches" uses: bimberlabinternal/DevOps/githubActions/branch-create@master with: - source_repo: "labkey/DiscvrLabkeyModules" + source_repo: "labkey/BimberLabKeyModules" source_branch_prefix: "release" - destination_repo: "BimberLab/DiscvrLabkeyModules" + destination_repo: "BimberLabInternal/BimberLabKeyModules" destination_branch_prefix: "discvr-" github_token: ${{ secrets.GITHUB_TOKEN }} From af65763242972c2f6531c7583b0cb660eaefb7f7 Mon Sep 17 00:00:00 2001 From: bbimber Date: Fri, 7 Feb 2020 17:53:45 -0800 Subject: [PATCH 25/25] Update name of branch sync task --- .github/workflows/sync-release-branches.yml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/.github/workflows/sync-release-branches.yml b/.github/workflows/sync-release-branches.yml index 2c414779e..b8474416f 100644 --- a/.github/workflows/sync-release-branches.yml +++ b/.github/workflows/sync-release-branches.yml @@ -3,7 +3,7 @@ on: schedule: - cron: "*/15 * * * *" jobs: - sync-develop: + sync-release-branches: runs-on: ubuntu-latest steps: - name: "Sync Release Branches"