# Load libraries
library(SingleCellExperiment)
## Loading required package: SummarizedExperiment
## Loading required package: MatrixGenerics
## Loading required package: matrixStats
## 
## Attaching package: 'MatrixGenerics'
## The following objects are masked from 'package:matrixStats':
## 
##     colAlls, colAnyNAs, colAnys, colAvgsPerRowSet, colCollapse,
##     colCounts, colCummaxs, colCummins, colCumprods, colCumsums,
##     colDiffs, colIQRDiffs, colIQRs, colLogSumExps, colMadDiffs,
##     colMads, colMaxs, colMeans2, colMedians, colMins, colOrderStats,
##     colProds, colQuantiles, colRanges, colRanks, colSdDiffs, colSds,
##     colSums2, colTabulates, colVarDiffs, colVars, colWeightedMads,
##     colWeightedMeans, colWeightedMedians, colWeightedSds,
##     colWeightedVars, rowAlls, rowAnyNAs, rowAnys, rowAvgsPerColSet,
##     rowCollapse, rowCounts, rowCummaxs, rowCummins, rowCumprods,
##     rowCumsums, rowDiffs, rowIQRDiffs, rowIQRs, rowLogSumExps,
##     rowMadDiffs, rowMads, rowMaxs, rowMeans2, rowMedians, rowMins,
##     rowOrderStats, rowProds, rowQuantiles, rowRanges, rowRanks,
##     rowSdDiffs, rowSds, rowSums2, rowTabulates, rowVarDiffs, rowVars,
##     rowWeightedMads, rowWeightedMeans, rowWeightedMedians,
##     rowWeightedSds, rowWeightedVars
## Loading required package: GenomicRanges
## Loading required package: stats4
## Loading required package: BiocGenerics
## 
## Attaching package: 'BiocGenerics'
## The following objects are masked from 'package:stats':
## 
##     IQR, mad, sd, var, xtabs
## The following objects are masked from 'package:base':
## 
##     anyDuplicated, append, as.data.frame, basename, cbind, colnames,
##     dirname, do.call, duplicated, eval, evalq, Filter, Find, get, grep,
##     grepl, intersect, is.unsorted, lapply, Map, mapply, match, mget,
##     order, paste, pmax, pmax.int, pmin, pmin.int, Position, rank,
##     rbind, Reduce, rownames, sapply, setdiff, sort, table, tapply,
##     union, unique, unsplit, which.max, which.min
## Loading required package: S4Vectors
## 
## Attaching package: 'S4Vectors'
## The following objects are masked from 'package:base':
## 
##     expand.grid, I, unname
## Loading required package: IRanges
## Loading required package: GenomeInfoDb
## Loading required package: Biobase
## Welcome to Bioconductor
## 
##     Vignettes contain introductory material; view with
##     'browseVignettes()'. To cite Bioconductor, see
##     'citation("Biobase")', and for packages 'citation("pkgname")'.
## 
## Attaching package: 'Biobase'
## The following object is masked from 'package:MatrixGenerics':
## 
##     rowMedians
## The following objects are masked from 'package:matrixStats':
## 
##     anyMissing, rowMedians
library(tidyverse)
## ── Attaching packages
## ───────────────────────────────────────
## tidyverse 1.3.2 ──
## ✔ ggplot2 3.3.6      ✔ purrr   0.3.4 
## ✔ tibble  3.1.8      ✔ dplyr   1.0.10
## ✔ tidyr   1.2.1      ✔ stringr 1.4.1 
## ✔ readr   2.1.3      ✔ forcats 0.5.2 
## ── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
## ✖ dplyr::collapse()   masks IRanges::collapse()
## ✖ dplyr::combine()    masks Biobase::combine(), BiocGenerics::combine()
## ✖ dplyr::count()      masks matrixStats::count()
## ✖ dplyr::desc()       masks IRanges::desc()
## ✖ tidyr::expand()     masks S4Vectors::expand()
## ✖ dplyr::filter()     masks stats::filter()
## ✖ dplyr::first()      masks S4Vectors::first()
## ✖ dplyr::lag()        masks stats::lag()
## ✖ ggplot2::Position() masks BiocGenerics::Position(), base::Position()
## ✖ purrr::reduce()     masks GenomicRanges::reduce(), IRanges::reduce()
## ✖ dplyr::rename()     masks S4Vectors::rename()
## ✖ dplyr::slice()      masks IRanges::slice()
library(Seurat)
## Attaching SeuratObject
## Attaching sp
## 
## Attaching package: 'Seurat'
## 
## The following object is masked from 'package:SummarizedExperiment':
## 
##     Assays
library(RColorBrewer)
library(scDblFinder)
library(magrittr)
## 
## Attaching package: 'magrittr'
## 
## The following object is masked from 'package:purrr':
## 
##     set_names
## 
## The following object is masked from 'package:tidyr':
## 
##     extract
## 
## The following object is masked from 'package:GenomicRanges':
## 
##     subtract
library(pheatmap)

