Seurat: warnings from SCTransform

[zbgaber]

Hello,
I am trying to use SCTransform() on some data I've previously successfully analysed with Seurat 2 and the standard log-transformation normalisation. When I run the command:

har72a <- SCTransform(object = har72a, vars.to.regress = c("percentMito","Sex"), verbose = TRUE)

I get the following output

Calculating cell attributes for input UMI matrix
Variance stabilizing transformation of count matrix of size 15614 by 4093
Model formula is y ~ log_umi
Get Negative Binomial regression parameters per gene
Using 2000 genes, 4093 cells
  |==================================================================================================================================================| 100%
Found 7 outliers - those will be ignored in fitting/regularization step

Second step: Get residuals using fitted parameters for 15614 genes
  |==================================================================================================================================================| 100%
Computing corrected count matrix for 15614 genes
  |==================================================================================================================================================| 100%
Calculating gene attributes
Wall clock passed: Time difference of 1.743396 mins
Placing corrected UMI matrix in counts slot
Determine variable features
Set 3000 variable features
Regressing out percentMito, Sex
  |==================================================================================================================================================| 100%
Centering data matrix
  |==================================================================================================================================================| 100%
Setting default assay to SCT
There were 50 or more warnings (use warnings() to see the first 50)

The warnings are:

1: In theta.ml(y = y, mu = fit$fitted) : iteration limit reached
2: In theta.ml(y = y, mu = fit$fitted) : iteration limit reached
3: In theta.ml(y = y, mu = fit$fitted) : iteration limit reached

Do you know what this problem means? Other people at my institute that I've spoken to are having the same issue with SCTransform and none of us are sure what is causing the issue or what it might mean. I don't think it has anything to do with the regression as I removed the vars.to.regress argument completely and still got the same error. The dataset comes from two 10x runs, combined with cellranger aggr, and imported into Seurat with your Read10X() function. I also added a some metadata but otherwsise the Seurat object comes straight from CreateSeuratObject using default parameters.

Thanks for any help you have and so far loving Seurat 3.0 - the cluster distinctions appear to be much sharper with the addition of anchors.

Zachary

[ChristophH]

Hi Zachary,
These warnings are related to the estimation of the theta parameter of the Negative Binomial regression model used to describe the molecule count distribution of each gene. When a gene has very few non-zero observations estimation of that parameter is not reliable and that's what the message indicates. This is not a problem, since the estimated theta parameters are not used directly for normalization but are regularized (sharing information across genes and thus lowering uncertainty).

[ahoffrichter]

Hello,
I have a similar problem using sctransform on my data.
Using this command:
SCTransform(ctrl, vars.to.regress = c("S.Score", "G2M.Score"), verbose = TRUE)

I get the following messages returned:

Calculating cell attributes for input UMI matrix
Variance stabilizing transformation of count matrix of size 18456 by 2230
Model formula is y ~ log_umi
Get Negative Binomial regression parameters per gene
Using 2000 genes, 2230 cells
|===========================================================================================================| 100%
iteration limit reached

"iteration limit reached" is repeated a couple of hundreds of times.

Found 2 outliers - those will be ignored in fitting/regularization step

Second step: Get residuals using fitted parameters for 18456 genes
|===========================================================================================================| 100%
Computing corrected count matrix for 18456 genes
|===========================================================================================================| 100%
Calculating gene attributes
Wall clock passed: Time difference of 36.10467 secs
Determine variable features
Set 3000 variable features
Place corrected count matrix in counts slot
Regressing out S.Score, G2M.Score
|===========================================================================================================| 100%
Centering data matrix
|===========================================================================================================| 100%
Set default assay to SCT
An object of class Seurat
38414 features across 2230 samples within 2 assays
Active assay: SCT (18456 features)
1 other assay present: RNA

I thought it worked anyways because @ChristophH said "This is not a problem" and because I got the message "Active assay: SCT...".
But when I run:
ctrl <- RunPCA(object = ctrl, npcs = 30, verbose = FALSE)
it returns:

Error in PrepDR(object = object, features = features, verbose = verbose) : Variable features haven't been set. Run FindVariableFeatures() or provide a vector of feature names.

And when I call my seurat object ctrl it returns:

An object of class Seurat
19958 features across 2230 samples within 1 assay
Active assay: RNA (19958 features)

indicating that the SCT data was not saved in my seurat object.

Does that mean using sctransform with this data won't work and I just have to use the usual linear regression (i.e. NormalizeData(), FindVariableFeatures(), ScaleData() )?

Edit: I'm using Seurat_3.0.2

[davisidarta]

Hello,
I have a similar problem using sctransform on my data.
Using this command:
SCTransform(ctrl, vars.to.regress = c("S.Score", "G2M.Score"), verbose = TRUE)

I get the following messages returned:

Calculating cell attributes for input UMI matrix
Variance stabilizing transformation of count matrix of size 18456 by 2230
Model formula is y ~ log_umi
Get Negative Binomial regression parameters per gene
Using 2000 genes, 2230 cells
|===========================================================================================================| 100%
iteration limit reached

"iteration limit reached" is repeated a couple of hundreds of times.

Found 2 outliers - those will be ignored in fitting/regularization step
Second step: Get residuals using fitted parameters for 18456 genes
|===========================================================================================================| 100%
Computing corrected count matrix for 18456 genes
|===========================================================================================================| 100%
Calculating gene attributes
Wall clock passed: Time difference of 36.10467 secs
Determine variable features
Set 3000 variable features
Place corrected count matrix in counts slot
Regressing out S.Score, G2M.Score
|===========================================================================================================| 100%
Centering data matrix
|===========================================================================================================| 100%
Set default assay to SCT
An object of class Seurat
38414 features across 2230 samples within 2 assays
Active assay: SCT (18456 features)
1 other assay present: RNA

I thought it worked anyways because @ChristophH said "This is not a problem" and because I got the message "Active assay: SCT...".
But when I run:
ctrl <- RunPCA(object = ctrl, npcs = 30, verbose = FALSE)
it returns:

Error in PrepDR(object = object, features = features, verbose = verbose) : Variable features haven't been set. Run FindVariableFeatures() or provide a vector of feature names.

And when I call my seurat object ctrl it returns:

An object of class Seurat
19958 features across 2230 samples within 1 assay
Active assay: RNA (19958 features)

indicating that the SCT data was not saved in my seurat object.

Does that mean using sctransform with this data won't work and I just have to use the usual linear regression (i.e. NormalizeData(), FindVariableFeatures(), ScaleData() )?

Edit: I'm using Seurat_3.0.2

I happen to have the same issue. I'm also using Seurat_3.02. I'm running on R version 3.5.3.

[andrewwbutler]

@ahoffrichter, did you assign the output of SCTransform to ctrl? The command you provided doesn't make any assignment so therefore ctrl won't have changed from before (and won't have the SCT assay).

[kmshort]

Has anyone managed to work out how to not have sctransform spit out heaps of "iteration limit reached" warnings in R markdown console and "## Warning in theta.ml......" in knitr output? A knitr rendered PDF of mine just output 487 pages of "## Warning in theta.ml...." messages.

[jagregory17]

I'm getting the same "iteration limit reached" just following the spatial vignette. Anybody have any luck with this?

nevermind, found this - https://github.com/satijalab/seurat/issues/2578