Last updated: 2026-01-27

Checks: 6 1

Knit directory: paed-airway-allTissues/

This reproducible R Markdown analysis was created with workflowr (version 1.7.1). The Checks tab describes the reproducibility checks that were applied when the results were created. The Past versions tab lists the development history.


The R Markdown is untracked by Git. To know which version of the R Markdown file created these results, you’ll want to first commit it to the Git repo. If you’re still working on the analysis, you can ignore this warning. When you’re finished, you can run wflow_publish to commit the R Markdown file and build the HTML.

Great job! The global environment was empty. Objects defined in the global environment can affect the analysis in your R Markdown file in unknown ways. For reproduciblity it’s best to always run the code in an empty environment.

The command set.seed(20230811) was run prior to running the code in the R Markdown file. Setting a seed ensures that any results that rely on randomness, e.g. subsampling or permutations, are reproducible.

Great job! Recording the operating system, R version, and package versions is critical for reproducibility.

Nice! There were no cached chunks for this analysis, so you can be confident that you successfully produced the results during this run.

Great job! Using relative paths to the files within your workflowr project makes it easier to run your code on other machines.

Great! You are using Git for version control. Tracking code development and connecting the code version to the results is critical for reproducibility.

The results in this page were generated with repository version 378760b. See the Past versions tab to see a history of the changes made to the R Markdown and HTML files.

Note that you need to be careful to ensure that all relevant files for the analysis have been committed to Git prior to generating the results (you can use wflow_publish or wflow_git_commit). workflowr only checks the R Markdown file, but you know if there are other scripts or data files that it depends on. Below is the status of the Git repository when the results were generated:


Ignored files:
    Ignored:    .DS_Store
    Ignored:    .RData
    Ignored:    .Rhistory
    Ignored:    .Rproj.user/
    Ignored:    analysis/.DS_Store
    Ignored:    analysis/naiveB_top10_markers_whitebg.png
    Ignored:    data/.DS_Store
    Ignored:    data/RDS/
    Ignored:    figure/.DS_Store
    Ignored:    output/.DS_Store
    Ignored:    output/CAM_results/
    Ignored:    output/CSV/.DS_Store
    Ignored:    output/CSV_v2/.DS_Store
    Ignored:    output/CSV_v2/Adenoids/.DS_Store
    Ignored:    output/CSV_v2/BAL/.DS_Store
    Ignored:    output/CSV_v2/Bronchial_brushings/.DS_Store
    Ignored:    output/CSV_v2/Nasal_brushings/.DS_Store
    Ignored:    output/CSV_v2/Tonsils/.DS_Store
    Ignored:    output/DGE/.DS_Store
    Ignored:    output/DGE/RUVseq_TonsilAtlas/.DS_Store
    Ignored:    output/DGE/RUVseq_TonsilAtlas_k1/
    Ignored:    output/DGE/RUVseq_earlyAIR_Tonsils/.DS_Store
    Ignored:    output/DGE/Tonsils/RUV_Tonsils_Double negative T/
    Ignored:    output/RDS/
    Ignored:    output/plots/

Untracked files:
    Untracked:  All_metadata_combinedplots.Rhistory
    Untracked:  analysis/03_Batch_Integration.Rmd
    Untracked:  analysis/Age_proportions.Rmd
    Untracked:  analysis/Age_proportions_AllBatches.Rmd
    Untracked:  analysis/All_Batches_QCExploratory_v2.Rmd
    Untracked:  analysis/All_Umaps_and_marker_plots.Rmd
    Untracked:  analysis/All_metadata.Rmd
    Untracked:  analysis/Annotation_BAL.Rmd
    Untracked:  analysis/Annotation_Bronchial_brushings.Rmd
    Untracked:  analysis/Annotation_Nasal_brushings.Rmd
    Untracked:  analysis/BatchCorrection_Adenoids.Rmd
    Untracked:  analysis/BatchCorrection_Nasal_brushings.Rmd
    Untracked:  analysis/BatchCorrection_Tonsils.Rmd
    Untracked:  analysis/Batch_Integration_&_Downstream_analysis.Rmd
    Untracked:  analysis/Batch_correction_&_Downstream.Rmd
    Untracked:  analysis/Cell_cycle_regression.Rmd
    Untracked:  analysis/Clustering_Tonsils_v2.Rmd
    Untracked:  analysis/Cross-tissue_analysis.Rmd
    Untracked:  analysis/DGE_RUVSeq_Adenoids.Rmd
    Untracked:  analysis/DGE_RUVSeq_Bronchial_brushings.Rmd
    Untracked:  analysis/DGE_RUVSeq_Nasal_brushings.Rmd
    Untracked:  analysis/DGE_earlyAIR_tAtlas_Naïve_Bcells_Adenoids.Rmd
    Untracked:  analysis/DGE_earlyAIR_tAtlas_Naïve_Bcells_activated.Rmd
    Untracked:  analysis/DGE_earlyAIR_tAtlas_Plasma_Bcells.Rmd
    Untracked:  analysis/DGE_earlyAIR_vs_TonsilAtlas.Rmd
    Untracked:  analysis/Explore_census.Rmd
    Untracked:  analysis/Lineage_subclustering_Bcells.Rmd
    Untracked:  analysis/Lineage_subclustering_Myeloid_cells.Rmd
    Untracked:  analysis/Lineage_subclustering_Tcells.Rmd
    Untracked:  analysis/Lineages_subclustering.Rmd
    Untracked:  analysis/Master_metadata.Rmd
    Untracked:  analysis/Pediatric_Vs_Adult_Atlases.Rmd
    Untracked:  analysis/Preprocessing_Batch1_Nasal_brushings.Rmd
    Untracked:  analysis/Preprocessing_Batch2_Tonsils.Rmd
    Untracked:  analysis/Preprocessing_Batch3_Adenoids.Rmd
    Untracked:  analysis/Preprocessing_Batch4_Bronchial_brushings.Rmd
    Untracked:  analysis/Preprocessing_Batch5_Nasal_brushings.Rmd
    Untracked:  analysis/Preprocessing_Batch6_BAL.Rmd
    Untracked:  analysis/Preprocessing_Batch7_Bronchial_brushings.Rmd
    Untracked:  analysis/Preprocessing_Batch8_Adenoids.Rmd
    Untracked:  analysis/Preprocessing_Batch9_Tonsils.Rmd
    Untracked:  analysis/TonsilsVsAdenoids.Rmd
    Untracked:  analysis/cell_cycle_regression.R
    Untracked:  analysis/testing_age_all.Rmd
    Untracked:  code/DGE_utils.R
    Untracked:  code/pseudobulk_analysis.R
    Untracked:  data/Cell_labels_Gunjan_v2/
    Untracked:  data/Cell_labels_Mel/
    Untracked:  data/Cell_labels_Mel_v2/
    Untracked:  data/Cell_labels_Mel_v3/
    Untracked:  data/Cell_labels_modified_Gunjan/
    Untracked:  data/Gene_sets/
    Untracked:  data/Hs.c2.cp.reactome.v7.1.entrez.rds
    Untracked:  data/Raw_feature_bc_matrix/
    Untracked:  data/cell_labels_Mel_v4_Dec2024/
    Untracked:  data/celltypes_Mel_GD_v3.xlsx
    Untracked:  data/celltypes_Mel_GD_v4_no_dups.xlsx
    Untracked:  data/celltypes_Mel_modified.xlsx
    Untracked:  data/celltypes_Mel_v2.csv
    Untracked:  data/celltypes_Mel_v2.xlsx
    Untracked:  data/celltypes_Mel_v2_MN.xlsx
    Untracked:  data/celltypes_for_mel_MN.xlsx
    Untracked:  data/col_palette.xlsx
    Untracked:  data/earlyAIR_sample_sheets_combined.xlsx
    Untracked:  figure/DGE_Tonsils_CD8.Rmd/
    Untracked:  metadata_export.rds
    Untracked:  metadata_export1.rds
    Untracked:  naiveB_top10_markers.png
    Untracked:  naiveB_top10_markers_whitebg.png
    Untracked:  output/CSV/All_tissues.propeller.xlsx
    Untracked:  output/CSV/Bronchial_brushings/
    Untracked:  output/CSV/Bronchial_brushings_Marker_gene_clusters.limmaTrendRNA_snn_res.0.4/
    Untracked:  output/CSV/G000231_Neeland_Adenoids.propeller.xlsx
    Untracked:  output/CSV/G000231_Neeland_Bronchial_brushings.propeller.xlsx
    Untracked:  output/CSV/G000231_Neeland_Nasal_brushings.propeller.xlsx
    Untracked:  output/CSV/G000231_Neeland_Tonsils.propeller.xlsx
    Untracked:  output/CSV/Nasal_brushings/
    Untracked:  output/CSV_v3/B_cell_lineage/
    Untracked:  output/DGE/Adenoids/
    Untracked:  output/DGE/Adenoids_B_memory/
    Untracked:  output/DGE/Adenoids_CD4/
    Untracked:  output/DGE/Adenoids_CD8/
    Untracked:  output/DGE/Adenoids_Cycling_GCB/
    Untracked:  output/DGE/Adenoids_IFN/
    Untracked:  output/DGE/Bronchial_brushings/
    Untracked:  output/DGE/DGE_Age_summary_BAL.xlsx
    Untracked:  output/DGE/Nasal_brushings/
    Untracked:  output/DGE/RUV/
    Untracked:  output/DGE/RUVseq_earlyAIR_Adenoids/
    Untracked:  output/DGE/RUVseq_earlyAIR_Bronchial_brushings/
    Untracked:  output/DGE/RUVseq_earlyAIR_Nasal_brushings/
    Untracked:  output/DGE/TonsilAtlas/
    Untracked:  output/DGE/Tonsils/DGE_Age_summary_Tonsils.xlsx
    Untracked:  output/DGE/Tonsils/RUVTonsils_Memory B cells/
    Untracked:  output/DGE/Tonsils/RUV_Tonsils_Memory B cells/
    Untracked:  output/DGE/Tonsils/Tonsils_Naïve B cells/
    Untracked:  output/DGE/Tonsils_CD4/
    Untracked:  output/DGE/Tonsils_CD8/
    Untracked:  output/DGE/Tonsils_IFN/
    Untracked:  paed_sub_L1_cells.rds