# Define color palettes
WtOrRd_pal <- colorRampPalette(c("#FFFFFF", brewer.pal(9, "OrRd")))(100)
  1. Creating the Seurat Object
# Data Wrangling Processes
VTcounts2 <- read.delim("raw_data_10x.txt")
VTcounts2[1:10, 1:10]

# We notice that the gene names are too lengthy
VTcounts2[,1] <- sub("_.*", "", VTcounts2[,1])

# check for duplicated gene names before setting row names
gene_dups <- VTcounts2[,1] %>% 
  duplicated()

# check the dimensions of the count matrix
dim(VTcounts2)

# filter out the duplicated genes
VTcounts2 <- VTcounts2[!gene_dups,]
dim(VTcounts2)

rownames(VTcounts2) <- NULL

# set the row names
VTcounts2 %<>% 
  column_to_rownames("Gene")

head(VTcounts2)

# Create the Seurat object
vtSeurat2 <- CreateSeuratObject(counts = VTcounts2)

# Let's save the Seurat object as an .rds file to facilitate analysis.
saveRDS(vtSeurat2, "vtSeurat2.rds")

For convenience, we can start running from this chunk of code, since we have already saved the Seurat object as a file in our directory.

Since we have the metadata for this dataset, let’s add the cell annotations first, and filter the dataset to only include cells from the placenta. This will prevent the Seurat object becoming too large and time-costly to process.

vtSeurat2 <- readRDS("vtSeurat2.rds")
# View(vtSeurat2@meta.data)

# Check the dimensions of this object
dim(vtSeurat2)

# Load the cell annotation data
annotations <- read_delim(file = "meta_10x.txt")

# Check if the number of samples (cells) in vtSeurat2 is identical to the number of rows of this annotation file.

dim(annotations)[1] == dim(vtSeurat2)[2]

# make the annotations object compatible for data analysis.
anno_revised <- annotations %>% 
  separate(annotation, 
           into = c("Cluster", "Annotation"),
           sep = "\t") %>% 
  mutate(Location = final_cluster,
         Sample = location) %>% 
  select(-location, -final_cluster)

# Extract the metadata
meta_data.copy <- vtSeurat2@meta.data

# make new column that copies the rownames
meta_data.copy %<>% 
  mutate(Fetus = rownames(meta_data.copy))

# add the annotations!!!
meta_data.copy <- 
  inner_join(meta_data.copy, anno_revised, by = 'Fetus')

# Reassign the metadata
vtSeurat2@meta.data$Cluster <- meta_data.copy$Cluster
vtSeurat2@meta.data$Annotation <- meta_data.copy$Annotation
vtSeurat2@meta.data$Location <- meta_data.copy$Location
vtSeurat2@meta.data$Sample <- meta_data.copy$Sample

# Set the identity (classifier) of the cells based on the location

# Set the column "Location" as characters
vtSeurat2@meta.data$Location <- as.character(vtSeurat2@meta.data$Location)

Idents(object = vtSeurat2) <- "Location"
levels(Idents(vtSeurat2))


vtSeurat2@meta.data %>% 
  filter(Location == "Placenta")
# Subset the vtSeurat2 object to a smaller one.
vtSeurat2 <- subset(vtSeurat2, idents = "Placenta")

# The Seurat object became smaller.
dim(vtSeurat2)

# Let's save this object as another RDS file
saveRDS(vtSeurat2, "vtSeurat2_placenta.rds")
vtSeurat2 <- readRDS("vtSeurat2_placenta.rds")

# Add in the Mitochondrial PCT% information
vtSeurat2$percent.mt <- PercentageFeatureSet(vtSeurat2, pattern = "^MT-")

# nCount_RNA is the number of UMI counts in a cell
hist(vtSeurat2$nCount_RNA)

