Skip to content

Commit

Permalink
Merge pull request #78 from JunxiangXu/main
Browse files Browse the repository at this point in the history
some typo and format error for spat proteomics and co registration
  • Loading branch information
jiajic authored Aug 5, 2024
2 parents 4565fd7 + dd0de23 commit 4afa19d
Show file tree
Hide file tree
Showing 2 changed files with 24 additions and 23 deletions.
9 changes: 5 additions & 4 deletions 02_session4.Rmd
Original file line number Diff line number Diff line change
Expand Up @@ -70,16 +70,16 @@ list.files(paste0(download_dir,'/Lunaphore'))
```
We provide a way to extract meta data information directly from ome.tiffs. Please note that different platforms may store the meta data such as channel information in a different format, we will probably need to change the node names of the ome-XML.
```{r eval=FALSE}
img_path <- paste0(download_dir,"Lunaphore/Lunaphore_example.ome.tiff")
img_path <- paste0(download_dir,"/Lunaphore/Lunaphore_example.ome.tiff")
img_meta <- ometif_metadata(img_path, node = "Channel", output = "data.frame")
img_meta
```

However, sometimes a cropping could result in a loss of ome-XML information from the ome.tiff file. That way, we can use a different strategy to parse the xml information seperately and get channel information from it.
However, sometimes a simple ometiff file manipulation like cropping could result in a loss of ome-XML information from the ome.tiff file. That way, we can use a different strategy to parse the xml information seperately and get channel information from it.
```{r eval=FALSE}
##Get channel information
Luna = paste0(download_dir,"Lunaphore/Lunaphore_example.ome.tiff")
xmldata = xml2::read_xml(paste0(download_dir,"Lunaphore/Lunaphore_sample_metadata.xml"))
Luna = paste0(download_dir,"/Lunaphore/Lunaphore_example.ome.tiff")
xmldata = xml2::read_xml(paste0(download_dir,"/Lunaphore/Lunaphore_sample_metadata.xml"))
node <- xml2::xml_find_all(xmldata, "//d1:Channel", ns = xml2::xml_ns(xmldata))
channel_df <- as.data.frame(Reduce("rbind", xml2::xml_attrs(node)))
channel_df
Expand Down Expand Up @@ -159,6 +159,7 @@ knitr::include_graphics("img/02_session4/mesmer_overlap.png")
```

### Using Giotto wrapper of Cellpose to perform segmentation
Here, we create a mini example by cropping the image to a smaller area. Note that crop() is probably easier to use to directly crop image, unless cropping the image when the image is inside of a giotto object.
```{r eval = F}
gimg_cropped <- cropGiottoLargeImage(giottoLargeImage = gimg_DAPI, crop_extent = terra::ext(zoom))
writeGiottoLargeImage(gimg_cropped, filename = paste0(download_dir,'/Lunaphore/DAPI_forcellpose.tiff'),overwrite = T)
Expand Down
38 changes: 19 additions & 19 deletions 03_session2.Rmd
Original file line number Diff line number Diff line change
@@ -1,6 +1,7 @@
# Spatial multi-modal analysis

George Chen

Junxiang Xu