Unstaged changes:
    Deleted:    02_QC_exploratoryPlots.Rmd
    Deleted:    02_QC_exploratoryPlots.html
    Modified:   analysis/00_AllBatches_overview.Rmd
    Modified:   analysis/01_QC_emptyDrops.Rmd
    Modified:   analysis/02_QC_exploratoryPlots.Rmd
    Modified:   analysis/Adenoids.Rmd
    Modified:   analysis/Adenoids_v2.Rmd
    Modified:   analysis/Age_modeling.Rmd
    Modified:   analysis/Age_modelling_Adenoids.Rmd
    Modified:   analysis/Age_modelling_Tonsils.Rmd
    Modified:   analysis/AllBatches_QCExploratory.Rmd
    Modified:   analysis/BAL.Rmd
    Modified:   analysis/BAL_v2.Rmd
    Modified:   analysis/Bronchial_brushings.Rmd
    Modified:   analysis/Bronchial_brushings_v2.Rmd
    Modified:   analysis/DGE_TonsilAtlas.Rmd
    Modified:   analysis/DGE_Tonsils_CD8.Rmd
    Modified:   analysis/DGE_analysis_George.Rmd
    Modified:   analysis/DGE_earlyAIR_tAtlas_CD4_TCM.Rmd
    Modified:   analysis/DGE_earlyAIR_tAtlas_CD4_TFH.Rmd
    Modified:   analysis/DGE_earlyAIR_tAtlas_CD4_TN.Rmd
    Modified:   analysis/DGE_earlyAIR_tAtlas_CD8_TF.Rmd
    Modified:   analysis/DGE_earlyAIR_tAtlas_CD8_TN.Rmd
    Modified:   analysis/DGE_earlyAIR_tAtlas_DZtoLZ_GCB.Rmd
    Modified:   analysis/DGE_earlyAIR_tAtlas_Memory_Bcells.Rmd
    Modified:   analysis/DGE_earlyAIR_tAtlas_Naïve_Bcells.Rmd
    Modified:   analysis/Nasal_brushings.Rmd
    Modified:   analysis/Nasal_brushings_v2.Rmd
    Modified:   analysis/Subclustering_Adenoids.Rmd
    Modified:   analysis/Subclustering_BAL.Rmd
    Modified:   analysis/Subclustering_Bronchial_brushings.Rmd
    Modified:   analysis/Subclustering_Nasal_brushings.Rmd
    Modified:   analysis/Subclustering_Tonsils.Rmd
    Modified:   analysis/Tonsils.Rmd
    Modified:   analysis/Tonsils_v2.Rmd
    Modified:   analysis/index.Rmd
    Modified:   output/CSV/BAL_Marker_gene_clusters.limmaTrendRNA_snn_res.0.4/REACTOME-cluster-limma-c0.csv
    Modified:   output/CSV/BAL_Marker_gene_clusters.limmaTrendRNA_snn_res.0.4/REACTOME-cluster-limma-c1.csv
    Modified:   output/CSV/BAL_Marker_gene_clusters.limmaTrendRNA_snn_res.0.4/REACTOME-cluster-limma-c10.csv
    Modified:   output/CSV/BAL_Marker_gene_clusters.limmaTrendRNA_snn_res.0.4/REACTOME-cluster-limma-c11.csv
    Modified:   output/CSV/BAL_Marker_gene_clusters.limmaTrendRNA_snn_res.0.4/REACTOME-cluster-limma-c12.csv
    Modified:   output/CSV/BAL_Marker_gene_clusters.limmaTrendRNA_snn_res.0.4/REACTOME-cluster-limma-c13.csv
    Modified:   output/CSV/BAL_Marker_gene_clusters.limmaTrendRNA_snn_res.0.4/REACTOME-cluster-limma-c14.csv
    Modified:   output/CSV/BAL_Marker_gene_clusters.limmaTrendRNA_snn_res.0.4/REACTOME-cluster-limma-c15.csv
    Modified:   output/CSV/BAL_Marker_gene_clusters.limmaTrendRNA_snn_res.0.4/REACTOME-cluster-limma-c16.csv
    Modified:   output/CSV/BAL_Marker_gene_clusters.limmaTrendRNA_snn_res.0.4/REACTOME-cluster-limma-c17.csv
    Modified:   output/CSV/BAL_Marker_gene_clusters.limmaTrendRNA_snn_res.0.4/REACTOME-cluster-limma-c2.csv
    Modified:   output/CSV/BAL_Marker_gene_clusters.limmaTrendRNA_snn_res.0.4/REACTOME-cluster-limma-c3.csv
    Modified:   output/CSV/BAL_Marker_gene_clusters.limmaTrendRNA_snn_res.0.4/REACTOME-cluster-limma-c4.csv
    Modified:   output/CSV/BAL_Marker_gene_clusters.limmaTrendRNA_snn_res.0.4/REACTOME-cluster-limma-c5.csv
    Modified:   output/CSV/BAL_Marker_gene_clusters.limmaTrendRNA_snn_res.0.4/REACTOME-cluster-limma-c6.csv
    Modified:   output/CSV/BAL_Marker_gene_clusters.limmaTrendRNA_snn_res.0.4/REACTOME-cluster-limma-c7.csv
    Modified:   output/CSV/BAL_Marker_gene_clusters.limmaTrendRNA_snn_res.0.4/REACTOME-cluster-limma-c8.csv
    Modified:   output/CSV/BAL_Marker_gene_clusters.limmaTrendRNA_snn_res.0.4/REACTOME-cluster-limma-c9.csv
    Modified:   output/CSV/BAL_Marker_gene_clusters.limmaTrendRNA_snn_res.0.4/up-cluster-limma-c0.csv
    Modified:   output/CSV/BAL_Marker_gene_clusters.limmaTrendRNA_snn_res.0.4/up-cluster-limma-c1.csv
    Modified:   output/CSV/BAL_Marker_gene_clusters.limmaTrendRNA_snn_res.0.4/up-cluster-limma-c10.csv
    Modified:   output/CSV/BAL_Marker_gene_clusters.limmaTrendRNA_snn_res.0.4/up-cluster-limma-c11.csv
    Modified:   output/CSV/BAL_Marker_gene_clusters.limmaTrendRNA_snn_res.0.4/up-cluster-limma-c12.csv
    Modified:   output/CSV/BAL_Marker_gene_clusters.limmaTrendRNA_snn_res.0.4/up-cluster-limma-c13.csv
    Modified:   output/CSV/BAL_Marker_gene_clusters.limmaTrendRNA_snn_res.0.4/up-cluster-limma-c14.csv
    Modified:   output/CSV/BAL_Marker_gene_clusters.limmaTrendRNA_snn_res.0.4/up-cluster-limma-c15.csv
    Modified:   output/CSV/BAL_Marker_gene_clusters.limmaTrendRNA_snn_res.0.4/up-cluster-limma-c16.csv
    Modified:   output/CSV/BAL_Marker_gene_clusters.limmaTrendRNA_snn_res.0.4/up-cluster-limma-c17.csv
    Modified:   output/CSV/BAL_Marker_gene_clusters.limmaTrendRNA_snn_res.0.4/up-cluster-limma-c2.csv
    Modified:   output/CSV/BAL_Marker_gene_clusters.limmaTrendRNA_snn_res.0.4/up-cluster-limma-c3.csv
    Modified:   output/CSV/BAL_Marker_gene_clusters.limmaTrendRNA_snn_res.0.4/up-cluster-limma-c4.csv
    Modified:   output/CSV/BAL_Marker_gene_clusters.limmaTrendRNA_snn_res.0.4/up-cluster-limma-c5.csv
    Modified:   output/CSV/BAL_Marker_gene_clusters.limmaTrendRNA_snn_res.0.4/up-cluster-limma-c6.csv
    Modified:   output/CSV/BAL_Marker_gene_clusters.limmaTrendRNA_snn_res.0.4/up-cluster-limma-c7.csv
    Modified:   output/CSV/BAL_Marker_gene_clusters.limmaTrendRNA_snn_res.0.4/up-cluster-limma-c8.csv
    Modified:   output/CSV/BAL_Marker_gene_clusters.limmaTrendRNA_snn_res.0.4/up-cluster-limma-c9.csv
    Modified:   output/CSV_v2/G000231_Neeland_Nasal_brushings.propeller.xlsx
    Modified:   output/DGE/RUVseq_TonsilAtlas/DGE_RUVSeq_Age_TonsilAtlas.xlsx
    Modified:   output/DGE/RUVseq_earlyAIR_Tonsils/CD4 TCM/all_res.csv
    Modified:   output/DGE/RUVseq_earlyAIR_Tonsils/CD4 TCM/down_res.csv
    Modified:   output/DGE/RUVseq_earlyAIR_Tonsils/CD4 TCM/up_res.csv
    Modified:   output/DGE/RUVseq_earlyAIR_Tonsils/CD4 TFH/all_res.csv
    Modified:   output/DGE/RUVseq_earlyAIR_Tonsils/CD4 TFH/down_res.csv
    Modified:   output/DGE/RUVseq_earlyAIR_Tonsils/CD4 TFH/up_res.csv
    Modified:   output/DGE/RUVseq_earlyAIR_Tonsils/CD4 TN/all_res.csv
    Modified:   output/DGE/RUVseq_earlyAIR_Tonsils/CD4 TN/down_res.csv
    Modified:   output/DGE/RUVseq_earlyAIR_Tonsils/CD4 TN/up_res.csv
    Modified:   output/DGE/RUVseq_earlyAIR_Tonsils/CD4 Treg-eff/all_res.csv
    Modified:   output/DGE/RUVseq_earlyAIR_Tonsils/CD4 Treg-eff/down_res.csv
    Modified:   output/DGE/RUVseq_earlyAIR_Tonsils/CD4 Treg-eff/up_res.csv
    Modified:   output/DGE/RUVseq_earlyAIR_Tonsils/CD4 Treg/all_res.csv
    Modified:   output/DGE/RUVseq_earlyAIR_Tonsils/CD4 Treg/down_res.csv
    Modified:   output/DGE/RUVseq_earlyAIR_Tonsils/CD4 Treg/up_res.csv
    Modified:   output/DGE/RUVseq_earlyAIR_Tonsils/CD4 effector/all_res.csv
    Modified:   output/DGE/RUVseq_earlyAIR_Tonsils/CD4 effector/down_res.csv
    Modified:   output/DGE/RUVseq_earlyAIR_Tonsils/CD4 effector/up_res.csv
    Modified:   output/DGE/RUVseq_earlyAIR_Tonsils/CD8 TF/all_res.csv
    Modified:   output/DGE/RUVseq_earlyAIR_Tonsils/CD8 TF/down_res.csv
    Modified:   output/DGE/RUVseq_earlyAIR_Tonsils/CD8 TF/up_res.csv
    Modified:   output/DGE/RUVseq_earlyAIR_Tonsils/CD8 TN/all_res.csv
    Modified:   output/DGE/RUVseq_earlyAIR_Tonsils/CD8 TN/down_res.csv
    Modified:   output/DGE/RUVseq_earlyAIR_Tonsils/CD8 TN/up_res.csv
    Modified:   output/DGE/RUVseq_earlyAIR_Tonsils/Cycling GCB/all_res.csv
    Modified:   output/DGE/RUVseq_earlyAIR_Tonsils/Cycling GCB/down_res.csv
    Modified:   output/DGE/RUVseq_earlyAIR_Tonsils/Cycling GCB/up_res.csv
    Modified:   output/DGE/RUVseq_earlyAIR_Tonsils/DGE_RUVSeq_Age_Tonsils.xlsx
    Modified:   output/DGE/RUVseq_earlyAIR_Tonsils/DGE_RUVSeq_Tonsils_all_results.rds
    Modified:   output/DGE/RUVseq_earlyAIR_Tonsils/DZ G2Mphase/all_res.csv
    Modified:   output/DGE/RUVseq_earlyAIR_Tonsils/DZ G2Mphase/down_res.csv
    Modified:   output/DGE/RUVseq_earlyAIR_Tonsils/DZ G2Mphase/up_res.csv
    Modified:   output/DGE/RUVseq_earlyAIR_Tonsils/DZ GCB/all_res.csv
    Modified:   output/DGE/RUVseq_earlyAIR_Tonsils/DZ GCB/down_res.csv
    Modified:   output/DGE/RUVseq_earlyAIR_Tonsils/DZ GCB/up_res.csv
    Modified:   output/DGE/RUVseq_earlyAIR_Tonsils/DZ early  Sphase/all_res.csv
    Modified:   output/DGE/RUVseq_earlyAIR_Tonsils/DZ early  Sphase/down_res.csv
    Modified:   output/DGE/RUVseq_earlyAIR_Tonsils/DZ early  Sphase/up_res.csv
    Modified:   output/DGE/RUVseq_earlyAIR_Tonsils/DZ late Sphase/all_res.csv
    Modified:   output/DGE/RUVseq_earlyAIR_Tonsils/DZ late Sphase/down_res.csv
    Modified:   output/DGE/RUVseq_earlyAIR_Tonsils/DZ late Sphase/up_res.csv
    Modified:   output/DGE/RUVseq_earlyAIR_Tonsils/DZtoLZ GCB transition/all_res.csv
    Modified:   output/DGE/RUVseq_earlyAIR_Tonsils/DZtoLZ GCB transition/down_res.csv
    Modified:   output/DGE/RUVseq_earlyAIR_Tonsils/DZtoLZ GCB transition/up_res.csv
    Modified:   output/DGE/RUVseq_earlyAIR_Tonsils/Double negative T/all_res.csv
    Modified:   output/DGE/RUVseq_earlyAIR_Tonsils/Double negative T/down_res.csv
    Modified:   output/DGE/RUVseq_earlyAIR_Tonsils/Double negative T/up_res.csv
    Modified:   output/DGE/RUVseq_earlyAIR_Tonsils/Early GC-committed NBC/all_res.csv
    Modified:   output/DGE/RUVseq_earlyAIR_Tonsils/Early GC-committed NBC/down_res.csv
    Modified:   output/DGE/RUVseq_earlyAIR_Tonsils/Early GC-committed NBC/up_res.csv
    Modified:   output/DGE/RUVseq_earlyAIR_Tonsils/Early MBC/all_res.csv
    Modified:   output/DGE/RUVseq_earlyAIR_Tonsils/Early MBC/down_res.csv
    Modified:   output/DGE/RUVseq_earlyAIR_Tonsils/Early MBC/up_res.csv
    Modified:   output/DGE/RUVseq_earlyAIR_Tonsils/Early PC precursor/all_res.csv
    Modified:   output/DGE/RUVseq_earlyAIR_Tonsils/Early PC precursor/down_res.csv
    Modified:   output/DGE/RUVseq_earlyAIR_Tonsils/Early PC precursor/up_res.csv
    Modified:   output/DGE/RUVseq_earlyAIR_Tonsils/Follicular dendritic cells/all_res.csv
    Modified:   output/DGE/RUVseq_earlyAIR_Tonsils/Follicular dendritic cells/down_res.csv
    Modified:   output/DGE/RUVseq_earlyAIR_Tonsils/Follicular dendritic cells/up_res.csv
    Modified:   output/DGE/RUVseq_earlyAIR_Tonsils/GC-commited metabolic activation/all_res.csv
    Modified:   output/DGE/RUVseq_earlyAIR_Tonsils/GC-commited metabolic activation/down_res.csv
    Modified:   output/DGE/RUVseq_earlyAIR_Tonsils/GC-commited metabolic activation/up_res.csv
    Modified:   output/DGE/RUVseq_earlyAIR_Tonsils/Gamma delta T/all_res.csv
    Modified:   output/DGE/RUVseq_earlyAIR_Tonsils/Gamma delta T/down_res.csv
    Modified:   output/DGE/RUVseq_earlyAIR_Tonsils/Gamma delta T/up_res.csv
    Modified:   output/DGE/RUVseq_earlyAIR_Tonsils/Memory B cells/all_res.csv
    Modified:   output/DGE/RUVseq_earlyAIR_Tonsils/Memory B cells/down_res.csv
    Modified:   output/DGE/RUVseq_earlyAIR_Tonsils/Memory B cells/up_res.csv
    Modified:   output/DGE/RUVseq_earlyAIR_Tonsils/Monocytes_macrophages/all_res.csv
    Modified:   output/DGE/RUVseq_earlyAIR_Tonsils/Monocytes_macrophages/down_res.csv
    Modified:   output/DGE/RUVseq_earlyAIR_Tonsils/Monocytes_macrophages/up_res.csv
    Modified:   output/DGE/RUVseq_earlyAIR_Tonsils/Naïve B cell-IFN/all_res.csv
    Modified:   output/DGE/RUVseq_earlyAIR_Tonsils/Naïve B cell-IFN/down_res.csv
    Modified:   output/DGE/RUVseq_earlyAIR_Tonsils/Naïve B cell-IFN/up_res.csv
    Modified:   output/DGE/RUVseq_earlyAIR_Tonsils/Naïve B cells activated/all_res.csv
    Modified:   output/DGE/RUVseq_earlyAIR_Tonsils/Naïve B cells activated/down_res.csv
    Modified:   output/DGE/RUVseq_earlyAIR_Tonsils/Naïve B cells activated/up_res.csv
    Modified:   output/DGE/RUVseq_earlyAIR_Tonsils/Naïve B cells/all_res.csv
    Modified:   output/DGE/RUVseq_earlyAIR_Tonsils/Naïve B cells/down_res.csv
    Modified:   output/DGE/RUVseq_earlyAIR_Tonsils/Naïve B cells/up_res.csv
    Modified:   output/DGE/RUVseq_earlyAIR_Tonsils/Plasma B cells/all_res.csv
    Modified:   output/DGE/RUVseq_earlyAIR_Tonsils/Plasma B cells/down_res.csv
    Modified:   output/DGE/RUVseq_earlyAIR_Tonsils/Plasma B cells/up_res.csv
    Modified:   output/DGE/RUVseq_earlyAIR_Tonsils/Pre-BCRi II/all_res.csv
    Modified:   output/DGE/RUVseq_earlyAIR_Tonsils/Pre-BCRi II/down_res.csv
    Modified:   output/DGE/RUVseq_earlyAIR_Tonsils/Pre-BCRi II/up_res.csv
    Modified:   output/DGE/RUVseq_earlyAIR_Tonsils/T-IFN/all_res.csv
    Modified:   output/DGE/RUVseq_earlyAIR_Tonsils/T-IFN/down_res.csv
    Modified:   output/DGE/RUVseq_earlyAIR_Tonsils/T-IFN/up_res.csv
    Modified:   output/DGE/RUVseq_earlyAIR_Tonsils/TFH-LZ-GC/all_res.csv
    Modified:   output/DGE/RUVseq_earlyAIR_Tonsils/TFH-LZ-GC/down_res.csv
    Modified:   output/DGE/RUVseq_earlyAIR_Tonsils/TFH-LZ-GC/up_res.csv
    Deleted:    test