summary(vtSeurat2@meta.data$nCount_RNA)
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##    2894    9097   15033   19055   24238  125759
# nFeature_RNA is the number of different genes that had any reads
hist(vtSeurat2$nFeature_RNA)

summary(vtSeurat2@meta.data$nFeature_RNA)
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##    1086    2700    3584    3655    4500    8477
# percent.mt is the percent mitochondrial reads
summary(vtSeurat2$percent.mt)
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##       0       0       0       0       0       0
# We see that the mitochondrial genes were removed.

# Check for the percentage of ribosomal genes
vtSeurat2[["percent.rb"]] <- PercentageFeatureSet(vtSeurat2, pattern = "^RP[SL]")
summary(vtSeurat2$percent.rb)
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##   3.671  22.224  27.515  27.348  32.371  59.986
dim(vtSeurat2)
## [1] 31743 18547

We now have a nice vtSeurat object that contains only cells from the placenta. Luckily, we still have 18547 cells, which would be sufficient for downstream analysis. Another important part for QC is removing doublets.

# Find for doublets
vtSeurat2 <- as.SingleCellExperiment(vtSeurat2)
vtSeurat2 <- scDblFinder(vtSeurat2)
## Creating ~14838 artificial doublets...
## Dimensional reduction
## Evaluating kNN...
## Training model...
## iter=0, 1579 cells excluded from training.
## iter=1, 2139 cells excluded from training.
## iter=2, 2222 cells excluded from training.
## Threshold found:0.366
## 1727 (9.3%) doublets called
# Convert back to Seurat object.
vtSeurat2 <- as.Seurat(vtSeurat2)

# Notice that our identifiers have changed.
levels(Idents(object = vtSeurat2))
## [1] "SingleCellExperiment"
# We want to specify the identifier whether the cell is a singlet or doublet.
table(vtSeurat2@meta.data$scDblFinder.class)
## 
## singlet doublet 
##   16820    1727
Idents(object = vtSeurat2) <- "scDblFinder.class"
vtSeurat2@meta.data$scDblFinder.class <- as.character(vtSeurat2@meta.data$scDblFinder.class)

# Subset the Seurat object to only include singlets.
vtSeurat2 <- subset(vtSeurat2, idents = "singlet")
levels(as.factor(vtSeurat2@meta.data$scDblFinder.class))
## [1] "singlet"
# Check the dimensions of the dataset again.
dim(vtSeurat2)
## [1] 31743 16820

Now, note that we have 16779 cells that are derived from the placenta.

colnames(vtSeurat2@meta.data)
##  [1] "orig.ident"             "nCount_RNA"             "nFeature_RNA"          
##  [4] "Cluster"                "Annotation"             "Location"              
##  [7] "Sample"                 "percent.mt"             "percent.rb"            
## [10] "ident"                  "scDblFinder.class"      "scDblFinder.score"     
## [13] "scDblFinder.weighted"   "scDblFinder.cxds_score"
# convert column names to make the protocol coherent
colnames(vtSeurat2@meta.data)[2] <- "nCount_RNA"
colnames(vtSeurat2@meta.data)[3] <- "nFeature_RNA"

# Make a violin plot of the QC columns
plt <- VlnPlot(vtSeurat2, 
               features = c("nFeature_RNA", "nCount_RNA" ,"percent.rb"),
               ncol = 3,
               pt.size = 0.1) & 
  theme(plot.title = element_text(size=10))

plt

ggsave(filename = "QC.png", plot = plt, width = 7, height = 3.5)

Looking at the violin plot, we do see some extreme values that should be filtered out. However, let’s put the subsetting part for later and proceed to next steps and come back later.

# Make scatter plots to compare
plot1 <- FeatureScatter(vtSeurat2, feature1 = "nCount_RNA", feature2 = "nFeature_RNA")
plot1

# Looking at the scatterplot, we should subset to remove some outliers
# This is because further data analysis steps requires the assumption that the data has relatively moderate, homogeneous technical variants.
dim(vtSeurat2)
## [1] 31743 16820
# Let's perform some filtering!
vtSeurat2 <- subset(vtSeurat2, subset = percent.rb > 5)

Next, we need to normalize the data.

# Next, we normalize the data
vtSeurat2 <- NormalizeData(vtSeurat2)

# Explore the most variable features, let's use the default version first.
vtSeurat2 <- FindVariableFeatures(vtSeurat2)
# vtSeurat2 <- FindVariableFeatures(vtSeurat2, selection.method = "vst", nfeatures = 2000)

