Reviewers Responses

A second question stemming from this is how good (realistic) these organoids are. The scRNA-seq is in the main figures and text mainly (only?) discussed and interpreted from the nephron point of view (which indeed looks impressive but as argued before the comparison to established methods is minimal). Only Sup Fig 3 shows there is a huge stromal component in these organoids which is (at this moment) not addressed in the text. What is this? Where is it coming from? Are these specific renal stromal cell types or more general (‘stromal’ after all means something different to different people). Is this also the case with the old methods, if so is this stromal compartment increasing or not?

Further analysis of existing data: This is asking for more detailed analysis of the stroma, which is a little beyond the scope of this paper. However, the “what is this?” question could be addressed by further scRNAseq analysis - isolation of stromal component and comparison to existing mouse datasets where stroma has been well-characterised, as well as immunofluorescence for stromal markers. The “where is this coming from” would require detailed lineage tracing, part of a separate question that is not the focus of this paper. We could argue this. We cannot be expected to study the origin of individual cell types within these organoids. This is unreasonable. Sean has offered to do some direct stromal comparisons, but I feel this still remains largely supplementary data.

Only further in the manuscript (Fig 5) we return to this stromal core in the organoids as an explanation for the radial morphology of the PE-organoids as a source of WNT inhibition. Here it turns out there is mainly cartilage developing at the stromal core of the organoids, and only now it is made clear that clusters 2 and 5 of the seq data are in fact ‘cartilage clusters’ (line 422). Sure, other protocols show cartilage development (as mentioned in the discussion) but again the comparison to the established methods becomes essential and is lacking (Fig 5B only shows PE-organoids). Further analysis of existing data potentially possible: Reviewer interpreted the data to show that ‘mainly cartilage’ develops in the centre of the organoid – I am not sure of exact proportion but it is not all of the core. Some further analysis of the core (eg- sox9 staining) could be beneficial here. In terms of comparison with existing protocols – this could be addressed through further scRNAseq analysis and could be rolled together with the stromal analysis above).
Add to the above analysis.

Strategy

  1. Re-evaluate the clustering and DKCC analysis from the original paper
  2. Compare to the England et al. 2020 stromal genes
  3. Perform an analysis in Azimuth using the developing fetal atlas

Part 1: Re-evaluate

## This version of bslib is designed to work with shiny version 1.6.0 or higher.
## Registered S3 method overwritten by 'cli':
##   method     from    
##   print.boxx spatstat
## ── Attaching packages ─────────────────────────────────────── tidyverse 1.3.0 ──
## ✓ ggplot2 3.3.5     ✓ purrr   0.3.4
## ✓ tibble  3.0.4     ✓ dplyr   1.0.6
## ✓ tidyr   1.1.4     ✓ stringr 1.4.0
## ✓ readr   1.4.0     ✓ forcats 0.5.1
## ── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
## x dplyr::filter() masks stats::filter()
## x dplyr::lag()    masks stats::lag()

get unassigned cells from devkidcc and recluster

We can grab both the “Stroma” and “unassigned” cells from the DevKidCC analysis. These represent the stroma and mesenchymal populations present, thus the cells we are interested in.

Can look at the gene expression within each subset of stroma.

How do these cells look when plotted out in the UMAP coords?

## Warning: Using `as.character()` on a quosure is deprecated as of rlang 0.3.0.
## Please use `as_label()` or `as_name()` instead.
## This warning is displayed once per session.

So using a resolution of 0.2, what is the breakdown of cells from DevKidCC analysis in these clusters?

##             
##                 0    1    2    3    4    5    6    7
##   CS            0 2538 1184 1374  129    0   18    0
##   MesS          0   90   97   66    0    0   17    4
##   MS            0   24   14    6    0    0    7    0
##   unassigned 4950  880 1433  310 1499 1148   63   72

clusters 0, 5 and 7 were predominantly unassigned, while clusters 1 and 3 were predominantly assigned as stroma. clusters 2, 3 and 4 were a mix.

The clusters are divided cleanly by the two ages of the samples, which is expected.

Some interesting and useful markers jump out here. OGN is expressed by developing cartilage. CXCL12 and CRABP1 are strongly coexpressed with OGN. These cells are in clusters 1,2,3 and 7. There are collagens present, COL2A1, COL3A1 and COL21A1 coming up as markers. Clusters 0 and 5 have a PAX8 signature, indicating these may be mesenchymal cells that go on to contribute to the nephron lineages.

## [[1]]

## 
## [[2]]

## 
## [[3]]

## 
## [[4]]

## 
## [[5]]

## 
## [[6]]

run through toppfun and stringdb

Can export these gene lists and run them through a geneset analysis process like StringDB or ToppFun.

The marker list is available at StromaUnassignedMarkers.xlsx.

I ran the top 50 DE genes through string db for cluster 1, which has OGN as it’s most DE gene. The result indicates there is an enrichment for biological processes such as Biomineral Tissue Development, Ossification, Skeletal System Developmet. Lots of collagens.

