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 new file mode 100644 index 000000000..b8474416f --- /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-release-branches: + runs-on: ubuntu-latest + steps: + - name: "Sync Release Branches" + uses: bimberlabinternal/DevOps/githubActions/branch-create@master + with: + source_repo: "labkey/BimberLabKeyModules" + source_branch_prefix: "release" + destination_repo: "BimberLabInternal/BimberLabKeyModules" + destination_branch_prefix: "discvr-" + github_token: ${{ secrets.GITHUB_TOKEN }} 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/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"); } 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()); + } + } } } } 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 diff --git a/tcrdb/resources/external/scRNAseq/Seurat3.rmd b/tcrdb/resources/external/scRNAseq/Seurat3.rmd new file mode 100644 index 000000000..48eb50371 --- /dev/null +++ b/tcrdb/resources/external/scRNAseq/Seurat3.rmd @@ -0,0 +1,172 @@ +--- +title: 'Seurat scRNA-seq Analysis' +--- + +```{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 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/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/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); diff --git a/tcrdb/src/org/labkey/tcrdb/pipeline/CellRangerCellHashingHandler.java b/tcrdb/src/org/labkey/tcrdb/pipeline/CellRangerCellHashingHandler.java index c07f6584a..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) @@ -94,7 +95,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 +198,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 +258,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()); @@ -270,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, 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); @@ -304,7 +311,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..52dc9016d --- /dev/null +++ b/tcrdb/src/org/labkey/tcrdb/pipeline/CellRangerSeuratHandler.java @@ -0,0 +1,603 @@ +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), + ToolParameterDescriptor.create("minCountPerCell", "Min Reads/Cell", null, "ldk-integerfield", null, 3) + )); + } + + @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) + { + ctx.getLogger().info("No hashing readset for: " + rs.getReadsetId() + ", this probably indicates either hashing is not used or the hashing data is not available."); + continue; + } + + 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[7])) + { + htosForReadset++; + bcWriter.writeNext(new String[]{line[6], line[5]}); + } + } + } + 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"); + } + } + else + { + ctx.getLogger().info("Cell hashing was not used"); + } + } + + 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("initialCells <- ncol(seuratObj)"); + 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("if (ncol(seuratObj) != initialCells) { stop('Cell count not equal after appending cell hashing calls!') }"); + 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 ca5250b53..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) @@ -105,7 +106,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 @@ -195,19 +196,18 @@ 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); - ctx.getFileManager().addStepOutputs(action, output); - - boolean useCellHashing = utils.useCellHashing(ctx.getSequenceSupport()); - if (useCellHashing) + 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) { - if (cellToHto == null) - { - throw new PipelineJobException("Missing cell to HTO file"); - } + throw new PipelineJobException("Missing cell to HTO file"); + } + + ctx.getFileManager().addStepOutputs(action, output); + } } @@ -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 dfc1a8d1a..a7402e002 100644 --- a/tcrdb/src/org/labkey/tcrdb/pipeline/CellRangerVDJUtils.java +++ b/tcrdb/src/org/labkey/tcrdb/pipeline/CellRangerVDJUtils.java @@ -75,15 +75,15 @@ 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"); - SequenceAnalysisService.get().writeAllCellHashingBarcodes(_sourceDir); + SequenceAnalysisService.get().writeAllCellHashingBarcodes(_sourceDir, job.getUser(), job.getContainer()); Map colMap = QueryService.get().getColumns(cDNAs, PageFlowUtil.set( FieldKey.fromString("rowid"), @@ -93,6 +93,7 @@ public void prepareVDJHashingFilesIfNeeded(PipelineJob job, SequenceAnalysisJobS FieldKey.fromString("sortId/hto"), FieldKey.fromString("sortId/hto/sequence"), FieldKey.fromString("hashingReadsetId"), + FieldKey.fromString("hashingReadsetId/totalFiles"), FieldKey.fromString("status")) ); @@ -101,15 +102,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", "HasHashingReads"}); 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"))); @@ -124,7 +125,8 @@ public void prepareVDJHashingFilesIfNeeded(PipelineJob job, SequenceAnalysisJobS 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,11 @@ 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()) + { + throw new PipelineJobException("There were no readsets found."); } } @@ -172,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"); @@ -208,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()) @@ -217,11 +223,27 @@ public File runRemoteCellHashingTasks(PipelineStepOutput output, String outputCa return null; } - _log.debug("total cached readset/HTO pairs: " + readsetToHashing.size()); + //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 HTO is used, will not use hashing"); + return null; + } + + _log.debug("total cached readset/hashing readset pairs: " + readsetToHashing.size()); + _log.debug("unique HTOs: " + lineCount); //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), ',')) { @@ -235,32 +257,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) @@ -268,11 +294,20 @@ public File runRemoteCellHashingTasks(PipelineStepOutput output, String outputCa throw new PipelineJobException(e); } - //prepare whitelist of barcodes, based on cDNA records - File htoBarcodeWhitelist = getValidHashingBarcodeFile(); - if (!htoBarcodeWhitelist.exists()) + if (uniqueBarcodes.size() < 500 && uniqueBarcodesIncludingNoCDR3.size() > uniqueBarcodes.size()) { - throw new PipelineJobException("Unable to find file: " + htoBarcodeWhitelist.getPath()); + _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())); @@ -283,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, minCountPerCell, sourceDir, editDistance, scanEditDistances, rs, genomeId); if (!hashtagCalls.exists()) { throw new PipelineJobException("Unable to find expected file: " + hashtagCalls.getPath()); @@ -381,7 +416,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; } @@ -433,6 +468,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 +482,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 +499,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 +522,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 +551,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,9 +584,13 @@ else if ("Negative".equals(hto)) { doubletSkipped++; } + else if (discordantBarcodes.contains(barcode)) + { + discordantSkipped++; + } else { - _log.info("skipping cell barcode without HTO call: " + barcode); + //_log.info("skipping cell barcode without HTO call: " + barcode); totalSkipped++; } continue; @@ -591,8 +640,9 @@ 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 because they are doublets: " + doubletSkipped); + _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); @@ -611,8 +661,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 +675,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 +703,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) { @@ -875,103 +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; - } - - 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); - - 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); - - 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 (getCachedReadsetMap(support).isEmpty()) + return false; - 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) + File htoBarcodeWhitelist = getValidHashingBarcodeFile(); + if (!htoBarcodeWhitelist.exists()) { - throw new PipelineJobException(e); + throw new PipelineJobException("Unable to find file: " + htoBarcodeWhitelist.getPath()); } + + return SequencePipelineService.get().getLineCount(htoBarcodeWhitelist) > 1; } public static class CDNA diff --git a/tcrdb/src/org/labkey/tcrdb/pipeline/CellRangerVDJWrapper.java b/tcrdb/src/org/labkey/tcrdb/pipeline/CellRangerVDJWrapper.java index 7a215554d..9bf7d9bcb 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() @@ -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/MiXCRAnalysis.java b/tcrdb/src/org/labkey/tcrdb/pipeline/MiXCRAnalysis.java index f28fd8d26..09b89b5c8 100644 --- a/tcrdb/src/org/labkey/tcrdb/pipeline/MiXCRAnalysis.java +++ b/tcrdb/src/org/labkey/tcrdb/pipeline/MiXCRAnalysis.java @@ -1450,9 +1450,11 @@ private void possiblyAddOrphans(RunData rd, File outDir) throws PipelineJobExcep try (CSVReader reader = new CSVReader(Readers.getReader(possibleNovels), '\t')) { String[] line; + int idx = 0; while ((line = reader.readNext()) != null) { - if (line[0].startsWith("Readset")) + idx++; + if (idx == 1) { continue; } diff --git a/tcrdb/src/org/labkey/tcrdb/pipeline/SeuratCellHashingHandler.java b/tcrdb/src/org/labkey/tcrdb/pipeline/SeuratCellHashingHandler.java index 09dc0b10c..a09d7f91a 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() { @@ -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 (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); @@ -79,7 +80,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