# Identify the 10 most highly variable genes
top10 <- head(VariableFeatures(vtSeurat2), 10) 
top10
##  [1] "PLA2G2A" "PSG3"    "AREG"    "MMP12"   "PRG2"    "CSH1"    "PSG5"   
##  [8] "PAPPA2"  "CHI3L1"  "PSG2"
# plot variable features with and without labels
plot2 <- VariableFeaturePlot(vtSeurat2)
plot2 <- LabelPoints(plot = plot2, points = top10, repel = TRUE)
## When using repel, set xnudge and ynudge to 0 for optimal results
plot2
## Warning: Transformation introduced infinite values in continuous x-axis
## Warning: Removed 6543 rows containing missing values (geom_point).

# Convert normalized gene expression to Z-score
all.genes <- rownames(vtSeurat2)
vtSeurat2 <- ScaleData(vtSeurat2, features = all.genes)
## Centering and scaling data matrix

Now, we are ready to run PCA.

# Run PCA analysis
vtSeurat2 <- RunPCA(vtSeurat2, features = VariableFeatures(object = vtSeurat2))
## PC_ 1 
## Positive:  KRT8, KRT7, GCSH, KRT19, GATA3, CDKN1C, LY6E, HMGB3, CLIC3, BZW2 
##     MORF4L2, FHL2, ISM2, EFNA1, S100P, PRSS8, GPX3, PTN, MFAP5, CCNE1 
##     GJA5, IGF2, CDK7, LIMA1, PGF, GATA2, CYP11A1, PTPRF, FLT1, OPN3 
## Negative:  TMSB4X, FTL, FCGRT, TYROBP, FCER1G, HLA-A, LAPTM5, AIF1, CD68, SH3BGRL 
##     ITGB2, CTSS, LY96, GMFG, CD36, TREM2, SPI1, CD53, HCST, HLA-B 
##     LST1, GPSM3, LGALS9, VIM, C1QC, UCP2, EMP3, SPP1, C1QA, RGS10 
## PC_ 2 
## Positive:  SAT1, HMGA1, CTSB, NPL, TYROBP, FTH1, FCER1G, S100A9, LAPTM5, AIF1 
##     ARHGDIB, CTSS, CREG1, ITGB2, CD53, TREM2, CD68, SPI1, HCST, KRT7 
##     LST1, LCP1, ENO1, CORO1A, PYCARD, C1QC, SPP1, LGALS9, C1QB, C1QA 
## Negative:  EGFL6, COL6A2, MEG3, DCN, COL6A1, DLK1, LUM, PCOLCE, HAPLN1, COL1A1 
##     COL3A1, PITX2, PLAC9, CXCL14, FRZB, CDH11, CTHRC1, COL1A2, WNT2, GPC3 
##     HGF, OLFML3, SPARC, VCAN, CCND2, SRPX, PLAGL1, ENPP1, MGP, HSD17B2 
## PC_ 3 
## Positive:  QSOX1, TGFB1, FSTL3, ITGA5, HAPLN3, C12orf75, ASCL2, DIO2, GLRX, FXYD5 
##     C1QTNF6, SLC16A3, MYCNUT, AIF1L, MT1X, PLAC8, GLIPR1, TPM1, ABCA7, COTL1 
##     HN1, HLA-C, HTRA4, MT1F, NOTUM, SLCO4A1, HLA-G, MT1E, AES, CA2 
## Negative:  SLC40A1, KRT23, BTG3, SLC26A2, SLC13A4, SLC52A1, GSTA3, INSL4, HSD17B1, HSPB8 
##     TMEM40, SOD3, CCR7, ARG2, CSTA, RP11-303E16.6, TBX3, SECISBP2L, CTNNB1, SEPP1 
##     SQSTM1, TACSTD2, TIMP3, ST3GAL4, GADD45G, MUC15, HMGA1, SLC1A2, PRKCZ, NPL 
## PC_ 4 
## Positive:  TOP2A, NUSAP1, BIRC5, CCNA2, CDK1, AURKB, TYMS, HMGB2, TUBA1B, STMN1 
##     PRC1, UBE2C, CENPF, MKI67, CDC20, HMMR, PTTG1, PBK, CDCA3, GTSE1 
##     HMGB1, KIAA0101, SMC4, CDKN3, SPC25, TPX2, TK1, MAD2L1, PLK1, H2AFZ 
## Negative:  S100P, PVRL4, GDF15, CLIC3, FLT1, CYP19A1, EPS8L1, TMEM40, PLIN2, INSL4 
##     SQSTM1, LYPD3, PRSS8, SAT1, MUC15, CYP11A1, INHA, PRKCZ, EBI3, CEBPB 
##     NUCB2, GADD45G, SERPINE1, GATA3, DHRS9, SGK1, CGA, HERPUD1, DNAJB9, LGALS16 
## PC_ 5 
## Positive:  LYVE1, F13A1, FOLR2, HPGDS, CD5L, EGFL7, LRRC39, IGFBP4, DAB2, IGF1 
##     MAF, SCN9A, MRC1, TMIGD3, C2, GAPLINC, FAM20A, VSIG4, SLC28A1, LILRB5 
##     C10orf128, CD14, METTL7B, GPR34, GAS6, ITM2B, SIGLEC1, FAM212A, ABCG2, FEZ1 
## Negative:  HLA-DRA, HLA-DPA1, HLA-DRB1, HLA-DPB1, LSP1, HLA-DQB1, HLA-DRB5, HLA-DQA1, CD52, CD74 
##     HLA-DMA, LILRB4, SDS, APOC1, FBP1, MSR1, LINC01272, CXCR4, LYZ, RGS1 
##     C1orf162, IL4I1, IFI30, HLA-DQA2, FPR3, SERPINA1, CD44, LPL, HAMP, CYP27A1
# Plot the PCA
DimPlot(vtSeurat2, reduction = "pca")