We can focus on completing this analysis if you think it’s useful. For now i’ll kick on with the single cell specific stuff.

Part 2: Comparisons to Stroma Publications

England et al. 2020 comparison

Here we can plot on the dataset where the genes that were identified as stromal population markers in the England et al., 2020 paper.

## Warning in FetchData(object = object, vars = features, slot = slot): The
## following requested variables were not found: CLCA3A1

Cluster 2, with an enrichment of FIBIN, SMOC2, DLK1 is reminiscent of a subset of the cortical stroma in E18.5 mice.

Tanigawa et al. 2022 comparison

Recent paper claiming to differentiate stromal progenitors from mouse ESC and then culture them with other mouse ESC diff pops (ub, np) to make higher order organoids.
They published some cool single cell data and lists of de genes so we can use those. 

## New names:
## * `` -> ...1
## New names:
## * `` -> ...1
## New names:
## * `` -> ...1
## New names:
## * `` -> ...1
## New names:
## * `` -> ...1
## New names:
## * `` -> ...1
## New names:
## * `` -> ...1
## New names:
## * `` -> ...1
## New names:
## * `` -> ...1
## * `ave.exp(P0)` -> `ave.exp(P0)...7`
## * `ave.exp(P0)` -> `ave.exp(P0)...8`
## New names:
## * `` -> ...1
## Warning in FeaturePlot(mes, features = .x[1:9], order = T): All cells have the
## same value (0) of MUSTN1.
## Warning in FeaturePlot(mes, features = .x[1:9], order = T): All cells have the
## same value (0) of MYOC.
## Warning in FeaturePlot(mes, features = .x[1:9], order = T): All cells have the
## same value (0) of PGR.
## $`REN and Mes(#22)`

## 
## $`Mes(#17)`

## 
## $`SP(#8)`

## 
## $`CS(#7)`

## 
## $`OM(#5)`

## 
## $`IM(#2)`

## 
## $`US(#15)`

## 
## $`iSP specific(#35)`

## 
## $`nd-iS specific(#0)`

## 
## $`nd-iS specific(#9)`

Can’t see much of a pattern with these, but it is a limited gene list. I’ll plot the heatmaps, using the top 30 genes for each segment.

## $`REN and Mes(#22)`

## 
## $`Mes(#17)`

## 
## $`SP(#8)`

## 
## $`CS(#7)`

## 
## $`OM(#5)`

## 
## $`IM(#2)`

## 
## $`US(#15)`

## 
## $`iSP specific(#35)`

## 
## $`nd-iS specific(#0)`

## 
## $`nd-iS specific(#9)`

IN general, some of the genes for each identity are expressed but there are not definite patterns that identify a match between gene populations.

try SCINA to predict identity

The last effort here would be using a predictive tool. SCINA takes lists of genes and predicts cell identity from these.

Well that is curious! while the majority of cells do come up as unknown, there are some clearly defined regions that have been identified as a particular stromal population. The D13 cells are split between a SP (stromal progenitor) population, an IM (inner medullary) population, and some few Mes (mesangial) cells. Moving into the older cells, D13p14, there is a small group of CS (cortical stroma), a larger group of OM (outer medullary) and some cells that associated with the induced stroma (iSP and nd-iS).

## 
##             CS(#7)             IM(#2)  iSP specific(#35)           Mes(#17) 
##                624                317                 48                543 
## nd-iS specific(#0)             OM(#5)   REN and Mes(#22)             SP(#8) 
##                240               1747                 89               2452 
##            unknown            US(#15) 
##               9826                 37

Part 3: Azimuth

run cells through azimuth

I can only run small subsets of the data through Azimuth (upload limits). As such I subset the data and took 10% of all the cells, evenly distributed between the two ages.

I ran this through Azimuth, using the “Human Fetal Development” reference.

