---
title: "Analysis of skin with Seurat- 2"
author: "B.Cheaib"
date: 'Compiled: `r Sys.Date()`'
output:
  html_document: 
    toc: yes
    fig_caption: yes
  df_print: default
  pdf_document: default
  word_document: default
editor_options: 
  chunk_output_type: inline
---

```{r setup, include=FALSE}
all_times <- list()  # store the time for each chunk
knitr::knit_hooks$set(time_it = local({
  now <- NULL
  function(before, options) {
    if (before) {
      now <<- Sys.time()
    } else {
      res <- difftime(Sys.time(), now, units = "secs")
      all_times[[options$label]] <<- res
    }
  }
}))
knitr::opts_chunk$set(
  tidy = TRUE,
  tidy.opts = list(width.cutoff = 95),
  message = FALSE,
  warning = FALSE,
  fig.width = 10,
  time_it = TRUE
)
```

# Overview

This tutorial demonstrates how to use Seurat (\>=3.2) to analyze spatially-resolved RNA-seq data. While the analytical pipelines are similar to the Seurat workflow for [single-cell RNA-seq analysis](pbmc3k_tutorial.html), we introduce updated interaction and visualization tools, with a particular emphasis on the integration of spatial and molecular information. This tutorial will cover the following tasks, which we believe will be common for many spatial analyses:

-   Normalization
-   Dimensional reduction and clustering
-   Detecting spatially-variable features
-   Interactive visualization
-   Integration with single-cell RNA-seq data
-   Working with multiple slices

