Differential expression and volcano

[antoine4ucsd]

Hello I am working on harmonized scRNA data integrated with Harmony and I want to look at DE between 2 groups of samples I followed this tutorial https://satijalab.org/seurat/articles/de_vignette and considered the default and pseudobulk strategy.

if I run the following

mydata_agg <- AggregateExpression(mydata, assays = "SCT", return.seurat = T, group.by = c( "pid", "hivstatus"))
mydata_agg$celltype.stim <- mydata_agg$status

mydata_agg.de <- FindMarkers(object = mydata_agg,
                                assay = "SCT",
                         ident.1 = "status1",
                         ident.2 = "status2",
                         test.use = "DESeq2")

the resulting Volcano plot seems shifted (i.e. log2 change not centered to zero). I assume it has to do with the scaling /log transformation but I am not sure what would be the optimal way to perform DE on these data. I also tried with RNA assay, default test (Wilcox) and without aggregation. all lead to different results.

any thoughts on the 'optimal' strategy to consider? sorry for the naive question and thank you in advance for your input! volcano

[pati-ni]

By default, Harmony is not touching your count data. It works with the embedding values and transforms them. The counts remain uncorrected and contain batch effects. Therefore, the phenotype you are reporting is related to the groups you are comparing and not the batch correction.

This is off-topic, but this behavior may be due to the different numbers of transcripts in the 2 groups you are comparing.

[antoine4ucsd]

thank you. really appreciated I did more investigation here https://github.com/satijalab/seurat/issues/8341

let me know if it makes sense

[pati-ni]

Yes, it makes sense. The fact that you see increased mitochondria counts in your control may be troublesome. I would definitely investigate more the QC parts of your analysis.

[antoine4ucsd]

thank you.