## Loading required package: SingleCellExperiment
## Loading required package: SummarizedExperiment
## Loading required package: GenomicRanges
## Loading required package: stats4
## Loading required package: BiocGenerics
## Loading required package: parallel
## 
## Attaching package: 'BiocGenerics'
## The following objects are masked from 'package:parallel':
## 
##     clusterApply, clusterApplyLB, clusterCall, clusterEvalQ,
##     clusterExport, clusterMap, parApply, parCapply, parLapply,
##     parLapplyLB, parRapply, parSapply, parSapplyLB
## The following objects are masked from 'package:dplyr':
## 
##     combine, intersect, setdiff, union
## The following objects are masked from 'package:stats':
## 
##     IQR, mad, sd, var, xtabs
## The following objects are masked from 'package:base':
## 
##     anyDuplicated, append, as.data.frame, basename, cbind, colnames,
##     dirname, do.call, duplicated, eval, evalq, Filter, Find, get, grep,
##     grepl, intersect, is.unsorted, lapply, Map, mapply, match, mget,
##     order, paste, pmax, pmax.int, pmin, pmin.int, Position, rank,
##     rbind, Reduce, rownames, sapply, setdiff, sort, table, tapply,
##     union, unique, unsplit, which, which.max, which.min
## Loading required package: S4Vectors
## 
## Attaching package: 'S4Vectors'
## The following objects are masked from 'package:dplyr':
## 
##     first, rename
## The following object is masked from 'package:tidyr':
## 
##     expand
## The following object is masked from 'package:base':
## 
##     expand.grid
## Loading required package: IRanges
## 
## Attaching package: 'IRanges'
## The following objects are masked from 'package:dplyr':
## 
##     collapse, desc, slice
## The following object is masked from 'package:purrr':
## 
##     reduce
## Loading required package: GenomeInfoDb
## Loading required package: Biobase
## Welcome to Bioconductor
## 
##     Vignettes contain introductory material; view with
##     'browseVignettes()'. To cite Bioconductor, see
##     'citation("Biobase")', and for packages 'citation("pkgname")'.
## Loading required package: DelayedArray
## Loading required package: matrixStats
## 
## Attaching package: 'matrixStats'
## The following objects are masked from 'package:Biobase':
## 
##     anyMissing, rowMedians
## The following object is masked from 'package:dplyr':
## 
##     count
## Loading required package: BiocParallel
## 
## Attaching package: 'DelayedArray'
## The following objects are masked from 'package:matrixStats':
## 
##     colMaxs, colMins, colRanges, rowMaxs, rowMins, rowRanges
## The following object is masked from 'package:purrr':
## 
##     simplify
## The following objects are masked from 'package:base':
## 
##     aperm, apply, rowsum
## 
## Attaching package: 'SummarizedExperiment'
## The following object is masked from 'package:Seurat':
## 
##     Assays
## Warning: Keys should be one or more alphanumeric characters followed by an
## underscore, setting key from PC__ to PC_
## Warning: All keys should be one or more alphanumeric characters followed by an
## underscore '_', setting key to PC_
## Warning: Keys should be one or more alphanumeric characters followed by an
## underscore, setting key from UMAP__ to UMAP_
## Warning: All keys should be one or more alphanumeric characters followed by an
## underscore '_', setting key to UMAP_
## Warning: Keys should be one or more alphanumeric characters followed by an
## underscore, setting key from UMAP__ to UMAP_
## Warning: All keys should be one or more alphanumeric characters followed by an
## underscore '_', setting key to UMAP_
## Warning: Cannot add objects with duplicate keys (offending key: UMAP_), setting
## key to 'azimuth_'

The above is some trickery to get the Version 4 Seurat object into Version 3.

THe top plot uses the coordinates from the HFD dataset UMAP, while the bottom is the stromal dataset UMAP. All the cells get assigned their most similar population within the reference.

We can ask the question, which of these have predictions above, say, 0.5?

All in all, not that many.

tabulate results

Ask the question, what cells are classified as and is that score greater than 0.5 similarity?

##                                            
##                                             FALSE TRUE
##   Astrocytes                                   26    0
##   Bronchiolar and alveolar epithelial cells     1    0
##   CCL19_CCL21 positive cells                    3    0
##   Corneal and conjunctival epithelial cells     1    0
##   Endocardial cells                             1    0
##   ENS glia                                     86    5
##   ENS neurons                                  72    3
##   Epicardial fat cells                         91   94
##   Intestinal epithelial cells                   5    0
##   Mesangial cells                               5    0
##   Mesothelial cells                             4    0
##   Metanephric cells                           103  695
##   Neuroendocrine cells                         32    3
##   PAEP_MECOM positive cells                     1    0
##   Satellite cells                             157    3
##   Smooth muscle cells                          10   29
##   Squamous epithelial cells                     1    0
##   Stromal cells                                65   27
##   Vascular endothelial cells                    1    0
##   Visceral neurons                             69    0
## Warning: Groups with fewer than two data points have been dropped.

## Warning: Groups with fewer than two data points have been dropped.

## Warning: Groups with fewer than two data points have been dropped.

## Warning: Groups with fewer than two data points have been dropped.

## Warning: Groups with fewer than two data points have been dropped.

## Warning: Groups with fewer than two data points have been dropped.

Summary

As expected, there is no “smoking gun” here, no crystal clear answer to what these cells are. There are genes that indicate a cell type is present similar to that of cartilage or pre-cartilage, and that many collagen genes are being expressed on top of that. What we could draw from this is that many of the non-kidney cells in the organoid are over-expressing ECM genes such as collagens, which may be leading to a cartilagenous identity crisis or a differentation pathway towards cartilage from a general mesenchymal precursor.

I don’t think there is much in the way of a defined stromal sub-population signature, at least when using the gene lists provided through England or Tanigawa. There are genes that are certainly enriched in the renal stroma, and in some cases many of those in a population.

Finally when we look at the projection onto a fetal development dataset only one population is reliably identified, and that is the “Metanephric” cells at day 13. The day 13+14 cells almost all have low similarity scores for their best correlation.