scRNAseq complementary

DoHeatmap(

  BC.combined,

  features = markers.to.plot,

  idents = "4",

  group.by = "condition",

  group.bar = TRUE,

  group.colors = NULL,

  disp.min = -2.5,

  disp.max = NULL,

  slot = "scale.data",

  assay = NULL,

  label = TRUE,

  size = 5.5,

  hjust = 0,

  angle = 45,

  raster = TRUE,

  draw.lines = TRUE,

  lines.width = NULL,

  group.bar.height = 0.02,

  combine = TRUE

)

library(dplyr)

library(Seurat)

library(patchwork)

library(ggplot2)

dat1 <- Read10X(“C:/Users/~/Desktop/data”) #data(matrix.mtx.gz, barcode.tsv.gz, features.tsv.gz)

dat1 <- CreateSeuratObject(counts = dat1, project = “NT”, min.cells = 3, min.features = 200)

dat1

dat1[[“percent.mt”]] <- PercentageFeatureSet(dat1, pattern = “^MT-“)

VlnPlot(dat1, features = c(“nFeature_RNA”, “nCount_RNA”, “percent.mt”), ncol = 3)

plot1 <- FeatureScatter(dat1, feature1 = “nCount_RNA”, feature2 = “percent.mt”)

plot2 <- FeatureScatter(dat1, feature1 = “nCount_RNA”, feature2 = “nFeature_RNA”)

plot1 + plot2

dat1 <- subset(dat1, subset = nFeature_RNA > 200 & nFeature_RNA < 8000 & percent.mt < 25)

dat1 <- NormalizeData(dat1)

dat1 <- FindVariableFeatures(dat1, selection.method = “vst”, nfeatures = 7000)

all.genes <- rownames(dat1)

dat1 <- ScaleData(dat1, features = all.genes)

dat1 <- RunPCA(dat1, features = VariableFeatures(object = dat1))

DimPlot(dat1, reduction = “pca”)

dat1 <- JackStraw(dat1, num.replicate = 100)

dat1 <- ScoreJackStraw(dat1, dims = 1:20)

JackStrawPlot(dat1, dims = 1:15)

ElbowPlot(dat1)

dat1 <- FindNeighbors(dat1, dims = 1:12)

dat1 <- FindClusters(dat1, resolution = 0.5)

head(Idents(dat1), 5)

dat1 <- RunUMAP(dat1, dims = 1:12)

DimPlot(dat1, reduction = “umap”)

dat1.markers <- FindAllMarkers(dat1, only.pos = TRUE, min.pct = 0.25, logfc.threshold = 0.25)

dat1.markers %>%

    group_by(cluster) %>%

    slice_max(n = 2, order_by = avg_log2FC)

VlnPlot(dat1, features = c(“EPCAM”, “PTPRC”))

VlnPlot(dat1, features = c(“COL1A1”, “COL1A2”, “PDGFRA”, “FAP”))

VlnPlot(dat1, features = c(“ACTA2”, “THY1”))

VlnPlot(dat1, features = c(“CD3D”, “CD8A”))

VlnPlot(dat1, features = c(“PECAM1”, “CD200”))

VlnPlot(dat1, features = c(“CD79A”, ” MS4A1″))

new.cluster.ids <- c(“1”, “2”, “3”, “4”, “5”, “6”, “7”, “8”, “9”, “10”, “11”, “12”, “13”, “14”)

names(new.cluster.ids) <- levels(dat1)

dat1 <- RenameIdents(dat1, new.cluster.ids)

DimPlot(dat1, reduction = “umap”, label = TRUE, pt.size = 0.5) + NoLegend()

library(Seurat)

library(SeuratData)

library(patchwork)

dat1 <- Read10X(“C:/Users/~/Desktop/dat1”)

NT <- CreateSeuratObject(counts = dat1, project = “NT”, min.cells = 200, min.features = 200)

dat2 <- Read10X(“C:/Users/~/Desktop/dat2”)

TT1 <- CreateSeuratObject(counts = dat2, project = “TT1”, min.cells = 200, min.features = 200)

dat3 <- Read10X(“C:/Users/~/Desktop/dat3”)

TT2 <- CreateSeuratObject(counts = dat3, project = “BC3”, min.cells = 200, min.features = 200)

NT[[“percent.mt”]] <- PercentageFeatureSet(NT, pattern = “^MT-“)

TT1[[“percent.mt”]] <- PercentageFeatureSet(TT1, pattern = “^MT-“)

TT2[[“percent.mt”]] <- PercentageFeatureSet(TT2, pattern = “^MT-“)

VlnPlot(NT, features = c(“nFeature_RNA”, “nCount_RNA”, “percent.mt”), ncol = 3)

plot1 <- FeatureScatter(NT, feature1 = “nCount_RNA”, feature2 = “percent.mt”)

