Thinking process templates for constructing data stories with SCDNEY

Selection of data stories to illustrate scdney

Here we showcase two data stories to illustrate scdney. These two data stories and the accompanying workshops were not derived from previous studies nor have they been published elsewhere. The first data story describes the use of scdney on cell level analysis through inferring and assessing the developmental trajectory of individual cells. The second data story details the use of scdney on patient level analysis by extracting and summarising information obtained from each cell. The code for both data stories are hosted on Github as reproducible Rmarkdown files, reported in the code availability section. The underlying data are reported in the data availability section.

Workshop for data story 1

A summary of the case study is provided below with detailed information including R code hosted on our Github (Lin, Kim and Chen, 2023).

To predict the gene pairs associated with the developmental course of the differentiation of mouse hippocampal cells, we downloaded the publicly available data (from GEO with accession number GSE104323) profiling eight cell types from neural lineages of the mouse hippocampus harvested from two post-natal timepoints (day 0 and 5) (La Manno et al., 2018). For speed, we removed the Nbl1, Nbl2, Granule and CA cell types from the dataset and reduced the dataset from 18,213 to 12,935 cells. To evaluate the accuracy of the original cell type labels, we applied scReClassify (Kim et al., 2019) from the scdney package. scReClassify generates cell-type-specific probabilities for each cell, where a probability of 1 denotes the highest accuracy in classification and 0 denotes lowest accuracy. Using the maximum probability assigned to each cell, we re-labelled the cell-type annotations of cells that have inconsistent labels and have a maximum probability greater than 0.9. Then, we used the re-labelled cell-type annotations to perform marker gene analysis using Cepo, a method to determine cell-type-specific differentially stable genes (Kim et al., 2021).

To build the trajectories, we applied two commonly used trajectory inference tools, Slingshot (Street et al., 2018) and destiny (Angerer et al., 2016). Finally, to predict gene-pairs that change over the trajectory course, we used our previously developed package scHOT (Ghazanfar et al., 2020), which is available on Bioconductor. scHOT enables detection of changes in higher-order interactions in single-cell gene expression data.

Workshop for data story 2

A summary of the case study is provided below with detailed information including R code hosted on our Github (Cao and Tran, 2023).

We predict patient disease outcome using COVID-19 datasets and the packages scFeatures (Cao et al., 2022), and ClassifyR (Strbenac et al., 2015). To build the prediction model on distinguishing mild and severe outcomes, we used the publicly available Schulte-Schrepping data (Schulte-Schrepping et al., 2020). We randomly sampled 20 mild and 20 severe patient samples for the purpose of demonstrating the workshop in a reasonable amount of time. Then, we applied scFeatures from the scdney package to generate patient representations from the single-cell data. scFeatures generates interpretable molecular representations from various feature types. By doing so, we were able to represent each patient with more information than a matrix of gene expressions. At the same time, it also transformed the scRNA-seq data into a matrix of samples by features, which is a standard form for machine learning models. We generated a total of 13 matrices, one for each feature type across the feature categories of (i) cell type proportions, (ii) cell type specific gene expressions, (iii) cell type specific pathway expressions, (iv) cell type specific CCI scores and (v) overall aggregated gene expressions. The details of the feature types can be found in the scFeatures publication (Cao et al., 2022).

