Marker genes disappeared in integrated dataset with SCTransform

[stanaka6]

Hi Seurat team,

Could you please explain possible reasons why the listed genes from FindConservedMarkers function are different between A. integrated data normalized without SCTransfrom and B. with SCTransform?

Shortly, most of the potential marker genes disappeared in my SCTransformed integrated data. For instance, while gh1 is a significant marker gene for pituitary somatotropes, gh1 was not listed as a conserved marker in integrated data with SCTransform. When I did analyses of each dataset separately (normalized data with SCTransform), I was able to obtain the potential marker genes including gh1 using FindAllMarkers function.

When I checked the expression level using AverageExpression function, in RNA assay, gh1's value was "4840.174" in cluster0 of A. integrated data without SCTransform, and "Inf" in that of B. SCTransform data. Also, cluster 1 to 11 also showed "Inf" in data B. I guess I might have missed something. I would like to use SCTransform to remove the batch effects and get reliable results for my experiments.

I would appreciate it if anyone could provide me any suggestions or thoughts.

Thank you!

---

The data description and scripts are below:

My Dataset:

1. WT: WT sample (15299 features across 2366 cells after QC)
2. KO: Mutant sample (15817 features across 5046 cells after QC)

#### A. Integration without SCtransform normalization ####

# Normalization and find variable features 
WT <- NormalizeData(WT, normalization.method = "LogNormalize", scale.factor = 10000)
KO <- NormalizeData(KO, normalization.method = "LogNormalize", scale.factor = 10000)
WT <- FindVariableFeatures(WT, selection.method = "vst", nfeatures = 10000)
KO <- FindVariableFeatures(KO, selection.method = "vst", nfeatures = 10000)

# Merge ----
x <- merge(x = WT, y = KO, add.cell.id = c("WT", "KO")) # Merge
[email protected]$cells <- rownames([email protected]) # Add cell IDs  and sample (WT or KO) column to metadata
[email protected]$sample <- NA
[email protected]$sample[which(str_detect(x$cells, "^WT_"))] <- "WT"
[email protected]$sample[which(str_detect(x$cells, "^KO_"))] <- "KO"

# Integrate ----
list <- SplitObject(x, split.by = "sample")
features <- SelectIntegrationFeatures(object.list = list)
anchors <- FindIntegrationAnchors(object.list = list, anchor.features = features)
integ <- IntegrateData(anchorset = anchors)

# Clustering ----
DefaultAssay(integ) <- "integrated"
integ <- ScaleData(integ)
integ <- RunPCA(integ)
integ <- FindNeighbors(integ, dims = 1:20)
integ <- FindClusters(integ , resolution = 0.2)

# Find markers ----
DefaultAssay(sct) <- "RNA"
C0.markers <- FindConservedMarkers(sct, ident.1 = 0, grouping.var = "sample", verbose = FALSE)

head(C0.markers)

## WT_p_val WT_avg_log2FC WT_pct.1 WT_pct.2  WT_p_val_adj KO_p_val KO_avg_log2FC KO_pct.1
## ak5l    4.203177e-287      2.172704    0.930    0.165 6.756606e-283        0      1.690238    0.861
## plp1b   5.330965e-284      2.513617    0.853    0.099 8.569526e-280        0      2.248159    0.809
## gh1     1.101822e-256      7.602656    0.987    0.877 1.771179e-252        0      5.636869    0.999
## fabp11a 1.117020e-237      1.961364    0.566    0.010 1.795610e-233        0      2.178432    0.593
## krt15   2.946687e-197      1.496218    0.476    0.008 4.736800e-193        0      1.778407    0.523
## ecrg4a  3.569782e-181      2.927163    0.493    0.023 5.738424e-177        0      3.375706    0.612
## KO_pct.2 KO_p_val_adj      max_pval minimump_p_val
## ak5l       0.164            0 4.203177e-287              0
## plp1b      0.118            0 5.330965e-284              0
## gh1        0.978            0 1.101822e-256              0
## fabp11a    0.018            0 1.117020e-237              0
## krt15      0.017            0 2.946687e-197              0
## ecrg4a     0.056            0 3.569782e-181              0

#### B. Integration with SCtransform normalization ####

# SCTransform ----
WT <- SCTransform(WT, method = "glmGamPoi", 
                  vars.to.regress = c("percent.mt", "S.Score", "G2M.Score"),
                  verbose = FALSE)

KO <- SCTransform(KO, method = "glmGamPoi", 
                  vars.to.regress = c("percent.mt", "S.Score", "G2M.Score"),
                  verbose = FALSE)

# Merge ----
x <- merge(x = WT, y = KO, add.cell.id = c("WT", "KO")) # Merge
[email protected]$cells <- rownames([email protected]) # Add cell IDs  and sample (WT or KO) column to metadata
[email protected]$sample <- NA
[email protected]$sample[which(str_detect(x$cells, "^WT_"))] <- "WT"
[email protected]$sample[which(str_detect(x$cells, "^KO_"))] <- "KO"

# Integrate ----
x <- SplitObject(x, split.by = "sample")
features <- SelectIntegrationFeatures(object.list = list, nfeatures = 10000)
list <- PrepSCTIntegration(object.list = list, anchor.features = features)
anchors <- FindIntegrationAnchors(object.list = list, normalization.method = "SCT", anchor.features = features)
sct <- IntegrateData(anchorset = anchors, normalization.method = "SCT")

# Clustering ----
sct <- RunPCA(sct)
sct <- FindNeighbors(sct, dims = 1:30)
sct <- FindClusters(sct, verbose = FALSE, resolution = 0.3)

# Find markers ----
DefaultAssay(sct) <- "RNA"
C0.markers <- FindConservedMarkers(sct, ident.1 = 0, grouping.var = "sample", verbose = FALSE)

head(C0.markers)

##            WT_p_val WT_avg_log2FC WT_pct.1 WT_pct.2  WT_p_val_adj KO_p_val KO_avg_log2FC KO_pct.1
## ak5l    2.572492e-285     19.542075    0.946    0.186 4.135281e-281        0      7.149989    0.868
## fabp11a 1.081340e-271    -86.041636    0.630    0.012 1.738253e-267        0     19.552099    0.613
## plp1b   2.542426e-251     21.972940    0.847    0.126 4.086950e-247        0     14.276752    0.813
## krt15   2.757512e-219     16.942512    0.525    0.011 4.432701e-215        0     12.981604    0.541
## kcnk2a  6.744695e-216      6.626783    0.820    0.148 1.084210e-211        0      3.279256    0.597
## ecrg4a  4.161518e-214    340.992522    0.554    0.023 6.689641e-210        0    182.093992    0.632
##        KO_pct.2 KO_p_val_adj      max_pval minimump_p_val
## ak5l       0.171            0 2.572492e-285              0
## fabp11a    0.019            0 1.081340e-271              0
## plp1b      0.126            0 2.542426e-251              0
## krt15      0.018            0 2.757512e-219              0
## kcnk2a     0.096            0 6.744695e-216              0
## ecrg4a     0.057            0 4.161518e-214              0

gh1 is gone. And similar things happened in all other clusters.

[saketkc]

In your second workflow, you switch back to RNA assay to perform DE, but you need to normalize the data first

DefaultAssay(sct) <- "RNA"
sct <- NormalizeData(sct)
C0.markers <- FindConservedMarkers(sct, ident.1 = 0, grouping.var = "sample", verbose = FALSE)

This will also explain the super high (or low) log fold changes that you see.

Alternatively, if the library sizes across WT and KO are comparable you can also use the data slot of SCT assay to perform DE.

[stanaka6]

Thank you very much for your reply. I was able to identify DEGs (including proper marker genes) by following your answer. So, if I use the data slot of the SCT assay, the values are already normalized and I don't need to perform NormalizeData followed by performing DE, is that correct? Also, I do not fully understand why SCT assay works better to perform DE when the library sizes between the groups are similar. I would appreciate it if you could explain the reason or provide me any resources to make it more clear.

Thanks!

[saketkc]

Sure. Does the reply in this comment help answer your question?

[stanaka6]

Very clear answer. Thank you so much!