# Plot the loadings
VizDimLoadings(vtSeurat2, dims = 1:2, reduction = "pca") & 
  theme(axis.text=element_text(size=8), axis.title=element_text(size = 8,face="bold"))

# Check the elbow plot to see how much variation is covered for each PC
ElbowPlot(vtSeurat2)

# FeaturePlot of PC1 top features
FeaturePlot(vtSeurat2, features = c("KRT8", 
                                   "KRT7"))

# FeaturePlot of PC2 top features
FeaturePlot(vtSeurat2, features = c("ENPP1", 
                                   "MGP"))

Next, we move on to derive clusters.

# Now perform clustering
vtSeurat2 <- FindNeighbors(vtSeurat2, dims = 1:10)
## Computing nearest neighbor graph
## Computing SNN
vtSeurat2 <- FindClusters(vtSeurat2, resolution = 0.5)
## Modularity Optimizer version 1.3.0 by Ludo Waltman and Nees Jan van Eck
## 
## Number of nodes: 16818
## Number of edges: 535612
## 
## Running Louvain algorithm...
## Maximum modularity in 10 random starts: 0.9225
## Number of communities: 18
## Elapsed time: 2 seconds
# For visualization purposes, we need to generate UMAP reduced dimensionality representation
vtSeurat2 <- RunUMAP(vtSeurat2, dims = 1:10, verbose = F)
## Warning: The default method for RunUMAP has changed from calling Python UMAP via reticulate to the R-native UWOT using the cosine metric
## To use Python UMAP via reticulate, set umap.method to 'umap-learn' and metric to 'correlation'
## This message will be shown once per session
# Plot the UMAP with clustering
DimPlot(vtSeurat2, reduction = "umap")

# Get the TSNE embedding
vtSeurat2 <- RunTSNE(vtSeurat2, dims = 1:10)

# Plot the t-SNE with clustering
DimPlot(vtSeurat2, reduction = "tsne")

# Let's look at cluster sizes
table(vtSeurat2@meta.data$seurat_clusters)
## 
##    0    1    2    3    4    5    6    7    8    9   10   11   12   13   14   15 
## 2679 2088 1977 1795 1514 1215 1067  949  841  767  553  455  344  231  143  108 
##   16   17 
##   56   36
# Let's plot those clusters!
DimPlot(vtSeurat2,label.size = 4,repel = T,label = T)

# Check expression patterns for some genes in top10
FeaturePlot(vtSeurat2, features = c("FN1",
                                   "SERPINE2",
                                   "IGF2",
                                   "RGS5"))

Let’s calculate the cell cycle scores at this point.

vtSeurat2 <- SetIdent(vtSeurat2, value = "Annotation")

# Calculate the cell cycle scores
s.genes <- cc.genes.updated.2019$s.genes
g2m.genes <- cc.genes.updated.2019$g2m.genes