To build a patient outcome classification model from the patient representations, we used our previously developed package ClassifyR (Strbenac et al., 2015), which is available on Bioconductor (https://bioconductor.org/packages/ClassifyR/). ClassifyR provides an implementation of cross-validated classifications, including implementation for a range of commonly used classifiers and evaluation metrics. For this case study, we ran SVM on each of the feature types using a repeated five-fold cross-validation framework with 20 repeats. The accuracy was measured using the metric ‘balanced accuracy’ that is implemented in ClassifyR.

To assess the generalisability of the constructed model, we used the Schulte-Schrepping data as training data and another publicly available COVID-19 scRNA-seq dataset, the Wilk data (Wilk et al., 2020), as an independent testing data. First, we processed the dataset in the same way using scFeatures to generate the patient representations. Given that different datasets generate slightly different sets of features, for example, due to the difference in the genes recorded, we subset the features derived from the Schulte-Schrepping dataset and the Wilk data by their common features. We then rebuilt the model using the Schulte-Schrepping dataset as the training dataset using the same cross-validation framework as above. The best model from the 100 models (i.e., from the 20 repeated five-fold cross-validation) was identified based on balanced accuracy and evaluated on the Wilk dataset.

Thinking Process Template

Typically, in scientific research papers involving cellular data technologies, there are three key components: (1) Data, (2) Narratives, and (3) Visuals (Figure 1). Through narratives, we explain the data; through visuals, we illuminate the data; through narratives and visuals we engage. At the intersection of the three components are the product: the data stories. However, what is hidden behind these components are the critical thinking questions such as evaluation and parameter choices that happen behind the decision-making process.

Figure 1. Critical thinking questions from the data drive our decision with the Narrative and illuminate with the Visuals.

Here, we present a Thinking Process Template, to uncover the thinking process behind the construction of data stories, guided by analytical decisions. We demonstrate this in two distinct data analytical scenarios, presented as scientific questions. First, we ask, what are the cell types present in our developmental single-cell dataset, and what are the correlated gene pairs in each trajectory? Second, what features are important for disease outcome classification? In both cases we illuminate the underlying thinking strategy taken by analysts/data scientists in extracting biological knowledge from the data and drawing from the vast compendium of prior knowledge to reveal novel scientific knowledge.

Scdney - Single cell data integrative analysis

As a vehicle to demonstrate the Thinking Process Template, we present scdney (Figure 2), a series of foundational methods for single cell data analysis, including

Figure 2. scdney workflow.

(A) Collection of Data - The data stories start with data. (B) Collection of methods of scdney - The collection of methods are used for the computational analysis of data. (C) Critical Thinking - Through critical thinking, we derive the final data story.

Building upon the collection of vignettes, the Thinking Process Template examines various critical thinking questions that analysts need to make, which drives the decision for the next step in the analysis workflow. Next, using the scdney workflow, we will illustrate the process of generating two data stories. The scdney workflow starts with data (Figure 2A), the series of methods are used for the analysis of data (Figure 2B), and through the critical thinking (Figure 2C), we derive the final data story.

Narrative for data story 1 - to identify key gene-pairs associated with the developmental course

In the first data story, the aim was to identify key gene-pairs associated with the developmental course of the differentiation of mouse hippocampal cells, enabling us to find key gene sets that distinguish hippocampal development in mice using scRNA-seq data (La Manno et al., 2018) (Figure 3). Box 1 lists some questions and our thought process during the development of the story.

Figure 3. Thinking process template for analysing a single-cell RNA-seq data with a lineage trajectory.

The thinking process begins from the processed data with cell type annotations and proceeds to constructing a trajectory and extracting biological insight through identification of correlated gene pairs. The orange diamonds highlight potential questions that help us quality check the data analysis, and the orange hexagonal shapes denote the specific computational tasks that are required to answer the questions above.

The dataset we use contains eight cell types from neural lineages of the mouse hippocampus harvested from two post-natal timepoints (day 0 and 5) (La Manno et al., 2018). Whilst the main goal in the original study was to demonstrate the RNA velocity fields that describe the fate decisions governing mouse hippocampal development, our data story aims to uncover novel gene-pairs associated with these neural lineages using scHOT (Ghazanfar et al., 2020).

We start by asking whether the cell type annotations in the original data are accurate. Here our expectation is that most of the labels are accurate, and by using scReClassify (Kim et al., 2019) we demonstrate that approximately 88.4% of cells show an original classification accuracy over 0.9 (Figure 4A). Among these cells, only 1.5% (177 cells) were re-classified, suggesting that a small proportion of cells may have been mislabelled. These findings were confirmed through marker analysis using Cepo (Figure 4B), and the cells with high confidence scores were re-labelled for subsequent analyses. Once we have confirmed with further quality control questions as shown in the box and ensured the quality of the cell type annotations, we can then use these labels to perform marker gene analysis and to construct the lineage trajectories (Figure 4C).

Figure 4. Assessment of cell type labels and re-annotation of sub-optimal labels with scReClassify.

(A) Shows the proportion of cells in each confidence level, defined by scReClassify, for each cell type group. (B) The distribution of gene expression of top five marker genes in Immature Granule 2 cells as per the original labels (bottom panel) and re-classified labels (top panel). (C) UMAP of mouse brain cells coloured by cell type and faceted by cells that maintain their original labels (left) and those that have been re-classified (right).

After performing the quality control of the original annotations, we then can ask questions relating to trajectory reconstruction. In the trajectory building stage, we ask questions (see Box 1) to ensure the stability and robustness of the trajectories by testing the concordance of the pseudo-times between various trajectory reconstruction tools (Cao et al., 2019; Street et al., 2018). In our Thinking Process Template, we indicate at various points at which one can use prior knowledge (indicated by the glasses icon) to guide the analysis. For example, we can use prior knowledge to ask whether the reconstructed trajectories show the correct branching expected in the underlying biology of the differentiation and whether key gene-pairs that are known to be co-regulated are identified by scHOT (Ghazanfar et al., 2020), as well as asking which genes are differentially expressed across a pseudotime using tradeSeq (Van den Berge et al., 2020) and performing functional annotation of these gene sets through clusterProfiler (Yu et al., 2012). Together, these analyses demonstrate that the final trajectories are in line with our expectations and provide more confidence in the new biological insights extracted from these trajectories. The story includes other downstream analysis of the data such as cell-cell communication using CellChat (Jin et al., 2021) and RNA velocity analysis using scVelo (Bergen et al., 2020), which the users can perform to further explore their data.

Narrative for data story 2 - develop a PBMC biomarker model to predict COVID-19 patient outcomes

In our second data story, we aim to predict COVID-19 patient outcomes (mild or severe) from scRNA-seq data of peripheral blood mononuclear cells (Schulte-Schrepping et al., 2020) (Figure 5). Below, we list some questions and our thought process during the development of the story. Here, we showcase the story we derived on the COVID-19 patient outcome prediction. The story begins with the question of what models and input format we will use to build a prediction model (see Box 2). Here, we decided to use classical machine learning instead of deep learning, given the small sample size of 20 mild and 20 severe patients. We utilise scFeatures, a package that generates interpretable multiscale features from scRNA-seq data, such as cell-type proportions, pathway expression, ligand-receptor interactions and more. These features can then be used as input to facilitate an interpretable classification model. Once we have asked quality control questions as shown in the box and ensured the quality of generated features, we then used these features to build models to predict mild or severe outcomes.

Figure 5. Thinking process template for analysing a single-cell RNA-seq data for disease outcome classification.

The thinking process begins from processed data with cell type annotations and branches into two questions, each with a different focus. The top part focuses on using the disease classification model to extract biological insights into the disease, such as what features are important towards disease classification. The bottom part focuses on examining the model properties, such as whether the model is generalisable.

When building machine learning models, it is crucial that we ask questions on the model performance on a variety of models and metrics. Therefore, we choose to use ClassifyR, as it provides a user-friendly implementation on a number of common machine learning models and evaluation metrics. We created models for each feature type, resulting in a total of 13 models and compared the utility of these feature types for patient classification. We found that a support vector machine classifier consistently achieves a cross-validation accuracy over 0.7 (Figure 6A), demonstrating the usefulness of these features to classify disease outcomes.

Figure 6. Assessment of disease outcome classification accuracy using scFeatures' generated features.

(A) Shows the balanced accuracy of each feature type on classifying the mild and severe patients in the Schulte-Schrepping dataset. Models were run using five fold cross-validation with 20 repeats. For each feature type, the best model from the cross-validation was then selected and used to predict on the mild and severe patients in the Wilk dataset, as shown in (B). (C) Rankings of each feature in the feature type “gene mean celltype” across all cross-validated models.

Once the final models are obtained, we ask questions on the robustness of the models. One approach on this involves assessing the generalisability of the model on an independent dataset. We examined the performance of the 13 models on a different COVID-19 dataset obtained from the Wilk study that also contain mild and severe patients (Wilk et al., 2020). We found while the 13 models have close balanced accuracy between 0.75 and 0.88 on the Schulte-Schrepping dataset, their performance varied greatly in the Wilk dataset and ranged between 0.49 to 0.78 (Figure 6A, B). It is noteworthy that models built from feature types such as “gene proportion cell type” that have high accuracy do not necessarily maintain good performance on the Wilk dataset. On the other hand, the feature type “CCI” achieved an accuracy of over 0.75 in both datasets, indicating potential for further examination. Finally, to extract biological insights from the fitted models, we guide users to interpret the fitted models to identify important features and reflect on whether the features make sense. Here, it is important not to select the top features based on a single model, but to ask about the stability of these top features. To illustrate this idea, we examined all features that appeared at least once as top 10 features in the cross-validated models. Figure 6C highlights that while the majority of the features were consistently ranked as top features across all models, a proportion of features were ranked in the hundreds and thousands position in some models. These two scenarios illustrate the importance of critical thinking to avoid heading down a wrong decision path.