Note that any generated files, e.g. HTML, png, CSS, etc., are not included in this status report because it is ok for generated content to have uncommitted changes.


There are no past versions. Publish this analysis with wflow_publish() to start tracking its development.


Load libraries

suppressPackageStartupMessages({
  library(BiocStyle)
  library(tidyverse)
  library(here)
  library(dplyr)
  library(Seurat)
  library(clustree)
  library(paletteer)
  library(viridis)
  library(ggforce)
  library(ggridges)
  library(kableExtra)
  library(RColorBrewer)
  library(data.table)
  library(dplyr)
  library(cowplot)
  library(ggplot2)
  library(paletteer)
  library(patchwork)
  library(harmony)
  library(BiocParallel)
  library(circlize)
  library(presto)
  library(gtools)
})

Myeloid cells

paed_sub <- readRDS(here("output/RDS/Lineages_RDS_combined/SEU_Myeloid_Cells.rds"))

if (!"umap.sub" %in% names(paed_sub@reductions)) {
  
  paed_sub <- paed_sub |>
    FindVariableFeatures() |>
    ScaleData() |>
    RunPCA()
  
  gc()
  
paed_sub[["RNA"]]
paed_sub[["RNA"]] <- JoinLayers(paed_sub[["RNA"]])
gc()

  paed_sub <- RunUMAP(paed_sub, dims = 1:30, reduction = "pca",
                      reduction.name = "umap.sub")
  
  paed_sub <- paed_sub |>
    FindNeighbors(dims = 1:30, reduction = "pca") |>
    FindClusters(resolution = 0.1)
  gc()
  
  saveRDS(paed_sub, here("output/RDS/Lineages_RDS_combined/SEU_Myeloid_Cells.rds")) 
  } else {
  message("skipping processing")
}
skipping processing
paed_sub$tissue <- factor(
  paed_sub$tissue,
  levels = c("Tonsils", "Adenoids", "Nasal_brushings",
             "Bronchial_brushings", "BAL")
)

