notes from first run

[LPerlaza]

• explain a bit the experiment we are analysing
• change parameters that don't make sense, like cutoff for doublets (nFeature_RNA < 2500)
• add the clustering code to the main code

rm(list=ls());
library(langevitour)
library(crosstalk)
library(ggplot2)
library(GGally)
library(plotly)
library(DT)
library(htmltools)
library(palmerpenguins)
library(dplyr)
library(Seurat)
library(patchwork)
library(kableExtra)
library(clustree)
library(SeuratDisk)

pbmc.data <- Read10X(data.dir = "/mnt/ceph/mbp/servers/bioinformatics-platform/home/lper0012/tasks/training/filtered_gene_bc_matrices/hg19/")
pbmc <- CreateSeuratObject(counts = pbmc.data, project = "pbmc3k", min.cells = 3, min.features = 200)
pbmc$percent.mt <- PercentageFeatureSet(pbmc, pattern = "^MT-")
pbmc <- subset(pbmc, subset = nFeature_RNA > 200 & nFeature_RNA < 2500 & percent.mt < 5)
pbmc <- NormalizeData(pbmc, normalization.method = "LogNormalize", scale.factor = 1e4)
pbmc <- FindVariableFeatures(pbmc, selection.method = 'vst', nfeatures = 2000)
all.genes <- rownames(pbmc)
pbmc <- ScaleData(pbmc, features = all.genes)
pbmc <- RunPCA(pbmc, features = VariableFeatures(object = pbmc))
pbmc <- FindNeighbors(pbmc, dims = 1:10)

resolution= 2

markers_list<-list()
for (res in seq(0.1,resolution,0.1)) {
  pbmc<- FindClusters(pbmc, resolution = res, algorithm = 3,
                      print.output = FALSE)
  clustername<-paste0("RNA_snn_res.",res)
  Idents(pbmc)<- pbmc[[clustername]] 
  markers_list[[clustername]] <- FindAllMarkers(pbmc, only.pos = TRUE, min.pct = 0.25, logfc.threshold = 0.25)
}

pbmc <- RunUMAP(pbmc, dims = 1:10)

resolution=2

pbmc_cluster=pbmc
metadata_cols<-c("nCount_RNA","nFeature_RNA","orig.ident","percent.mt",paste0("RNA_snn_res.",seq(0,resolution,0.1)))
clustername<-paste0("RNA_snn_res.",resolution)
Idents(pbmc_cluster)<-pbmc_cluster@meta.data[clustername]
pbmc_cluster@meta.data <- pbmc_cluster@meta.data[, which(colnames(pbmc_cluster@meta.data) %in% metadata_cols)]

# clustree plot 
# https://cran.r-project.org/web/packages/clustree/vignettes/clustree.html
clustree(pbmc_cluster, prefix = "RNA_snn_res.")

#UMAP plot
# https://satijalab.org/seurat/reference/findclusters
# https://satijalab.org/seurat/reference/dimplot
DimPlot(pbmc_cluster,label = TRUE, repel = TRUE,label.box=TRUE)+ NoLegend()

# heatmap plot
#https://satijalab.org/seurat/reference/doheatmap
top10 <-markers_list[[clustername]] %>% group_by(cluster)  %>% top_n(n = 10,wt = avg_log2FC)
DoHeatmap(pbmc_cluster, features = top10$gene) + NoLegend()

#table of markers
markers_list[[clustername]] %>% group_by(cluster) %>% slice_max(n = 10, order_by = avg_log2FC) 

• how to decide in number of dimensions for UMAP. Shiny app?

• add "why you need to do this?" and "how to make a decision (interpreting plots/numbres)" to each session

• rotate Dotplot x-axis

• improve genes labels for heatmaps. Like: DoHeatmap(subset_obj, features = top10$gene, size = 3) + NoLegend() + theme(axis.text.x =element_text(size = 5),text =element_text(size = 5) ,axis.text.y =element_text(size = 5))

• addModuleScore to label cell types using markers

#this helps you to annotate your cluster using a list of markers
library(RColorBrewer)
#find markers
cluster5.markers <- FindMarkers(pbmc, ident.1 = 5, logfc.threshold = 0.25, test.use = "roc", only.pos = TRUE)
#pick the top 10
genes_markers<-list(cluster5=rownames(cluster5.markers)[1:10])
# find cells with markers
pbmc<- AddModuleScore( object =pbmc,features = genes_markers,ctrl = 5,name = "cluster_5", search=TRUE)

#color scale for better visualization
plotCol = rev(brewer.pal(n = 7, name = "RdYlBu"))

#notice the name of the cluster has a 1 at the end
FeaturePlot(pbmc,
                  features = "cluster_51", label=TRUE , repel = TRUE, ) +
  scale_color_gradientn(colors = plotCol)

 # label that cell type 
 Ident(pbmc[pbmc$cluster_51>1])="CellTypeX"

[aabarug]

More challenges/break out in the later sections

SingleR challenge:

## Challenge

# See if you can annotate the data with the Monoco reference dataset and whether it improves the cell type annotation resolution`
# Remember you can view the list of references with ls('package:celldex')

Or alternatively - look at the fine labels for the human cell atlas

[LPerlaza]

change all tutorial to KANG data!! we like it better

[aabarug]

Idle thought - might be worth mentioning that Seurat has its own integration strategy? And demo how that works as well and compare with harmony