scRNAseq

library(Seurat)

library(tibble)

library(magrittr)

library(patchwork)

 

# Read in data

data <- data.table::fread("~/Downloads/~.txt.gz", data.table = F)

 

# Gene names are in the first column so we need to move them to rownames

data <- data %>% 

  column_to_rownames("GENE")

 

# Create Seurat

NT1<-CreateSeuratObject(counts = data, min.cells = 3, min.features = 350, project = "NT1")

data2 <- data.table::fread("C:/Users/~/Downloads/~.txt.gz", data.table = F)

data2 <- data2 %>% column_to_rownames("GENE")

NT2<-CreateSeuratObject(counts = data2, min.cells = 3, min.features = 350, project = "PS")

data3 <- data.table::fread("C:/Users/~/Downloads/~.txt.gz", data.table = F)

data3 <- data3 %>% column_to_rownames("GENE")

TT1<-CreateSeuratObject(counts = data3, min.cells = 3, min.features = 350, project = "PS")

data4 <- data.table::fread("C:/Users/~/Downloads/~.txt.gz", data.table = F)

data4 <- data4 %>% column_to_rownames("GENE")

TT2<-CreateSeuratObject(counts = data4, min.cells = 3, min.features = 350, project = "PS")

NT1[["percent.mt"]] <- PercentageFeatureSet(NT1, pattern = "^MT-")

NT2[["percent.mt"]] <- PercentageFeatureSet(NT2, pattern = "^MT-")

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

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

VlnPlot(NT1, features = c("nFeature_RNA", "nCount_RNA", "percent.mt"), ncol = 3)

VlnPlot(NT2, features = c("nFeature_RNA", "nCount_RNA", "percent.mt"), ncol = 3)

VlnPlot(TT1, features = c("nFeature_RNA", "nCount_RNA", "percent.mt"), ncol = 3) VlnPlot(TT2, features = c("nFeature_RNA", "nCount_RNA", "percent.mt"), ncol = 3)

NT1 <- subset(NT1, subset = nFeature_RNA > 100 & nFeature_RNA < 4000 & percent.mt < 20)

NT2 <- subset(NT2, subset = nFeature_RNA > 100 & nFeature_RNA < 4000 & percent.mt < 20)

TT1 <- subset(TT1, subset = nFeature_RNA > 100 & nFeature_RNA < 4000 & percent.mt < 20)

TT2 <- subset(TT2, subset = nFeature_RNA > 100 & nFeature_RNA < 4000 & percent.mt < 20)

[e-mail]$condition <- c('NT1')

[e-mail]$condition <- c('NT2')

 

[e-mail]$condition <- c('TT1')

[e-mail]$condition <- c('TT2')

dat <- merge(NT1, y = c(NT2, TT1, TT2), add.cell.ids = c("NT1", "NT2", "TT1", "TT2"), project = "PS")

PS.list<-SplitObject(dat, split.by = "condition")

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

    x <- NormalizeData(x)

    x <- FindVariableFeatures(x, selection.method = "vst", nfeatures = 4000)

})

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

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

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

DefaultAssay(PS.combined) <- "integrated"

 

# Run the standard workflow for visualization and clustering

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

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

PS.combined<- RunUMAP(PS.combined, reduction = "pca", dims = 1:30)

PS.combined<- FindNeighbors(PS.combined, reduction = "pca", dims = 1:30)

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

p1 <- DimPlot(PS.combined, reduction = "umap", group.by = "condition")

p2 <- DimPlot(PS.combined, reduction = "umap", label = TRUE, repel = TRUE)

p1 + p2

DimPlot(PS.combined, reduction = "umap", split.by = "condition")

DefaultAssay(PS.combined) <- "RNA"

seurat <- NormalizeData(object = PS.combined, normalization.method = "LogNormalize", scale.factor = 1e4)

seurat <- FindVariableFeatures(seurat, selection.method = "vst", nfeatures = 4000)

top10 <- head(VariableFeatures(seurat), 10)

[1] "MYH11"   "LTF"     "OLFM4"   "SPARCL1" "APOE"    "SPP1"    "DCN"     "GPNMB"   "RGS5"    "ACTA2"

[1] "MMP7"    "OLFM4"   "MYH11"   "CPA3"    "LTF"     "IL1B"    "SPARCL1" "GPNMB"   "KIT"     "SCGB1A1"

nk.markers <- FindConservedMarkers(PS.combined, ident.1 = 6, grouping.var = "condition", verbose = FALSE)

head(nk.markers)

library(ggplot2)

library(cowplot)

theme_set(theme_cowplot())

plots <- VlnPlot(PS.combined, features = c("KRT8", "KRT18", "CD24"), split.by = "condition", group.by = "celltype",

    pt.size = 0, combine = FALSE)

wrap_plots(plots = plots, ncol = 1)

plots <- VlnPlot(PS.combined, features = c("ACTA2", "COL1A1", "KRT19", "PIK3C2A"), split.by = "condition", idents = "9",

                             pt.size = 0, combine = FALSE)

 

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

> 

FeaturePlot(BC.combined, features = c("COMMD1", "NFKB1", "HIF1A", "SLC16A1", "SLC16A3", "TIGAR", "ZEB1", "CDH1", "CAV1"), min.cutoff = "q9")

markers.to.plot <- c("COL1A1", "COL1A2", "PDPN", "ACTA2", "EPCAM", "CDH1", "MSMB", "CAV1", "FTH1", "RPLP1", "HSPB1", "KLK3", "KRT5", "LCN2", "PSCA", "PTPRC", "CD68", "CD163", "CD14", "SLC16A3", "ITGAX", "CD3D", "CD8A", "MS4A1", "CD79A", "VWF", "CD200")

DotPlot(PS.combined, features = markers.to.plot, cols = c("blue", "red", "green", "#00BFC4"), dot.scale = 8, split.by = "condition") + RotatedAxis() 

markers.to.plot <- c("COL1A1", "COL1A2", "PDPN", "ACTA2", "EPCAM", "CDH1", "MSMB", "CAV1", "FTH1", "RPLP1", "HSPB1", "KLK3", "KRT19", "PKM", "OS9", "KLF4", "KLF9", "PIK3C2A", "PTPRC", "CD68", "VWF", "CD200")

> DotPlot(PS.combined, features = markers.to.plot, cols = c("blue", "red", "green", "#00BFC4"), dot.scale = 8, split.by = "condition") + RotatedAxis()

C4 <- read.csv("C:/Users/~/Desktop/gsm5353232TC/nk.markers4.csv")

CC4 <- merge(NTF, C4)

write.csv(CC4, "cluster4.csv")