Previous labels

Excluding counts less than 100 to avoid too many labels in the plot

#sort(table(paed_sub$cell_labels_v2), decreasing = T)

k <- names(which(table(paed_sub$cell_labels_v2) > 50))
paed_sub$cell_labels_v22 <- ifelse(
    as.character(paed_sub$cell_labels_v2) %in% k,
    as.character(paed_sub$cell_labels_v2),
    NA
)

DimPlot(
    paed_sub,
    group.by = "cell_labels_v22",
    reduction = "umap.sub",
    label = TRUE,
    repel = TRUE,
    label.size = 4.5,
    shuffle = TRUE, raster = F)

#table( paed_sub$cell_labels_v2,paed_sub$RNA_snn_res.0.5)
opt_res <- "RNA_snn_res.0.1"  
n <- nlevels(paed_sub$RNA_snn_res.0.1)
paed_sub$RNA_snn_res.0.1 <- factor(paed_sub$RNA_snn_res.0.1, levels = seq(0,n-1))
paed_sub$seurat_clusters <- NULL
paed_sub$cluster <- paed_sub$RNA_snn_res.0.1
Idents(paed_sub) <- paed_sub$cluster
DimPlot(paed_sub, reduction = "umap.sub", raster =F, repel=T, label = T)

DimPlot(paed_sub, reduction = "umap.sub", raster =F, repel=T, label = T, split.by = "tissue")

Adding CellTypist info

files  <- c(
  paeds  = here("../EarlyAir_paper/output/CellTypist_results/cellTtypist_Myeloid_Paed_Covid19.csv"),
  immune = here("../EarlyAir_paper/output/CellTypist_results/cellTtypist_Myeloid_Immune_low.csv"),
  lung_airway   = here("../EarlyAir_paper/output/CellTypist_results/cellTtypist_Myeloid_Lung_Airway.csv"),
   lung_atlas   = here("../EarlyAir_paper/output/CellTypist_results/cellTtypist_Myeloid_Lung_atlas.csv"),
  tonsil = here("../EarlyAir_paper/output/CellTypist_results/cellTtypist_Myeloid_Human_Tonsil.csv")
)

for (nm in names(files)) {
  df <- read.csv(files[nm])
  paed_sub[[paste0("celltypist_", nm)]] <- df$predicted_labels
}
plot_celltype_heatmap <- function(
  seu,
  cluster_col = "cluster",
  celltype_col,
  min_n = 100,
  title,
  scale = "row",
  fontsize = 10
) {
  seu@meta.data %>%
    dplyr::count(.data[[cluster_col]], .data[[celltype_col]]) %>%
    dplyr::filter(n >= min_n) %>%
    tidyr::pivot_wider(
      names_from = .data[[celltype_col]],
      values_from = n,
      values_fill = 0
    ) %>%
    as.data.frame() %>%
    tibble::column_to_rownames(cluster_col) %>%
    t() %>%
    pheatmap::pheatmap(
      scale = scale,
      main = title,
      fontsize = fontsize
    )
}

plot_celltype_heatmap(paed_sub,celltype_col = "celltypist_paeds",title = "CellTypist Predictions (Paed Covid Study) by Cluster")

plot_celltype_heatmap(paed_sub,celltype_col = "celltypist_immune",title = "CellTypist Predictions (Immune Low) by Cluster")

plot_celltype_heatmap(paed_sub,celltype_col = "celltypist_lung_airway",title = "CellTypist Predictions (Lung Atlas) by Cluster")

plot_celltype_heatmap( paed_sub,celltype_col = "celltypist_lung_atlas",title = "CellTypist Predictions (Lung Atlas) by Cluster")

plot_celltype_heatmap(paed_sub,celltype_col = "celltypist_tonsil", title = "Azimuth Predictions (Tonsil Atlas) by Cluster")

plot_ct <- function(obj, col) {
  k <- names(which(table(obj[[col]][,1]) > 100))
  DimPlot(
    obj,
    group.by = col,
    reduction = "umap.sub",
    cells = colnames(obj)[obj[[col]][,1] %in% k],
    label = TRUE, repel = TRUE, label.size = 2.5,
    shuffle = TRUE, raster = FALSE
  )
}

plot_ct(paed_sub, "celltypist_paeds") + ggtitle("Paed Covid Study")

plot_ct(paed_sub, "celltypist_immune") + ggtitle("Immune low dataset")