For our first vignette, we analyze a dataset generated with the [Visium technology](https://www.10xgenomics.com/spatial-transcriptomics/) from 10x Genomics. We will be extending Seurat to work with additional data types in the near-future, including [SLIDE-Seq](https://science.sciencemag.org/content/363/6434/1463), [STARmap](https://science.sciencemag.org/content/361/6400/eaat5691), and [MERFISH](https://science.sciencemag.org/content/362/6416/eaau5324).

First, we load Seurat and the other packages necessary for this vignette.

```{r install}
library(Seurat)
# library(SeuratData)
library(ggplot2)
library(patchwork)
library(dplyr)
library("devtools")

```

```{r libraries.for.rmd, echo = FALSE}
library("htmltools")
library("vembedr")
```

# 10x Visium

## Dataset

Here, we will be using a recently released dataset of sagital mouse skin slices generated using the Visium v1 chemistry. There are two serial anterior sections, and two (matched) serial posterior sections.

You can download the data [here](https://support.10xgenomics.com/spatial-gene-expression/datasets), and load it into Seurat using the `Load10X_Spatial()` function. This reads in the output of the [spaceranger](https://support.10xgenomics.com/spatial-gene-expression/software/pipelines/latest/what-is-space-ranger) pipeline, and returns a Seurat object that contains both the spot-level expression data along with the associated image of the tissue slice. You can also use our [SeuratData package](https://github.com/satijalab/seurat-data) for easy data access, as demonstrated below. After installing the dataset, you can type `?stxskin` to learn more.

```{r data.import, include=TRUE}
Datadir<-"/mnt/Data/Bach_asterix_workspace/Spaceranger_work/SRanger_outpouts_brighter_images_not_Edited/104984-001-001-002/outs"

list.files(Datadir)

# Read10X_Image("Spaceranger_work/SRanger_outputs_Edited_LoupeBrowser_Cropped_figures/104984-001-001-002/outs/spatial",image.name="tissue_lowres_image.png")


skin= Load10X_Spatial(data.dir = Datadir)
setwd("/mnt/Data/Bach_asterix_workspace/Seurat_Analysis/V2_Seurat_Analysis")


```

<details>

<summary>**How is the spatial data stored within Seurat?**</summary>

The visium data from 10x consists of the following data types:

-   A spot by gene expression matrix
-   An image of the tissue slice (obtained from H&E staining during data acquisition)
-   Scaling factors that relate the original high resolution image to the lower resolution image used here for visualization.

In the Seurat object, the spot by gene expression matrix is similar to a typical "RNA" `Assay` but contains spot level, not single-cell level data. The image itself is stored in a new `images` slot in the Seurat object. The `images` slot also stores the information necessary to associate spots with their physical position on the tissue image.

</details>

\## Data preprocessing

The initial preprocessing steps that we perform on the spot by gene expression data are similar to a typical scRNA-seq experiment. We first need to normalize the data in order to account for variance in sequencing depth across data points. We note that the variance in molecular counts / spot can be substantial for spatial datasets, particularly if there are differences in cell density across the tissue. We see substantial heterogeneity here, which requires effective normalization.

```{r qc, warning=FALSE, include=TRUE}

plot1 <- VlnPlot(skin, features = 'nCount_Spatial', pt.size = 0.1) + NoLegend()
plot2 <- SpatialFeaturePlot(skin, features = 'nCount_Spatial') + theme(legend.position = "right")
wrap_plots(plot1, plot2)
```

These plots demonstrate that the variance in molecular counts across spots is not just technical in nature, but also is dependent on the tissue anatomy. For example, regions of the tissue that are depleted for neurons (such as the cortical white matter), reproducibly exhibit lower molecular counts. As a result, standard approaches (such as the `LogNormalize()` function), which force each data point to have the same underlying 'size' after normalization, can be problematic.

As an alternative, we recommend using sctransform (Hafemeister and Satija, Genome Biology 2019), which which builds regularized negative binomial models of gene expression in order to account for technical artifacts while preserving biological variance. For more details on sctransform, please see the paper [here](https://doi.org/10.1186/s13059-019-1874-1) and the Seurat vignette [here](sctransform_vignette.html). sctransform normalizes the data, detects high-variance features, and stores the data in the `SCT` assay.

```{r preprocess}
# This trial could provide Seurat object for Spatial not RNA

# GetAssayData(skin, slot="counts")
# counts<-GetAssayData(skin, slot="counts")
# nonzero <- counts > 0
# keep_genes<- Matrix::rowSums(nonzero) >=5
# skin_filtred <-counts[keep_genes,]
# skin_fil_new <-CreateAssayObject(skin_filtred, assays="Spatial")
# skin_fil_new@assay.orig<-"Spatial"
# skinfake <- SCTransform(skin_fil_new, assay = "Spatial", verbose = FALSE)

# Create a personalized Seurat Object
# devtools::install_github("jbergenstrahle/STUtility")
# library(STutility)
# library(hdf5r)
# infoTable<- data.frame (samples="/mnt/Data/Bach_asterix_workspace/Spaceranger_work/SRanger_outpouts_brighter_images_not_Edited/104984-001-001-001/outs/filtered_feature_bc_matrix.h5" ,
#                         spotfiles="/mnt/Data/Bach_asterix_workspace/Spaceranger_work/SRanger_outpouts_brighter_images_not_Edited/104984-001-001-001/outs/spatial/tissue_positions_list.csv",
#                         imgs="/mnt/Data/Bach_asterix_workspace/Spaceranger_work/SRanger_outpouts_brighter_images_not_Edited/104984-001-001-001/outs/spatial/tissue_hires_image.png",
#                         json="/mnt/Data/Bach_asterix_workspace/Spaceranger_work/SRanger_outpouts_brighter_images_not_Edited/104984-001-001-001/outs/spatial/scalefactors_json.json")
# 
# se <- InputFromTable(infotable = infoTable, 
#                      min.gene.count = 10, 
#                      min.gene.spots = 5,
#                      min.spot.count = 50,
#                      platform =  "Visium",
#                      assay="Spatial")
# 
# # Plot the number of unique genes and the number of UMIs per spot
# ST.FeaturePlot(se, features = c("nFeature_RNA", "nCount_RNA"), cols = c("lightgray", "mistyrose", "red", "darkred", "black"), ncol = 2, pt.size = 1.5, show.sb = TRUE)
# 
# # Plot selected genes
# ST.FeaturePlot(se, features = c("Cck", "Dcn"))
# 
# # Plot normalized values
# se <- SCTransform(se)


```

<details>

<summary>**How do results compare to log-normalization?**</summary>

To explore the differences in normalization methods, we examine how both the sctransform and log normalization results correlate with the number of UMIs. For this comparison, we first rerun sctransform to store values for all genes and run a log-normalization procedure via `NormalizeData()`.

```{r norm.test, include=TRUE}
# check if cells have no UMIs:
min(skin$nCount_Spatial)
skin <- subset(skin, nCount_Spatial>3000)

# rerun normalization to store sctransform residuals for all genes
skin <- SCTransform(skin, assay = "Spatial", return.only.var.genes = FALSE, verbose = FALSE)
# also run standard log normalization for comparison
skin <- NormalizeData(skin, verbose = FALSE, assay = "Spatial")
```

```{r norm.test2, include=TRUE}
# Computes the correlation of the log normalized data and sctransform residuals with the number of UMIs
skin <- GroupCorrelation(skin, group.assay = "Spatial", assay = "Spatial", slot = "data", do.plot = FALSE)
skin <- GroupCorrelation(skin, group.assay = "Spatial", assay = "SCT", slot = "scale.data", do.plot = FALSE)
```

```{r norm.test3, include=TRUE}
p1 <- GroupCorrelationPlot(skin, assay = "Spatial", cor = "nCount_Spatial_cor") + ggtitle("Log Normalization") + theme(plot.title = element_text(hjust = 0.5))
p2 <- GroupCorrelationPlot(skin, assay = "SCT", cor = "nCount_Spatial_cor") + ggtitle("SCTransform Normalization") + theme(plot.title = element_text(hjust = 0.5))
wrap_plots(p1 + p2)
```

For the boxplots above, we calculate the correlation of each feature (gene) with the number of UMIs (the `nCount_Spatial` variable here). We then place genes into groups based on their mean expression, and generate boxplots of these correlations. You can see that log-normalization fails to adequately normalize genes in the first three groups, suggesting that technical factors continue to influence normalized expression estimates for highly expressed genes. In contrast, sctransform normalization substantially mitigates this effect.

</details>

\## Gene expression visualization

In Seurat, we have functionality to explore and interact with the inherently visual nature of spatial data. The `SpatialFeaturePlot()` function in Seurat extends `FeaturePlot()`, and can overlay molecular data on top of tissue histology. For example, in this data set of the mouse skin, the gene Hpca is a strong hippocampus marker and Ttr is a marker of the choroid plexus.

```{r featureplot, include=TRUE}

## for fib1
f1<-SpatialFeaturePlot(skin, features = c("Ifi27l2a", "Mfap5", "Col14a1"))
## for fib2
f2<-SpatialFeaturePlot(skin, features = c("Bgn", "Dcn", "Col1a1", "Col1a2", "Col12a1","Col15a1", "Col16a1", "Cyp26b1"))
## for fib3
f3<-SpatialFeaturePlot(skin, features = c("Mgp", "Sparcl1", "Tmem176b",  "Tmem176a"))
## for fib4
f4<-SpatialFeaturePlot(skin, features = c("lamc3", "Crabp1","Col11a1", "Co23a1"))

pdf("Fibroblasts.pdf")
f1
f2
f3
f4
dev.off()

# for KC markers per cluster
# # for Basal-outer bulge-upper hair follic
# kc1<-SpatialFeaturePlot(skin, features = c("Krt35",	"Krt1", "Cst6", "Krt77"))
# # for inner hair follic
# kc2<-SpatialFeaturePlot(skin, features = c("Fabp5", "Krt10", 	"Krt14", 	"Aldoa",	"Krt79", "Serpinb2","Krt10"))
# # Basal2 may be inflammatory
# kc3<-SpatialFeaturePlot(skin, features = c("Krt25", 	"Krtdap", 	"Defb6"	, "Krt16", 	"Krt16", 	"Krt5"))

# for Hair Follic 1
kc1<-SpatialFeaturePlot(skin, features = c("Krt35", "Fabp5", "Krt25", "Hmgb2", "Ube2c","Krt71", "Lef1", "Krt27", "Ccnb1", "Dcn", "Kitl","Krt75", "Selenbp1", "Fabp4", "Krt73", "Tchh"))
# for Superbasal1
kc2<-SpatialFeaturePlot(skin, features = c("Krt1","Krt10","Krtdap","Krt77"))
# for Basal1
kc3<-SpatialFeaturePlot(skin, features = c("Krt14", "Krt5", "Bhlhe40", "Fth1", "Ly6a", "Mt2", "Avpi1"))
# for Hair follic 2
kc4<-SpatialFeaturePlot(skin, features = c("Postn", "Krt17", "Cd34", "Lgr5", "Krt6a"))
# for Superbasal2
kc5<-SpatialFeaturePlot(skin, features = c("Aldoa", "Ly6a", "Krt77", "Mt4"))
# for Hair follic 3
kc6<-SpatialFeaturePlot(skin, features = c("Cst6","Krt79","Defb6","Krt17", "Rbp1","Lrig1", "Cryab"))
# for Basal 2
kc7<-SpatialFeaturePlot(skin, features = c("Serpinb2", "Krt16", "Krt14"))
# for Superbasal 3
kc8<-SpatialFeaturePlot(skin, features = c("Krt5", "Ly6a", "Krt14", "Fth1"))
# for Basal3
kc9<-SpatialFeaturePlot(skin, features = c("Krt77", "Krt10", "Krtdap"))


pdf("Keratynocytes.pdf")
kc1
kc2
kc3
kc4
kc5
kc6
kc7
kc8
kc9
dev.off()

# # Tcells : TH1	vγ4 T cells 1	vγ4 T cells 2

# tc3<-SpatialFeaturePlot(skin, features = c("Ms4a4b", "Tcrg-C1", "Ctla2a"))
# TH1
tc1<-SpatialFeaturePlot(skin, features = c("Ms4a4b", "Ccr2", "Cd28", "Icos", "Cd4", "Ifng"))
# vγ4-T-cells-1
tc2<-SpatialFeaturePlot(skin, features = c("Tcrg-C1", "Trdc", "Trdv4", "Cd3e", "Cd3g", "Cd3d"))
# vγ4-T-cells-2
tc3<-SpatialFeaturePlot(skin, features = c("Cd3g", "Trdc", "Cd3e", "Tcrg-C1", "Cd3d", "Trdv4"))


pdf("TCells.pdf")
tc1
tc2
tc3
dev.off()


# CDS cells
## for Mac1
Mac1<-SpatialFeaturePlot(skin, features = c("S100a4", "H2-DMb1", "Pla2g7", "Cd300c2", "Cd14", "H2-Ab1", "H2-DMa"))
## for Mac2
Mac2<-SpatialFeaturePlot(skin, features = c("Cd209f", "Cd209g", "Lyve1", "Pla2g2d", "Clec10a","Adgre1"))
## for Mac3
Mac3<-SpatialFeaturePlot(skin, features = c("Hdc", "Gata2", "Ms4a2"))
## For DCs
dcs<-SpatialFeaturePlot(skin, features = c("Clec9a", "Xcr1", "Btla", "Itgae", "H2-Oa"))

pdf("MacDcs.pdf")
Mac1
Mac2
Mac3
dcs
dev.off()

```

```{r save.img, include=TRUE}
library(ggplot2)
plot <- SpatialFeaturePlot(skin, features = c("Fth1")) + 
  theme(legend.text = element_text(size = 0), legend.title = element_text(size = 20), legend.key.size = unit(1, "cm")) 

jpeg("spatial_vignette_ttr_mouse_skin.jpeg", height = 700, width = 1200, quality=50)
print(plot)

```

The default parameters in Seurat emphasize the visualization of molecular data. However, you can also adjust the size of the spots (and their transparency) to improve the visualization of the histology image, by changing the following parameters:

-   `pt.size.factor`- This will scale the size of the spots. Default is 1.6
-   `alpha` - minimum and maximum transparency. Default is c(1, 1).
-   Try setting to `alpha` c(0.1, 1), to downweight the transparency of points with lower expression

```{r fpe1 , include=TRUE}
p1 <- SpatialFeaturePlot(skin, features = "Fth1", pt.size.factor = 1)
p2 <- SpatialFeaturePlot(skin, features = "Fth1", alpha = c(0.1, 1))
p1 + p2
```

## Dimensionality reduction, clustering, and visualization

We can then proceed to run dimensionality reduction and clustering on the RNA expression data, using the same workflow as we use for scRNA-seq analysis.

```{r dim.cluster, include=TRUE}
skin <- RunPCA(skin, assay = "SCT", verbose = FALSE)
skin <- FindNeighbors(skin, reduction = "pca", dims = 1:30)
skin <- FindClusters(skin, verbose = FALSE)
skin <- RunUMAP(skin, reduction = "pca", dims = 1:30)
```

We can then visualize the results of the clustering either in UMAP space (with `DimPlot()`) or overlaid on the image with `SpatialDimPlot()`.

```{r dim.plots,fig.height=5, include=TRUE}
p1 <- DimPlot(skin, reduction = "umap", label = TRUE) 
p2 <- SpatialDimPlot(skin, label = TRUE, label.size = 3) 
p1 + p2
```

As there are many colors, it can be challenging to visualize which voxel belongs to which cluster. We have a few strategies to help with this. Setting the `label` parameter places a colored box at the median of each cluster (see the plot above).

You can also use the `cells.highlight` parameter to demarcate particular cells of interest on a `SpatialDimPlot()`. This can be very useful for distinguishing the spatial localization of individual clusters, as we show below:

```{r facetdim, include=TRUE}
SpatialDimPlot(skin, cells.highlight = CellsByIdentities(object = skin,idents = c(0,1,2,3,4,5,6,7)), facet.highlight = TRUE, ncol = 3)
```

## Interactive plotting

We have also built in a number of interactive plotting capabilities. Both `SpatialDimPlot()` and `SpatialFeaturePlot()` now have an `interactive` parameter, that when set to `TRUE`, will open up the Rstudio viewer pane with an interactive Shiny plot. The example below demonstrates an interactive `SpatialDimPlot()` in which you can hover over spots and view the cell name and current identity class (analogous to the previous `do.hover` behavior).

{r ispatialdimplot, include=TRUE}

SpatialDimPlot(skin, interactive = TRUE)

{r, echo = FALSE, include=FALSE} embed_url("<https://youtu.be/E1aZjmG1neQ>")

For `SpatialFeaturePlot()`, setting interactive to `TRUE` brings up an interactive pane in which you can adjust the transparency of the spots, the point size, as well as the `Assay` and feature being plotted. After exploring the data, selecting the done button will return the last active plot as a ggplot object.

{r ispatialfeatureplot, include=TRUE} SpatialFeaturePlot(skin, features = "Ttr", interactive = TRUE)

```{r, echo = FALSE, eval=FALSE, include=FALSE}
embed_url("https://youtu.be/ILmb8XNlgEM")
```

The `LinkedDimPlot()` function links the UMAP representation to the tissue image representation and allows for interactive selection. For example, you can select a region in the UMAP plot and the corresponding spots in the image representation will be highlighted.

{r linkedplot, include=TRUE} LinkedDimPlot(skin)

```{r, echo = FALSE, eval=FALSE, include=FALSE}
embed_url("https://youtu.be/10PZqjcSKrg")
```

## Identification of Spatially Variable Features

Seurat offers two workflows to identify molecular features that correlate with spatial location within a tissue. The first is to perform differential expression based on pre-annotated anatomical regions within the tissue, which may be determined either from unsupervised clustering or prior knowledge. This strategy works will in this case, as the clusters above exhibit clear spatial restriction.

```{r de, fig.height = 4, include=TRUE}
# de_markers <- FindMarkers(skin, ident.1=5, ident.2=4)
# SpatialFeaturePlot(object = skin, features = rownames(de_markers)[1:3], alpha = c(0.1, 1), ncol = 3)

Combined.filt.markers <- FindAllMarkers(skin, only.pos = TRUE, min.pct = 0.05, logfc.threshold = 0.25)
top8<-Combined.filt.markers %>% group_by(cluster) %>% top_n(n = 8, wt = avg_log2FC)
## SpatialFeaturePlot(object = skin, features = rownames(Combined.filt.markers)[1:8], alpha = c(0.1, 1), ncol = 3)

p1<-DoHeatmap(skin, features= top8$gene, slot="data", label = TRUE) + theme(axis.text.y = element_text(size = 10)) + scale_fill_gradientn(colors = c("black", "red")) 

p2 <- SpatialDimPlot(skin, label = TRUE, label.size = 3) 

wrap_plots(p1+p2)

pdf("heatmap-spatial.pdf")
p1
p2
dev.off()

write.csv(Combined.filt.markers, "allmarkers_minpct0.25_Logfc0.25.csv", row.names = TRUE)

 
```

An alternative approach, implemented in `FindSpatiallyVariables()`, is to search for features exhibiting spatial patterning in the absence of pre-annotation. The default method (`method = 'markvariogram`), is inspired by the [Trendsceek](https://www.nature.com/articles/nmeth.4634), which models spatial transcriptomics data as a mark point process and computes a 'variogram', which identifies genes whose expression level is dependent on their spatial location. More specifically, this process calculates gamma(r) values measuring the dependence between two spots a certain "r" distance apart. By default, we use an r-value of '5' in these analyses, and only compute these values for variable genes (where variation is calculated independently of spatial location) to save time.

We note that there are multiple methods in the literature to accomplish this task, including [SpatialDE](https://www.nature.com/articles/nmeth.4636), and [Splotch](https://www.biorxiv.org/content/10.1101/757096v1.article-metrics). We encourage interested users to explore these methods, and hope to add support for them in the near future.

```{r spatial.vf, include=TRUE}
skin <- FindSpatiallyVariableFeatures(skin, assay = 'SCT', features = VariableFeatures(skin)[1:1000], selection.method = 'markvariogram')
```

Now we visualize the expression of the top 6 features identified by this measure.

```{r spatial.vf.plot, fig.height=8, include=TRUE}
top.features <- head(SpatiallyVariableFeatures(skin, selection.method = 'markvariogram'),6)
SpatialFeaturePlot(skin, features = top.features, ncol = 3, alpha = c(0.1, 1))
```

```{r violon.plot.top, fig.height=8, include=TRUE}
VlnPlot(skin, features=top.features, ncol = 3,  pt.size = 0.2)
```
```{r Feature.plot.top, fig.height=8, include=TRUE}
FeaturePlot(skin, features=top.features, ncol = 3,  pt.size = 0.2)
```

```{r featureplotmarkers, include=TRUE}

## for fib1
f1<-FeaturePlot(skin, features = c("Ifi27l2a", "Mfap5", "Col14a1"))
## for fib2
f2<-FeaturePlot(skin, features = c("Bgn", "Dcn", "Col1a1", "Col1a2", "Col12a1","Col15a1", "Col16a1", "Cyp26b1"))
## for fib3
f3<-FeaturePlot(skin, features = c("Mgp", "Sparcl1", "Tmem176b",  "Tmem176a"))
## for fib4
f4<-FeaturePlot(skin, features = c("lamc3", "Crabp1","Col11a1", "Co23a1"))

pdf("UMAP-Fibroblasts.pdf")
f1
f2
f3
f4
dev.off()

# for KC markers per cluster
# # for Basal-outer bulge-upper hair follic
# kc1<-SpatialFeaturePlot(skin, features = c("Krt35",	"Krt1", "Cst6", "Krt77"))
# # for inner hair follic
# kc2<-SpatialFeaturePlot(skin, features = c("Fabp5", "Krt10", 	"Krt14", 	"Aldoa",	"Krt79", "Serpinb2","Krt10"))
# # Basal2 may be inflammatory
# kc3<-SpatialFeaturePlot(skin, features = c("Krt25", 	"Krtdap", 	"Defb6"	, "Krt16", 	"Krt16", 	"Krt5"))

# for Hair Follic 1
kc1<-FeaturePlot(skin, features = c("Krt35", "Fabp5", "Krt25", "Hmgb2", "Ube2c","Krt71", "Lef1", "Krt27", "Ccnb1", "Dcn", "Kitl","Krt75", "Selenbp1", "Fabp4", "Krt73", "Tchh"))
# for Superbasal1
kc2<-FeaturePlot(skin, features = c("Krt1","Krt10","Krtdap","Krt77"))
# for Basal1
kc3<-FeaturePlot(skin, features = c("Krt14", "Krt5", "Bhlhe40", "Fth1", "Ly6a", "Mt2", "Avpi1"))
# for Hair follic 2
kc4<-FeaturePlot(skin, features = c("Postn", "Krt17", "Cd34", "Lgr5", "Krt6a"))
# for Superbasal2
kc5<-FeaturePlot(skin, features = c("Aldoa", "Ly6a", "Krt77", "Mt4"))
# for Hair follic 3
kc6<-FeaturePlot(skin, features = c("Cst6","Krt79","Defb6","Krt17", "Rbp1","Lrig1", "Cryab"))
# for Basal 2
kc7<-FeaturePlot(skin, features = c("Serpinb2", "Krt16", "Krt14"))
# for Superbasal 3
kc8<-FeaturePlot(skin, features = c("Krt5", "Ly6a", "Krt14", "Fth1"))
# for Basal3
kc9<-FeaturePlot(skin, features = c("Krt77", "Krt10", "Krtdap"))


pdf("UMAP-Keratynocytes.pdf")
kc1
kc2
kc3
kc4
kc5
kc6
kc7
kc8
kc9
dev.off()

# # Tcells : TH1	vγ4 T cells 1	vγ4 T cells 2

# tc3<-SpatialFeaturePlot(skin, features = c("Ms4a4b", "Tcrg-C1", "Ctla2a"))
# TH1
tc1<-FeaturePlot(skin, features = c("Ms4a4b", "Ccr2", "Cd28", "Icos", "Cd4", "Ifng"))
# vγ4-T-cells-1
tc2<-FeaturePlot(skin, features = c("Tcrg-C1", "Trdc", "Trdv4", "Cd3e", "Cd3g", "Cd3d"))
# vγ4-T-cells-2
tc3<-FeaturePlot(skin, features = c("Cd3g", "Trdc", "Cd3e", "Tcrg-C1", "Cd3d", "Trdv4"))


pdf("UMAP-TCells.pdf")
tc1
tc2
tc3
dev.off()


# CDS cells
## for Mac1
Mac1<-FeaturePlot(skin, features = c("S100a4", "H2-DMb1", "Pla2g7", "Cd300c2", "Cd14", "H2-Ab1", "H2-DMa"))
## for Mac2
Mac2<-FeaturePlot(skin, features = c("Cd209f", "Cd209g", "Lyve1", "Pla2g2d", "Clec10a","Adgre1"))
## for Mac3
Mac3<-FeaturePlot(skin, features = c("Hdc", "Gata2", "Ms4a2"))
## For DCs
dcs<-FeaturePlot(skin, features = c("Clec9a", "Xcr1", "Btla", "Itgae", "H2-Oa"))

pdf("UMAP-MacDcs.pdf")
Mac1
Mac2
Mac3
dcs
dev.off()

```



```{r violonplotmarkers, include=TRUE}

## for fib1
f1<-VlnPlot(skin, features = c("Ifi27l2a", "Mfap5", "Col14a1"))
## for fib2
f2<-VlnPlot(skin, features = c("Bgn", "Dcn", "Col1a1", "Col1a2", "Col12a1","Col15a1", "Col16a1", "Cyp26b1"))
## for fib3
f3<-VlnPlot(skin, features = c("Mgp", "Sparcl1", "Tmem176b",  "Tmem176a"))
## for fib4
f4<-VlnPlot(skin, features = c("lamc3", "Crabp1","Col11a1", "Co23a1"))

pdf("Plots_Fibroblasts.pdf")
f1
f2
f3
f4
dev.off()

# for KC markers
# for Hair Follic 1
kc1<-VlnPlot(skin, features = c("Krt35", "Fabp5", "Krt25", "Hmgb2", "Ube2c","Krt71", "Lef1", "Krt27", "Ccnb1", "Dcn", "Kitl","Krt75", "Selenbp1", "Fabp4", "Krt73", "Tchh"))
# for Superbasal1
kc2<-VlnPlot(skin, features = c("Krt1","Krt10","Krtdap","Krt77"))
# for Basal1
kc3<-VlnPlot(skin, features = c("Krt14", "Krt5", "Bhlhe40", "Fth1", "Ly6a", "Mt2", "Avpi1"))
# for Hair follic 2
kc4<-VlnPlot(skin, features = c("Postn", "Krt17", "Cd34", "Lgr5", "Krt6a"))
# for Superbasal2
kc5<-VlnPlot(skin, features = c("Aldoa", "Ly6a", "Krt77", "Mt4"))
# for Hair follic 3
kc6<-VlnPlot(skin, features = c("Cst6","Krt79","Defb6","Krt17", "Rbp1","Lrig1", "Cryab"))
# for Basal 2
kc7<-VlnPlot(skin, features = c("Serpinb2", "Krt16", "Krt14"))
# for Superbasal 3
kc8<-VlnPlot(skin, features = c("Krt5", "Ly6a", "Krt14", "Fth1"))
# for Basal3
kc9<-VlnPlot(skin, features = c("Krt77", "Krt10", "Krtdap"))


pdf("Plots_Keratynocytes.pdf")
kc1
kc2
kc3
kc4
kc5
kc6
kc7
kc8
kc9
dev.off()

# # Tcells : TH1	vγ4 T cells 1	vγ4 T cells 2

tc1<-VlnPlot(skin, features = c("Ms4a4b", "Ccr2", "Cd28", "Icos", "Cd4", "Ifng"))
# vγ4-T-cells-1
tc2<-VlnPlot(skin, features = c("Tcrg-C1", "Trdc", "Trdv4", "Cd3e", "Cd3g", "Cd3d"))
# vγ4-T-cells-2
tc3<-VlnPlot(skin, features = c("Cd3g", "Trdc", "Cd3e", "Tcrg-C1", "Cd3d", "Trdv4"))


pdf("Plots_TCells.pdf")
tc1
tc2
tc3
dev.off()


# CDS cells
## for Mac1
Mac1<-VlnPlot(skin, features = c("S100a4", "H2-DMb1", "Pla2g7", "Cd300c2", "Cd14", "H2-Ab1", "H2-DMa"))
## for Mac2
Mac2<-VlnPlot(skin, features = c("Cd209f", "Cd209g", "Lyve1", "Pla2g2d", "Clec10a","Adgre1"))
## for Mac3
Mac3<-VlnPlot(skin, features = c("Hdc", "Gata2", "Ms4a2"))
## For DCs
dcs<-VlnPlot(skin, features = c("Clec9a", "Xcr1", "Btla", "Itgae", "H2-Oa"))

pdf("Plots_MacDcs.pdf")
Mac1
Mac2
Mac3
dcs
dev.off()


```

```{r integrationSingleCell , include=TRUE}
sc_reference <- readRDS("combined_mcs.rds")
getwd()


```

 note that setting ncells=3000 normalizes the full dataset but learns noise models on 3k
 cells this speeds up SCTransform dramatically with no loss in performance


```{r integrationSCTtransformUMAP , include=TRUE}
sc_reference <- SCTransform(sc_reference, ncells = 3000, verbose = FALSE) %>%
    RunPCA(verbose = FALSE) %>%
    RunUMAP(dims = 1:30)

```



```{r integrationSCNorm , include=TRUE}
# After subsetting, we renormalize 
skin<- SCTransform(skin, assay = "Spatial", verbose = FALSE) %>%
    RunPCA(verbose = FALSE)
# the annotation is stored in the 'subclass' column of object metadata
DimPlot(sc_reference, group.by = "CellType", label = TRUE)
```


```{r integrationSCTtransform , include=TRUE}
anchors <- FindTransferAnchors(reference = sc_reference, query = skin, normalization.method = "SCT")
predictions.assay <- TransferData(anchorset = anchors, refdata = sc_reference$CellType, prediction.assay = TRUE,
    weight.reduction = skin[["pca"]], dims = 1:30)
skin[["predictions"]] <- predictions.assay

```

```{r SpatialFeaturePlot1, include=TRUE}
DefaultAssay(skin) <- "predictions"
## select cells using the following command 
row.names(predictions.assay@data)
SpatialFeaturePlot(skin, features=c("pro KC", "KC 1", "KC 2",  "KC 3", "KC 4", "KC 5"), pt.size.factor = 1.2, ncol = 3, crop =TRUE,   combine = TRUE)
```

```{r SpatialFeaturePlot2 , include=TRUE}
SpatialFeaturePlot(skin, features=c("FB", "Melanocyte", "Langerhan's",  "Hair follicle KC"), pt.size.factor = 1.2, ncol = 2, crop =TRUE, combine = TRUE)
```

```{r SpatialFeaturePlot3 , include=TRUE}
SpatialFeaturePlot(skin, features=c("T/NK", "Endo", "Exclude 1",  "Exclude 2", "Exclude 3" ), pt.size.factor = 1.2, ncol = 3, crop =TRUE,   combine = TRUE)
```

```{r SpatialFeaturePlot4 , include=TRUE}
SpatialFeaturePlot(skin, features=c("DC/Mac", "max"), pt.size.factor = 1.2, ncol = 3, crop =TRUE,   combine = TRUE)
```


```{r integrationSCTtransformpredictions , include=TRUE}
skin <- FindSpatiallyVariableFeatures(skin, assay = "predictions", selection.method = "markvariogram",
    features = rownames(skin), r.metric = 5, slot = "data")
top.clusters <- head(SpatiallyVariableFeatures(skin), 4)
SpatialPlot(object = skin, features = top.clusters, ncol = 2)

```

saveRDS(skin, file="skin-B1.rds")



