Run your own analysis

Configure your data once, then work through the eight checkpoints of the primary module in order. At each checkpoint, run its targets, review its plots and revise its settings until you accept the result. Start with the defaults and revise as the plots suggest: each plot’s subtitle and caption say what to look for, and the output gallery shows an example of every plot.

Run the R commands from the repository root, as in the demo. Replace my_aggregation and my_GEM_well with your own names. <store> is the targets store set in _targets.yaml, normally outputs/ in the repository root. Edit the files in the selected configuration directory; linked parameters belong in your aggregation’s entry in cfg_aggregations.yaml.

Each command also builds the earlier results it depends on. After changing a setting, rerun the checkpoint: targets rebuilds only the results the change affects, and later checkpoints pick up the change when they run.

Configure your data

  1. GEM wells: add one row per GEM well to cfg_GEM_wells.tsv, with GEM_well_QC_exclude_list set to NA for the first run; see GEM well table.
  2. Donors: prepare the donor metadata TSV, with one row for each donor in these GEM wells; see Donor metadata table.
  3. Aggregation: in cfg_aggregations.yaml, copy the template_aggregation entry and rename the copy my_aggregation. Set is_active: true, list your GEM wells in aggregation_GEM_well_IDs, give the path to the donor table in aggregation_donor_id_metadata_tsv, and replace the placeholder aggregation_GEX_marker_genes with markers for the cell types you expect. Optionally, list the donor or GEM well columns to show in the plots in aggregation_categorical_vars and aggregation_continuous_vars. See Aggregation configuration.
  4. Demo: set is_active: false for the immune_human_2x aggregation and GEM_well_is_active to FALSE for its GEM wells, healthy_PBMC_human and controller_PBMC_human. Change both together, because an active aggregation cannot use an inactive GEM well. A later targets::tar_make() without names then builds only your data.

Checkpoint 1: Pre-aggregation QC

Computes QC metrics and donor assignments for each GEM well, applies the GEM well QC filters and compares the GEM wells of the aggregation. Method details

Run pipeline:

R
targets::tar_make(
  names = tidyselect::ends_with("1_pre_aggregation_QC.my_aggregation")
)

The command also builds every per-GEM-well target the aggregation needs: Cell Ranger input conversion, genotyping of pooled wells with cellsnp-lite and Vireo, and per-well QC. After new or re-run sequencing, these dominate its run time, so list the outdated targets first:

R
targets::tar_outdated(
  names = tidyselect::ends_with("1_pre_aggregation_QC.my_aggregation"),
  callr_function = NULL
)

Review plots (examples):

<store>/plots/my_aggregation/1_pre_aggregation_QC/
├── per_aggregation_GEM_well_QC_comparisons/
├── demultiplexing_assignment_bars.png
├── nuclei_per_donor_id_bars.png
├── aggregation_excluded_barcodes_by_type_upset.png
├── aggregation_excluded_cellranger_only_barcodes_by_type_upset.png
└── cell_retention_flow_plot.png

For the exclusion overlaps within one GEM well, run its checkpoint 1 targets:

R
targets::tar_make(
  names = tidyselect::ends_with("1_pre_aggregation_QC.my_GEM_well")
)
<store>/plots/my_GEM_well/1_pre_aggregation_QC/
├── excluded_barcodes_by_type_upset.png
└── excluded_cellranger_only_barcodes_by_type_upset.png

Revise: choose cutoffs from the per_aggregation_GEM_well_QC_comparisons/ distributions and enter them in GEM_well_QC_exclude_list for each GEM well; see Set QC filters after the first run. After the rerun, these plots mark numeric cutoffs, and the UpSet and retention plots show the nuclei each filter removes. Remove a GEM well that fails QC from aggregation_GEM_well_IDs.

Checkpoint 2: GEX dimension reduction

Normalizes the combined GEX data, selects variable genes and computes PCA, with Harmony correction when covariates are configured. Method details

Run pipeline:

R
targets::tar_make(
  names = tidyselect::ends_with("2_GEX_PCA_QC.my_aggregation")
)

Review plots (examples):