# Somehow this code does not work.
vtSeurat2 <- CellCycleScoring(vtSeurat2,
                            search = TRUE,
                            s.features = s.genes, 
                            g2m.features = g2m.genes)
## Warning: The following features are not present in the object: PIMREG, JPT1,
## attempting to find updated synonyms
## Warning: No updated symbols found
## Warning: The following features are still not present in the object: PIMREG,
## JPT1
table(vtSeurat2[[]]$Phase)
## 
##    G1   G2M     S 
## 10451  2767  3600

Around 1/5 of the cells are in the S phase.

# Check if ribosomal genes show dependency on cluster
FeaturePlot(vtSeurat2,features = "percent.rb",
            label.size = 4,
            repel = T,
            label = F) & 
  theme(plot.title = element_text(size=10))

# Check using the violin plot
VlnPlot(vtSeurat2,
        features = "percent.rb") & 
  theme(plot.title = element_text(size=10),
        legend.position = 'none')
## Warning: Groups with fewer than two data points have been dropped.
## Groups with fewer than two data points have been dropped.
## Groups with fewer than two data points have been dropped.
## Groups with fewer than two data points have been dropped.
## Groups with fewer than two data points have been dropped.
## Groups with fewer than two data points have been dropped.

# Check # of genes and # of transcripts per cluster
VlnPlot(vtSeurat2,features = "nCount_RNA") & 
  theme(plot.title = element_text(size=10),
        legend.position = 'none')
## Warning: Groups with fewer than two data points have been dropped.
## Groups with fewer than two data points have been dropped.
## Groups with fewer than two data points have been dropped.
## Groups with fewer than two data points have been dropped.
## Groups with fewer than two data points have been dropped.
## Groups with fewer than two data points have been dropped.

VlnPlot(vtSeurat2,features = "nFeature_RNA") & 
  theme(plot.title = element_text(size=10),
        legend.position = 'none')
## Warning: Groups with fewer than two data points have been dropped.
## Groups with fewer than two data points have been dropped.
## Groups with fewer than two data points have been dropped.
## Groups with fewer than two data points have been dropped.
## Groups with fewer than two data points have been dropped.
## Groups with fewer than two data points have been dropped.

# Cell cycle score does not seem to depend on the cell type much
FeaturePlot(vtSeurat2,
            features = c("S.Score","G2M.Score"),
            label.size = 4,
            repel = T,
            label = T) & 
  theme(plot.title = element_text(size = 10))

VlnPlot(vtSeurat2, features = c("S.Score","G2M.Score")) & 
  theme(plot.title = element_text(size = 10),
        legend.position = 'none')
## Warning: Groups with fewer than two data points have been dropped.
## Groups with fewer than two data points have been dropped.
## Groups with fewer than two data points have been dropped.
## Groups with fewer than two data points have been dropped.
## Groups with fewer than two data points have been dropped.
## Groups with fewer than two data points have been dropped.
## Groups with fewer than two data points have been dropped.
## Groups with fewer than two data points have been dropped.
## Groups with fewer than two data points have been dropped.
## Groups with fewer than two data points have been dropped.
## Groups with fewer than two data points have been dropped.
## Groups with fewer than two data points have been dropped.

In this case, it looks like some cells do show some cell cycle scores, which should be examined more carefully.

# Plot in UMAP
DimPlot(vtSeurat2, 
        group.by = "Phase")

# Plot with clusters
compTable.int <- table(vtSeurat2$Annotation, vtSeurat2$Phase)
compTable.int <- (compTable.int / rowSums(compTable.int)) * 100
pheatmap(compTable.int, color = WtOrRd_pal, breaks = 1:100)

# FeaturePlot for showing mod score
FeaturePlot(vtSeurat2, 
            features = c("S.Score", "G2M.Score"))

# This part is replaced when using the SCTransform
DefaultAssay(vtSeurat2) <- "RNA"

# Perform log-normalization
vtSeurat2 <- NormalizeData(vtSeurat2)

# Find variable features
vtSeurat2 <- FindVariableFeatures(vtSeurat2, 
                                 selection.method = "vst", 
                                 nfeatures = 2000)
all.genes <- rownames(vtSeurat2)
vtSeurat2 <- ScaleData(vtSeurat2, features = all.genes)
## Centering and scaling data matrix
# Designate the markers
all.markers <- FindAllMarkers(vtSeurat2,
                              only.pos = T, 
                              min.pct = 0.5, 
                              logfc.threshold = 0.5)
