ABSTRACT (original article): We presented a method to find potential cancer attractors using single-cell RNA sequencing (scRNA-seq) data. We tested our method in a Glioblastoma Multiforme (GBM) dataset, an aggressive brain tumor presenting high heterogeneity. Using the cancer attractor concept, we argued that the GBM’s underlying dynamics could partially explain the observed heterogeneity, with the dataset covering a representative region around the attractor. Exploratory data analysis revealed promising GBM’s cellular clusters within a 3-dimensional marker space. We approximated the clusters’ centroid as stable states and each cluster covariance matrix as defining confidence regions. To investigate the presence of attractors inside the confidence regions, we constructed a GBM gene regulatory network, defined a model for the dynamics, and prepared a framework for parameter estimation. An exploration of hyperparameter space allowed us to sample time series intending to simulate myriad variations of the tumor microenvironment. We obtained different densities of stable states across gene expression space and parameters displaying multistability across different clusters. Although we used our methodological approach in studying GBM, we would like to highlight its generality to other types of cancer. Therefore, this report contributes to an advance in the simulation of cancer dynamics and opens avenues to investigate potential therapeutic targets. CITATION (original article): Vieira Junior, M.G., de Almeida Côrtes, A.M., Gonçalves Carneiro, F.R. et al. A method for in silico exploration of potential glioblastoma multiforme attractors using single-cell RNA sequencing. Scientific Reports 14, 26003 (2024). https://doi.org/10.1038/s41598-024-74985-2​
GitHub: https://github.com/marcosgvjunior/Biomarker-Guided-scRNA-Seq-Cancer-Attractor-Analysis

General description

This notebook presents a methodology for exploring potential cancer attractors using single-cell RNA sequencing (scRNA-seq) data, with glioblastoma multiforme (GBM) as a case study. The approach combines biomarker-guided clustering and elliptical statistics to identify high-density regions in gene expression space by delineating confidence regions around clusters. To refine the identification of attractors, the pipeline estimates parameters supporting stability in these regions considering the stochastic dynamics of gene regulatory networks (GRNs). It also investigates multistability, revealing the potential for transitions between cellular states closely tied to cancer recurrence. While the methodology was developed for GBM, it is broadly applicable to other cancers, offering a versatile framework for characterizing tumor heterogeneity, probing dynamic cancer progression, and identifying biomarkers for therapeutic intervention.

Code Overview

The code is divided into two main components. The first component focuses on statistical analysis, where gene expression datasets are loaded. It employs clustering techniques such as k-means (centroid-based) and neighborhood contraction (density-based) to identify patterns in gene expression data. Using visualization functions, users can explore the outcomes of dimensionality reduction methods, inspecting gene combinations across varying dimensions and observing the resulting cluster formations. These visualizations also allow for examining specific gene combinations and their respective clusters. Additionally, confidence regions are defined to help identify potential attractors, providing a deeper understanding of the gene expression landscape. The second component addresses dynamic analysis. It simulates GRN dynamics using stochastic differential equations (SDEs) and estimates parameters via linear programming to ensure (multi)stability for various cluster combinations (centroids) under different regulatory function models. This section also includes functions for visualizing time series, histograms, and other relevant outputs of the dynamical simulations, offering insights into the stability and behavior of the identified attractors.

Input Data

The input includes gene expression datasets from multiple patients (datapoints_seurat_BT_ALL.xlsx, datapoints_seurat_BT_S1.xlsx, datapoints_seurat_BT_S2.xlsx, datapoints_seurat_BT_S4.xlsx, datapoints_seurat_BT_S6.xlsx, and datapoints_seurat_Regular.xlsx) and GRN adjacency matrices (act_BT_ALL_seurat.xlsx and sup_BT_ALL_seurat.xlsx).

Workspace Setup


Input GRN files, definition of marker genes, and confidence level


Defining Functions for Analysis


Part 1 - Clusters’ investigation using ellipsoid statistics