plot_ct(paed_sub, "celltypist_lung_airway") + ggtitle("Lung Airway Study")

plot_ct(paed_sub, "celltypist_lung_atlas") + ggtitle("Lung Atlas")

plot_ct(paed_sub, "celltypist_tonsil") + ggtitle("Tonsil Atlas")

plot_ct_bar <- function(obj, col) {
  k <- names(which(table(obj[[col]][,1]) > 100))
  
  df <- obj@meta.data %>%
    filter(get(col) %in% k) %>%   
    count(cluster, !!sym(col)) %>%
    group_by(cluster) %>%
    mutate(frac = n / sum(n)) %>%
    ungroup()
  
  ggplot(df, aes(x = cluster, y = frac, fill = !!sym(col))) +
    geom_col(position = "stack") +
    scale_y_continuous(labels = scales::percent) +
    labs(x = "Cluster", y = "Fraction of cells", fill = col) +
    theme_minimal() +
    theme(axis.text.x = element_text(angle = 45, hjust = 1))
}

plot_ct_bar(paed_sub, "celltypist_paeds") + ggtitle("Paed Covid Study")

plot_ct_bar(paed_sub, "celltypist_immune") + ggtitle("Immune low dataset")

plot_ct_bar(paed_sub, "celltypist_lung_airway") + ggtitle("Lung Airway study")

plot_ct_bar(paed_sub, "celltypist_lung_atlas") + ggtitle("Lung Atlas")

plot_ct_bar(paed_sub, "celltypist_tonsil") + ggtitle("Tonsil Atlas")

features <- list(
  Myeloid_core = c("LYZ", "CTSD", "LST1"),
  Mono_classical = c("CD14", "S100A8", "S100A9", "LGALS3"),
  Mono_nonclassical = c("FCGR3A", "LILRB1", "MS4A7"),
  Macrophage_core = c("CD68", "CSF1R", "APOE", "CTSB"),
  Macro_M2_like = c("CD163", "MRC1", "MSR1"),
  Macro_CCL_inflam = c("CCL2", "CCL3", "CCL4", "IL1B"),
  Macro_prolif = c("MKI67", "TOP2A", "PCNA"),
  cDC2 = c("CD1C", "FCER1A", "CLEC10A"),
  pDC = c("IL3RA", "GZMB", "TCF4", "IRF7"),
  Neutrophil = c("S100A8", "S100A9", "FCGR3B", "CSF3R"),
  Mast = c("TPSAB1", "CPA3", "KIT"),
  FDC = c("CXCL13", "CR2", "FDCSP")
)

for (g in features) {
  print(
    FeaturePlot(
      paed_sub,
      features  = g,
      reduction = "umap.sub",
      raster    = FALSE,
      label     = TRUE
      #split.by  = "tissue"
    ) #
  )
}
Warning: The `slot` argument of `FetchData()` is deprecated as of SeuratObject 5.0.0.
ℹ Please use the `layer` argument instead.
ℹ The deprecated feature was likely used in the Seurat package.
  Please report the issue at <https://github.com/satijalab/seurat/issues>.
This warning is displayed once every 8 hours.
Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
generated.

Warning: The following requested variables were not found: TPSAB1

for (g in features) {
  print(
    FeaturePlot(
      paed_sub,
      features  = g,
      reduction = "umap.sub",
      raster    = FALSE,
      label     = TRUE,
      split.by  = "tissue",
      min.cutoff = "q10",
      max.cutoff = "q90"
    )
  )
}

Warning: The following requested variables were not found: TPSAB1

paed_sub <- FindSubCluster(paed_sub, cluster = 1, graph.name = "RNA_snn", resolution = 0.2, subcluster.name = "M_sub")
Modularity Optimizer version 1.3.0 by Ludo Waltman and Nees Jan van Eck

Number of nodes: 15619
Number of edges: 500867

Running Louvain algorithm...
Maximum modularity in 10 random starts: 0.9065
Number of communities: 7
Elapsed time: 1 seconds
Idents(paed_sub) <- paed_sub$M_sub
#levels(Idents(paed_sub)) <- mixedsort(levels(Idents(paed_sub)))
#levels(paed_sub$T_sub) <- mixedsort(levels(paed_sub$T_sub))
DimPlot(paed_sub, reduction = "umap.sub", raster =F, repel=T, label = T, group.by = "M_sub")

new_levels <- c(
  "0" = "Macrophages",
  "1_0" = "Monocytes",
    "1_1" = "Macrophages",
    "1_2" = "Monocytes",
    "1_3" = "DC2",
    "1_4" = "Monocytes",
    "1_5" = "Activated DCs",
    "1_6" = "Monocytes",
  "2" = "Macrophages",
  "3" = "Macrophages",
  "4" = "Neutrophils",
  "5" = "Plasmacytoid DCs",
  "6" = "Macrophages",
  "7" = "Macrophages",
  "8" = "Mast cells",
  "9" = "Macrophages",
  "10" = "Follicular DCs",
  "11" = "DC1",
  "12" = "Ionocytes",
  "13" = "Melanocyte"
)
paed_sub <- RenameIdents(paed_sub, new_levels)
paed_sub$cell_labels <- NULL

Correcting previous labels and defining level 1

paed_sub$cell_labels_l1 <- Idents(paed_sub)
paed_sub$cell_labels_l1 <- factor(
  Idents(paed_sub),
  levels = c(
    "Monocytes",
    "Macrophages",
    "Neutrophils",
    "DC1",
    "DC2",
    "Activated DCs",
    "Plasmacytoid DCs",
    "Follicular DCs",
    "Mast cells",
    "Ionocytes",
    "Melanocyte"
  )
)

Idents(paed_sub) <- paed_sub$cell_labels_l1
DimPlot(
  paed_sub,
  reduction = "umap.sub",
  label = TRUE,
  repel = TRUE,
  label.size = 2.5, raster = F
)

Macrophages subclustering

macro <-  subset(paed_sub, idents = "Macrophages")
keep_cells <- paed_sub$cell_labels_l1 %in% "Macrophages"
macro@assays$RNA@cells@.Data <- paed_sub@assays$RNA@cells@.Data[keep_cells, ]
macro
An object of class Seurat 
18076 features across 50832 samples within 1 assay 
Active assay: RNA (18076 features, 2000 variable features)
 3 layers present: data, counts, scale.data
 4 dimensional reductions calculated: pca, umap.unintegrated, umap.l1, umap.sub
out <- here("output/RDS/Lineages_RDS_combined/SEU_Macrophages_Cells.rds")

if (!file.exists(out)) {
  macro <- macro %>%
  FindVariableFeatures() %>%
  ScaleData() %>%
  RunPCA() 
  
  gc() 
  macro <- RunUMAP(macro, dims = 1:30, reduction = "pca", reduction.name = "umap.mm")

  meta_data_columns <- colnames(macro@meta.data)
  columns_to_remove <- grep("^RNA_snn_res", meta_data_columns, value = TRUE)
  macro@meta.data <- macro@meta.data[, !(colnames(macro@meta.data) %in% columns_to_remove)]
  resolutions <- seq(0.1, 0.5, by = 0.1)
  macro <- FindNeighbors(macro, dims = 1:30, reduction = "pca")
  macro <- FindClusters(macro, resolution = resolutions )
  
  saveRDS(macro, out)

} else {
  macro <- readRDS(out)
}
clustree(macro, prefix = "RNA_snn_res.")

DimPlot(macro, group.by = "RNA_snn_res.0.2", reduction = "umap.mm", label = TRUE, label.size = 2.5, repel = TRUE, raster = FALSE )

DimPlot(macro, reduction = "umap.mm", label = TRUE, label.size = 2.5, repel = TRUE, raster = FALSE, split.by = "tissue" )

DimPlot(macro, group.by = "donor", reduction = "umap.mm", label = TRUE, label.size = 2.5, repel = TRUE, raster = FALSE, split.by = "tissue" )
Warning: ggrepel: 30 unlabeled data points (too many overlaps). Consider
increasing max.overlaps
Warning: ggrepel: 32 unlabeled data points (too many overlaps). Consider
increasing max.overlaps
Warning: ggrepel: 21 unlabeled data points (too many overlaps). Consider
increasing max.overlaps

opt_res <- "RNA_snn_res.0.2"  
n <- nlevels(macro$RNA_snn_res.0.2)
macro$RNA_snn_res.0.2 <- factor(macro$RNA_snn_res.0.2, levels = seq(0,n-1))
macro$seurat_clusters <- NULL
macro$cluster <- macro$RNA_snn_res.0.2
Idents(macro) <- macro$cluster
DimPlot(macro, group.by = "donor", reduction = "umap.mm", label = TRUE, label.size = 2.5, repel = TRUE, raster = FALSE, split.by = "tissue" )
Warning: ggrepel: 30 unlabeled data points (too many overlaps). Consider
increasing max.overlaps
Warning: ggrepel: 32 unlabeled data points (too many overlaps). Consider
increasing max.overlaps
Warning: ggrepel: 21 unlabeled data points (too many overlaps). Consider
increasing max.overlaps

