STOmics / Stereopy

A toolkit of spatial transcriptomic analysis.
MIT License
197 stars 65 forks source link

SpatialFeaturePlot() fails to generate the plot of Seurat object converted from Stereopy h5ad #126

Closed limin321 closed 11 months ago

limin321 commented 1 year ago

I first read XX_tissue.GEF file using StereoPy and save it to seurat compatible h5ad like this adata = st.io.stereo_to_anndata(data, flavor="seurat", output = "seurat_out.h5ad"). Then I ran Rscript annh5ad2rds.R --infile <h5ad file> --outfile <rds file> to generate the RDS file. It actually generates 2 files, one is xx.RDS which includes an image, the other is the h5seurat object which doesn't have an image. Am I right?

When I read this RDS into R using seurat, it doesn't have nFeature_spatial and nCount_spatial columns. Instead, it has x and y columns, which I think is coordinate for each spot? I have no issue running Vlnplot(obj, features = c("nFeature_spatial", "nCount_spatial")). However, the SpatialFeaturePlot() is not working. It only says writing error.

Could you please help me troubleshoot the SpatialFeaturePlot() issue? Thanks in advance.

I am waiting for this issue to be solved to move forward with the analysis.

Your help is greatly appreciated.

Best,

LyonLen commented 1 year ago

Hi~ Before you ran stereo_to_anndata, what else method did you run? nFeature_spatial and nCount_spatial is produced by data.tl.cal_qc.

limin321 commented 1 year ago

Hi~ Before you ran stereo_to_anndata, what else method did you run? nFeature_spatial and nCount_spatial is produced by data.tl.cal_qc.

Hi LyonLen,

I ran data.tl.cal_qc() before running st.io.stereo_to_anndata(data,flavor='seurat',output='seurat_out.h5ad'). And then converted to Seurat obj.

Then in R I ran the following 2 codes:

VlnPlot(wt, features = "nCount_Spatial", pt.size = 0.1) + NoLegend()
SpatialFeaturePlot(wt, features = "nCount_Spatial") + theme(legend.position = "right")

Here is the output:

Screen Shot 2023-06-02 at 11 23 33 AM

Why the heatmap doesn't show? Same as nFeature_Spatial plot.

I really appreciate your help.

Best, LC

LyonLen commented 1 year ago

Sorry that I am not so familiar with SpatialFeaturePlot in Seurat. Could you print the nCount_Spatial out? I would like to check whether the data is valid.

limin321 commented 1 year ago

Sorry that I am not so familiar with SpatialFeaturePlot in Seurat. Could you print the nCount_Spatial out? I would like to check whether the data is valid.

Can you see my screenshot? The nCount_Spatial VlnPlot is good. I don't think it is the issue because of Seurat. It has something to do with the RDS data converted from Stereopy.

Screen Shot 2023-06-04 at 8 45 24 PM

Any suggestion? Or could you try converting the stereo GEF data to RDS and test in Seurat? I think that will help you understand the issue more.

Looking forwarding to your updates.

Thank you so much. Best,

TheSallyGardens commented 1 year ago

@limin321 Hi,please try increasing the pt.size.factor.

SpatialFeaturePlot(brain, features = "nFeature_Spatial", pt.size.factor = 20) + theme(legend.position = "right")
LyonLen commented 1 year ago

@limin321 You can try as @TheSallyGardens informed. I am thinking this may not be a transformation bug. Thx~

limin321 commented 1 year ago

@limin321 Hi,please try increasing the pt.size.factor.

SpatialFeaturePlot(brain, features = "nFeature_Spatial", pt.size.factor = 20) + theme(legend.position = "right")

Thank you so much. This work. @LyonLen . Really appreciate for both of you.

May I have a follow-up question? During the conversion of stereoSeq data to Seurat RDS. It automatically ran the NormalizeData(). Why did this normalization step need to happen during conversion? Where is the normalized data stored in Seurat Obj? I saw there is only one assays called "Spatial". My understanding is that the raw counts is stored under the "Spatial" assay data slot, am I right?

Will this already-run NormalizeData() step affect later data integration in Seurat? Because Seurat suggests using SCTransform for normalization. Here is the Seurat Object when I first read-in the converted RDS in Seurat, where you can see the "Commands" list is not empty.

Screen Shot 2023-06-05 at 2 43 15 PM

Sorry, I just throw out so many questions at a time. I just try to understand how the data is structured in Seurat which will help me understand the data better.

Your help is greatly appreciated.

Best,

LyonLen commented 1 year ago

haha~ It's ok with that. Unfolding the spatial assay, you could see count, data, and scaled.data. count means the raw count of the expression matrix. data means the normalized count of the expression matrix. Hope this can help~

limin321 commented 1 year ago

