Simulation results

[alex-d13]

Hi all,

I spent the last days generating some results for the simulation. I used only quantiseq for deconvolution of the simulated datasets and for now mainly focused on the effect of scaling factors. One example can be seen here: I used the Travaglini-Dataset (the one with spike-in data) and compared the correlations with and without scaling factors. Basically one bar is the correlation of one scatter plot from here (the empty scatter-plots are cell-types which are not in the dataset; also, I compared Macrophages.M1 from quantiseq with the annotated Macrophages from Travaglini):

I feel like its pretty hard to analyze these results, since there is no obvious trend when using scaling-factors.. I am planing on doing this same analysis now for the 3 other datasets i was already using for the scaling-factor analysis (maynard, hao & vento-tormo).

I was wondering if you feel like I should also use different deconvolution tools like EPIC, or if that would be more part of the full benchmark?

[grst]

I wouldn't expect the correlation to vary a lot between the different scaling methods (In fact, I'm surprised it varies at all). What should change though, is the slope of the linear model fit. Could you make something like your first figure, but plot the slope instead of the correlation?

Trying out EPIC should be straightforward and can't harm. Would be great if you could try it out and you can then still decide if you include both in your thesis.

[FFinotello]

Hi Alex,

Very interesting!

Usually it is easier to see the bias or ots absence when you also plot all cell types together (coulored with the cell type-specific colours).

And do you have the scatterplots for all the three options? I would expect to see some differences for the DCs in terms of scaling rather than correlation.

A small note for quanTIseq: total CD4 T cells to be compared with the gold standard should be computed as CD4 + Tregs.

Also regarding quanTIseq, which parameter settings are you using?

Cheers, Francesca

On Sun, 7 Nov 2021, 12:32 Alexander Dietrich, @.***> wrote:

Hi all,

I spent the last days generating some results for the simulation. I used only quantiseq for deconvolution of the simulated datasets and for now mainly focused on the effect of scaling factors. One example can be seen here: [image: image] https://user-images.githubusercontent.com/40594032/140642763-2e07cf81-3ada-4db7-a3b9-193023d2605f.png I used the Travaglini-Dataset (the one with spike-in data) and compared the correlations with and without scaling factors. Basically one bar is the correlation of one scatter plot from here (the empty scatter-plots are cell-types which are not in the dataset; also, I compared Macrophages.M1 from quantiseq with the annotated Macrophages from Travaglini): [image: image] https://user-images.githubusercontent.com/40594032/140642861-3017b591-2cec-4c62-bd3a-3187c9daa294.png

I feel like its pretty hard to analyze these results, since there is no obvious trend when using scaling-factors.. I am planing on doing this same analysis now for the 3 other datasets i was already using for the scaling-factor analysis (maynard, hao & vento-tormo).

I was wondering if you feel like I should also use different deconvolution tools like EPIC, or if that would be more part of the full benchmark?

— You are receiving this because you are subscribed to this thread. Reply to this email directly, view it on GitHub https://github.com/omnideconv/simulator/issues/10, or unsubscribe https://github.com/notifications/unsubscribe-auth/AF6C62UDHLL6J53X6K6JGXLUKZPUJANCNFSM5HQW5DUA . Triage notifications on the go with GitHub Mobile for iOS https://apps.apple.com/app/apple-store/id1477376905?ct=notification-email&mt=8&pt=524675 or Android https://play.google.com/store/apps/details?id=com.github.android&referrer=utm_campaign%3Dnotification-email%26utm_medium%3Demail%26utm_source%3Dgithub.

[alex-d13]

I wouldn't expect the correlation to vary a lot between the different scaling methods (In fact, I'm surprised it varies at all). What should change though, is the slope of the linear model fit. Could you make something like your first figure, but plot the slope instead of the correlation?

Thats true, here's the same plot with slope values (its all 0 for Monocytes):

And do you have the scatterplots for all the three options?

A small note for quanTIseq: total CD4 T cells to be compared with the gold standard should be computed as CD4 + Tregs.

I will keep it in mind for the other datasets, but here there were no Trges annotated for this dataset.

Also regarding quanTIseq, which parameter settings are you using?

I am using immunedeconv to run it, like this: immunedeconv::deconvolute_quantiseq(simulation$pseudo_bulk, scale_mrna = F, tumor = F, arrays = F)