<store>/plots/my_aggregation/2_GEX_PCA_QC/
├── variable_feature_plot.png
├── VizDimLoadings_plots/
├── PCA_singular_values_elbow_plot.png
├── PCA_embedding_sdev_plot.png
└── PCA_metadata_association_barplots.png

Revise: choose the PCs used for clustering and the UMAP with aggregation_GEX_data_PCs. The PCA computes at least 50 components, so the plots show components beyond your choice, and a choice within them does not recompute the PCA. When components track a batch variable such as GEM_well_ID, list it in aggregation_harmony_correction_metadata_col_names; Harmony then corrects both modalities, and the embedding and association plots add Harmony panels. Do not correct for GEM_well_ID when the GEM wells hold different tissues; see Combine several tissues.

Checkpoint 3: GEX clusters and cell types

Clusters the GEX data, labels the clusters from your marker genes and removes GEX doublets. Method details

Run pipeline:

R
targets::tar_make(
  names = tidyselect::ends_with("3_GEX_QC.my_aggregation")
)

The command also builds GEX_Seurat_object.3_GEX_QC.my_aggregation, a Seurat export that is costly for large aggregations and outdated whenever the labels change. While you revise the labels, leave it out, and run the command above once you accept them:

R
targets::tar_make(
  names = tidyselect::ends_with("3_GEX_QC.my_aggregation") &
    !tidyselect::starts_with("GEX_Seurat_object")
)

Review plots (examples):

<store>/plots/my_aggregation/3_GEX_QC/
├── UMAPs/
│   ├── categorical/{harmony,non_harmony}/   # non_harmony only with Harmony covariates
│   ├── continuous/{harmony,non_harmony}/
│   └── cross/
├── markers_by_cluster_dot_plot.png
├── markers_by_cell_type_dot_plot.png
├── module_scores_by_cluster_dot_plot.png
├── module_scores_by_cell_type_dot_plot.png
├── cluster_UCell_advantage_plots/
├── cluster_marker_volcano_plots/
├── cell_type_marker_volcano_plots.png
├── categorical_by_cluster_bars_plots/
├── categorical_by_cell_type_bars_plots/
├── continuous_by_cluster_violin_plot/
├── continuous_by_cell_type_violin_plot/
├── resolution_sweep_plot.png
└── cell_retention_flow_plot.png
<store>/files/my_aggregation/3_GEX_QC/
├── cluster_composition_table.tsv
└── resolution_sweep_table.tsv

Revise:

Checkpoint 4: Peak QC

Calls peaks within groups of nuclei, by default the GEX clusters, merges them into one consensus peak set and computes peak-based ATAC QC metrics. Method details

Run pipeline:

R
targets::tar_make(
  names = tidyselect::ends_with("4_peak_QC.my_aggregation")
)

Review plots (examples):

<store>/plots/my_aggregation/4_peak_QC/
├── peaks_QC_violins_plot/
└── peaks_similarity_tiles_plot.png

Revise: choose ATAC filters from the peaks_QC_violins_plot/ distributions and enter them in aggregation_QC_exclude_list_combined_object; see QC filters after peak calling. Checkpoint 5 applies them.

Checkpoint 5: ATAC filtering

Applies the ATAC filters from checkpoint 4 to the GEX-retained nuclei. Method details

Run pipeline:

R
targets::tar_make(
  names = tidyselect::ends_with("5_pre_LSI_QC.my_aggregation")
)

Review plots (examples):

<store>/plots/my_aggregation/5_pre_LSI_QC/
├── QC_excluded_upset_plot.png
└── cell_retention_flow_plot.png

Revise: if the filters remove more nuclei than you can justify, or remove one GEM well disproportionately, revise aggregation_QC_exclude_list_combined_object.

Checkpoint 6: ATAC dimension reduction

Computes LSI embeddings of the filtered peak matrix, with Harmony correction when covariates are configured. Method details

Run pipeline:

R
targets::tar_make(
  names = tidyselect::ends_with("6_ATAC_LSI_QC.my_aggregation")
)

Review plots (examples):

<store>/plots/my_aggregation/6_ATAC_LSI_QC/
├── LSI_singular_values_elbow_plot.png
├── LSI_embedding_sdev_plot.png
├── VizDimLoadings_plots/
└── LSI_metadata_association_barplots.png