k <- names(which(table(macro$celltypist_lung_airway) > 150))
macro$celltypist_lung2 <- ifelse(macro$celltypist_lung_airway %in% k,
                                  macro$celltypist_lung_airway, NA)

DimPlot(
    macro,
    group.by = "celltypist_lung2",
    reduction = "umap.mm",
    label = TRUE,
    repel = TRUE,
    label.size = 2.5,
    shuffle = TRUE
)
plot_ct_bar(macro, "celltypist_paeds")

plot_ct_bar(macro, "celltypist_lung_airway")

plot_ct_bar(macro, "celltypist_lung_atlas")

plot_ct_bar(macro, "celltypist_immune")

plot_ct_bar(macro, "celltypist_tonsil")

#plot_ct_bar(macro, "predicted.celltype.l2") # Tonsil Atlas level 2 
plot_ct_bar(macro, "cell_labels_v2") # previous labels

Checking for known markers-

features <-c("C1QB", "C1QC", "APOE", "FABP4", "MARCO", "INHBA", "MCEMP1")

    FeaturePlot(
      macro,                     
      features  = features,
      reduction = "umap.mm",     
      raster    = FALSE,
      label     = TRUE,
      ncol = 1,
      split.by  = "tissue"
    )
Warning: All cells have the same value (0) of "FABP4"
All cells have the same value (0) of "FABP4"

Renaming idents

new_levels <- c(
  "0" = "Alveolar Macrophages",
  "1" = "Alveolar Macrophages",
  "2" = "Alveolar Macrophages",
   "3" = "Alveolar Macrophages",
  "4" =    "Macro CCL",
  "5" =  "Intermediate Macrophages",
    "6" ="Macro_tonsils_adenoids",
"7" = "Monocyte-derived Macrophages",
"8" = "Proliferating Macrophages",
"9" = "Alveolar Macrophages",
"10" =  "Interstitial Macrophages"
)
macro <- RenameIdents(macro, new_levels)

macro$cell_labels_l2 <- Idents(macro)
macro$cell_labels_l3 <- Idents(macro)
table(macro$cell_labels_l2, macro$tissue)
                              
                               Tonsils Adenoids Nasal_brushings
  Alveolar Macrophages               0        0               1
  Macro CCL                          0        0               3
  Intermediate Macrophages          48       31             368
  Macro_tonsils_adenoids          1506      801              27
  Monocyte-derived Macrophages       0        0               1
  Proliferating Macrophages          0        0              27
  Interstitial Macrophages           0        0               0
                              
                               Bronchial_brushings   BAL
  Alveolar Macrophages                        6468 29523
  Macro CCL                                   1739  2092
  Intermediate Macrophages                    1789  1231
  Macro_tonsils_adenoids                        13    22
  Monocyte-derived Macrophages                  11  2138
  Proliferating Macrophages                    160  1742
  Interstitial Macrophages                     972   119
df_counts <- macro@meta.data |>
  count(cell_labels_l2, tissue, name = "n")
df_totals <- df_counts |>
  group_by(cell_labels_l2) |>
  summarise(total = sum(n), .groups = "drop")

ggplot(df_counts, aes(x = cell_labels_l2, y = n, fill = tissue)) +
  geom_bar(stat = "identity", position = "stack") +
  geom_text(data = df_totals,
            aes(x = cell_labels_l2, y = total, label = total),
            vjust = -0.3, size = 3, inherit.aes = FALSE) +
  theme_minimal() +
  theme(axis.text.x = element_text(angle = 90, hjust = 1)) +
  labs(x = "cluster at l2",
       y = "Number of cells",
       fill = "Tissue")

macro@meta.data %>%
    count(cell_labels_l2, donor) %>%
    group_by(cell_labels_l2) %>%
    mutate(prop = n / sum(n)) %>%  # proportion per L2
    ungroup() %>%
    ggplot(aes(x = cell_labels_l2, y = prop, fill = donor)) +
    geom_col(width = 0.9) +
    scale_y_continuous(labels = scales::percent) +
    scale_fill_viridis_d(option = "turbo", guide = guide_legend(ncol = 2)) +
    theme_classic() +
    theme(axis.text.x = element_text(angle = 45, hjust = 1)) +
    labs(
        x = NULL,
        y = "Proportion",
        fill = "Donor",
        title = "Donor composition per cluster"
    )

DE markers

The marker genes for this subclustering can be found here-

Macrophages_population_subclusters

features.no.mt <- rownames(paed_sub)[!grepl("^MT-", rownames(paed_sub))]
gc()
             used   (Mb) gc trigger    (Mb) limit (Mb)   max used    (Mb)
Ncells    4436866  237.0   10214208   545.5         NA   10214208   545.5
Vcells 1272521453 9708.6 2092132540 15961.8     163840 2092111009 15961.6
paed_sub.markers <- FindAllMarkers(
  macro,
  only.pos = TRUE,
  min.pct = 0.25,
  logfc.threshold = 0.25,
  features = features.no.mt
)
gc()
             used   (Mb) gc trigger    (Mb) limit (Mb)   max used    (Mb)
Ncells    4501622  240.5   10214208   545.5         NA   10214208   545.5
Vcells 1272706287 9710.0 3012846857 22986.2     163840 2991650449 22824.5
top10 <- paed_sub.markers %>%
    group_by(cluster) %>%
    top_n(n = 10, wt = avg_log2FC) %>%
    ungroup() %>%
    distinct(gene, .keep_all = TRUE) %>%
    arrange(cluster, desc(avg_log2FC))


best.wilcox.gene.per.cluster <- paed_sub.markers |>
  group_by(cluster) |>
  slice_max(avg_log2FC, n = 1) |>
  select(cluster, gene)
best.wilcox.gene.per.cluster
# A tibble: 7 × 2
# Groups:   cluster [7]
  cluster                      gene   
  <fct>                        <chr>  
1 Alveolar Macrophages         SLC19A3
2 Macro CCL                    IL6    
3 Intermediate Macrophages     SDS    
4 Macro_tonsils_adenoids       PLA2G2D
5 Monocyte-derived Macrophages CCR2   
6 Proliferating Macrophages    KIF18B 
7 Interstitial Macrophages     CXCL10 
FeaturePlot(
  macro,
  features = best.wilcox.gene.per.cluster$gene,
  reduction = "umap.mm",
  raster = FALSE,
  repel = TRUE,
  ncol = 1,
  label = TRUE,
  split.by = "tissue"
) 

## Seurat top markers
cluster_colors <- paletteer::paletteer_d("pals::glasbey")[factor(top10$cluster)]

DotPlot(macro,    
        features = unique(top10$gene),
        group.by = "cell_labels_l3",
        cols = c("azure1", "blueviolet"),
        dot.scale = 3, assay = "RNA") +
    RotatedAxis() +
    FontSize(y.text = 8, x.text = 12) +
    labs(y = element_blank(), x = element_blank()) +
    coord_flip() +
    theme(axis.text.y = element_text(color = cluster_colors)) +
    ggtitle("Top 10 marker genes per cluster (Seurat)")

marker_counts <- paed_sub.markers %>%
    filter(avg_log2FC > 1 & p_val_adj < 0.05) %>%  
    count(cluster, name = "n_markers") %>%
    arrange(desc(n_markers))

print(marker_counts)
                       cluster n_markers
1     Intermediate Macrophages       443
2    Proliferating Macrophages       395
3       Macro_tonsils_adenoids       321
4     Interstitial Macrophages       261
5 Monocyte-derived Macrophages       152
6         Alveolar Macrophages       151
7                    Macro CCL       130
ggplot(marker_counts, aes(x = reorder(cluster, n_markers), y = n_markers)) +
    geom_col(fill = "steelblue", alpha = 0.8) +
    coord_flip() +
    labs(title = "High-confidence markers per subcluster (logFC > 1)",
         x = "Subcluster", y = "Number of markers") +
    theme_minimal()

out_markers <- here("output",
            "CSV_v3","Myeloid_lineage",
            paste("Marker_genes_Reclustered_Macrophages.",opt_res, sep = ""))

dir.create(out_markers, recursive = TRUE, showWarnings = FALSE)

for (cl in unique(paed_sub.markers$cluster)) {
  cluster_data <- paed_sub.markers %>% dplyr::filter(cluster == cl)
  file_name <- here(out_markers, paste0("G000231_Neeland_", cl, ".csv"))
  if (!file.exists(file_name)) {
  write.csv(cluster_data, file = file_name)
  }
}

Pairwise markers within Macrophages

ids <- unique(macro$cell_labels_l3)
res_list <- list()