plot2 <- FeatureScatter(NT, feature1 = “nCount_RNA”, feature2 = “nFeature_RNA”)

plot1 + plot2

NT <- subset(NT, subset = nFeature_RNA > 200 & nFeature_RNA < 8000 & percent.mt < 20)

VlnPlot(TT1, features = c(“nFeature_RNA”, “nCount_RNA”, “percent.mt”), ncol = 3)

plot1 <- FeatureScatter(TT1, feature1 = “nCount_RNA”, feature2 = “percent.mt”)

plot2 <- FeatureScatter(TT1, feature1 = “nCount_RNA”, feature2 = “nFeature_RNA”)

plot1 + plot2

TT1 <- subset(TT1, subset = nFeature_RNA > 200 & nFeature_RNA < 2500 & percent.mt < 20)

VlnPlot(TT2, features = c(“nFeature_RNA”, “nCount_RNA”, “percent.mt”), ncol = 3)

plot1 <- FeatureScatter(TT2, feature1 = “nCount_RNA”, feature2 = “percent.mt”)

plot2 <- FeatureScatter(TT2, feature1 = “nCount_RNA”, feature2 = “nFeature_RNA”)

plot1 + plot2

TT2 <- subset(TT2, subset = nFeature_RNA > 200 & nFeature_RNA < 15000 & percent.mt < 20)

#Add condition column to metadata

[email protected]$condition <- c(‘NT’)

[email protected]$condition <- c(‘TT1’)

[email protected]$condition <- c(‘TT2’)

dat <- merge(NT, y = TT1, add.cell.ids = c(“NT”, “TT1”), project = “BC”)

dat <- merge(NT, y = c(TT1, TT2), add.cell.ids = c(“NT”, “TT1”, “TT2”), project = “BC”)

BC.list<-SplitObject(dat, split.by = “condition”)

BC.list <- lapply(X = BC.list, FUN = function(x) {

    x <- NormalizeData(x)

    x <- FindVariableFeatures(x, selection.method = “vst”, nfeatures = 2000)

})

features <- SelectIntegrationFeatures(object.list = BC.list)

BC.anchors <- FindIntegrationAnchors(object.list = BC.list, anchor.features = features)

BC.combined <- IntegrateData(anchorset = BC.anchors)

DefaultAssay(BC.combined) <- “integrated”

# Run the standard workflow for visualization and clustering

BC.combined<- ScaleData(BC.combined, verbose = FALSE)

BC.combined<- RunPCA(BC.combined, npcs = 30, verbose = FALSE)

BC.combined<- RunUMAP(BC.combined, reduction = “pca”, dims = 1:30)

BC.combined<- FindNeighbors(BC.combined, reduction = “pca”, dims = 1:30)

BC.combined<- FindClusters(BC.combined, resolution = 0.5)

p1 <- DimPlot(BC.combined, reduction = “umap”, group.by = “condition”)

p2 <- DimPlot(BC.combined, reduction = “umap”, label = TRUE, repel = TRUE)

p1 + p2

DimPlot(BC.combined, reduction = “umap”, split.by = “condition”)

DefaultAssay(BC.combined) <- “RNA”

nk.markers <- FindConservedMarkers(BC.combined, ident.1 = 6, grouping.var = “condition”, verbose = FALSE)

head(nk.markers)

library(ggplot2)

library(cowplot)

theme_set(theme_cowplot())

plots <- VlnPlot(BC.combined, features = c(“COL1A1″, “COL1A2”), split.by = “condition”, idents = “4”,

    pt.size = 0, combine = FALSE)

wrap_plots(plots = plots, ncol = 1)

groups <- sample(c(“NT”, “TT1”), replace = TRUE)

names(groups) <- colnames(data)

data <- AddMetaData(object = data, metadata = groups, col.name = “group”)

BC.list <- SplitObject(data, split.by = “group”)

NT <- PercentageFeatureSet(NT, pattern = “^MT-“, col.name = “percent.mt”)

TT1 <- PercentageFeatureSet(TT1, pattern = “^MT-“, col.name = “percent.mt”)

TT2 <- PercentageFeatureSet(TT2, pattern = “^MT-“, col.name = “percent.mt”)

#find markers for every cluster compared to all remaining cells, report only the positive ones

Idents(object = NT) <- [email protected]$clusters

NT.markers <- FindAllMarkers(NT,min.pct = 0.25, logfc.threshold = 0.25)

clus7.markers <- FindMarkers(liver.combined.sct, ident.1 = ‘7’, min.pct = 0.25)

clus15markers <- FindMarkers(liver.combined.sct, ident.1 = ’15’, min.pct = 0.25)

#Selecting markers from Fibroblast sub populations

VlnPlot(liver.combined.sct, features = “COL1A1”, idents = c(“7″,”15”))