[FFinotello]

Hi Alex,

Many thanks! Two quick notes:

If you want to correct the mRNA bias with quanTIseq, the corresponding option should be set to TRUE.

As for the CD4, I meant you should add up conventional CD4 and Tregs to get total CD4 to be compared to the gold standard.

Cheers, Francesca

On Sun, 7 Nov 2021, 19:12 Alexander Dietrich, @.***> wrote:

I wouldn't expect the correlation to vary a lot between the different scaling methods (In fact, I'm surprised it varies at all). What should change though, is the slope of the linear model fit. Could you make something like your first figure, but plot the slope instead of the correlation?

Thats true, here's the same plot with slope values (its all 0 for Monocytes): [image: image] https://user-images.githubusercontent.com/40594032/140655631-76ef8ff4-c140-4ed5-869e-15cd48fb59c9.png

And do you have the scatterplots for all the three options?

[image: image] https://user-images.githubusercontent.com/40594032/140656651-9d95aaa0-83e2-4266-9417-1114d47f480b.png [image: image] https://user-images.githubusercontent.com/40594032/140656656-fcd9afe8-e21a-4dcb-8ef4-8d95f48592f0.png [image: image] https://user-images.githubusercontent.com/40594032/140656660-c7218a52-df98-4e36-8170-b076a6242dbb.png

A small note for quanTIseq: total CD4 T cells to be compared with the gold standard should be computed as CD4 + Tregs.

I will keep it in mind for the other datasets, but here there were no Trges annotated for this dataset.

Also regarding quanTIseq, which parameter settings are you using?

I am using immunedeconv to run it, like this: immunedeconv::deconvolute_quantiseq(simulation$pseudo_bulk, scale_mrna = F, tumor = F, arrays = F)

— You are receiving this because you commented. Reply to this email directly, view it on GitHub https://github.com/omnideconv/simulator/issues/10#issuecomment-962656771, or unsubscribe https://github.com/notifications/unsubscribe-auth/AF6C62XQOEOAW5QY6YVIBGTUK26QZANCNFSM5HQW5DUA . Triage notifications on the go with GitHub Mobile for iOS https://apps.apple.com/app/apple-store/id1477376905?ct=notification-email&mt=8&pt=524675 or Android https://play.google.com/store/apps/details?id=com.github.android&referrer=utm_campaign%3Dnotification-email%26utm_medium%3Demail%26utm_source%3Dgithub.

[alex-d13]

If you want to correct the mRNA bias with quanTIseq, the corresponding option should be set to TRUE.

But then I am correcting twice, right? Once with the scaling factor on the simulated dataset and then again when I use this dataset in quantiseq and it is corrected again by the quantiseq approach.

As for the CD4, I meant you should add up conventional CD4 and Tregs to get total CD4 to be compared to the gold standard.

Ah no I get it; ok I did, but it changed nothing here, as Trges is always 0 (except for 2-3 cases, where its at 0.002 or something)

[FFinotello]

In the simulation you should ADD the bias (not remove it).

Whereas a deconvolution methods might correct it (e.g. quanTIseq, EPIC, Monaco) or not (e.g. CIBERSORT).

On Sun, 7 Nov 2021, 19:39 Alexander Dietrich, @.***> wrote:

If you want to correct the mRNA bias with quanTIseq, the corresponding option should be set to TRUE.

But then I am correcting twice, right? Once with the scaling factor on the simulated dataset and then again when I use this dataset in quantiseq and it is corrected again by the quantiseq approach.

As for the CD4, I meant you should add up conventional CD4 and Tregs to get total CD4 to be compared to the gold standard.

Ah no I get it; ok I did, but it changed nothing here, as Trges is always 0 (except for 2-3 cases, where its at 0.002 or something)

— You are receiving this because you commented. Reply to this email directly, view it on GitHub https://github.com/omnideconv/simulator/issues/10#issuecomment-962660503, or unsubscribe https://github.com/notifications/unsubscribe-auth/AF6C62XBJH45H2QRLKKTK43UK3BW5ANCNFSM5HQW5DUA . Triage notifications on the go with GitHub Mobile for iOS https://apps.apple.com/app/apple-store/id1477376905?ct=notification-email&mt=8&pt=524675 or Android https://play.google.com/store/apps/details?id=com.github.android&referrer=utm_campaign%3Dnotification-email%26utm_medium%3Demail%26utm_source%3Dgithub.

