fixup! fixup! WIP

Author: arteymix

README.md

@@ -102,13 +102,15 @@ The output is organized as follow:
 
 ```
 pipeline-output/
-    genomes/<reference_id>/                 # Genomic references
-    references/<reference_id>/              # RSEM/STAR indexes
-    data/<source>                           # FASTQs (note that GEO source uses SRA)
-    data-qc/<experiment_id>/<sample_id>/    # FastQC reports
-    aligned/<reference_id>/<experiment_id>/ # alignments and quantification results
-    quantified/<reference_id>               # quantification matrices for isoforms and genes
-    report/<reference_id>/<experiment_id>/  # MultiQC reports for reads and alignments
+    genomes/<reference_id>/                       # Genomic references
+    references/<reference_id>/                    # RSEM/STAR indexes
+    data/<source>/                                # FASTQs (organization is source-specific; note that GEO source uses SRA)
+    data-qc/<experiment_id>/<sample_id>/          # FastQC reports
+    data-single-cell/<experiment_id>/<sample_id>/ # Single-cell data (hard links to files from data/)
+    aligned/<reference_id>/<experiment_id>/       # alignments and quantification results
+    quantified/<reference_id>                     # quantification matrices for isoforms and genes
+    quantified-single-cell/<reference_id>         # quantified single-cell data (Cell Ranger outputs)
+    report/<reference_id>/<experiment_id>/        # MultiQC reports for reads and alignments
 ```
 
 You can adjust the pipeline output directory by setting `OUTPUT_DIR` under

example.luigi.cfg

@@ -30,6 +30,9 @@ submit_batch_info_jobs=2
 scheduler=local
 scheduler_partition=
 scheduler_extra_args=[]
+# Default tools, override as needed
+#cutadapt_bin=cutadapt
+#cell_ranger_bin=cellranger
 
 #
 # This section contains the necessary variables for the pipeline execution
@@ -40,7 +43,7 @@ scheduler_extra_args=[]
 OUTPUT_DIR=pipeline-output
 GENOMES=genomes
 REFERENCES=references
-REFERENCES_CELL_RANGER=references-cell-ranger
+SINGLE_CELL_REFERENCES=references-single-cell
 METADATA=metadata
 DATA=data
 DATAQCDIR=data-qc
@@ -52,9 +55,6 @@ BATCHINFODIR=batch-info
 # RSEM
 RSEM_DIR=contrib/RSEM
 
-# Cell Ranger
-cell_ranger_bin=cellranger
-
 SLACK_WEBHOOK_URL=
 
 [rnaseq_pipeline.sources.sra]
@@ -72,3 +72,6 @@ appdata_dir=/space/gemmaData
 human_reference_id=hg38_ncbi
 mouse_reference_id=mm10_ncbi
 rat_reference_id=rn7_ncbi
+human_single_cell_reference_id=refdata-gex-GRCh38-2024-A
+mouse_single_cell_reference_id=refdata-gex-GRCm39-2024-A
+rat_single_cell_reference_id=refdata-gex-mRatBN7-2-2024-A

luigi.cfg

@@ -0,0 +1,77 @@
+#
+# This is a configuration example for Luigi and the RNA-Seq pipeline.
+#
+
+#
+# This section contains scheduler resources dedicated to the pipeline
+# execution.
+#
+
+[core]
+autoload_range=true
+
+[resources]
+# in number of available CPUs
+cpus=16
+# in gigabytes
+memory=32
+geo_http_connections=4
+edirect_http_connections=4
+array_express_http_connections=4
+sra_connections=4
+# If you specify the 'slurm' scheduler in Bioluigi, you must set this resource
+slurm_jobs=384
+prefetch_jobs=2
+fastq_dump_jobs=40
+submit_data_jobs=1
+submit_batch_info_jobs=2
+
+[bioluigi]
+scheduler=local
+scheduler_partition=
+scheduler_extra_args=[]
+
+#
+# This section contains the necessary variables for the pipeline execution
+#
+
+[rnaseq_pipeline]
+# pipeline output directories (relative to OUTPUT_DIR)
+OUTPUT_DIR=pipeline-output
+GENOMES=genomes
+REFERENCES=references
+SINGLE_CELL_REFERENCES=references-single-cell
+METADATA=metadata
+DATA=data
+DATAQCDIR=data-qc
+ALIGNDIR=aligned
+ALIGNQCDIR=aligned-qc
+QUANTDIR=quantified
+BATCHINFODIR=batch-info
+
+# RSEM
+RSEM_DIR=contrib/RSEM
+
+# Cell Ranger
+cell_ranger_bin=cellranger
+
+SLACK_WEBHOOK_URL=
+
+[rnaseq_pipeline.sources.sra]
+# location where tools like prefetch and fastq-dump will store downloaded SRA files
+# you can get this value with vdb-config -p
+ncbi_public_dir=/home/guillaume/Projets/Gemma/gemma-data/cache/ncbi/public
+
+[rnaseq_pipeline.gemma]
+cli_bin=gemma-cli
+# values for $JAVA_HOME and $JAVA_OPTS environment variables
+cli_JAVA_HOME=
+cli_JAVA_OPTS=
+baseurl=https://gemma.msl.ubc.ca
+appdata_dir=/space/gemmaData
+human_reference_id=hg38_ncbi
+mouse_reference_id=mm10_ncbi
+rat_reference_id=rn7_ncbi
+human_single_cell_reference_id=hg38_ncbi
+mouse_single_cell_reference_id=mm10_ncbi
+rat_single_cell_reference_id=rn7_ncbi

rnaseq_pipeline/gemma.py

@@ -18,6 +18,9 @@ class gemma(luigi.Config):
     human_reference_id: str = luigi.Parameter()
     mouse_reference_id: str = luigi.Parameter()
     rat_reference_id: str = luigi.Parameter()