DimPlot(BC.combined, reduction = “umap”, split.by = “condition”, label = TRUE, pt.size = 0.5)

DefaultAssay(BC.combined) <- “RNA”

nk.markers <- FindConservedMarkers(BC.combined, ident.1 = 6, grouping.var = “condition”, verbose = FALSE)

head(nk.markers)

FeaturePlot(BC.combined, features = c(“COL1A1”, “COL1A2”, “ACTA2”, “THY1”, “EPCAM”, “KRT8”, “KRT18”, “CD24”, “PTPRC”), min.cutoff = “q9”)

DimPlot(BC.combined, label = TRUE)

Idents(BC.combined) <- factor(Idents(BC.combined), levels = c(“0”, “1”, “2”, “3”, “4”, “5”, “6”, “7”, “8”, “9”, “10”, “11”, “12”, “13”, “14”))

markers.to.plot <- c(“COL1A1”, “COL1A2”, “ACTA2”, “THY1”, “EPCAM”, “KRT8”, “PTPRC”, “PECAM1”, “CD3D”, “CD8A”, “MS4A1”, “CD79A” , “CD163”, “CD68″”, “JCHAIN”)

DotPlot(BC.combined, features = markers.to.plot, cols = c(“blue”, “red”, “green”), dot.scale = 8, split.by = “condition”) + RotatedAxis()

DotPlot(BC.combined, features = markers.to.plot, cols = c(“blue”, “red”, “green”), dot.scale = 8, split.by = “condition”, idents = c(“5″, “8”, “11”, “13”)) + RotatedAxis()

VlnPlot(BC.combined, features = “COL1A1”, idents = c(“7″,”15”))

FC <- subset(BC.combined, idents = c(“5”, “8”, “13”))

Idents(FC) <- “condition”

avg.FC <- as.data.frame(log1p(AverageExpression(FC, verbose = FALSE)$RNA))

avg.FC$gene <- rownames(avg.FC)

genes.to.label = c(“COL1A1”, “COL1A2”, “ACTA2”, “THY1”, “KRT8”, “KRT18”, “CD24”, “CD44”, “ZEB1”)

p1 <- ggplot(avg.FC, aes(NT, TT1, TT2)) + geom_point() + ggtitle(“Fibroblasts”)

p1 <- LabelPoints(plot = p1, points = genes.to.label, repel = TRUE)

p1 

plots <- VlnPlot(BC.combined, features = c(“COL1A1”, “COL1A2”, “ACTA2”, “THY1”, “KRT8”, “KRT18”, “CD24”, “CD44”, “ZEB1”), split.by = “condition”, idents = “13”, pt.size = 0, combine = FALSE)

wrap_plots(plots = plots, ncol = 1)

plots <- VlnPlot(BC.combined, features = c(“KRT8”, “KRT18”, “CD24”, “CD44”), split.by = “condition”, group.by = “ident”,

    pt.size = 0, combine = FALSE)

wrap_plots(plots = plots, ncol = 1)

plots <- VlnPlot(BC.combined, features = c(“COL1A1”, “COL1A2”), split.by = “condition”, idents = ” Fibroblasts”, pt.size = 0, combine = FALSE)

wrap_plots(plots = plots, ncol = 1)

plots <- VlnPlot(BC.combined, features = c(“ACTA2”, “THY1”), split.by = “condition”, idents = “Fibroblasts”, pt.size = 0, combine = FALSE)

wrap_plots(plots = plots, ncol = 1)

plots <- VlnPlot(BC.combined, features = c(“ZEB1”, “CDH1”), split.by = “condition”, idents = “Fibroblasts”, pt.size = 0, combine = FALSE)

wrap_plots(plots = plots, ncol = 1)

BC.combined <- RenameIdents(BC.combined, `0` = “Epithelial”, `1` = “Epithelial”, `2` = “Epithelial”, `3` = “Epithelial”, `4` = “Fibroblasts”, `5` = “T-cells”, `6` = “Epithelial”, `7` = “Myeloid”, `8` = “Fibroblasts”, `9` = “Endothelial”, `10` = “Fibroblasts”, `11` = “Plasma cells”, `12` = “T-cells”)

DimPlot(BC.combined, label = TRUE)

markers.to.plot <- c(“COL1A1”, “COL1A2”, “ACTA2”, “THY1”, “KRT8”, “KRT18”, “ZEB1”, “COMMD1”)

DotPlot(BC.combined, features = markers.to.plot, cols = c(“blue”, “red”, “green”), dot.scale = 8, split.by = “condition”, idents = “Fibroblasts”) + RotatedAxis()

Fibroblasts1<- subset(BC.combined, idents = “Fibroblasts”)

avg. Fibroblasts<- as.data.frame(log1p(AverageExpression(Fibroblasts1, verbose = FALSE)$RNA))

avg. Fibroblasts$gene <- rownames(avg. Fibroblasts)

Source