[FFinotello]

I forgot to reply about Tregs.

Tregs are never changed.

What you change are the CD4 fractions (you have to add also the Tregs). So the CD4 T cell plot will change.

On Sun, 7 Nov 2021, 19:39 Alexander Dietrich, @.***> wrote:

If you want to correct the mRNA bias with quanTIseq, the corresponding option should be set to TRUE.

But then I am correcting twice, right? Once with the scaling factor on the simulated dataset and then again when I use this dataset in quantiseq and it is corrected again by the quantiseq approach.

As for the CD4, I meant you should add up conventional CD4 and Tregs to get total CD4 to be compared to the gold standard.

Ah no I get it; ok I did, but it changed nothing here, as Trges is always 0 (except for 2-3 cases, where its at 0.002 or something)

— You are receiving this because you commented. Reply to this email directly, view it on GitHub https://github.com/omnideconv/simulator/issues/10#issuecomment-962660503, or unsubscribe https://github.com/notifications/unsubscribe-auth/AF6C62XBJH45H2QRLKKTK43UK3BW5ANCNFSM5HQW5DUA . Triage notifications on the go with GitHub Mobile for iOS https://apps.apple.com/app/apple-store/id1477376905?ct=notification-email&mt=8&pt=524675 or Android https://play.google.com/store/apps/details?id=com.github.android&referrer=utm_campaign%3Dnotification-email%26utm_medium%3Demail%26utm_source%3Dgithub.

[alex-d13]

In the simulation you should ADD the bias (not remove it).

Yes, thats true, but I am doing that. Sorry, I probably just didn't think clearly yesterday :D I updated the plots with using the mRNA correction of quantiseq.

What you change are the CD4 fractions (you have to add also the Tregs).

Yes, thats also how I understood it. But I just wanted to say, that the deconvolution result for Trges is always 0 for this dataset (as there are no Trges annotated in the Travaglini dataset, so thats a good thing I guess). So I just add 0 to the CD4 fractions here.

[FFinotello]

Hey @alex-d13

Thanks for the explanations. Things are definitely clearer now and the Treg results are indeed good as you say.

I do see an mRNA bias in the "all cell" plots you shared (some cell types are clearly off from the diagonal), but I am surprised to see it in all settings.

Could you maybe make a violin plot of the number of expressed genes and total counts per cell type in the three scenarios (Census, spike-in, none)?

Thanks a lot, Francesca

[alex-d13]

I will do that. In the meantime, I also added the deconvolution results of epic & CIBERSORT:

[FFinotello]

Thanks for sharing, Alex. Maybe you could add also a method that does not correct for mRNA bias like CIBERSORT (absolute)?

[alex-d13]

Could you maybe make a violin plot of the number of expressed genes and total counts per cell type in the three scenarios (Census, spike-in, none)?

The number of expressed genes stays exactly the same; no gene will be set to 0 due to a scaling factor, since I only multiply a value with the cell vector. But I could try a different cutoff for the genes to be considered expressed?

Maybe you could add also a method that does not correct for mRNA bias like CIBERSORT (absolute)?

Updated above.

[grst]

Regarding the slopes, I'm a bit surprised that the differences are so small between the different scaling methods. Two more quick questions:

• How do you calculate the simulated TPM values (based on scaling factor and the single cell expression)
• Could you plot the actual scaling factors for the cell-types shown in the slope plot and the three scaling methods?

[alex-d13]

How do you calculate the simulated TPM values (based on scaling factor and the single cell expression)

First i get the scaling factor. Thats a vector with one value for each cell. I multiply my count matrix with this scaling vector. Then I generate a pseudo-bulk sample by sampling the cells i want based on the scenario; this gives me a new matrix with only these cells. I calculate the mean count value for each gene to get a vector again. This is now one pseudo-bulk sample. The sample is finally normalized by scaling the counts to 1e6. So basically i only do relative expression, not TPM normalization. But i do not really know how to get the gene length information into there, since I cannot expect that a dataset provides the gene lengths. Or could I just work with a pre-defined set of gene lengths here?