August 7th 2024
Expand Down Expand Up @@ -227,13 +228,13 @@ Currently _giotto_ image objects are not fully compatible with .ome.tif files. _
10X Genomics Released a comprehensive dataset on 2022. To capture spatial structure by complementing different spatial resolutions and modalities across different assays, they provided a dataset with Xenium in situ transcriptomics data, together with Visium on closely adjacent sections. Additional IF staining was also performed on the Xenium slides. For more information, please refer to the pre-release dataset [page](https://www.10xgenomics.com/products/xenium-in-situ/preview-dataset-human-breast) as well as the [publication](https://www.nature.com/articles/s41467-023-43458-x).

- Visium
-- H&E Histology
-- 55um spot level expression with transcriptome coverage
- H&E Histology
- 55um spot level expression with transcriptome coverage
- Xenium
-- H&E Histology
-- IF image staining DAPI, HER2 and CD20
-- in situ transcripts
-- cooresponding centroid locations
- H&E Histology
- IF image staining DAPI, HER2 and CD20
- in situ transcripts
- cooresponding centroid locations

The goal of creating this multi-modal dataset is to register all the modalities listed above to the same coordinate system as Xenium in situ transcripts as the coordinate represents a certain micron distance.
```{r eval=FALSE}
Expand All @@ -247,14 +248,14 @@ destfile <- file.path(download_dir,'Multimodal_registration.zip')
if (!dir.exists(download_dir)) { dir.create(download_dir, recursive = TRUE) }
download.file('https://zenodo.org/records/13208139/files/Multimodal_registration.zip?download=1', destfile = destfile)
unzip(paste0(download_dir,'/Multimodal_registration.zip'), exdir = download_dir)
Xenium_dir <- paste0(download_dir,'/Multimodal_registration/Xenium/')
Visium_dir <- paste0(download_dir,'/Multimodal_registration/Visium/')
Xenium_dir <- paste0(download_dir,'/Xenium/')
Visium_dir <- paste0(download_dir,'/Visium/')
```

### Target Coordinate system
Xenium transcripts, polygon information and corresponding centroids are output from the Xenium instrument and are in the same coordinate system from the raw output. We can start with checking the centroid information as a representation of the target coordinate system.
```{r eval=FALSE}
xen_cell_df <- read.csv(paste0(Xenium_dir,"cells.csv.gz"))
xen_cell_df <- read.csv(paste0(Xenium_dir,"/cells.csv.gz"))
xen_cell_pl <- ggplot2::ggplot() + ggplot2::geom_point(data = xen_cell_df, ggplot2::aes(x = x_centroid , y = y_centroid),size = 1e-150,,color = 'orange') + ggplot2::theme_classic()
xen_cell_pl
```
Expand Down Expand Up @@ -329,7 +330,7 @@ knitr::include_graphics("img/03_session2/Visium_resgister_check.png")
#### Create Pseudo Visium dataset for comparison
Giotto provides a way to create different shapes on certain locations, we can use that to create a pseudo-visium polygons to aggregate transcripts or image intensities. To do that, we will need the centroid locations, which can be get using getSpatialLocations(). and also the radius information to create circles. We know that Visium certer to center distance is 100um and spot diameter is 55um, thus we can estimate the radius from certer to center distance. And we can use a spatial network created by nearest neighbor = 2 to capture the distance.
```{r eval=FALSE}
V_final <- createSpatialNetwork(V_final, k = 1,method = 'kNN')
V_final <- createSpatialNetwork(V_final, k = 1)
spat_network <- getSpatialNetwork(V_final,output = 'networkDT')
spatPlot2D(V_final,
show_network = T,
Expand Down Expand Up @@ -358,7 +359,7 @@ knitr::include_graphics("img/03_session2/Visium_polygon.png")

Create Xenium object with pseudo Visium polygon. To save run time, example shown here only have MS4A1 and ERBB2 genes to create Giotto points
```{r eval = FALSE}
xen_transcripts <- data.table::fread(paste0(Xenium_dir,'Xen_2_genes.csv.gz'))
xen_transcripts <- data.table::fread(paste0(Xenium_dir,'/Xen_2_genes.csv.gz'))
gpoints <- createGiottoPoints(xen_transcripts)
Xen_obj <-createGiottoObjectSubcellular(gpoints = list('rna' = gpoints),
gpolygons = list('visium' = pseudo_visium_poly))
Expand Down Expand Up @@ -432,8 +433,8 @@ knitr::include_graphics("img/03_session2/Visium_spatfeatplot.png")
### Register post-Xenium H&E and IF image
For Xenium instrument output, Giotto provide a convenience function to load the output from the Xenium ranger output. Note that 10X created the affine image alignment file by applying rotation, scale at (0,0) of the top left corner and translation last. Thus, it will look different than the affine matrix created from landmarks above. In this example, we used a 0.05X compressed ometiff and the alignment file is also create by first rescale at 20X, then apply the affine matrix provided by 10X Genomics.
```{r eval = F}
HE_xen <- read10xAffineImage(file = paste0(Xenium_dir, "HE_ome_compressed.tiff"),
imagealignment_path = paste0(Xenium_dir,"Xenium_he_imagealignment.csv"),
HE_xen <- read10xAffineImage(file = paste0(Xenium_dir, "/HE_ome_compressed.tiff"),
imagealignment_path = paste0(Xenium_dir,"/Xenium_he_imagealignment.csv"),
micron = 0.2125)
plot(HE_xen)
```
Expand Down Expand Up @@ -479,7 +480,7 @@ Xen_obj <- addGiottoLargeImage(gobject = Xen_obj,

Get the cell polygons, as Xenium and IF are both subcellular resolution
```{r eval=FALSE}
cellpoly_dt <- data.table::fread(paste0(Xenium_dir,'cell_boundaries.csv.gz'))
cellpoly_dt <- data.table::fread(paste0(Xenium_dir,'/cell_boundaries.csv.gz'))
colnames(cellpoly_dt) <- c('poly_ID','x','y')
cellpoly <- createGiottoPolygonsFromDfr(cellpoly_dt)
Xen_obj <- addGiottoPolygons(Xen_obj,gpolygons = list('cell' = cellpoly))
Expand Down Expand Up @@ -595,7 +596,7 @@ Pin landmarks or use compounded affine transforms to register image usually prov

Here, we provide an example of two compressed images to show the automatic alignment pipeline.
```{r eval = F}
HE <- createGiottoLargeImage(paste0(Xenium_dir,'mini_HE.png'),negative_y = F)
HE <- createGiottoLargeImage(paste0(Xenium_dir,'/mini_HE.png'),negative_y = F)
plot(HE)
```

Expand All @@ -604,7 +605,7 @@ knitr::include_graphics("img/03_session2/mini_HE.png")
```

```{r eval = F}
IF <- createGiottoLargeImage(paste0(Xenium_dir,'mini_IF.tif'),negative_y = F,flip_horizontal = T)
IF <- createGiottoLargeImage(paste0(Xenium_dir,'/mini_IF.tif'),negative_y = F,flip_horizontal = T)
terra::plotRGB(IF@raster_object,r=1, g=2, b=3,, stretch="lin")
```

Expand All @@ -614,11 +615,11 @@ knitr::include_graphics("img/03_session2/mini_IF.png")

Now, we can use the automated transformation pipeline. Note that we will start with a path to the images, run the preprocessImageToMatrix() first to meet the requirement of estimateAutomatedImageRegistrationWithSIFT() function. The images will be preprocessed to gray scale. And for that purpose, we use the DAPI channel from the miniIF, and set invert = T for mini HE as HE image so that grayscle HE will have higher value for high intensity pixels. The function will output an estimation of the transform.
```{r eval = F}
estimation <- estimateAutomatedImageRegistrationWithSIFT(x = preprocessImageToMatrix(paste0(Xenium_dir,'mini_IF.tif'),
estimation <- estimateAutomatedImageRegistrationWithSIFT(x = preprocessImageToMatrix(paste0(Xenium_dir,'/mini_IF.tif'),
flip_horizontal = T,
use_single_channel = T,
single_channel_number = 3),
y = preprocessImageToMatrix(paste0(Xenium_dir,'mini_HE.png'),
y = preprocessImageToMatrix(paste0(Xenium_dir,'/mini_HE.png'),
invert = T),
plot_match = T,
max_ratio = 0.5,estimate_fun = 'Projective')
Expand All @@ -635,7 +636,6 @@ transformed <- affine(IF, mtx)
To_see_overlay <- transformed@funs$realize_magick(size = 2e6)
plot(HE)
plot(To_see_overlay@raster_object[[2]], add=TRUE, alpha=0.5)
SIFT_registration_overlay
```

```{r, echo=FALSE, out.width="80%", fig.align='center'}
Expand Down

0 comments on commit 4afa19d

Please sign in to comment.