Scientist Stories: Samuel Alber, Random Subsampling Meets Confounders

My name is Samuel Albrecht. So unfortunately, Alvaro Mani couldn't make it, so I'm going to be presenting on his behalf. And today I'm going to be discussing a novel statistical framework based on random subsampling that allows for replicable single-cell analysis. So let's say we have some sort of single-cell dataset measuring the gene expression of each individual cell in our set of samples. So here we have a UMAP visualization of such a dataset. Here we can see the cells are clustered based on their cell type. And from here we can do a variety of downstream tasks. And one of the more simple ones that we can do is just take two groups of cells and try to determine which genes are differentially expressed between these two groups. So one example of this is we could just take one of these clusters. We can compare them to every other cluster to find the cell type markers of that cluster. But when doing this kind of analysis, one natural question that arises is how consistent are the signals that we detect with the true biological variation in our dataset? So to help quantify this, we made this graph on the left where each point is representing a statistically significant gene. On the x-axis is the log fold change of that gene in a large 962-sample dataset that we can effectively think as a ground truth. And then on the y-axis, we're taking 5 samples from 2 separate single-cell experiments respectively, and we're comparing the log fold change seen in those 5 samples with our ground truth. And we can see that even though we're just using 5 samples, we see very nice correlation. A different type of differential expression analysis that we can do is we can zoom in on one particular cluster. So let's say CD4 naive T cells., and we can group the cells by something like the age of the individuals that they come from. So the question is, do we still see this nice consistency between separate experiments? And the answer is no. And here we're actually increasing the number of samples from 5 to 50, and yet still we see very poor correlation for both of these separate single-cell experiments. And so why is this happening? Well, on the left, we're working in a regime where the signal is particularly potent. So we speculate that it dominates any sort of inter-individual variation. And so the signal that we detect is fairly homogeneous. However, on the right, this is a regime where inter-individual variation seems to dominate the signal from age. And so the signal that we see is fairly heterogeneous. And so small datasets are not as reliable. However, this becomes problematic when you consider that if you look at all the single-cell datasets published in the past 5 years and take the 90th quantile of their size, that is just 45 samples. And yet here we're demonstrating that for 50 samples, if we're examining a heterogeneous condition, we see very poor consistency. So to examine this further and why inter-individual variation matters, let's say that we have this 10-sample dataset. Where we have 5 cases, 5 controls. Each patient in the study has different types of individual-level types of variation. This could be their age, ethnicity, sex, and so on. For the purposes of this illustration, we're assuming that each patient just have 3 types of variation, which are colored in green, blue, and purple. Problems arise when there are imbalances in these types of variation. So for example, in our controls, we have 4 patients with this purple type of variation, whereas in our cases, we have just 1 patient with that type of variation. And this leads to a potential confounder in our analysis. So how might we solve this? Well, if we know what the confounders are and who's affected by it, then a simple thing to do would just be to subsample our controls and our cases in a manner that is balanced with respect to all these potential confounders. However, if we want to apply this to a real dataset, this becomes impractical because often our confounders are hidden, meaning we don't know what they are or much less who they're affecting. So our solution is to repeatedly randomly subsample our dataset, in this case m times, with expectation that at least some of these subsamples will be balanced with respect to all these potential confounders. In practice, we don't know what these m these particular subsamples will be. So we aggregate all the signals together. So for example, if we're examining something like differential expression or differential cell composition, where we're trying to see if there's some biological signal between our condition Y and our feature X, then we can take all the p-values from these m subsamples, take the maximum p-value, and note that if the maximum p-value is low, that means even in the balanced subsamples, we still see statistical evidence that our biological signal is real. Therefore, we can conclude that it is real. However, if the maximum p-value is high, then that means we cannot confidently conclude that this biological signal is real. In practice, a maximum function leads to low power. So to combat this, we use a quantile. The way this quantile is selected is that we adjust it until it achieves calibration under permutation of the case-control labels. We do this separately for each feature since these types of variation might affect different features differently. Based on this distribution of p-values from these different subsamples, we can calculate a final p-value with basically the assumption that when we select a particular quantile, under the null, everything at that quantile and above should be uniformly distributed from 0 to 1 if there's no signal.. And so the p-value at that specific quantile therefore will be the minimum of that set of uniform random variables, which can be modeled with a beta distribution. So first we apply this to differential expression analysis. We use just a simple t-test, nothing more sophisticated, since it has high sensitivity and low runtime. We compare 3 separate single-cell experiments to this large 962-sample dataset that we're using as a ground truth. For our heterogeneous condition, we're again using age, and then as a more positive control where we expect more methods to do well since there's a more dominant signal, we use sex and cell type markers. So if we look at, again, a dataset size of 50 and we compare the ROC curves of our method to all the other methods, we can see that our method highlighted in pink for age exceeds state-of-the-art performance.. And then for sex and cell type markers, it equals it. However, if we look at the false discovery rates of all these separate methods, we can see that our method in pink, which is on the far right of each one of these separate subplots, is the only method that is consistently calibrated across all three conditions. And to ensure that we weren't overfitting to just age as a heterogeneous condition, we also repeated this analysis with COVID-19 and lupus where we found a similar increase in performance. Performance, if not more. So one question that arises is, we demonstrated this for small datasets, but what happens if, as sequencing advances, datasets get larger? Does our method still become applicable? And so to test this, we took this 962-sample dataset and we split it 4 ways. And then between each pair of these 4 subdatasets, we use one as discovery and one as validation. And then calculated a Matthews Correlation Coefficient, or MCC, which you can think of as an improved version of an F1 score. And then across all of these different pairs, we just took the average MCC. And so we repeated this analysis for all 3 conditions and saw that for all 3 conditions, our method performs the best when compared to the next 3 best benchmark methods. Um, okay. So then we extended this to differential cell composition. So here we're doing a similar kind of analysis except using the sample-level cell type counts rather than the single-cell level gene counts. Here we're using a simple just chi-square test as our input into our procedure. And then since the number of features here is much smaller, we can just look up each feature in the literature and make a literature-curated ground truth. And then we took about 100 or so sub-datasets from 4 separate single-cell experiments and then calculated their empirical FDR versus their sensitivity. And we found that our method, this robustified chi-squared, which is kind of circled in red over there, has— is the only method that for a dataset size of greater than 10 is consistently calibrated. Furthermore, among the methods that have a similarly low FDR but aren't really calibrated, we still have a higher sensitivity. Then repeat this analysis for COVID-19 where we saw a much stronger signal, so more methods tended to be consistently calibrated. However, amongst those methods, we can see that our method circled in red is still the only method— sorry, it's the method with the highest sensitivity in this regime. So once we apply this to differential expression and differential cell composition, the question arose of what other downstream tasks can we apply this to? And one task that we focused on was contrastive learning, where we're trying to find the parameters of our model theta that maximize the ratio of the variance of our cases with respect to the variance of our controls. And if you just kind of optimize this directly, then you might see this nice enrichment or even separation. Of cases versus controls, seen in the red versus the gray. However, in reality, this might be driven by some sort of hidden confounder that has an incidental correlation with our condition of interest. And so to robustify this, you can first minimize with respect to all these subsamples to ensure that the parameters that we find lead to a high variance ratio even in the subsamples balanced with respect to all these confounders.. And in this manner, we can help ensure that these hidden confounders with incidental correlations are not driving the signal that we see. So we applied this to a dataset for COVID-19 where we applied it to CD4 naive T cells. Each row here is representing a separate contrastive learning method. Our method is on the top row and we benchmark it with 3 other methods shown below it. In the leftmost column, we're just coloring the cells based on cases versus controls. We can see that for our method, we have this nice enrichment or even separation between the cells. And then while doing this analysis, we also saw that some of the genes were correlated with the directions that we found, and these led to the identification of potential cell populations. And so the two columns next to it are showing the activation of those 2 genes which are attributed to CD4 naive T cell activation. And although some of the other methods also seem to capture some structure in our dataset, if we take the representations learned from each of these methods and compare them to a completely different dataset, also for COVID-19, then we can see that although the other methods tend to not really replicate the structure that was seen in the training dataset, our method does. In particular, the case control structure as well as the structure of these subpopulations. We also repeated this analysis for other cell types, which I didn't include here just for simplicity. But in this manner, we're essentially allowing our method to find the variation of interest in this heterogeneous condition. So in conclusion, we developed this framework based on random subsampling that allows for more robust analysis for separate types of single-cell analyses. One big kind of advantage of using this procedure is that you can just input an existing statistical model into our procedure. In fact, we used something as simple as a t-test to show that under our procedure, it exceeds state-of-the-art performance, and also something as simple as just a chi-square test, which also exceeded state-of-the-art performance for differential cell composition. And then, as you just saw in the previous slide, if you apply this same type of idea to contrastive learning, you can robustify it so that it's more replicable across separate experiments. And then in that manner, you can isolate the variation of interest that you're interested in. So in conclusion, I would like to acknowledge the rest of my team, in particular Professor Ali Rahmani, as well as Professor Nir Yosef, who mentored me throughout the course of this project. And I'd like to thank the audience for taking the time to listen to my presentation.