Could you plot the actual scaling factors for the cell-types shown in the slope plot and the three scaling methods?

[grst]

I have an idea what the problem might be: The scaling factor is somewhat proportional to the number of counts per cell. If you multiply the counts with the scaling factor the effects cancel each other out and you get something similar to just library-size normalized values.

I think you should multiply the normalized counts (or even better, TPM values, if you have them) with the scaling factors.

The number of counts per cell could then be treated as an additional scaling factor which should give you the same result as what is currently "NONE".

Also it could be nice to normalize the scaling factors such that they are better comparable.

@alex-d13, @FFinotello, does that make any sense?

[FFinotello]

Hey @alex-d13

Following up our intense discussion on zoom we could:

1. Select two single-cell datasets we can reasonably trust;

2. Consider the following simulation scenarios:

   • No added mRNA bias;
   • Multiplication by Census scaling factor;
   • Multiplication by spike-in scaling factor (when available);
   • Multiplication by extreme scaling factors (even applied just to a single cell type to simplify things);
   • Division by Census scaling factor (to test a real "no-bias" scenario.

3. Apply the following deconvolution approaches:

   • CIBERSORT
   • quanTIseq with no mRNA bias correction
   • quanTIseq with mRNA bias correction
   • EPIC with no mRNA bias correction
   • EPIC with mRNA bias correction.

Let's see if the results align with our expectations :D

[alex-d13]

Thank you Francesca for that summary!

Additionally what I will do is to rescale all the used scaling factors (census & spike-in) to 0-1 range. Also for the Travaglini dataset (since it has spike-ins I would keep on using it), I will normalize the counts from now on.

I'll try to get that done before the weekend so we can discuss the results soon :)

[alex-d13]

So, I found the bug that messed up all the results: i was applying the scaling factor on the genes, not cells (forgot to transpose the matrix). Here are the results for Travaglini [whitelisted cell-types which are also in quantiseq results: 2131 cells; SS2; lung] For the extreme scenarios i applied a scaling factor of 10 to Monocytes, Dendritic cells and B cells in 3 separate runs.

Will do the same for the CITE-seq dataset as well tomorrow, where we also might have more cells per cell-type.

[FFinotello]

Hi @alex-d13!

Thanks for sharing these results. Very nice work!

I will have a closer look tomorrow, but from the scatterplots I recognize the expected patterns now! I think we can remove the Census_div that was a sort of diagnostic test.

Curious to see the results on the other dataset as well!

Cheers, Francesca

[FFinotello]

Which exact cell types are you using to represent DCs?

[alex-d13]

Only "Dendritic_P1", which are only ~10 cells. But I could also add "Plasmacytoid Dendritic", "Myeloid Dendritic Type 2" and "IGSF21+ Dendritic", then I would have 138 cells in that group.

Also for Monocytes, I used "Classical Monocytes", but they would also have annotations for "Nonclassical Monocyte" & "Intermediate Monocyte", which only add ~30 cells though.

[FFinotello]

If you want to evaluate fairly quanTIseq, you should use only myeloid dendritic cells (definitely not plasmacytoid, that are very different). How many myeloid DCs do we have? EPIC has no DCs instead.

For monocytes, I would consider all of them together.

Of course, for the final benchmarking we have to discuss this and plan it wisely. But for now this should work well.

[alex-d13]

How many myeloid DCs do we have?

Only 7 unfortunately. So now 17 Dendritic cells in total.

But for the Hao dataset we have much more cells, it just takes forever to process them :D

[FFinotello]

Only 7 unfortunately. So now 17 Dendritic cells in total.

Is there a UMAP/t-SNE plot showing whether these cells cluster with other similar DC subtypes that are more abundant?

[alex-d13]

https://www.nature.com/articles/s41586-020-2922-4/figures/5 Here figure b is the only thing I found in their paper. But thats for the 10x data, I am only using the SS2 samples.

[alex-d13]

More results are ready! I decided to replace the census algorithm with the number of expressed genes per cell now. For this dataset census calculated "NaN" so many times for the cells (I also checked that again with the monocle implementation, there the same thing happens).

Hao [whitelisted cell-types: "Monocytes","T cells CD8", "T cells CD4","NK cells","T regulatory cells","B cells","Dendritic cells": 108k cells; CITE-seq; PBMC]

[FFinotello]

Very interesting, Alex. Thanks for sharing. Do you have also CIBERSORT results?

[alex-d13]

Yes, I updated the both results above.

[FFinotello]

Hey @alex-d13 would it be possible to make the same scatter plots (without extreme cases) considering only three cell types: DC, neutrophils, CD8?

[alex-d13]

Should I just remove the other cell-types from the plot or should I rather only simulate these 3 mentioned cell-types?

[FFinotello]

I meant doing just a separate simulation where only there three cell types are mixed. It's an oversimplification, but the plots should be then easy to read

[alex-d13]

Wow, I wish I had done this sooner, makes it much easier to see something :D

I think it's quite interesting how census (sim_census) for epic & quantiseq (with scaling) really improves the predictions, compare to using no scaling (sim_none)! Also overall I feel like census works quite well here. I guess for DCs this dataset just is too small (only 17 cells), thats why they dont really get detected. Also, as we expected, using the quantiseq scaling factors in the simulation (sim_qtseq) and then applying quantiseq deconvolution with mRNA scaling cancle each other out.

Here we have no Neutrophils in the dataset unfortunately. I think here its nice to see how all scaling factors improve the deconvolution for samples with low amounts of DCs/CD8. Also it looks like quantiseq is better(?) without mRNA scaling here?

Overall really cool to look at I feel like, would you suggest to maybe include this plot with one or two more cell-types? Or maybe do another one with B, NK and Monocytes?

[FFinotello]

Hello Alex,

Thank you so much for sharing this so quickly. I will get back to you ASAP. Meanwhile, could you maybe summarize briefly here how the different scaling factors were computed and how were the CPM or TPM scaled, so we are all on the same page?

Thanks a lot, Francesca

[FFinotello]

Also, may be worth reducing CIBERSORT y-axis to improve visibility

[alex-d13]

Scaling factors are always applied on the (normalized) expression matrix. Then I sample cells and take the mean of the expression values to get one pseudo-bulk sample.

Two datasets: Travaglini (manually TPM normalized)[508 cells with these 3 cell-types] Hao (manually CPM normalized)[14281 cells with these 2 cell-types]

• sim_none: multiply by 1 (no effect)
• sim_reads: multiply by number of reads from cell
• sim_census: multiply by census/1e6 (while census == number of expresses features)
• sim_spike: multiply by % of reads NOT mapped to spike-ins
• sim_epic: multiply all cells from cell-type X with scaling value of EPIC for cell-type X
• sim_qtseq: same as EPIC, but with quantiseq values

Also, may be worth reducing CIBERSORT y-axis to improve visibility

I already did that for the negative y-axis, somehow CIBERSORT had ~10 cases were it produces negative estimates? I will also decrease the positive y-axis.

[FFinotello]

Thanks!

For CIBERSORT, nu-SVR can produce negative estimates but CIBERSORT algorithm, to my understanding, should automatically set negative estimates to 0. Strange this is not happening here...

One more question: which DC types are you currently considering for the two studies?

[FFinotello]

Meanwhile, quick feedback from my side:

I am still not sure that considering simple mixes of few cell types is the best option as quanTIseq, for instance, is not well-designed for this. Maybe better considering more complex simulations and, then, only plotting results for a few cells, as Alex proposed.

Also, the type of DCs is very important (quanTIseq, for instance, is not able to disentangle pDC as its signature is based on mDC). I do not want to favor it, it is just that in that case, we would not see the mRNA bias for DCs.

Interesting patterns I see on the Travaglini data with quanTIseq (qs):

• none: qs works well, its scaled version corrects to much as there is no mRNA bias -> as expected
• census (i.e. mRNA bias): qs works badly, its scaled version works well -> as expected
• EPIC scores bias: qs bad, scaled version good -> expected considering that qs and EPIC scores correlate
• quanTIseq scores bias: qs bad, scaled version good -> expected
• spike-in: here the behavior is a bit strange... any idea why, considering that these scores should correlate with the census ones?

Hao data and EPIC results are somehow less clear, and CIBERSORT results not clearly visible.

A comment also for @grst: when I say one method works well, I do not expect that the slope is 1 but only that the cell-type-specific slopes are comparable (i.e. the points are aligned along the same line). To get a slope of 1 we should not only correct for the mRNA bias, but also get perfect deconvolution estimates which might not be always the case.

[alex-d13]

For Travaglini I use "regular" DCs(n=10) and myeolid DCs(n=7). They also have Plasmacytoid(n=13) and IGSF21+ DCs(n=8). For Hao I only use "regular" DCs(n=2652). They also have plasmacytoid DCs(n=861).

[alex-d13]

census (i.e. mRNA bias): qs works badly, its scaled version works well -> as expected

Also, when introducing the bias with the reads (sim_reads), we get almost the opposite of sim_none in terms of Neutrophils and T cells. So maybe the reads add too much bias? Same for EPIC as well.

[alex-d13]

I am still not sure that considering simple mixes of few cell types is the best option as quanTIseq, for instance, is not well-designed for this. Maybe better considering more complex simulations and, then, only plotting results for a few cells, as Alex proposed.

I think simulating only 2-3 cell-types gives clearer results for the effects of our scaling factors here, even though these simulations are not something a tools would see in practice (and so as you said, they are not designed for it).

This is the simulation with 8 cell-types (same as here: https://github.com/omnideconv/simulator/issues/10#issuecomment-966600691) but I only colored the three we are comparing right now (DC, Neutro, T CD8). Basically we see most of the effects you mentioned earlier, but not as pronounced.

[I changed the color palette, the linear model is now the red line and it is only based on the 3 colored cell-types]

[FFinotello]

A comment also for @grst: when I say one method works well, I do not expect that the slope is 1 but only that the cell-type-specific slopes are comparable (i.e. the points are aligned along the same line). To get a slope of 1 we should not only correct for the mRNA bias, but also get perfect deconvolution estimates which might not be always the case.

Regarding this, maybe the "slope plots" would be more informative as boxplots:

• Still one panel per method
• One "box" per simulation type (i.e., Census, none, etc.) on the x-axis
• Each "point" is a cell type (i.e., for each simulation you consider all cell types together)

[alex-d13]

Did you think of something like this? One box now contains the slopes for all cell-types for one simulation+deconvolution combination. I guess we are looking for boxes, which are as "narrow" as possible, so that the slopes do not differ as much between the cell-types.

[FFinotello]

Did you think of something like this? One box now contains the slopes for all cell-types for one simulation+deconvolution combination. I guess we are looking for boxes, which are as "narrow" as possible, so that the slopes do not differ as much between the cell-types.

Exactly! Maybe having also comparable x-axes and showing the single points (jitter) would facilitate the comparison

[alex-d13]

Maybe having also comparable x-axes

Then I would have to exclude cibersort, there the slopes are ~2-3 times higher then the rest

[FFinotello]

Maybe having also comparable x-axes

Then I would have to exclude cibersort, there the slopes are ~2-3 times higher then the rest

How would plots look if we consider CIBERSORT axis?

[alex-d13]

How would plots look if we consider CIBERSORT axis?

[FFinotello]

It's not bad at all. We could add a line at x=1 representing optimal results

[alex-d13]

Then we would have this one in the end. I tried to show the points in different shapes for the cell-types, but that gets a little bit too cluttered in my opinion:

[mlist]

This plot gives a really nice summary of the results. I'm wondering if it makes sense to suggest better correction factors from the simulated data, i.e. such that the slope 1 is exactly 1. We would need independent data to validate this then of course...

[alex-d13]

I'm wondering if it makes sense to suggest better correction factors from the simulated data, i.e. such that the slope 1 is exactly 1. We would need independent data to validate this then of course...

Mh what I could try out is to for example get scaling factors from the Hao dataset so that the slopes are close to 1 and then apply them on the Travaglini dataset. I think its not possible to get the slopes to 1 for all 5 deconv approaches, but maybe we pick one?

Btw, for Hao the slopes are much closer together than for Travaglini. I think that might be because of more cells in the Hao dataset: Also, I think for cibersort, the "perfect" slope would not be 1, right? Because cibersort gives estimates not in fractions but their own values. Would it make sense to scale them to the 0-1 range as well in order to have comparable slopes?