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)
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()
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
NT@meta.data$condition <- c('NT')
TT1@meta.data$condition <- c('TT1')
TT2@meta.data$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)
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) <- NT@meta.data$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)
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
)