## Calculating cluster dS1
## Calculating cluster VCT
## Calculating cluster dNK2
## Calculating cluster dNK1
## Calculating cluster dM2
## Calculating cluster dM1
## Calculating cluster Tcells
## Calculating cluster fFB1
## Calculating cluster EVT
## Calculating cluster NK CD16+
## Calculating cluster HB
## Calculating cluster SCT
## Calculating cluster dM3
## Calculating cluster MO
## Calculating cluster dNK p
## Calculating cluster Plasma
## Calculating cluster DC1
## Calculating cluster DC2
## Calculating cluster Granulocytes
## Calculating cluster ILC3
## Calculating cluster fFB2
## Calculating cluster Endo (m)
## Calculating cluster Endo L
## Calculating cluster Endo (f)
## Calculating cluster Epi1
## Calculating cluster Epi2
dim(all.markers)
## [1] 6896    7
table(all.markers$cluster)
## 
##          dS1          VCT         dNK2         dNK1          dM2          dM1 
##            0          450            0          337          380          302 
##       Tcells         fFB1          EVT     NK CD16+           HB          SCT 
##          235          349          364            0          469          332 
##          dM3           MO        dNK p       Plasma          DC1          DC2 
##          427          311            0          126            0          345 
## Granulocytes         ILC3         fFB2     Endo (m)       Endo L     Endo (f) 
##          290            0          405          434          419          420 
##         Epi1         Epi2 
##            0          501
top3_markers <- as.data.frame(all.markers %>%
                              group_by(cluster) %>% 
                              top_n(n = 3, wt = avg_log2FC))
head(top3_markers)
##           p_val avg_log2FC pct.1 pct.2     p_val_adj cluster  gene
## 1  0.000000e+00   3.734920 0.999 0.894  0.000000e+00     VCT PAGE4
## 2  0.000000e+00   2.692421 0.985 0.480  0.000000e+00     VCT SMAGP
## 3  0.000000e+00   2.625993 0.996 0.542  0.000000e+00     VCT PEG10
## 4 4.881180e-184   6.012708 1.000 0.006 1.549433e-179    dNK1  CTSW
## 5  1.959920e-61   6.039312 1.000 0.018  6.221374e-57    dNK1  NKG7
## 6  4.210265e-36  10.512587 1.000 0.033  1.336464e-31    dNK1  GNLY

Clustering based on Cell Type annotation

DimPlot(vtSeurat2, label = T , repel = T, label.size = 3) + NoLegend()
## Warning: ggrepel: 9 unlabeled data points (too many overlaps). Consider
## increasing max.overlaps

DimPlot(vtSeurat2, reduction = "pca")

DimPlot(vtSeurat2, reduction = "tsne")

DimPlot(vtSeurat2, reduction = "umap")

# ggsave(filename = "Check.png", plot = p1, width = 30, height = 10)
saveRDS(vtSeurat2, "vtSeurat2_annotated.rds")

Still Work in Progress

# Pseudotime
# Convert to monocle3 cell_data_set
vtSeurat2.cds <- as.cell_data_set(vtSeurat2)
class(vtSeurat2.cds)

# Re-cluster (required for trajectories in monocle3)
vtSeurat2.cds <- cluster_cells(vtSeurat2.cds)

# Learn the trajectory
vtSeurat2.cds <- learn_graph(vtSeurat2.cds)

# Plot the trajectory
plot_cells(vtSeurat2.cds, color_cells_by = "seurat_clusters", show_trajectory_graph = TRUE)

# Find the pseudotime order
DimPlot(vtSeurat2, group.by = "Phase")
vtSeurat2.cds <- order_cells(vtSeurat2.cds)

# Plot pseudotime
plot_cells(vtSeurat2.cds, color_cells_by = "pseudotime")

# Add pseudotime score back into the dataset
vtSeurat2$pseudotime_score <- vtSeurat2.cds@principal_graph_aux@listData$UMAP$pseudotime
vtSeurat2$pseudotime_score[vtSeurat2$pseudotime_score == Inf] <- NA

# FeaturePlot with Pseudotime
FeaturePlot(vtSeurat2, features = c("pseudotime_score")) 
DimPlot(vtSeurat2, group.by = "Annotation")
FeaturePlot(vtSeurat2, features = c("pseudotime_score"))