+    human_single_cell_reference_id: str = luigi.Parameter()
+    mouse_single_cell_reference_id: str = luigi.Parameter()
+    rat_single_cell_reference_id: str = luigi.Parameter()
 
 cfg = gemma()
 
@@ -98,6 +101,13 @@ def reference_id(self):
         except KeyError:
             raise ValueError('Unsupported Gemma taxon {}.'.format(self.taxon))
 
+    def single_cell_reference_id(self):
+        try:
+            return {'human': cfg.human_single_cell_reference_id, 'mouse': cfg.mouse_single_cell_reference_id, 'rat': cfg.rat_single_cell_reference_id}[
+                self.taxon]
+        except KeyError:
+            raise ValueError('Unsupported Gemma taxon {}.'.format(self.taxon))
+
     @property
     def platform_short_name(self):
         return f'Generic_{self.taxon}_ncbiIds'

rnaseq_pipeline/sources/arrayexpress.py

@@ -50,6 +50,7 @@ class DownloadArrayExpressRun(luigi.Task):
 
     @property
     def platform(self):
+        # TODO: detect platforms from ArrayExpress metadata
         return IlluminaPlatform('HiSeq 2500')
 
     def run(self):

rnaseq_pipeline/targets.py

@@ -85,7 +85,7 @@ def is_stale(self):
     def exists(self):
         return super().exists() and not self.is_stale()
 
-class DownloadRunTarget(luigi.LocalTarget):
+class DownloadRunTarget(luigi.Target):
     run_id: str
     files: list[str]
     layout: list[str]

rnaseq_pipeline/tasks.py

@@ -4,7 +4,6 @@
 import uuid
 from glob import glob
 from os.path import join, dirname
-from typing import Optional
 
 import luigi
 import luigi.task
@@ -456,8 +455,8 @@ class OrganizeSingleCellSample(luigi.Task):
     def run(self):
         runs = self.input()
         os.makedirs(self.output().path)
-        for run in runs:
-            for lane, (f, read_type) in enumerate(zip(run.files, run.layout)):
+        for lane, run in enumerate(runs):
+            for f, read_type in zip(run.files, run.layout):
                 dest = join(self.output().path, f'{self.sample_id}_S1_L{lane + 1:03}_{read_type}_001.fastq.gz')
                 if os.path.exists(dest):
                     os.unlink(dest)
@@ -470,10 +469,6 @@ def output(self):
 class AlignSingleCellSample(DynamicWrapperTask):
     experiment_id: str
     sample_id: str
-    expect_cells: Optional[int] = luigi.OptionalIntParameter(default=None, positional=False)
-    force_cells: Optional[int] = luigi.OptionalIntParameter(default=None, positional=False)
-    chemistry: Optional[str] = luigi.OptionalParameter(default=None, positional=False,
-        description='Chemistry to use for Cell Ranger (default is to auto-detect)')
 
     def run(self):
         fastqs_dir, transcriptome_dir = self.input()
@@ -482,15 +477,13 @@ def run(self):
             transcriptome_dir=transcriptome_dir,
             fastqs_dir=fastqs_dir,
             output_dir=self.output().path,
-            expect_cells=self.expect_cells,
-            force_cells=self.force_cells,
-            chemistry=self.chemistry,
-            # FIXME: request a node with AVX512
+            # TODO: add an avx feature on slurm
             scheduler_extra_args=['--constraint', 'thrd64']
     )
 
     def output(self):
-        return luigi.LocalTarget(join(cfg.OUTPUT_DIR, 'quantified-single-cell', self.experiment_id, self.sample_id))
+        return luigi.LocalTarget(
+            join(cfg.OUTPUT_DIR, 'quantified-single-cell', self.reference_id, self.experiment_id, self.sample_id))
 
 class AlignSingleCellExperiment(DynamicTaskWithOutputMixin, DynamicWrapperTask):
     experiment_id: str = luigi.Parameter()
@@ -591,6 +584,24 @@ def output(self):
         return luigi.LocalTarget(
             join(gemma_cfg.appdata_dir, 'metadata', self.experiment_id, 'MultiQCReports/multiqc_report.html'))
 
+class SubmitSingleCellExperimentDataToGemma(RerunnableTaskMixin, GemmaCliTask):
+    experiment_id: str = luigi.Parameter()
+    subcommand = 'loadSingleCellData'
+
+    def requires(self):
+        return AlignSingleCellExperiment(experiment_id=self.experiment_id,
+                                         reference_id=self.single_cell_reference_id(),
+                                         source='gemma')
+
+    def subcommand_args(self):
+        return ['-e', self.experiment_id, '-a', self.platform_short_name,
+                '--data-path', self.input().path,
+                '--quantitation-type-recomputed-from-raw-data',
+                '--preferred-quantitation-type',
+                # TODO: add sequencing metadata
+                # FIXME: add --replace
+                ]
+
 @requires(SubmitExperimentDataToGemma, SubmitExperimentBatchInfoToGemma, SubmitExperimentReportToGemma)
 class SubmitExperimentToGemma(TaskWithOutputMixin, WrapperTask):
     """

rnaseq_pipeline/webviewer/__init__.py

@@ -1,4 +1,3 @@
-from os.path import basename, getctime, join, dirname
 import datetime
 from glob import glob
 from os.path import basename, getctime, join, dirname

tests/test_webviewer.py

@@ -11,6 +11,7 @@ def test_experiment_summary(client):
     res = client.get('/experiment/GSE87750')
     assert res.status == '200 OK'
 
+@pytest.mark.skip()
 def test_experiment_batch_info(client):
     res = client.get('/experiment/GSE87750/batch-info')
     assert res.status == '200 OK'