Revise: choose the LSI components used for clustering and the UMAP with aggregation_ATAC_data_PCs. As for the PCA, at least 50 are computed, and a choice within them does not recompute the LSI. They start at component 2 by default; check the association plots for components that track sequencing depth or other technical metrics. Add ATAC-only Harmony covariates with aggregation_extra_harmony_covars_ATAC.

Checkpoint 7: ATAC clusters and motifs

Clusters the ATAC data, names the clusters from GEX marker scores, removes ATAC doublets, and computes motif accessibility. Method details

Configure: optionally, list expected transcription factors per cell type in aggregation_ATAC_marker_TFs to organize the motif-accessibility plots. Set aggregation_ATAC_gene_activity_scores to true to also compute gene activity scores and marker_gene_activity_dot_plot.png; they are off by default because they take hours of worker time in large aggregations.

Run pipeline:

R
targets::tar_make(
  names = tidyselect::ends_with("7_ATAC_QC.my_aggregation")
)

Review plots (examples):

<store>/plots/my_aggregation/7_ATAC_QC/
├── UMAPs/{categorical,continuous,cross}/
├── marker_gene_activity_dot_plot.png
├── motif_family_accessibility_by_ATAC_cluster_heatmap.png
├── motif_family_accessibility_by_GEX_cluster_heatmap.png
├── motif_family_accessibility_by_GEX_cell_type_heatmap.png
├── motif_family_accessibility_marker_volcano_plots.png
├── coverage_tracks_plots/
├── cluster_UCell_advantage_plots/
├── confusion_matrices_plots.png
├── categorical_by_cluster_bars_plots/
├── categorical_by_cell_type_bars_plots/
├── continuous_by_cluster_violin_plot/
├── continuous_by_cell_type_violin_plot/
├── resolution_sweep_plot.png
└── cell_retention_flow_plot.png
<store>/files/my_aggregation/7_ATAC_QC/
├── cluster_composition_table.tsv
└── resolution_sweep_table.tsv

Revise: choose aggregation_ATAC_cluster_res from the motif patterns, the agreement with GEX labels in confusion_matrices_plots.png and the cluster counts at half and twice the resolution in resolution_sweep_plot.png. aggregation_ATAC_cluster_annotation_min_advantage sets the assignment margin of the ATAC labels.

Checkpoint 8: WNN integration

Integrates the GEX and ATAC embeddings with weighted nearest neighbors (WNN), clusters and labels the integrated graph, and builds the final multimodal Seurat/Signac object. Method details

Run pipeline:

R
targets::tar_make(
  names = tidyselect::ends_with("8_multimodal_QC.my_aggregation")
)

Review plots (examples):

<store>/plots/my_aggregation/8_multimodal_QC/
├── UMAPs/
│   ├── categorical/
│   ├── continuous/
│   ├── cross/
│   ├── cluster_named_dim_tri_plot.png
│   └── cluster_cell_type_dim_tri_plot.png
├── markers_by_cluster_dot_plot.png
├── module_scores_by_cluster_dot_plot.png
├── cluster_UCell_advantage_plots/
├── confusion_matrices_plots/
├── WNN_weight_metadata_associations_plot/
├── categorical_by_cluster_bars_plots/
├── categorical_by_cell_type_bars_plots/
├── continuous_by_cluster_violin_plot/
├── continuous_by_cell_type_violin_plot/
└── cell_retention_flow_plot.png
<store>/files/my_aggregation/8_multimodal_QC/
└── cluster_composition_table.tsv

Revise: choose aggregation_WNN_cluster_res. To refine the WNN labels without rebuilding the earlier checkpoints, override marker sets in aggregation_WNN_marker_genes or change aggregation_WNN_cluster_annotation_min_advantage. The WNN clusters and cell types are the populations used in downstream summaries and comparisons.

The final object is multimodal_Seurat_object.8_multimodal_QC.my_aggregation; read it as in Inspect the demo results. Like the checkpoint 3 export, it is costly and outdated whenever the labels change: until you accept them, add & !tidyselect::starts_with("multimodal_Seurat_object") to names in the command above.

Choose your next analysis

To build every remaining output, run targets::tar_make() without names. Then continue with an optional module: