Thank you, really interesting thinking ! I dug a bit more into SCTransform and seems like the HVGs are not always the same, see below.
Warnings are the same in both cases, the same line is repeated except for warning 31 :
#Warnings
29: In theta.ml(y = y, mu = fit$fitted) : iteration limit reached
30: In theta.ml(y = y, mu = fit$fitted) : iteration limit reached
31: In sqrt(1/i) : production de NaN
32: In theta.ml(y = y, mu = fit$fitted) : iteration limit reached
33: In theta.ml(y = y, mu = fit$fitted) : iteration limit reached
Checking towards SCtransform
I have 7 different samples in SCRNARR.list
SCRNARR.list.LAB <- readRDS(file = "/mnt/raid1/Data/SCRNARR/SCRNARR_AfterSCT_LAB.rds")
object.size(SCRNARR.list.LAB)
2131213616 bytes
SCRNARR.list.HOME <- readRDS(file = "/mnt/raid1/Data/SCRNARR/SCRNARR_AfterSCT_HOME.rds")
object.size(SCRNARR.list.HOME)
2131212976 bytes
#For LAB
HVGs.LAB = list()
for (sample in 1:length(SCRNARR.list.LAB)) {
HVGs.LAB[[sample]] <- rownames(SCRNARR.list.LAB[[sample]]@assays$SCT@meta.features[SCRNARR.list.LAB[[sample]]@assays$SCT@meta.features$sct.variable,])
}
mean.list.LAB=list()
for (sample in 1:length(SCRNARR.list.LAB )) {
mean.list.LAB[[sample]] <- SCRNARR.list.LAB[[sample]]@assays$SCT@meta.features$sct.residual_mean
}
variance.list.LAB=list()
for (sample in 1:length(SCRNARR.list.LAB)) {
variance.list.LAB[[sample]] <- SCRNARR.list.LAB[[sample]]@assays$SCT@meta.features$sct.residual_variance
}
#For HOME
HVGs.HOME = list()
for (sample in 1:length(SCRNARR.list.HOME)) {
HVGs.HOME[[sample]] <- rownames(SCRNARR.list.HOME[[sample]]@assays$SCT@meta.features[SCRNARR.list.HOME[[sample]]@assays$SCT@meta.features$sct.variable,])
}
mean.list.HOME=list()
for (sample in 1:length(SCRNARR.list.HOME)) {
mean.list.HOME[[sample]] <- SCRNARR.list.HOME[[sample]]@assays$SCT@meta.features$sct.residual_mean
}
variance.list.HOME=list()
for (sample in 1:length(SCRNARR.list.HOME)) {
variance.list.HOME[[sample]] <- SCRNARR.list.HOME[[sample]]@assays$SCT@meta.features$sct.residual_variance
}
#Present in HOME but not in LAB
for (sample in 1:length(HVGs.HOME)) {
print(setdiff(HVGs.HOME[[sample]], HVGs.LAB[[sample]]))
}
character(0)
character(0)
[1] "Sntg1" "Tgfbr3l" "Apoa1"
character(0)
character(0)
character(0)
[1] "Arhgap29" "Gm36888" "Helz" "Nrg1" "Zfp369" "Zmym6"
#Present in LAB but not in HOME
for (sample in 1:length(HVGs.HOME)) {
print(setdiff(HVGs.LAB[[sample]], HVGs.HOME[[sample]]))
}
character(0)
character(0)
[1] "Cass4" "Galnt17" "Grid2"
character(0)
character(0)
character(0)
[1] "Aldh5a1" "Gm27162" "Gm36969" "Gm37452" "Ung" "Wfdc18"
for (sample in 1:length(mean.list.LAB)) {
print(sample)
print(sum(mean.list.HOME[[sample]]))
print(sum(mean.list.LAB[[sample]]))
}
[1] 1
[1] 302.9477
[1] 302.9477
[1] 2
[1] 373.7752
[1] 373.7752
[1] 3
[1] 203.9835
[1] 204.2463
[1] 4
[1] 270.4948
[1] 270.4948
[1] 5
[1] 154.8764
[1] 154.8764
[1] 6
[1] 247.2674
[1] 247.2674
[1] 7
[1] 196.836
[1] 198.1857
for (sample in 1:length(variance.list.LAB)) {
print(sample)
print(sum(variance.list.HOME[[sample]]))
print(sum(variance.list.LAB[[sample]]))
}
[1] 1
[1] 25587.02
[1] 25587.02
[1] 2
[1] 30130.43
[1] 30130.43
[1] 3
[1] 20651.4
[1] 20649.48
[1] 4
[1] 25024.74
[1] 25024.74
[1] 5
[1] 19100.73
[1] 19100.73
[1] 6
[1] 18098.05
[1] 18098.05
[1] 7
[1] 16157.47
[1] 16170.22
I only plot UMAPs for sample from 1 to 4 and the screenshot in the main thread is coming from the sample 3...
I finished the analysis for both LAB and HOME objects after SCTransform on my office workstation.
For the sample 3, LAB object outputted Africa-shaped cluster in the good way and for the HOME object I got Africa-shaped cluster in reverse
How can the SCTransform function choose different HVGs from the exact same input objects, parameters ?
• 0 views
•
link
written
by
Bastien Hervé