for (i in 1:(length(ids) - 1)) {
  for (j in (i + 1):length(ids)) {

    id1 <- ids[i]
    id2 <- ids[j]

    res <- FindMarkers(
      macro,
      ident.1 = id1,
      ident.2 = id2,
      only.pos = TRUE,
      min.pct = 0.25,
      logfc.threshold = 0.5,
      features = features.no.mt
    )

    top_genes <- res %>%
      filter(pct.1 > 0.1 & pct.2 > 0.1) %>%
      arrange(p_val_adj, desc(avg_log2FC)) %>%
      head(10) %>%
      rownames()

    res_list[[paste(id1, "vs", id2, sep = "_")]] <- res

    cat("\n\n###", paste(id1, "vs", id2), "\n\n")

    print(
      VlnPlot(
        macro,
        features = top_genes,
        idents = c(id1, id2),
        pt.size = 0
      ) + patchwork::plot_annotation(title = paste(id1, "vs", id2))
    )

    kable(
      head(res, 50),
      caption = paste("Top 50 markers for", id1, "vs", id2)
    )
  }
}

Intermediate Macrophages vs Proliferating Macrophages

Intermediate Macrophages vs Macro_tonsils_adenoids

Intermediate Macrophages vs Macro CCL

Intermediate Macrophages vs Alveolar Macrophages

Intermediate Macrophages vs Monocyte-derived Macrophages

Intermediate Macrophages vs Interstitial Macrophages

Proliferating Macrophages vs Macro_tonsils_adenoids

Proliferating Macrophages vs Macro CCL

Proliferating Macrophages vs Alveolar Macrophages

Proliferating Macrophages vs Monocyte-derived Macrophages

Proliferating Macrophages vs Interstitial Macrophages

Macro_tonsils_adenoids vs Macro CCL

Macro_tonsils_adenoids vs Alveolar Macrophages

Macro_tonsils_adenoids vs Monocyte-derived Macrophages

Macro_tonsils_adenoids vs Interstitial Macrophages

Macro CCL vs Alveolar Macrophages

Macro CCL vs Monocyte-derived Macrophages

Macro CCL vs Interstitial Macrophages

Alveolar Macrophages vs Monocyte-derived Macrophages

Alveolar Macrophages vs Interstitial Macrophages

Monocyte-derived Macrophages vs Interstitial Macrophages

Overview of l1, l2 and l3

paed_sub$cell_labels_l2 <- as.character(paed_sub$cell_labels_l1)
paed_sub$cell_labels_l3 <- as.character(paed_sub$cell_labels_l2)

for (obj in list(macro)) {
  cells <- intersect(colnames(paed_sub), colnames(obj))
  paed_sub$cell_labels_l2[cells] <- as.character(obj$cell_labels_l2[cells])
  paed_sub$cell_labels_l3[cells] <- as.character(obj$cell_labels_l3[cells])
}

paed_sub$cell_labels_l2 <- factor(paed_sub$cell_labels_l2)
paed_sub$cell_labels_l3 <- factor(paed_sub$cell_labels_l3)
DimPlot(paed_sub, reduction = "umap.sub", raster = F, repel = T, label = T, group.by = "lineage")

DimPlot(paed_sub, reduction = "umap.sub", raster = F, repel = T, label = T, group.by = "cell_labels_l1")

paed_sub$cell_labels_l2 <- factor(
  paed_sub$cell_labels_l2,
  levels = c(
    "Monocytes",
    "Monocyte-derived Macrophages",
    "Intermediate Macrophages",
    "Proliferating Macrophages",
    "Macro CCL",
    "Macro_tonsils_adenoids",
    "Interstitial Macrophages",
    "Alveolar Macrophages",
    "DC1",
    "DC2",
    "Activated DCs",
    "Plasmacytoid DCs",
    "Follicular DCs",
    "Neutrophils",
    "Mast cells",
    "Ionocytes",
    "Melanocyte"
  )
)

Idents(paed_sub) <- paed_sub$cell_labels_l2
DimPlot(paed_sub, reduction = "umap.sub", raster = F, repel = T, label = T, group.by = "cell_labels_l2")

DimPlot(paed_sub, reduction = "umap.sub", raster = F, repel = T, label = T, group.by = "cell_labels_l2", split.by = "tissue")

paed_sub$cell_labels_l3 <- paed_sub$cell_labels_l2
#DimPlot(paed_sub, reduction = "umap.sub", raster = F, repel = T, label = T, group.by = "cell_labels_l3", split.by = "tissue")

Save RDS

#saveRDS(paed_sub, "../EarlyAir_paper/output/RDS/Annotated_lineages/SEU_Myeloid_annotated.rds")
paed_sub@meta.data %>%
  count(cell_labels_l3, tissue) %>%
  group_by(cell_labels_l3) %>%
  mutate(total = sum(n)) %>%
  ungroup() %>%
  #mutate(cell_labels_l3 = reorder(cell_labels_l3, total)) %>%
  ggplot(aes(x = cell_labels_l3, y = n, fill = tissue)) +
  geom_col() +
  geom_text(
    data = function(df) distinct(df, cell_labels_l3, total),
    aes(x = cell_labels_l3, y = total, label = total),
    inherit.aes = FALSE,
    hjust = -0.1,
    size = 3
  ) +
  coord_flip() +
  theme_classic() +
  labs(
    x = NULL,
    y = "Cell count",
    fill = "Tissue",
    title = "Cell counts by L3 cell type (stacked by tissue)"
  )

paed_sub@meta.data %>%
    count(cell_labels_l3, donor) %>%
    group_by(cell_labels_l3) %>%
    mutate(prop = n / sum(n)) %>%  # proportion per L3
    ungroup() %>%
    ggplot(aes(x = cell_labels_l3, y = prop, fill = donor)) +
    geom_col(width = 0.9) +
    scale_y_continuous(labels = scales::percent) +
    scale_fill_viridis_d(option = "turbo", guide = guide_legend(ncol = 2)) +
    theme_classic() +
    theme(axis.text.x = element_text(angle = 45, hjust = 1)) +
    labs(
        x = NULL,
        y = "Proportion",
        fill = "Donor",
        title = "Donor composition per L3 cell type"
    )

Some notes/observations from Myeloid cell lineage-

  1. Macrophages mostly present in BAL and Bronchial brushings
  2. Macrophages found in Tonsils/Adenoids sit in a separate cluster, I’ve named it as Macro_tonsil_adenoids. It shows SLAN+ (slan-like signature)
  3. Macrophages have donor confounding effect

Session Info

sessionInfo()
R version 4.3.3 (2024-02-29)
Platform: aarch64-apple-darwin20 (64-bit)
Running under: macOS 15.7.3

Matrix products: default
BLAS:   /Library/Frameworks/R.framework/Versions/4.3-arm64/Resources/lib/libRblas.0.dylib 
LAPACK: /Library/Frameworks/R.framework/Versions/4.3-arm64/Resources/lib/libRlapack.dylib;  LAPACK version 3.11.0

locale:
[1] en_US.UTF-8/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8

time zone: Australia/Melbourne
tzcode source: internal

attached base packages:
[1] stats     graphics  grDevices utils     datasets  methods   base     

other attached packages:
 [1] gtools_3.9.5        presto_1.0.0        circlize_0.4.16    
 [4] BiocParallel_1.36.0 harmony_1.2.3       Rcpp_1.0.14        
 [7] patchwork_1.3.0     cowplot_1.1.3       data.table_1.17.2  
[10] RColorBrewer_1.1-3  kableExtra_1.4.0    ggridges_0.5.6     
[13] ggforce_0.4.2       viridis_0.6.5       viridisLite_0.4.2  
[16] paletteer_1.6.0     clustree_0.5.1      ggraph_2.2.1       
[19] Seurat_5.0.3        SeuratObject_5.1.0  sp_2.2-0           
[22] here_1.0.1          lubridate_1.9.4     forcats_1.0.0      
[25] stringr_1.5.1       dplyr_1.1.4         purrr_1.0.4        
[28] readr_2.1.5         tidyr_1.3.1         tibble_3.2.1       
[31] ggplot2_3.5.2       tidyverse_2.0.0     BiocStyle_2.30.0   
[34] workflowr_1.7.1    

