I was hoping someone might be able to help me conceptually understand something I find bizarre. I have scRNA data from 6 healthy controls and 12 IBD patients. I ran PySCENIC in two different ways, (1) on the entire set of plasma cells from both healthy and diseased donors; and (2) separately on the IBD donor set and on the healthy set. In each case there were 50 replicates, so I have 50 replicates for the former case and 100 total (50 each) in the latter case.
I have a particular transcription factor of interest in mind, so I focus on this TF. So in the separated case, as I wanted to see how important the TF is I ranked the median AUC value of regulons where the TF is the regulator as compared with the median AUCs I get for other regulons. I find that this TF is super important in IBD (ranks range from like 1-6, with the vast majority of replicates having it rank 3 or higher), while in healthy controls the ranks are way lower down and vary a lot more, say 15-30. So I thought, ok, thats a pretty strong signal.
But then, I go to my combined HC+IBD SCENIC results, and subset the AUCell values I get based on whether these are from IBD or healthy donors (in other words, it is only now that I divide the cells into the two cohorts). I still find in terms of ranks that the TF regulon ranks higher (i.e. more important) for IBD compared to control like 95% of the time, so still a strong signal in that sense, but the signal is way more diluted. Now there are no replicates where it is ranked #1-3, but almost all will be #12-20. And the differential in ranking compared to healthy will be way smaller as well.
My question is, can someone help me wrap my head around what is going on, and in particular, why does my TF of interest appear to be less active compared to other regulons when I divide cells into cohorts only after running SCENIC as opposed to before? Because crucially, we are talking about the exact same cells in either case, with the exact same transcriptomes.
It may be relevant to note, in the separated case there were 106 regulons total for IBD and 487 for healthy, while in combined there were about 280 regulons. So of course the denominator will affect the rankings of a particular regulon. But still, the signal is so so strong run one way and not the other that I dont think this explains it.
Thanks so much in advance.
1 answer
The regulons aren't the same objects in the two runs, and that's most of the story. GRNBoost2 learns co-expression modules from whatever matrix you hand it, so running on IBD cells alone builds your TF's regulon out of genes that co-vary with it within IBD - and then scoring IBD cells highly with that regulon is partly guaranteed. It's the same circularity as selecting features on your data and then testing on the same data. The joint run has to build one regulon that works across both groups, so it comes out a broader, blurrier gene set, and the weaker signal you get from it is the more honest number.
I'd also drop the ranking metric. AUCell scales with regulon size, so a regulon's rank against the others depends on that run's regulon size distribution - and yours differ a lot (106 vs 487 vs 280). Ranks aren't comparable across runs even in principle, so that comparison was never going to behave.
What you want lives inside the joint run: test directly whether that regulon's AUCell is higher in IBD cells than in HC. With 6 and 12 donors, aggregate to per-donor means first, otherwise your n is cells rather than people.
Log in to answer this question.