In[]:=
(*Loaddata*)​​loadAllData[EuclideanDistance,postotal];​​(*loadAllData[ManhattanDistance,postotal];*)​​(*Tablewithmeanandstandarddeviationsofeachcluster*)​​(*clustersAndMarkersStatistics[methodNBCdatabasinscentersTotal,methodNBCdatabasinssigmasTotal,postotal]*)​​​​listaDeDatasets={"methodNBCdataclustTotal","methodNBCdataclustBTS1","methodNBCdataclustBTS2","methodNBCdataclustBTS4","methodNBCdataclustBTS6","methodNBCdataclustAllRegularCells"};​​Print["\n"]​​Map[Print["Number of clusters - "<>#<>" :",Dimensions[ToExpression[#]]]&,listaDeDatasets];
Number of clusters - methodNBCdataclustTotal :{7}
Number of clusters - methodNBCdataclustBTS1 :{6}
Number of clusters - methodNBCdataclustBTS2 :{8}
Number of clusters - methodNBCdataclustBTS4 :{9}
Number of clusters - methodNBCdataclustBTS6 :{7}
Number of clusters - methodNBCdataclustAllRegularCells :{6}

Data and ellipsoids; 95%, 68%, and 20%; 40 genes and 3 genes


Verifying the unclustered points within ellipsoids


Ellipsoids visualization

elipsoidPlotsEuc={​​plotPointsAndElipsoid[orderedDataCluster[methodNBCdataclustBTS1,methodNBCdataclustTotal,postotal],methodNBCdataclustTotal,postotal,confidenceLevelDefaultValue],​​plotPointsAndElipsoid[orderedDataCluster[methodNBCdataclustBTS2,methodNBCdataclustTotal,postotal],methodNBCdataclustTotal,postotal,confidenceLevelDefaultValue],​​plotPointsAndElipsoid[orderedDataCluster[methodNBCdataclustBTS4,methodNBCdataclustTotal,postotal],methodNBCdataclustTotal,postotal,confidenceLevelDefaultValue],​​plotPointsAndElipsoid[orderedDataCluster[methodNBCdataclustBTS6,methodNBCdataclustTotal,postotal],methodNBCdataclustTotal,postotal,confidenceLevelDefaultValue],​​plotPointsAndElipsoid[methodNBCdataclustTotal,methodNBCdataclustTotal,postotal,confidenceLevelDefaultValue],​​plotPointsAndElipsoid[orderedDataCluster[methodNBCdataclustAllRegularCells,methodNBCdataclustTotal,postotal],methodNBCdataclustTotal,postotal,confidenceLevelDefaultValue]​​};​​​​plotsPartitionedEuc=Partition[elipsoidPlotsEuc,2];​​plotsPartitionedEucExport=Grid[plotsPartitionedEuc,Frame->None,Spacings->{0,0}]
Out[]=
A
B
C
D
E
F
G
H
A
B
C
D
E
F
G
H
A
B
C
D
E
F
G
H
I
A
B
C
D
E
F
G
A
B
C
D
E
F
G
A
B
C
D
E
F
G
H

Part 2 - GRN stochastic dynamics (timestep 0.05, Euler e SimTime 200)

Run

load data
loadAllData[EuclideanDistance,postotal];​​(*loadAllData[ManhattanDistance,postotal];*)​​(*---------------------------------------------*)
Par limits
parSinputs={0.5,4,0.5};​​parNinputs={1,4,0.5};​​parsLimits={1/100,10};​​​​selectedNoiseFactor=0.001;​​selectedTimeForSimulation=200;​​​​samplingsToRun=1;​​uniformSampling=0;​​pertubationAmplitude=0;​​combinationsToTest=1;(*allcombinationsmode*)​​​​AbsoluteTiming[aLLCOMBeuclideanListOfClutersCombinationsResult95=clustersCombinationsMoreThanTwo[methodNBCdataclustTotal,methodNBCdatabasinscentersTotal,methodNBCdatabasinssigmasTotal,methodNBCdataweightspercTotal,parSinputs,parNinputs,parsLimits,1,{1,2,3,4,5,6,7},samplingsToRun,uniformSampling,False,0.95,pertubationAmplitude,postotal,combinationsToTest,selectedNoiseFactor,selectedTimeForSimulation];]​​Dimensions[aLLCOMBeuclideanListOfClutersCombinationsResult95]
Total number of clusters combinations: 127
Out[]=
{20439.7,Null}
Out[]=
{3,127}
In[]:=
AbsoluteTiming[aLLCOMBeuclideanListOfClutersCombinationsResult68=clustersCombinationsMoreThanTwo[methodNBCdataclustTotal,methodNBCdatabasinscentersTotal,methodNBCdatabasinssigmasTotal,methodNBCdataweightspercTotal,parSinputs,parNinputs,parsLimits,1,{1,2,3,4,5,6,7},samplingsToRun,uniformSampling,False,0.68,pertubationAmplitude,postotal,combinationsToTest,selectedNoiseFactor,selectedTimeForSimulation];]​​Dimensions[aLLCOMBeuclideanListOfClutersCombinationsResult68]

BarChart

Grid of paired barcharts

Euclidean for parameters

Euclidean for combinations

Final totals bar chart

Part 2.1 - Testing sensitivity to initial conditions:

BarChart

Final totals bar chart

CITE THIS NOTEBOOK

A method for in silico exploration of potential glioblastoma multiforme attractors using single-cell RNA sequencing​
by Marcos Guilherme Vieira Junior, Adriano Maurício de Almeida Côrtes, Flávia Raquel Gonçalves Carneiro, Nicolas Carels & Fabrício Alves Barbosa da Silva
Wolfram Community, STAFF PICKS, December 20, 2024
​https://community.wolfram.com/groups/-/m/t/3342742