Questions for Sam. Hey, thanks for a cool talk. So you've got this really nice method to account for sort of finite sample effects in a much better way, but you're comparing it to methods like DESeq and LUMA, which really get their power from pooling information across features. So is there a good way to combine both of those good things?

Yeah, I think— is the microphone on, by the way? I think I can't hear Joshua. Yeah, actually, yeah, can you repeat the first half of that question?

Yeah, so I was saying that you have this cool method to account for finite sample effects in a better way, but the comparisons you're making are to methods like DESeq and LIMMA where the power is really coming from pooling information across features. So can you do both of those together at the same time?

Yeah, I guess we haven't experimented with that too much. I mean, you could. Apply a method like that within this subsampling framework. And I think we did try that, although we found that due to power reasons, it didn't perform very well. So that's like one of the reasons why we use something like a t-test, since it's overly sensitive. But yeah, it seems like pulling cost features for now doesn't seem to kind of lead to increasing performance. But this is still unpublished, so that might change in the future.

Thanks, great talk. How does this compare with an explicit confounder correction, like if you just included, started to include PCs into your analysis?

Right, yeah, so I guess for some of these methods, we're kind of explicitly modeling for these different confounders that are present in the dataset. And as we can see, the performance in those cases tends to not do as well. And also, as kind of mentioned in kind of like the introduction to this problem, A big issue arises when a lot of these confounders are hidden, and so you can't really explicitly model for them.

Well, but what if you just include the principal components? Oh, oh, I see. Sorry. So you do some kind of naive or complicated hidden confounder method.

Right. Yeah. So I guess we— I don't think we've tried that yet. So, yeah, perhaps kind of incorporating those principal confounders might also lead to kind of a similar type of performance. But I guess for now, no method that we know of has been made to kind of demonstrate this. So we, for now, haven't benchmarked it with our method.

Thank you. Thanks for the talk. I have a quick question on how do you select the quantiles? Like, it seems for extremely imbalanced data, you need to select more stringent quantile, the other word. But is there a data-driven way that given the dataset you have, select automatically select quantile?

Yeah, so the way we select quantiles is actually based on the data. What we do is we essentially permute the case control labels where we're assuming that under permutation, we shouldn't see some sort of signal. And we're driving up the quantile until the kind of p-value that we're calculating based on the beta distribution that I described is below 0.05 or whatever metric that we're using. And so this does tend to lead to pretty stringent quantiles that are near the maximum function, but they tend to be at least slightly lower, which leads to an increase in power. Okay, let's give Sam another round of applause. Okay, so next up we have Gong Li, who will be talking about pseudotime analysis with single-cell sequencing and imaging data.

He triggers.