Cluster Analysis of Single-Cell Gene Data

Question

I have a dataset of 480 cells (observations) by 57 variables (48 are gene expression levels). Multiple cells (approximately 24) are sampled from within each patient. Therefore, I have have a variable for patient ID, 48 genes, and clinical information for each patient. My main exposure of interest is location by which the cell was sampled. I have run a series of mixed models with each gene expressivity as outcome to determine if location is significant while controlling for other factors. However, due to the single-cell nature of the data, I am interested in the 'landscape' of expressivity. Therefore, I would like to do unsupervised clustering to determine possible clusters of cells based on gene expressivity. Of course, cells sampled from the same individual are going to naturally cluster together. I have visualized the data using MDS: (colors identify patient ID). There appears to be 4 or 5 major clusters despite patient ID. In order to identify which method and number of clusters would be best for the data, I looked at the internal validity and stability using the clValid package in R: . Lower stability scores are best. The best clustering algorithms and number of clusters were determined to be: 2 clusters using hierarchical agglomerative, and 9-10 clusters with SOM. So I decided to do both and visualize the results in MDS space. First I used SOM on the data and cut the distance tree between nodes at 5 (also k-means WSS plot for each cluster number also has an elbow at 5), 9, and 10 (5 because of the biggest increase in stability from 4 to 5 in the clValid plots).

I can't post further images of the results do to reputation points, but SOM 5 and 9 do not cluster into the intuitive groups seen in the MDS plot colored by Patient ID. Hierarchical clustering at 2 also doesn't look good (sorry wish I could post the plots). So my questions are:

1. Given grouped data such as this, do unsupervised clustering models exist which compensate for random effects for patients? I.E. Is there a way of looking at clustering due to gene expression levels while compensating for correlated data units within patients?
2. SOM provides results which do not cluster based on patient (MDS colorized by SOM looks very different to MDS colorized by patient); however, MDS colorized by k-means looks like it clusters more on patients.
3. In the clValid literature, it states that scores which are consistently low are the best, and this is how clustering methods are chosen. However, in this case - there seems to be a 'lock-in-effect' for stability at 5 groups. I imagine if more clusters are provided than that which exists in the dataspace, this would decrease stability. Can anyone comment on this?

Overall, my question is regarding the best way of finding clusters in this type of data where single cells are sampled from within patients and are thus correlated.

Thank you for your time,

Answer 1

Given grouped data such as this, do unsupervised clustering models exist which compensate for random effects for patients? I.E. Is there a way of looking at clustering due to gene expression levels while compensating for correlated data units within patients?

You could try Rahul Satija's recent batch-effect removal method, which is based on CCA. It's intended for RNA-seq data, not qPCR, but it might still work. In the linked paper below, he uses it to cluster cells by type while avoiding what would ordinarily be a massive effect due to human vs mouse samples.

https://www.biorxiv.org/content/early/2017/07/18/164889

SOM provides results which do not cluster based on patient (MDS colorized by SOM looks very different to MDS colorized by patient); however, MDS colorized by k-means looks like it clusters more on patients.

No comment. Sorry!

In the clValid literature, it states that scores which are consistently low are the best, and this is how clustering methods are chosen. However, in this case - there seems to be a 'lock-in-effect' for stability at 5 groups. I imagine if more clusters are provided than that which exists in the dataspace, this would decrease stability. Can anyone comment on this?

Typically, stability will decrease as you add clusters. It's best to adjust for this somehow, for example by using a uniform distribution as a reference. One easy alternative is to look for a "notch" in the stability curve: a place where it should go down, but doesn't. Here's a great review on stability for cluster model selection.

https://arxiv.org/abs/1007.1075

Answer 2

You can separate the gene expression matrix from the metadata variables (such as patients). And use patient IDs as batch ID to do batch correction. There are multiple algorithms to do that such as harmony, Seurat3 integration, bbknn, etc.

Of course you can also do a 2D hierarchical clustering and visualize the heatmap of gene expression to spot the nature of patient-specific gene expression patterns since you have a small matrix.