haha~ It's ok with that. Unfolding the spatial assay, you could see count, data, and scaled.data. count means the raw count of the expression matrix. data means the normalized count of the expression matrix. Hope this can help~

Thank you so much! That is pretty helpful! Could you please also explain why you run NormalizeData during conversion? I am worried how this may affect later integration from different batches. Because Seurat suggests using SCTranform for integration.

LyonLen commented 1 year ago

Another assay will be produced by SCTransform. I think it will not affact anything. It is just pre-processing in order to fit our workflow in other application.

limin321 commented 1 year ago

Another assay will be produced by SCTransform. I think it will not affact anything. It is just pre-processing in order to fit our workflow in other application.

@LyonLen , Thanks. You are right. I tested by running NormalizeData again, it will overwrite the original data. So it shouldn't be any problem.

By the way, Could you please elaborate how the coordinates work in the Stereopy? After converting to Seurat obj, there are 2 columns called "x" and "y". Are these x and y the original coordinate from the Stereo-seq chip? I tried to use these coordinates to subset the seurat object as Seurat has provide xx_imagerow and xx_imagecol function. I couldn't figure out how the coordinates work.

For example, I run stereopy to export a subset of tissue.gef following this tutorial (https://stereopy.readthedocs.io/en/latest/Tutorials/High_Resolution_Export.html). The region I subset is in the middle of the tissue. Then, I converted it to .rds following this tutorial (https://stereopy.readthedocs.io/en/latest/Tutorials/Format_Conversion.html). After that, I read the rds file as an seurat object with the following codes and check the coordinates. The subset gef file is bin40. Interestingly, the coordinates for both row and column start with 1. That means the converted seurat obj didn't inherit the original coordinates in the original tissue gef file. Here are my codes and results:

wtROI <- readRDS("./final/wtROI_seurat.RDS")
wtROI_coord <- GetTissueCoordinates(wtROI, scale = "lowres", cols = c("imagerow", "imagecol"))
sort(unique(wtROI_coord$imagerow))
sort(unique(wtROI_coord$imagecol))

[1] 1 41 81 121 161 201 241 281 321 361 401 441 481 521 561 601 641 681 721 [20] 761 801 841 881 921 961 1001 1041 1081 1121 1161 1201 1241 1281 1321 1361 1401 1441 1481 [39] 1521 1561 1601 1641 1681 1721 1761 1801 1841 1881 1921 [1] 1 41 81 121 161 201 241 281 321 361 401 441 481 521 561 601 641 681 721 [20] 761 801 841 881 921 961 1001 1041 1081 1121 1161 1201 1241 1281 1321 1361

Next, I tested by runnint subset() as below

wt1s <- subset(wtROI, slice1_imagerow < 121 & slice1_imagecol < 121, invert = TRUE) # filter out spots < 121, NOT keep them?
SpatialFeaturePlot(wtROI, features = "nFeature_Spatial", pt.size.factor = 40)
SpatialFeaturePlot(wt1s, features = "nFeature_Spatial", pt.size.factor = 40)

After plot them, here is the output. The left is the wtROI I subset from original tissue.gef using stereopy tool, the right is subset based on coordinate in Seurat. You will see the top left spots are filtered out, which I thought it should keep those spots.

Screen Shot 2023-06-06 at 10 48 58 PM Screen Shot 2023-06-07 at 9 28 38 AM

opy/assets/42760669/fc7f2aed-7e6e-40e0-b53a-8f2f103d021b">

Then I did a second test by only subset 3 consecutive rows based on the output of this code sort(unique(wtROI_coord$imagerow)), and plot out the results, the 3 rows selected are so spread out.

wtROI1 <- subset(wtROI, slice1_imagerow > 40 & slice1_imagerow < 161, invert = FALSE) 
wtROI1df4 <- GetTissueCoordinates(wtROI1, scale = "lowres", cols = c("imagerow", "imagecol"))
SpatialFeaturePlot(wtROI1, features = "nFeature_Spatial", pt.size.factor = 60)
Screen Shot 2023-06-07 at 9 44 27 AM

Not sure why in the first test, spots < 121 rows were filtered out, while in the second test, spots < 161 rows were kept. Any explanation?

Could you please explain how the coordinates were assigned during conversion to Seurat obj? How is this different from the original coordinates in the tissue.gef?

Sorry, if I have too many questions. Try to understand the data more. You are really appreciated.

Best, LC

zhongll19 commented 1 year ago

我遇见过同样的问题,通过修改 annh5ad2rds的代码解决的,将 "Fiducial Diameter Fullres=1" 改为 "Fiducial Diameter Fullres=50",或者其他数值就好了。转化过程中,这个值会导致seurat对象中管理point 大小的值非常小,以至于输出的图像是黑点。