loaded via a namespace (and not attached):
  [1] RcppAnnoy_0.0.22       splines_4.3.3          later_1.4.2           
  [4] prismatic_1.1.2        polyclip_1.10-7        fastDummies_1.7.5     
  [7] lifecycle_1.0.4        rprojroot_2.0.4        globals_0.17.0        
 [10] processx_3.8.6         lattice_0.22-5         MASS_7.3-60.0.1       
 [13] backports_1.5.0        magrittr_2.0.3         limma_3.58.1          
 [16] plotly_4.10.4          sass_0.4.10            rmarkdown_2.29        
 [19] jquerylib_0.1.4        yaml_2.3.10            httpuv_1.6.16         
 [22] sctransform_0.4.2      spam_2.11-1            spatstat.sparse_3.1-0 
 [25] reticulate_1.42.0      pbapply_1.7-2          abind_1.4-8           
 [28] Rtsne_0.17             tweenr_2.0.3           git2r_0.36.2          
 [31] ggrepel_0.9.6          irlba_2.3.5.1          listenv_0.9.1         
 [34] spatstat.utils_3.1-4   pheatmap_1.0.12        goftest_1.2-3         
 [37] RSpectra_0.16-2        spatstat.random_3.4-1  fitdistrplus_1.2-2    
 [40] parallelly_1.44.0      svglite_2.2.1          leiden_0.4.3.1        
 [43] codetools_0.2-19       xml2_1.3.8             tidyselect_1.2.1      
 [46] shape_1.4.6.1          farver_2.1.2           matrixStats_1.5.0     
 [49] spatstat.explore_3.4-3 jsonlite_2.0.0         tidygraph_1.3.1       
 [52] progressr_0.15.1       survival_3.5-8         systemfonts_1.2.3     
 [55] tools_4.3.3            ica_1.0-3              glue_1.8.0            
 [58] gridExtra_2.3          xfun_0.52              withr_3.0.2           
 [61] BiocManager_1.30.25    fastmap_1.2.0          callr_3.7.6           
 [64] digest_0.6.37          timechange_0.3.0       R6_2.6.1              
 [67] mime_0.13              textshaping_1.0.1      colorspace_2.1-1      
 [70] scattermore_1.2        tensor_1.5             spatstat.data_3.1-6   
 [73] utf8_1.2.5             generics_0.1.4         graphlayouts_1.2.2    
 [76] httr_1.4.7             htmlwidgets_1.6.4      whisker_0.4.1         
 [79] uwot_0.2.3             pkgconfig_2.0.3        gtable_0.3.6          
 [82] lmtest_0.9-40          htmltools_0.5.8.1      dotCall64_1.2         
 [85] scales_1.4.0           png_0.1-8              spatstat.univar_3.1-3 
 [88] knitr_1.50             rstudioapi_0.17.1      tzdb_0.5.0            
 [91] reshape2_1.4.4         checkmate_2.3.2        nlme_3.1-164          
 [94] cachem_1.1.0           zoo_1.8-14             GlobalOptions_0.1.2   
 [97] KernSmooth_2.23-22     vipor_0.4.7            parallel_4.3.3        
[100] miniUI_0.1.2           ggrastr_1.0.2          pillar_1.10.2         
[103] grid_4.3.3             vctrs_0.6.5            RANN_2.6.2            
[106] promises_1.3.2         xtable_1.8-4           cluster_2.1.6         
[109] beeswarm_0.4.0         evaluate_1.0.3         cli_3.6.5             
[112] compiler_4.3.3         rlang_1.1.6            crayon_1.5.3          
[115] future.apply_1.11.3    labeling_0.4.3         rematch2_2.1.2        
[118] ps_1.9.1               ggbeeswarm_0.7.2       getPass_0.2-4         
[121] plyr_1.8.9             fs_1.6.6               stringi_1.8.7         
[124] deldir_2.0-4           lazyeval_0.2.2         spatstat.geom_3.4-1   
[127] Matrix_1.6-5           RcppHNSW_0.6.0         hms_1.1.3             
[130] future_1.40.0          statmod_1.5.0          shiny_1.10.0          
[133] ROCR_1.0-11            igraph_2.1.4           memoise_2.0.1         
[136] bslib_0.9.0           

sessionInfo()
R version 4.3.3 (2024-02-29)
Platform: aarch64-apple-darwin20 (64-bit)
Running under: macOS 15.7.3

Matrix products: default
BLAS:   /Library/Frameworks/R.framework/Versions/4.3-arm64/Resources/lib/libRblas.0.dylib 
LAPACK: /Library/Frameworks/R.framework/Versions/4.3-arm64/Resources/lib/libRlapack.dylib;  LAPACK version 3.11.0

locale:
[1] en_US.UTF-8/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8

time zone: Australia/Melbourne
tzcode source: internal

attached base packages:
[1] stats     graphics  grDevices utils     datasets  methods   base     

other attached packages:
 [1] gtools_3.9.5        presto_1.0.0        circlize_0.4.16    
 [4] BiocParallel_1.36.0 harmony_1.2.3       Rcpp_1.0.14        
 [7] patchwork_1.3.0     cowplot_1.1.3       data.table_1.17.2  
[10] RColorBrewer_1.1-3  kableExtra_1.4.0    ggridges_0.5.6     
[13] ggforce_0.4.2       viridis_0.6.5       viridisLite_0.4.2  
[16] paletteer_1.6.0     clustree_0.5.1      ggraph_2.2.1       
[19] Seurat_5.0.3        SeuratObject_5.1.0  sp_2.2-0           
[22] here_1.0.1          lubridate_1.9.4     forcats_1.0.0      
[25] stringr_1.5.1       dplyr_1.1.4         purrr_1.0.4        
[28] readr_2.1.5         tidyr_1.3.1         tibble_3.2.1       
[31] ggplot2_3.5.2       tidyverse_2.0.0     BiocStyle_2.30.0   
[34] workflowr_1.7.1    

loaded via a namespace (and not attached):
  [1] RcppAnnoy_0.0.22       splines_4.3.3          later_1.4.2           
  [4] prismatic_1.1.2        polyclip_1.10-7        fastDummies_1.7.5     
  [7] lifecycle_1.0.4        rprojroot_2.0.4        globals_0.17.0        
 [10] processx_3.8.6         lattice_0.22-5         MASS_7.3-60.0.1       
 [13] backports_1.5.0        magrittr_2.0.3         limma_3.58.1          
 [16] plotly_4.10.4          sass_0.4.10            rmarkdown_2.29        
 [19] jquerylib_0.1.4        yaml_2.3.10            httpuv_1.6.16         
 [22] sctransform_0.4.2      spam_2.11-1            spatstat.sparse_3.1-0 
 [25] reticulate_1.42.0      pbapply_1.7-2          abind_1.4-8           
 [28] Rtsne_0.17             tweenr_2.0.3           git2r_0.36.2          
 [31] ggrepel_0.9.6          irlba_2.3.5.1          listenv_0.9.1         
 [34] spatstat.utils_3.1-4   pheatmap_1.0.12        goftest_1.2-3         
 [37] RSpectra_0.16-2        spatstat.random_3.4-1  fitdistrplus_1.2-2    
 [40] parallelly_1.44.0      svglite_2.2.1          leiden_0.4.3.1        
 [43] codetools_0.2-19       xml2_1.3.8             tidyselect_1.2.1      
 [46] shape_1.4.6.1          farver_2.1.2           matrixStats_1.5.0     
 [49] spatstat.explore_3.4-3 jsonlite_2.0.0         tidygraph_1.3.1       
 [52] progressr_0.15.1       survival_3.5-8         systemfonts_1.2.3     
 [55] tools_4.3.3            ica_1.0-3              glue_1.8.0            
 [58] gridExtra_2.3          xfun_0.52              withr_3.0.2           
 [61] BiocManager_1.30.25    fastmap_1.2.0          callr_3.7.6           
 [64] digest_0.6.37          timechange_0.3.0       R6_2.6.1              
 [67] mime_0.13              textshaping_1.0.1      colorspace_2.1-1      
 [70] scattermore_1.2        tensor_1.5             spatstat.data_3.1-6   
 [73] utf8_1.2.5             generics_0.1.4         graphlayouts_1.2.2    
 [76] httr_1.4.7             htmlwidgets_1.6.4      whisker_0.4.1         
 [79] uwot_0.2.3             pkgconfig_2.0.3        gtable_0.3.6          
 [82] lmtest_0.9-40          htmltools_0.5.8.1      dotCall64_1.2         
 [85] scales_1.4.0           png_0.1-8              spatstat.univar_3.1-3 
 [88] knitr_1.50             rstudioapi_0.17.1      tzdb_0.5.0            
 [91] reshape2_1.4.4         checkmate_2.3.2        nlme_3.1-164          
 [94] cachem_1.1.0           zoo_1.8-14             GlobalOptions_0.1.2   
 [97] KernSmooth_2.23-22     vipor_0.4.7            parallel_4.3.3        
[100] miniUI_0.1.2           ggrastr_1.0.2          pillar_1.10.2         
[103] grid_4.3.3             vctrs_0.6.5            RANN_2.6.2            
[106] promises_1.3.2         xtable_1.8-4           cluster_2.1.6         
[109] beeswarm_0.4.0         evaluate_1.0.3         cli_3.6.5             
[112] compiler_4.3.3         rlang_1.1.6            crayon_1.5.3          
[115] future.apply_1.11.3    labeling_0.4.3         rematch2_2.1.2        
[118] ps_1.9.1               ggbeeswarm_0.7.2       getPass_0.2-4         
[121] plyr_1.8.9             fs_1.6.6               stringi_1.8.7         
[124] deldir_2.0-4           lazyeval_0.2.2         spatstat.geom_3.4-1   
[127] Matrix_1.6-5           RcppHNSW_0.6.0         hms_1.1.3             
[130] future_1.40.0          statmod_1.5.0          shiny_1.10.0          
[133] ROCR_1.0-11            igraph_2.1.4           memoise_2.0.1         
[136] bslib_0.9.0