--- title: "Liver Droplet- Raw data Notebook" output: html_notebook --- ```{r} tissue_of_interest = "Liver" library(here) source("/restricted/projectnb/waxmanlab/kkarri/scRNAseq_data_integration/boilerplate.R") #tiss = load_tissue_droplet(tissue_of_interest) #library(scater) library(dplyr) library(Seurat) library(cowplot) #library(MAST) ########## function load_tissue_droplet############ droplet_metadata_G171B <- read.csv("/net/waxman-server/mnt/data/waxmanlabvm_home/kkarri/G171/Analysis/G171_metadata_droplet_liver.csv", sep=",", header = TRUE) colnames(droplet_metadata_G171B)[1] <- "channel" tissue_metadata_G171B = filter(droplet_metadata_G171B, tissue == tissue_of_interest)[,c('channel','tissue','subtissue','mouse.sex', 'mouse.id')] raw.data <- Read10X("/net/waxman-server/mnt/data/waxmanlabvm_home/kkarri/G171/Analysis/Transcript_Refined/Liver-10X_G171B/") #raw.data <- Read10X("/net/waxman-server/mnt/data/waxmanlabvm_home/kkarri/G171/Analysis/Transcript_downsampled_0.5/Liver-10X_G171B") colnames(raw.data) <- lapply(colnames(raw.data), function(x) paste0(tissue_metadata_G171B$channel[1],'-',x)) meta.data1 = data.frame(row.names = colnames(raw.data)) meta.data1['channel'] = tissue_metadata_G171B$channel[1] rnames = row.names(meta.data1) meta.data1 <- merge(meta.data1, tissue_metadata_G171B, sort = F) row.names(meta.data1) <- rnames # Order the cells alphabetically to ensure consistency. ordered_cell_names = order(colnames(raw.data)) raw.data = raw.data[,ordered_cell_names] meta.data1 = meta.data1[ordered_cell_names,] # Find ERCC's, compute the percent ERCC, and drop them from the raw data. erccs <- grep(pattern = "^ERCC-", x = rownames(x = raw.data), value = TRUE) percent.ercc <- Matrix::colSums(raw.data[erccs, ])/Matrix::colSums(raw.data) ercc.index <- grep(pattern = "^ERCC-", x = rownames(x = raw.data), value = FALSE) raw.data <- raw.data[-ercc.index,] ncRNA.genes <- grep(pattern = "^ncRNA", x = rownames(x = raw.data), value = TRUE) percent.ncRNA <- Matrix::colSums(raw.data[ncRNA.genes, ])/Matrix::colSums(raw.data) KRAB.genes <- grep(pattern = "^KRAB", x = rownames(x = raw.data), value = TRUE) cherry.genes <- grep(pattern = "^mcherry", x = rownames(x = raw.data), value = TRUE) # Create the Seurat object with all the data droplet <- CreateSeuratObject(raw.data) # dropseq droplet <- AddMetaData(object = droplet, meta.data1) droplet@meta.data$tech <- "G171B" #droplet <- SubsetData(droplet,subset.names = c("nGene", "nUMI"), low.thresholds = c(500, 1000)) # old version of seurat #droplet <- subset(droplet, subset = nFeature_RNA > 200 & nCount_RNA > 500) droplet <- NormalizeData(droplet, verbose = FALSE) droplet <- FindVariableFeatures(droplet, selection.method = "vst", nfeatures = 2000) droplet$stim <- "G171B" #lnc5998 <- subset(droplet, subset = `ncRNA-inter-chr7-5998` >1) #lnc5998 <- NormalizeData(lnc5998, verbose = FALSE) #lnc5998 <- FindVariableFeatures(lnc5998, selection.method = "vst", nfeatures = 2000) #lnc5998$stim <- "lnc5998G171B" ### this is tranformation option for develpment SCtransform program ####### droplet <- SCTransform(droplet,verbose =TRUE) droplet <- ScaleData(droplet, verbose = FALSE) droplet <- RunPCA(droplet, npcs = 30, verbose = FALSE) # t-SNE and Clustering droplet <- RunUMAP(droplet, reduction = "pca", dims = 1:25) droplet <- FindNeighbors(droplet, reduction = "pca", dims = 1:25) droplet <- FindClusters(droplet, resolution = 1) droplet <- RunTSNE(droplet, reduction = "pca", dims = 1:20) pd1 <- UMAPPlot(droplet, reduction = "umap", group.by = "stim", label=TRUE, label.size=5) pd1.1 <- UMAPPlot(droplet, reduction = "umap", split.by = "stim",label=TRUE, label.size=5) pd2 <- UMAPPlot(droplet, reduction = "umap", group.by = "mouse.sex") pd3 <- UMAPPlot(droplet, reduction = "umap", label = TRUE, label.size=3) pd4 <- UMAPPlot(droplet, label=TRUE, label.size=6) pd1 <- TSNEPlot(droplet, reduction = "tsne", label=TRUE, label.size=5) DefaultAssay(droplet) <- "RNA" droplet <- NormalizeData(droplet) d1 <- DotPlot(droplet, features = all_genes) lnc5998 <- SCTransform(lnc5998, verbose = TRUE) ####################### droplet G171C ####################################### droplet_metadata_G171C <- read.csv("/net/waxman-server/mnt/data/waxmanlabvm_home/kkarri/G171/Analysis/G171C_metadata_droplet_liver.csv", sep=",", header = TRUE) colnames(droplet_metadata_G171C)[1] <- "channel" tissue_metadata_G171C = filter(droplet_metadata_G171C, tissue == tissue_of_interest)[,c('channel','tissue','subtissue','mouse.sex', 'mouse.id')] raw.data1 <- Read10X("/net/waxman-server/mnt/data/waxmanlabvm_home/kkarri/G171/Analysis/Transcript_Refined/Liver-10X_G171C/") colnames(raw.data1) <- lapply(colnames(raw.data1), function(x) paste0(tissue_metadata_G171C$channel[1],'-',x)) meta.data = data.frame(row.names = colnames(raw.data1)) meta.data['channel'] = tissue_metadata_G171C$channel[1] rnames = row.names(meta.data) meta.data <- merge(meta.data, tissue_metadata_G171C, sort = F) row.names(meta.data) <- rnames # Order the cells alphabetically to ensure consistency. ordered_cell_names1 = order(colnames(raw.data1)) raw.data1 = raw.data1[,ordered_cell_names1] meta.data = meta.data[ordered_cell_names1,] # Find ERCC's, compute the percent ERCC, and drop them from the raw data. erccs1 <- grep(pattern = "^ERCC-", x = rownames(x = raw.data1), value = TRUE) percent.ercc <- Matrix::colSums(raw.data1[erccs1, ])/Matrix::colSums(raw.data1) ercc.index1 <- grep(pattern = "^ERCC-", x = rownames(x = raw.data1), value = FALSE) raw.data1 <- raw.data1[-ercc.index1,] # Create the Seurat object with all the data droplet1 <- CreateSeuratObject(raw.data1) # dropseq droplet1 <- AddMetaData(object = droplet1, meta.data) droplet1@meta.data$tech <- "G171C" #droplet <- SubsetData(droplet,subset.names = c("nGene", "nUMI"), low.thresholds = c(500, 1000)) # old version of seurat #droplet1 <- subset(droplet1, subset = nFeature_RNA > 200 & nCount_RNA > 500) droplet1 <- NormalizeData(droplet1, verbose = FALSE) droplet1 <- FindVariableFeatures(droplet1, selection.method = "vst", nfeatures = 2000) droplet1$stim <- "G171C" droplet1 <- SCTransform(droplet1,verbose =TRUE) ############################################################################### combined.features <- SelectIntegrationFeatures(object.list = list(droplet, droplet1), nfeatures = 2000) combined.list <- PrepSCTIntegration(object.list = list(droplet, droplet1), anchor.features = combined.features, verbose = TRUE) combined.anchors <- FindIntegrationAnchors(object.list = combined.list, normalization.method = "SCT",anchor.features = combined.features, verbose = TRUE) combined.integrated <- IntegrateData(anchorset = combined.anchors, normalization.method = "SCT", verbose = TRUE) combined.integrated <- RunPCA(combined.integrated, verbose = FALSE) combined.integrated <- RunUMAP(combined.integrated, dims = 1:30) plots <- DimPlot(combined.integrated, group.by = c("stim", "mouse.sex"), combine = FALSE) plots <- lapply(X = plots, FUN = function(x) x + theme(legend.position = "top") + guides(color = guide_legend(nrow = 3, byrow = TRUE, override.aes = list(size = 3)))) CombinePlots(plots) #######################3 option 2 standard workflow merging ################################ anchors <- FindIntegrationAnchors(object.list = list(droplet, droplet1), dims = 1:50, anchor.features = 3000) #anchors <- FindIntegrationAnchors(object.list = list(lnc5998, droplet1), dims = 1:50, anchor.features = 3000) combined <- IntegrateData(anchorset = anchors, dims = 1:50) DefaultAssay(combined) <- "integrated" # Run the standard workflow for visualization and clustering combined <- ScaleData(combined, verbose = FALSE) combined <- RunPCA(combined, npcs = 30, verbose = FALSE) # t-SNE and Clustering combined <- RunUMAP(combined, reduction = "pca", dims = 1:25) combined <- FindNeighbors(combined, reduction = "pca", dims = 1:10) combined <- FindClusters(combined, resolution = 0.5 ) combined <- RunTSNE(combined, reduction = "pca", dims = 1:20) # Visualization p1 <- UMAPPlot(combined, reduction = "umap", group.by = "stim", label=TRUE, label.size=5) p1.1 <- UMAPPlot(combined, reduction = "umap", split.by = "stim",label=TRUE, label.size=5) p2 <- UMAPPlot(combined, reduction = "umap", group.by = "mouse.sex") p3 <- UMAPPlot(combined, reduction = "umap", label = TRUE, label.size=3) p4 <- UMAPPlot(combined, label=TRUE, label.size=6) plot_grid(p1, p4) DimPlot(combined, reduction = "umap", split.by = "stim") #cluster.averages <- AverageExpression(combined,use.counts = TRUE,) p5 <- DimPlot(combined, reduction = "tsne", group.by = "stim") p6 <- DimPlot(combined, reduction = "tsne", group.by = "mouse.sex") p7 <- DimPlot(combined, reduction = "tsne", label = TRUE) p8 <- TSNEPlot(combined, label =T) plot_grid(p6, p7,p8) DimPlot(combined, reduction = "tsne", split.by = "stim") ``` diffusion plt code ```{r} # Before running MDS, we first calculate a distance matrix between all pairs of cells. Here we # use a simple euclidean distance metric on all genes, using scale.data as input d <- dist(t(GetAssayData(combined, slot = "scale.data"))) # Run the MDS procedure, k determines the number of dimensions mds <- cmdscale(d = d, k = 2) # cmdscale returns the cell embeddings, we first label the columns to ensure downstream # consistency colnames(mds) <- paste0("MDS_", 1:2) # We will now store this as a custom dimensional reduction called 'mds' combined[["mds"]] <- CreateDimReducObject(embeddings = mds, key = "MDS_", assay = DefaultAssay(combined)) # We can now use this as you would any other dimensional reduction in all downstream functions DimPlot(combined, reduction = "mds", pt.size = 0.5) ``` ```{r} droplet_metadata <- read.csv("/restricted/projectnb/waxmanlab/kkarri/scRNAseq_data_integration/metadata_droplet_liver.csv", sep=",", header = TRUE) colnames(droplet_metadata)[1] <- "channel" tissue_metadata1 = filter(droplet_metadata, tissue == tissue_of_interest)[,c('channel','tissue','subtissue','mouse.sex', 'mouse.id')] raw.data2 <- Read10X("/restricted/projectnb/waxmanlab/kkarri/scRNAseq_data_integration/Refined_cellmatrices/Liver-10X_P4_2/") colnames(raw.data2) <- lapply(colnames(raw.data2), function(x) paste0(tissue_metadata1$channel[1],'_',x)) meta.data2 = data.frame(row.names = colnames(raw.data2)) meta.data2['channel'] = tissue_metadata1$channel[1] if (length(tissue_metadata1$channel) > 1){ # Some tissues, like Thymus and Heart had only one channel for(i in 2:nrow(tissue_metadata1)){ subfolder = paste0("/restricted/projectnb/waxmanlab/kkarri/scRNAseq_data_integration/Refined_cellmatrices/",tissue_of_interest, '-', tissue_metadata1$channel[i]) new.data1 <- Read10X(data.dir = subfolder) colnames(new.data1) <- lapply(colnames(new.data1), function(x) paste0(tissue_metadata1$channel[i],'_', x)) new.metadata1 = data.frame(row.names = colnames(new.data1)) new.metadata1['channel'] = tissue_metadata1$channel[i] raw.data2 = cbind(raw.data2, new.data1) meta.data2 = rbind(meta.data2, new.metadata1) }} rnames = row.names(meta.data2) meta.data2 <- merge(meta.data2, tissue_metadata1, sort = F) row.names(meta.data2) <- rnames # Order the cells alphabetically to ensure consistency. ordered_cell_names = order(colnames(raw.data2)) raw.data2 = raw.data2[,ordered_cell_names] meta.data2 = meta.data2[ordered_cell_names,] # Find ERCC's, compute the percent ERCC, and drop them from the raw data. erccs <- grep(pattern = "^ERCC-", x = rownames(x = raw.data2), value = TRUE) percent.ercc <- Matrix::colSums(raw.data2[erccs, ])/Matrix::colSums(raw.data2) ercc.index <- grep(pattern = "^ERCC-", x = rownames(x = raw.data2), value = FALSE) raw.data2 <- raw.data2[-ercc.index,] # Create the Seurat object with all the data droplet2 <- CreateSeuratObject(raw.data2) # dropseq droplet2 <- AddMetaData(object = droplet2, meta.data2) droplet2@meta.data$tech <- "droplet" #droplet <- SubsetData(droplet,subset.names = c("nGene", "nUMI"), low.thresholds = c(500, 1000)) # old version of seurat droplet2 <- subset(droplet2, subset = nFeature_RNA > 500 & nCount_RNA > 1000) droplet2 <- NormalizeData(droplet2, verbose = FALSE) droplet2 <- FindVariableFeatures(droplet2, selection.method = "vst", nfeatures = 2000) droplet2$stim <- "control" ``` ```{r} DefaultAssay(combined) <- "RNA" #combined <- NormalizeData(combined, verbose = TRUE, normalization.method = "RC", scale.factor = 1e6) combined <- NormalizeData(combined, verbose = TRUE) # Normalize RNA data for visualization purposes if using sctranform #combined.integrated <- NormalizeData(combined.integrated, verbose = FALSE) DotPlot(combined, features = all_genes) FeaturePlot(combined, features = genes_hep_main, min.cutoff = "q9") #hepatocytes subtissplot <- DotPlot(combined, features = c(genes_hep_main, genes_endo, genes_bec_b_immune, genes_kuppfer, genes_nk)) PC <- DotPlot(combined, features = c(genes_hep_main,genes_zones)) NPC <- DotPlot(combined, features = c(genes_endo,genes_kuppfer, genes_nk)) all <- DotPlot(combined, features=c(all_genes)) ### coexpression plots#### f1 <- FeaturePlot(KO.cells, features = c('Cyp2b10','ncRNA-inter-chr7-5998'), reduction = "mds", order = TRUE,split.by = "stim", blend = TRUE,sort.cell = TRUE, max.cutoff = 0.5) ########## this is exact averaging formula ###############33 x <- (AverageExpression(KO.cells, verbose = TRUE, assays = "RNA" ,slot="counts")$RNA) x["ncRNA-inter-chr7-5998",] # G171B G171C #ncRNA-inter-chr7-5998 1.871795 1.091463 ######## #Idents(combined) <- factor(Idents(combined), levels = c(0,1,12)) markers.to.plot <- c("Alb","ncRNA-inter-chr7-5998") DotPlot(combined, features = rev(markers.to.plot), cols = c("blue", "red"), dot.scale = 8, split.by = "stim") + RotatedAxis() FeaturePlot(combined, features = c("Alb", "ncRNA-inter-chr7-5998","Cyp2b10","dSaCas9","KRAB","AAV8-mCherry"), split.by = "stim", max.cutoff = 3, cols = c("grey", "red")) ######################### vlnplot ########################## plots <- VlnPlot(combined, features = c("Alb", "ncRNA-inter-chr7-5998","Cyp2b10","dSaCas9"), split.by = "stim", group.by = "seurat_clusters", pt.size = 0, combine = FALSE) CombinePlots(plots = plots, ncol = 1) plots <- VlnPlot(combined, features = c("Lhx4","Dtna","Fam189a1","Galnt16","Kalrn"), split.by = "stim", group.by = "seurat_clusters", pt.size = 0, combine = FALSE) CombinePlots(plots = plots, ncol = 1) #endothelial DotPlot(combined, features = genes_endo) #zones zones <- DotPlot(combined, features = genes_zones) f1 <- FeaturePlot(combined, features = c('Cyp2e1','Cyp2f2','Ass1'), min.cutoff = "q9", reduction = "tsne") DimPlot(combined, label = TRUE) save(combined, file="Seurat_smart-drop_integrated.Robj") ################# save raw counts from cluster ##################### Idents(combined) <- "stim" ### to avergae out the matrix from KO cells combined.raw.data.0.1 <- as.matrix(GetAssayData(combined, slot = "data")[, WhichCells(combined, ident = 0,idents = "stim")]) combined.raw.data.1 <- as.matrix(GetAssayData(combined, slot = "counts")[, WhichCells(combined, ident = 0)]) combined.raw.data.2 <- as.matrix(GetAssayData(combined, slot = "counts")[, WhichCells(combined, ident = 0)]) #combined.raw.data.[i] <- as.matrix(GetAssayData(combined, slot = "counts")[, WhichCells(combined, ident = 1)]) #combined.raw.data.12 <-as.matrix(GetAssayData(combined, slot = "counts")[, WhichCells(combined, ident = 12)]) #combined.raw.data.1 <- as.matrix(GetAssayData(combined, slot = "counts")) x <- AverageExpression(test.combined,assays = "RNA",add.ident = "stim", slot = "data",use.scale = FALSE, use.counts = FALSE)$RNA #}######## CAR data avg.combined.cells <- (AverageExpression(combined, verbose = FALSE)$RNA) avg.combined.cells$gene <- rownames(avg.combined.cells) CAR_FP <- FeaturePlot(combined, features = c('Cyp2b10','Nr1i3'), reduction = "umap", order = TRUE,split.by = "stim", blend = TRUE,sort.cell = TRUE, max.cutoff = 1, min.cutoff = 0, pt.size = 0.5, repel = TRUE) CAR_DOT_NR <- DotPlot(combined, features = 'Nr1i3', col.min = 0) Cyp2b10_FP <- FeaturePlot(combined, features = 'Cyp2b10', reduction = "umap", min.cutoff = 0) ########### tSNE ################################# combined <- NormalizeData(object = combined) combined <- FindVariableFeatures(combined, selection.method = "vst", nfeatures = 2000) ``` ```{r} KO.cells <- subset(combined, idents = c("0","1","12")) Idents(KO.cells) <- "stim" ### to avergae out the matrix from KO cells raw.data.0 <- as.matrix(GetAssayData(combined, slot = c("counts","data"))[, WhichCells(combined, ident = 0)]) raw.data.1 <- as.matrix(GetAssayData(combined, slot = "counts")[, WhichCells(combined, ident = 1)]) raw.data.2 <-as.matrix(GetAssayData(combined, slot = "counts")[, WhichCells(combined, ident = 2)]) raw.data.3 <-as.matrix(GetAssayData(combined, slot = "counts")[, WhichCells(combined, ident = 3)]) raw.data.4 <-as.matrix(GetAssayData(combined, slot = "counts")[, WhichCells(combined, ident = 4)]) raw.data.5 <-as.matrix(GetAssayData(combined, slot = "counts")[, WhichCells(combined, ident = 5)]) raw.data.6 <-as.matrix(GetAssayData(combined, slot = "counts")[, WhichCells(combined, ident = 6)]) raw.data.7 <-as.matrix(GetAssayData(combined, slot = "counts")[, WhichCells(combined, ident = 7)]) raw.data.8 <-as.matrix(GetAssayData(combined, slot = "counts")[, WhichCells(combined, ident = 8)]) raw.data.9 <-as.matrix(GetAssayData(combined, slot = "counts")[, WhichCells(combined, ident = 9)]) raw.data.10 <-as.matrix(GetAssayData(combined, slot = "counts")[, WhichCells(combined, ident = 10)]) raw.data.11 <-as.matrix(GetAssayData(combined, slot = "counts")[, WhichCells(combined, ident = 11)]) raw.data.12 <-as.matrix(GetAssayData(combined, slot = "counts")[, WhichCells(combined, ident = 12)]) TPMcount0<- cbind(as.matrix(GetAssayData(combined, slot = "counts")[, WhichCells(combined, ident = 0)]), as.matrix(GetAssayData(combined, slot = "data")[, WhichCells(combined, ident = 0)])) TPMcount1<- cbind(as.matrix(GetAssayData(combined, slot = "counts")[, WhichCells(combined, ident = 1)]), as.matrix(GetAssayData(combined, slot = "data")[, WhichCells(combined, ident = 1)])) TPMcount2<- cbind(as.matrix(GetAssayData(combined, slot = "counts")[, WhichCells(combined, ident = 2)]), as.matrix(GetAssayData(combined, slot = "data")[, WhichCells(combined, ident = 2)])) TPMcount3<- cbind(as.matrix(GetAssayData(combined, slot = "counts")[, WhichCells(combined, ident = 3)]), as.matrix(GetAssayData(combined, slot = "data")[, WhichCells(combined, ident = 3)])) TPMcount4<- cbind(as.matrix(GetAssayData(combined, slot = "counts")[, WhichCells(combined, ident = 4)]), as.matrix(GetAssayData(combined, slot = "data")[, WhichCells(combined, ident = 4)])) TPMcount5<- cbind(as.matrix(GetAssayData(combined, slot = "counts")[, WhichCells(combined, ident = 5)]), as.matrix(GetAssayData(combined, slot = "data")[, WhichCells(combined, ident = 5)])) TPMcount6<- cbind(as.matrix(GetAssayData(combined, slot = "counts")[, WhichCells(combined, ident = 6)]), as.matrix(GetAssayData(combined, slot = "data")[, WhichCells(combined, ident = 6)])) TPMcount7<- cbind(as.matrix(GetAssayData(combined, slot = "counts")[, WhichCells(combined, ident = 7)]), as.matrix(GetAssayData(combined, slot = "data")[, WhichCells(combined, ident = 7)])) TPMcount8<- cbind(as.matrix(GetAssayData(combined, slot = "counts")[, WhichCells(combined, ident = 8)]), as.matrix(GetAssayData(combined, slot = "data")[, WhichCells(combined, ident = 8)])) TPMcount9<- cbind(as.matrix(GetAssayData(combined, slot = "counts")[, WhichCells(combined, ident = 9)]), as.matrix(GetAssayData(combined, slot = "data")[, WhichCells(combined, ident = 9)])) TPMcount10<- cbind(as.matrix(GetAssayData(combined, slot = "counts")[, WhichCells(combined, ident = 10)]), as.matrix(GetAssayData(combined, slot = "data")[, WhichCells(combined, ident = 10)])) TPMcount11<- cbind(as.matrix(GetAssayData(combined, slot = "counts")[, WhichCells(combined, ident = 11)]), as.matrix(GetAssayData(combined, slot = "data")[, WhichCells(combined, ident = 11)])) TPMcount12<- cbind(as.matrix(GetAssayData(combined, slot = "counts")[, WhichCells(combined, ident = 12)]), as.matrix(GetAssayData(combined, slot = "data")[, WhichCells(combined, ident = 12)])) write.csv(TPMcount0, "CountResult/counts.TPMcount0.csv") write.csv(TPMcount1, "CountResult/counts.TPMcount1.csv") write.csv(TPMcount2, "CountResult/counts.TPMcount2.csv") write.csv(TPMcount3, "CountResult/counts.TPMcount3.csv") write.csv(TPMcount4, "CountResult/counts.TPMcount4.csv") write.csv(TPMcount5, "CountResult/counts.TPMcount5.csv") write.csv(TPMcount6, "CountResult/counts.TPMcount6.csv") write.csv(TPMcount7, "CountResult/counts.TPMcount7.csv") write.csv(TPMcount8, "CountResult/counts.TPMcount8.csv") write.csv(TPMcount9, "CountResult/counts.TPMcount9.csv") write.csv(TPMcount10, "CountResult/counts.TPMcount10.csv") write.csv(TPMcount11, "CountResult/counts.TPMcount11.csv") write.csv(TPMcount12, "CountResult/counts.TPMcount12.csv") #avg.KO.cells <- log1p(AverageExpression(KO.cells, verbose = FALSE)$RNA) #original code log transformed avg.KO.cells <- (AverageExpression(KO.cells, verbose = FALSE)$RNA) avg.KO.cells$gene <- rownames(avg.KO.cells) genes.to.label= ("ncRNA-inter-chr7-5998") #genes.to.label = c("ISG15", "LY6E", "IFI6", "ISG20", "MX1", "IFIT2", "IFIT1", "CXCL10", "CCL8") p1 <- ggplot(avg.KO.cells, aes(CTRL, STIM)) + geom_point() + ggtitle("CD4 Naive T Cells") p1 <- LabelPoints(plot = p1, points = genes.to.label, repel = TRUE) p2 <- ggplot(avg.cd14.mono, aes(CTRL, STIM)) + geom_point() + ggtitle("CD14 Monocytes") p2 <- LabelPoints(plot = p2, points = genes.to.label, repel = TRUE) plot_grid(p1, p2) ``` ```{r} DefaultAssay(combined) <- "RNA" KO.cells.0 <- subset(combined, idents = c("0")) Idents(KO.cells.0) <- "All0" DefaultAssay(KO.cells.0) <- "RNA" KO.cells.0.lnc5998 <- subset(KO.cells.0 , subset = `ncRNA-inter-chr7-5998` >0) KO.cells.0.lnc5998 <- NormalizeData(KO.cells.0.lnc5998, verbose = FALSE) KO.cells.0.lnc5998 <- FindVariableFeatures(KO.cells.0.lnc5998, selection.method = "vst", nfeatures = 2000) KO.cells.0.lnc5998 <- ScaleData(KO.cells.0.lnc5998, verbose = FALSE) KO.cells.0.lnc5998 <- RunPCA(KO.cells.0.lnc5998, npcs = 30, verbose = FALSE) # t-SNE and Clustering KO.cells.0.lnc5998 <- RunUMAP(KO.cells.0.lnc5998, reduction = "pca", dims = 1:20) KO.cells.0.lnc5998 <- FindNeighbors(KO.cells.0.lnc5998, reduction = "pca", dims = 1:20) KO.cells.0.lnc5998 <- FindClusters(KO.cells.0.lnc5998, resolution = 0.5) KO.cells.0.lnc5998 <- RunTSNE(KO.cells.0.lnc5998, reduction = "pca", dims = 1:20) # Visualization p1 <- UMAPPlot(KO.cells.0.lnc5998, reduction = "umap", split.by = "tech", label=TRUE, label.size=5) p1.1 <- UMAPPlot(KO.cells.0.lnc5998, reduction = "umap", group.by = "tech",label=TRUE, label.size=5) p2 <- UMAPPlot(KO.cells.0.lnc5998, reduction = "umap", group.by = "mouse.sex") p3 <- UMAPPlot(KO.cells.0.lnc5998, reduction = "umap", label = TRUE, label.size=3) p4 <- UMAPPlot(KO.cells.0.lnc5998, label=TRUE, label.size=6) KO.cell.0.lnc5998_DE1 <- FindMarkers(KO.cells.0.lnc5998, ident.1 = c(1,2), ident.2 = 0,verbose = TRUE, logfc.threshold = FALSE,min.pct = FALSE) KO.cell.0.lnc5998_AllMarkers <- FindAllMarkers(KO.cells.0.lnc5998, only.pos = TRUE, min.pct = 0.25, thresh.use = 0.25) KO.cell.0.lnc5998_AllMarkers %>% group_by(cluster) %>% top_n(2, avg_logFC) write.csv(KO.cell.0.lnc5998_DE1, "Analysis/CountResult/Markers/KO.cell.0.lnc5998_DE1") write.csv(KO.cell.0.lnc5998_AllMarkers, "Analysis/CountResult/Markers/KO.cell.0.lnc5998_AllMarkers") KO.cells.0.lnc5998$celltype.stim <- paste(Idents(KO.cells.0.lnc5998), KO.cells.0.lnc5998$tech, sep = "_") KO.cells.0.lnc5998$celltype <- Idents(KO.cells.0.lnc5998) Idents(KO.cells.0.lnc5998) <- "celltype.stim" #############3 KO cells that express lnc5998 ################# DefaultAssay(combined) <- "RNA" KO.cells.1 <- subset(combined, idents = c("1")) Idents(KO.cells.1) <- "All0" DefaultAssay(KO.cells.0) <- "RNA" KO.cells.1.lnc5998 <- subset(KO.cells.1 , subset = `ncRNA-inter-chr7-5998` >0) KO.cells.1.lnc5998 <- NormalizeData(KO.cells.1.lnc5998, verbose = FALSE) KO.cells.1.lnc5998 <- FindVariableFeatures(KO.cells.1.lnc5998, selection.method = "vst", nfeatures = 2000) KO.cells.1.lnc5998 <- ScaleData(KO.cells.1.lnc5998, verbose = FALSE) KO.cells.1.lnc5998 <- RunPCA(KO.cells.1.lnc5998, npcs = 30, verbose = FALSE) # t-SNE and Clustering KO.cells.1.lnc5998 <- RunUMAP(KO.cells.1.lnc5998, reduction = "pca", dims = 1:20) KO.cells.1.lnc5998 <- FindNeighbors(KO.cells.1.lnc5998, reduction = "pca", dims = 1:20) KO.cells.1.lnc5998 <- FindClusters(KO.cells.1.lnc5998, resolution = 0.5) KO.cells.1.lnc5998 <- RunTSNE(KO.cells.1.lnc5998, reduction = "pca", dims = 1:20) # Visualization p11 <- UMAPPlot(KO.cells.1.lnc5998, reduction = "umap", split.by = "tech", label=TRUE, label.size=5) p11.1 <- UMAPPlot(KO.cells.1.lnc5998, reduction = "umap", group.by = "tech",label=TRUE, label.size=5) p11.2 <- UMAPPlot(KO.cells.1.lnc5998, reduction = "umap", label = TRUE, label.size=3) KO.cells.1.lnc5998$celltype.stim <- paste(Idents(KO.cells.1.lnc5998), KO.cells.1.lnc5998$stim, sep = "_") KO.cells.1.lnc5998$celltype <- Idents(KO.cells.1.lnc5998) Idents(KO.cells.1.lnc5998) <- "celltype.stim" KO.cell.1.lnc5998_DE1 <- FindMarkers(KO.cells.1.lnc5998, ident.1 = c('0_G171B','1_G171B'), ident.2 = c('0_G171C','1_G171C'),verbose = TRUE, logfc.threshold = FALSE,min.pct = FALSE) KO.cell.1.lnc5998_DE2 <- FindMarkers(KO.cells.1.lnc5998, ident.1 = c('1_G171B'), ident.2 = c('0_G171C'),verbose = TRUE, logfc.threshold = FALSE,min.pct = FALSE) KO.cell.1.lnc5998_DE3 <- FindMarkers(KO.cells.1.lnc5998, ident.1 = c('0_G171B'), ident.2 = c('0_G171C'),verbose = TRUE, logfc.threshold = FALSE,min.pct = FALSE) KO.cell.1.lnc5998_AllMarkers <- FindAllMarkers(KO.cells.1.lnc5998, only.pos = TRUE, min.pct = 0.25, thresh.use = 0.25) KO.cell.1.lnc5998_AllMarkers %>% group_by(cluster) %>% top_n(2, avg_logFC) write.csv(KO.cell.1.lnc5998_DE1, "Analysis/CountResult/Markers/KO.cell.1.lnc5998_DE1") write.csv(KO.cell.1.lnc5998_AllMarkers, "Analysis/CountResult/Markers/KO.cell.1.lnc5998_AllMarkers") write.csv(KO.cell.1.lnc5998_DE2, "Analysis/CountResult/Markers/KO.cell.1.lnc5998_DE2_1-G171B_0-G171C") write.csv(KO.cell.1.lnc5998_DE3, "Analysis/CountResult/Markers/KO.cell.1.lnc5998_DE3_0-G171B_0-G171C") ``` ######### cell clusters 0 and cluster 1 from main UMAP that do not epxress lnc5998 in G171B and G171C ################################## ```{r} DefaultAssay(combined) <- "RNA" KO.cells.0 <- subset(combined, idents = c("0")) Idents(KO.cells.0) <- "All0" DefaultAssay(KO.cells.0) <- "RNA" KO.cells.0.null <- subset(KO.cells.0 , subset = `ncRNA-inter-chr7-5998`== 0) KO.cells.0.null <- NormalizeData(KO.cells.0.null, verbose = FALSE) KO.cells.0.null <- FindVariableFeatures(KO.cells.0.null, selection.method = "vst", nfeatures = 2000) KO.cells.0.null <- ScaleData(KO.cells.0.null, verbose = FALSE) KO.cells.0.null <- RunPCA(KO.cells.0.null, npcs = 30, verbose = FALSE) # t-SNE and Clustering KO.cells.0.null <- RunUMAP(KO.cells.0.null, reduction = "pca", dims = 1:20) KO.cells.0.null <- FindNeighbors(KO.cells.0.null, reduction = "pca", dims = 1:20) KO.cells.0.null <- FindClusters(KO.cells.0.null, resolution = 0.5) KO.cells.0.null <- RunTSNE(KO.cells.0.null, reduction = "pca", dims = 1:20) # Visualization pnull1 <- UMAPPlot(KO.cells.0.null, reduction = "umap", split.by = "tech", label=TRUE, label.size=5) pnull1.1 <- UMAPPlot(KO.cells.0.null, reduction = "umap", group.by = "tech",label=TRUE, label.size=5) pnull2 <- UMAPPlot(KO.cells.0.null, reduction = "umap", group.by = "mouse.sex") pnull3 <- UMAPPlot(KO.cells.0.null, reduction = "umap", label = TRUE, label.size=3) pnull4 <- UMAPPlot(KO.cells.0.null, label=TRUE, label.size=6) KO.cell.0.null_DE1 <- FindMarkers(KO.cells.0.null, ident.1 = c(1,2,3), ident.2 = 0,verbose = TRUE, logfc.threshold = FALSE,min.pct = FALSE) KO.cell.0.null_AllMarkers <- FindAllMarkers(KO.cells.0.null, only.pos = TRUE, min.pct = 0.25, thresh.use = 0.25) KO.cell.0.null_AllMarkers %>% group_by(cluster) %>% top_n(2, avg_logFC) write.csv(KO.cell.0.null_DE1, "CountResult/Markers/KO_cell_0_null_DE1_All_G171B_vs_G171C") write.csv(KO.cell.0.null_AllMarkers, "CountResult/Markers/KO_cell_0_null_AllMarkers") KO.cells.0.null$celltype.stim <- paste(Idents(KO.cells.0.null), KO.cells.0.null$tech, sep = "_") KO.cells.0.null$celltype <- Idents(KO.cells.0.null) Idents(KO.cells.0.null) <- "celltype.stim" #############3 KO cells that express lnc5998 ################# DefaultAssay(combined) <- "RNA" KO.cells.1 <- subset(combined, idents = c("1")) Idents(KO.cells.1) <- "All0" DefaultAssay(KO.cells.1) <- "RNA" KO.cells.1.null <- subset(KO.cells.1 , subset = `ncRNA-inter-chr7-5998` ==0) KO.cells.1.null <- NormalizeData(KO.cells.1.null, verbose = FALSE) KO.cells.1.null <- FindVariableFeatures(KO.cells.1.null, selection.method = "vst", nfeatures = 2000) KO.cells.1.null <- ScaleData(KO.cells.1.null, verbose = FALSE) KO.cells.1.null <- RunPCA(KO.cells.1.null, npcs = 30, verbose = FALSE) # t-SNE and Clustering KO.cells.1.null <- RunUMAP(KO.cells.1.null, reduction = "pca", dims = 1:20) KO.cells.1.null <- FindNeighbors(KO.cells.1.null, reduction = "pca", dims = 1:20) KO.cells.1.null <- FindClusters(KO.cells.1.null, resolution = 0.5) KO.cells.1.null <- RunTSNE(KO.cells.1.null, reduction = "pca", dims = 1:20) # Visualization p11 <- UMAPPlot(KO.cells.1.null, reduction = "umap", split.by = "tech", label=TRUE, label.size=5) p11.1 <- UMAPPlot(KO.cells.1.null, reduction = "umap", group.by = "tech",label=TRUE, label.size=5) p11.2 <- UMAPPlot(KO.cells.1.null, reduction = "umap", label = TRUE, label.size=3) KO.cells.1.null$celltype.stim <- paste(Idents(KO.cells.1.null), KO.cells.1.null$stim, sep = "_") KO.cells.1.null$celltype <- Idents(KO.cells.1.null) Idents(KO.cells.1.null) <- "celltype.stim" KO.cell.1.null_DE1 <- FindMarkers(KO.cells.1.null, ident.1 = c('0_G171B','1_G171B'), ident.2 = c('0_G171C','1_G171C'),verbose = TRUE, logfc.threshold = FALSE,min.pct = FALSE) KO.cell.1.null_DE2 <- FindMarkers(KO.cells.1.null, ident.1 = c('1_G171B'), ident.2 = c('1_G171C'),verbose = TRUE, logfc.threshold = FALSE,min.pct = FALSE) KO.cell.1.null_DE3 <- FindMarkers(KO.cells.1.null, ident.1 = c('0_G171B'), ident.2 = c('0_G171C'),verbose = TRUE, logfc.threshold = FALSE,min.pct = FALSE) KO.cell.1.null_DE4 <- FindMarkers(KO.cells.1.null, ident.1 = c('1_G171B'), ident.2 = c('0_G171C'),verbose = TRUE, logfc.threshold = FALSE,min.pct = FALSE) KO.cell.1.null_DE5 <- FindMarkers(KO.cells.1.null, ident.1 = c('0_G171B'), ident.2 = c('1_G171C'),verbose = TRUE, logfc.threshold = FALSE,min.pct = FALSE) KO.cell.1.null_AllMarkers <- FindAllMarkers(KO.cells.1.null, only.pos = TRUE, min.pct = 0.25, thresh.use = 0.25) KO.cell.1.null_AllMarkers %>% group_by(cluster) %>% top_n(2, avg_logFC) write.csv(KO.cell.1.null_DE1, "CountResult/Markers/KO_cell_1_null_DE1_all_G171B_vs_G171C") write.csv(KO.cell.1.null_AllMarkers, "CountResult/Markers/KO.cell.1.null_AllMarkers") write.csv(KO.cell.1.null_DE2, "CountResult/Markers/KO_cell_1_null_DE2_1-G171B_1-G171C") write.csv(KO.cell.1.null_DE3, "CountResult/Markers/KO_cell_1_null_DE3_0-G171B_0-G171C") write.csv(KO.cell.1.null_DE4, "CountResult/Markers/KO_cell_1_null_DE4_1-G171B_0-G171C") write.csv(KO.cell.1.null_DE5, "CountResult/Markers/KO_cell_1_null_DE5_0-G171B_1-G171C") ``` All cell clusters that express lnc5998 ```{r} DefaultAssay(combined) <- "RNA" KO.cells.1.all <- subset(combined, idents = c("1")) Idents(KO.cells.1) <- "All1" DefaultAssay(KO.cells.1.all) <- "RNA" KO.cells.1.all <- NormalizeData(KO.cells.1.all, verbose = FALSE) KO.cells.1.all <- FindVariableFeatures(KO.cells.1.all, selection.method = "vst", nfeatures = 2000) KO.cells.1.all <- ScaleData(KO.cells.1.all, verbose = FALSE) KO.cells.1.all <- RunPCA(KO.cells.1.all, npcs = 30, verbose = FALSE) # t-SNE and Clustering KO.cells.1.all <- RunUMAP(KO.cells.1.all, reduction = "pca", dims = 1:20) KO.cells.1.all <- FindNeighbors(KO.cells.1.all, reduction = "pca", dims = 1:20) KO.cells.1.all <- FindClusters(KO.cells.1.all, resolution = 0.5) KO.cells.1.all <- RunTSNE(KO.cells.1.all, reduction = "pca", dims = 1:20) # Visualization pall1 <- UMAPPlot(KO.cells.1.all, reduction = "umap", split.by = "tech", label=TRUE, label.size=5) pall1.1 <- UMAPPlot(KO.cells.1.all, reduction = "umap", group.by = "tech",label=TRUE, label.size=5) pall1.2 <- UMAPPlot(KO.cells.1.all, reduction = "umap", label = TRUE, label.size=3) KO.cells.1.all$celltype.stim <- paste(Idents(KO.cells.1.all), KO.cells.1.all$stim, sep = "_") KO.cells.1.all$celltype <- Idents(KO.cells.1.all) Idents(KO.cells.1.all) <- "celltype.stim" KO.cell.1.all_DE1 <- FindMarkers(KO.cells.1.all, ident.1 = c('0_G171B','1_G171B', '2_G171B','3_G171B'), ident.2 = c('0_G171C','1_G171C','2_G171C','3_G171C'),verbose = TRUE, logfc.threshold = FALSE,min.pct = FALSE) KO.cell.1.all_DE2 <- FindMarkers(KO.cells.1.all, ident.1 = c('1_G171B'), ident.2 = c('1_G171C'),verbose = TRUE, logfc.threshold = FALSE,min.pct = FALSE) KO.cell.1.all_DE3 <- FindMarkers(KO.cells.1.all, ident.1 = c('0_G171B'), ident.2 = c('0_G171C'),verbose = TRUE, logfc.threshold = FALSE,min.pct = FALSE) KO.cell.1.all_DE4 <- FindMarkers(KO.cells.1.all, ident.1 = c('1_G171B'), ident.2 = c('0_G171C'),verbose = TRUE, logfc.threshold = FALSE,min.pct = FALSE) KO.cell.1.all_DE5 <- FindMarkers(KO.cells.1.all, ident.1 = c('0_G171B'), ident.2 = c('1_G171C'),verbose = TRUE, logfc.threshold = FALSE,min.pct = FALSE) KO.cell.1.all_AllMarkers <- FindAllMarkers(KO.cells.1.all, only.pos = TRUE, min.pct = 0.25, thresh.use = 0.25) KO.cell.1.all_AllMarkers %>% group_by(cluster) %>% top_n(2, avg_logFC) write.csv(KO.cell.1.all_DE1, "CountResult/Markers/KO_cell_1_all_DE1_all_G171B_vs_G171C") write.csv(KO.cell.1.all_AllMarkers, "CountResult/Markers/KO.cell.1.all_AllMarkers") write.csv(KO.cell.1.all_DE2, "CountResult/Markers/KO_cell_1_all_DE2_1-G171B_1-G171C") write.csv(KO.cell.1.all_DE3, "CountResult/Markers/KO_cell_1_all_DE3_0-G171B_0-G171C") write.csv(KO.cell.1.all_DE4, "CountResult/Markers/KO_cell_1_all_DE4_1-G171B_0-G171C") write.csv(KO.cell.1.all_DE5, "CountResult/Markers/KO_cell_1_all_DE5_0-G171B_1-G171C") ###################3 cluster 00 #################### DefaultAssay(combined) <- "RNA" KO.cells.0.all <- subset(combined, idents = c("0")) Idents(KO.cells.0.all) <- "All0" DefaultAssay(KO.cells.0.all) <- "RNA" KO.cells.0.all <- NormalizeData(KO.cells.0.all, verbose = FALSE) KO.cells.0.all <- FindVariableFeatures(KO.cells.0.all, selection.method = "vst", nfeatures = 2000) KO.cells.0.all <- ScaleData(KO.cells.0.all, verbose = FALSE) KO.cells.0.all <- RunPCA(KO.cells.0.all, npcs = 30, verbose = FALSE) # t-SNE and Clustering KO.cells.0.all <- RunUMAP(KO.cells.0.all, reduction = "pca", dims = 1:20) KO.cells.0.all <- FindNeighbors(KO.cells.0.all, reduction = "pca", dims = 1:20) KO.cells.0.all <- FindClusters(KO.cells.0.all, resolution = 0.5) KO.cells.0.all <- RunTSNE(KO.cells.0.all, reduction = "pca", dims = 1:20) # Visualization pall01 <- UMAPPlot(KO.cells.0.all, reduction = "umap", split.by = "tech", label=TRUE, label.size=5) pall0.1 <- UMAPPlot(KO.cells.0.all, reduction = "umap", group.by = "tech",label=TRUE, label.size=5) pall0.2 <- UMAPPlot(KO.cells.0.all, reduction = "umap", group.by = "mouse.sex") pall0.3 <- UMAPPlot(KO.cells.0.all, reduction = "umap", label = TRUE, label.size=3) pall0.4 <- UMAPPlot(KO.cells.0.all, label=TRUE, label.size=6) KO.cell.0.all_DE1 <- FindMarkers(KO.cells.0.all, ident.1 = c('0_G171B','1_G171B','2_G171B','3_G171B','4_G171B'), ident.2 = c('0_G171C','2_G171C'),verbose = TRUE, logfc.threshold = FALSE,min.pct = FALSE) KO.cell.0.all_AllMarkers <- FindAllMarkers(KO.cells.0.all, only.pos = TRUE, min.pct = 0.25, thresh.use = 0.25) KO.cell.0.all_AllMarkers %>% group_by(cluster) %>% top_n(2, avg_logFC) write.csv(KO.cell.0.all_DE1, "CountResult/Markers/KO_cell_0_all_DE1_All_G171B_vs_G171C") write.csv(KO.cell.0.all_AllMarkers, "CountResult/Markers/KO_cell_0_all_AllMarkers") KO.cells.0.all$celltype.stim <- paste(Idents(KO.cells.0.all), KO.cells.0.all$tech, sep = "_") KO.cells.0.all$celltype <- Idents(KO.cells.0.all) Idents(KO.cells.0.all) <- "celltype.stim" ######################### cluster 2 ############################## DefaultAssay(combined) <- "RNA" KO.cells.2.all <- subset(combined, idents = c("2")) Idents(KO.cells.2.all) <- "All2" DefaultAssay(KO.cells.2.all) <- "RNA" KO.cells.2.all <- NormalizeData(KO.cells.2.all, verbose = FALSE) KO.cells.2.all <- FindVariableFeatures(KO.cells.2.all, selection.method = "vst", nfeatures = 2000) KO.cells.2.all <- ScaleData(KO.cells.2.all, verbose = FALSE) KO.cells.2.all <- RunPCA(KO.cells.2.all, npcs = 30, verbose = FALSE) # t-SNE and Clustering KO.cells.2.all <- RunUMAP(KO.cells.2.all, reduction = "pca", dims = 1:20) KO.cells.2.all <- FindNeighbors(KO.cells.2.all, reduction = "pca", dims = 1:20) KO.cells.2.all <- FindClusters(KO.cells.2.all, resolution = 0.5) KO.cells.2.all <- RunTSNE(KO.cells.2.all, reduction = "pca", dims = 1:20) # Visualization pall21 <- UMAPPlot(KO.cells.2.all, reduction = "umap", split.by = "tech", label=TRUE, label.size=5) pall2.1 <- UMAPPlot(KO.cells.2.all, reduction = "umap", group.by = "tech",label=TRUE, label.size=5) pall2.2 <- UMAPPlot(KO.cells.2.all, reduction = "umap", group.by = "mouse.sex") pall2.3 <- UMAPPlot(KO.cells.2.all, reduction = "umap", label = TRUE, label.size=3) pall2.4 <- UMAPPlot(KO.cells.2.all, label=TRUE, label.size=6) KO.cell.2.all_DE1 <- FindMarkers(KO.cells.2.all, ident.1 = c(1,2,3), ident.2 = 0,verbose = TRUE, logfc.threshold = FALSE,min.pct = FALSE) KO.cell.2.all_AllMarkers <- FindAllMarkers(KO.cells.2.all, only.pos = TRUE, min.pct = 0.25, thresh.use = 0.25) KO.cell.2.all_AllMarkers %>% group_by(cluster) %>% top_n(2, avg_logFC) write.csv(KO.cell.2.all_DE1, "CountResult/Markers/KO_cell_0_all_DE1_All_G171B_vs_G171C") write.csv(KO.cell.2.all_AllMarkers, "CountResult/Markers/KO_cell_0_all_AllMarkers") KO.cells.2.all$celltype.stim <- paste(Idents(KO.cells.2.all), KO.cells.2.all$tech, sep = "_") KO.cells.2.all$celltype <- Idents(KO.cells.2.all) Idents(KO.cells.2.all) <- "celltype.stim" ``` ```{r} droplet_metadata <- read.csv("/restricted/projectnb/waxmanlab/kkarri/scRNAseq_data_integration/metadata_droplet_liver.csv", sep=",", header = TRUE) colnames(droplet_metadata)[1] <- "channel" tissue_metadata = filter(droplet_metadata, tissue == tissue_of_interest)[,c('channel','tissue','subtissue','mouse.sex', 'mouse.id')] raw.data <- Read10X("/restricted/projectnb/waxmanlab/kkarri/scRNAseq_data_integration/Refined_cellmatrices/Liver-10X_P4_2/") colnames(raw.data) <- lapply(colnames(raw.data), function(x) paste0(tissue_metadata$channel[1],'_',x)) meta.data1 = data.frame(row.names = colnames(raw.data)) meta.data1['channel'] = tissue_metadata$channel[1] if (length(tissue_metadata$channel) > 1){ # Some tissues, like Thymus and Heart had only one channel for(i in 2:nrow(tissue_metadata)){ subfolder = paste0("/restricted/projectnb/waxmanlab/kkarri/scRNAseq_data_integration/Refined_cellmatrices/",tissue_of_interest, '-', tissue_metadata$channel[i]) new.data1 <- Read10X(data.dir = subfolder) colnames(new.data1) <- lapply(colnames(new.data1), function(x) paste0(tissue_metadata$channel[i],'_', x)) new.metadata1 = data.frame(row.names = colnames(new.data1)) new.metadata1['channel'] = tissue_metadata$channel[i] raw.data = cbind(raw.data, new.data1) meta.data1 = rbind(meta.data1, new.metadata1) } } rnames = row.names(meta.data1) meta.data1 <- merge(meta.data1, tissue_metadata, sort = F) row.names(meta.data1) <- rnames # Order the cells alphabetically to ensure consistency. ordered_cell_names = order(colnames(raw.data)) raw.data = raw.data[,ordered_cell_names] meta.data1 = meta.data1[ordered_cell_names,] # Find ERCC's, compute the percent ERCC, and drop them from the raw data. erccs <- grep(pattern = "^ERCC-", x = rownames(x = raw.data), value = TRUE) percent.ercc <- Matrix::colSums(raw.data[erccs, ])/Matrix::colSums(raw.data) ercc.index <- grep(pattern = "^ERCC-", x = rownames(x = raw.data), value = FALSE) raw.data <- raw.data[-ercc.index,] # Create the Seurat object with all the data droplet <- CreateSeuratObject(raw.data) # dropseq droplet <- AddMetaData(object = droplet, meta.data1) droplet@meta.data$tech <- "droplet" #n.pcs = 10 #droplet <- SubsetData(droplet,subset.names = c("nGene", "nUMI"), low.thresholds = c(500, 1000)) # old version of seurat droplet <- subset(droplet, subset = nFeature_RNA > 500 & nCount_RNA > 1000) droplet <- NormalizeData(droplet, verbose = FALSE) droplet <- FindVariableFeatures(droplet, selection.method = "vst", nfeatures = 2000) droplet <- ScaleData(droplet, verbose = FALSE) #droplet <- RunPCA(droplet, npcs = 10, verbose = FALSE) droplet$stim <- "droplet" # droplet <- ScaleData(droplet, verbose = FALSE) droplet <- RunPCA(droplet, npcs = 30, verbose = FALSE) droplet <- RunUMAP(droplet, reduction = "pca", dims = 1:25) droplet <- FindNeighbors(droplet, reduction = "pca", dims = 1:10) droplet <- FindClusters(droplet, resolution = 0.5 ) p1<- UMAPPlot(droplet, reduction = "umap", group.by = "channel", label=TRUE, label.size=5) p2 <- UMAPPlot(droplet, label=TRUE, label.size=6) p3<- UMAPPlot(droplet, reduction = "umap", group.by = "mouse.sex", label=TRUE, label.size=5) res.used <- 1 droplet <- FindClusters(object = droplet, reduction.type = "pca", dims.use = 1:n.pcs, resolution = res.used, print.output = 0, save.SNN = TRUE, force.recalc = TRUE) droplet <- RunTSNE(object = droplet, dims.use = 1:n.pcs, seed.use = 10, perplexity=30) TSNEPlot(object = droplet, do.label = T, pt.size = 1.2, label.size = 4) ``` ```{r} KO.cells$celltype.stim <- paste(Idents(KO.cells), KO.cells$stim, sep = "_") KO.cells$celltype <- Idents(KO.cells) Idents(KO.cells) <- "celltype.stim" response3 <- FindMarkers(KO.cells, ident.1 = c("1_G171B","0_G171B","2_G171B"), ident.2 = c("1_G171C", "0_G171C","2_G171C"), verbose = TRUE, test.use = "MAST", logfc.threshold = FALSE,min.pct = FALSE) head(response3, n = 15) KO.cell.0.lnc5998_DE1 <- FindMarkers(KO.cells.0.lnc5998, ident.1 = "4_G171C", ident.2 = "4_G171B",verbose = TRUE, logfc.threshold = FALSE,min.pct = FALSE) KO.cells.1 <- subset(combined, idents = c("1")) Idents(KO.cells.1) <- "All1" ``` ``` ```{r} KO.cells <- RunUMAP(KO.cells, reduction = "pca", dims = 1:20 ) KO.cells <- FindNeighbors(KO.cells, reduction = "pca", dims = 1:20) KO.cells <- FindClusters(KO.cells, resolution = 0.5 ) KO.cells <- RunTSNE(KO.cells, reduction = "pca", dims = 1:20) Hep.cells <- RunUMAP(Hep.cells, reduction = "pca", dims = 1:20 ) Hep.cells <- FindNeighbors(Hep.cells, reduction = "pca", dims = 1:20) Hep.cells <- FindClusters(Hep.cells, resolution = 0.5 ) Hep.cells <- RunTSNE(Hep.cells, reduction = "pca", dims = 1:20) # Visualization p1 <- UMAPPlot(KO.cells, reduction = "umap", group.by = "stim") p2 <- UMAPPlot(KO.cells, reduction = "umap", group.by = "mouse.sex") p3 <- UMAPPlot(KO.cells, reduction = "umap", label = TRUE) p4 <- UMAPPlot(KO.cells, label=TRUE) plot_grid(p1,p4) DimPlot(KO.cells, reduction = "umap", split.by = "stim") #hepatocyte cells p5 <- UMAPPlot(Hep.cells, reduction = "umap", group.by = "stim") p6 <- UMAPPlot(Hep.cells, reduction = "umap", group.by = "mouse.sex") p7 <- UMAPPlot(Hep.cells, reduction = "umap", label = TRUE) p8 <- UMAPPlot(Hep.cells, label=TRUE) plot_grid(p5,p8) DimPlot(KO.cells, reduction = "umap", split.by = "stim") raw.data.KO.0 <- as.matrix(GetAssayData(KO.cells, slot = "counts")[, WhichCells(KO.cells, ident = 0)]) raw.data.KO.1 <- as.matrix(GetAssayData(KO.cells, slot = "counts")[, WhichCells(KO.cells, ident = 1)]) raw.data.KO.2 <-as.matrix(GetAssayData(KO.cells, slot = "counts")[, WhichCells(KO.cells, ident = 2)]) raw.data.KO.3 <-as.matrix(GetAssayData(KO.cells, slot = "counts")[, WhichCells(KO.cells, ident = 3)]) raw.data.KO.4 <-as.matrix(GetAssayData(KO.cells, slot = "counts")[, WhichCells(KO.cells, ident = 4)]) raw.data.KO.5 <-as.matrix(GetAssayData(KO.cells, slot = "counts")[, WhichCells(KO.cells, ident = 5)]) write.csv(raw.data.KO.0, "CountResult/Markers/raw.data.KO.0") write.csv(raw.data.KO.1, "CountResult/Markers/raw.data.KO.1") write.csv(raw.data.KO.2, "CountResult/Markers/raw.data.KO.2") write.csv(raw.data.KO.3, "CountResult/Markers/raw.data.KO.3") write.csv(raw.data.KO.4, "CountResult/Markers/raw.data.KO.4") write.csv(raw.data.KO.5, "CountResult/Markers/raw.data.KO.5") f1 <- FeaturePlot(KO.cells, features = c('ncRNA-inter-chr7-5998'), reduction = "umap", split.by = "stim") plot_grid(f1,p1,p4) DefaultAssay(KO.cells) <- "RNA" KO.cells <- NormalizeData(KO.cells, verbose = FALSE) plots <- VlnPlot(KO.cells, features = c("Alb", "ncRNA-inter-chr7-5998","Cyp2b10","Cyp2e1","Cyp2f2"), split.by = "stim", group.by = "seurat_clusters", pt.size = 0, combine = FALSE) CombinePlots(plots = plots, ncol = 1) Three_five_six <- subset(KO.cells, idents = c("5","6")) Three_five_six <- RunUMAP(Three_five_six, reduction = "pca", dims = 1:20 ) Three_five_six <- FindNeighbors(Three_five_six, reduction = "pca", dims = 1:20) Three_five_six <- FindClusters(Three_five_six, resolution = 1 ) Three_five_six <- RunTSNE(Three_five_six, reduction = "pca", dims = 1:20) p1 <- UMAPPlot(Three_five_six, reduction = "umap", group.by = "stim") p2 <- UMAPPlot(Three_five_six, reduction = "umap", group.by = "mouse.sex") p3 <- UMAPPlot(Three_five_six, reduction = "umap", label = TRUE) p4 <- UMAPPlot(Three_five_six, label=TRUE) plot_grid(p1,p4) DimPlot(Three_five_six, reduction = "umap", split.by = "stim") f1 <- FeaturePlot(Three_five_six, features = c('ncRNA-inter-chr7-5998'), reduction = "umap", split.by = "stim") lnc5998 <- subset(combined, cells = lnc5998.cells, idents = "1") lnc5998 <- RunUMAP(lnc5998, reduction = "pca", dims = 1:20 ) lnc5998 <- FindNeighbors(lnc5998, reduction = "pca", dims = 1:20) lnc5998 <- FindClusters(lnc5998, resolution = 1 ) lnc5998 <- RunTSNE(lnc5998, reduction = "pca", dims = 1:20) # Visualization p1 <- UMAPPlot(lnc5998, reduction = "umap", group.by = "stim") p2 <- UMAPPlot(lnc5998, reduction = "umap", group.by = "mouse.sex") p3 <- UMAPPlot(lnc5998, reduction = "umap", label = TRUE) p4 <- UMAPPlot(lnc5998, label=TRUE) plot_grid(p1,p4) DimPlot(lnc5998, reduction = "umap", split.by = "stim") lnc5998.cells <- WhichCells(object = combined, expression = "ncRNA-inter-chr7-5998" > 1) FeaturePlot(lnc5998, features = c("ncRNA-inter-chr7-5998"), split.by = "stim", + cols = c("grey", "red"), cells = lnc5998.cells,min.cutoff = 0.5) DefaultAssay(lnc5998) <- "RNA" lnc5998 <- NormalizeData(lnc5998, verbose = FALSE) plots <- VlnPlot(lnc5998, features = c("Alb", "ncRNA-inter-chr7-5998","Cyp2b10","Cyp2e1","Cyp2f2"), split.by = "stim", group.by = "seurat_clusters", pt.size = 0, combine = FALSE) CombinePlots(plots = plots, ncol = 1) ``` ```{r} d <- dist(t(GetAssayData(KO.cells, slot = "scale.data"))) # Run the MDS procedure, k determines the number of dimensions mds <- cmdscale(d = d, k = 2) # cmdscale returns the cell embeddings, we first label the columns to ensure downstream # consistency colnames(mds) <- paste0("MDS_", 1:2) # We will now store this as a custom dimensional reduction called 'mds' KO.cells[["mds"]] <- CreateDimReducObject(embeddings = mds, key = "MDS_", assay = DefaultAssay(KO.cells)) # We can now use this as you would any other dimensional reduction in all downstream functions DimPlot(KO.cells, reduction = "mds", pt.size = 0.5) ``` Find differential markers ```{r} KO.cells$celltype.stim <- paste(Idents(KO.cells), KO.cells$stim, sep = "_") KO.cells$celltype <- Idents(KO.cells) Idents(KO.cells) <- "celltype.stim" response3 <- FindMarkers(KO.cells, ident.1 = c("1_G171B","0_G171B","2_G171B"), ident.2 = c("1_G171C", "0_G171C","2_G171C"), verbose = TRUE, test.use = "MAST", logfc.threshold = FALSE,min.pct = FALSE) head(response3, n = 15) #DE1: Compare cluster 0+1+12 (that expressed lnc5998) with Other hepatocyte clusters (2+5+8) #DE2: compare cluster 1 (showed major effects in the KD) vs Cluster 0 (that showed little KD) #DE3: For KO.cells that formed five subcluster, compare clusters 3+4+5+1 vs 2+0 #DE4: for KO.cells that formed five clusters. Comapre cluster 4 vs. 2 DE0112.258 <- FindMarkers(combined, ident.1 = c("0","1","12" ), ident.2 = c("2","5","8"),verbose = TRUE, logfc.threshold = FALSE,min.pct = FALSE) DE0112.All <- FindMarkers(combined, ident.1 = c("0","1","12" ), verbose = TRUE, logfc.threshold = FALSE,min.pct = FALSE) DE1.2 <- FindMarkers(combined, ident.1 = c("1" ), ident.2 = c("2"),verbose = TRUE, logfc.threshold = FALSE,min.pct = FALSE) DE1.5 <- FindMarkers(combined, ident.1 = c("1" ), ident.2 = c("5"),verbose = TRUE, logfc.threshold = FALSE,min.pct = FALSE) DE1.8 <- FindMarkers(combined, ident.1 = c("1" ), ident.2 = c("8"),verbose = TRUE, logfc.threshold = FALSE,min.pct = FALSE) DE2.5 <- FindMarkers(combined, ident.1 = c("2" ), ident.2 = c("5"),verbose = TRUE, logfc.threshold = FALSE,min.pct = FALSE) DE2.8 <- FindMarkers(combined, ident.1 = c("2" ), ident.2 = c("8"),verbose = TRUE, logfc.threshold = FALSE,min.pct = FALSE) DE5.8 <- FindMarkers(combined, ident.1 = c("5" ), ident.2 = c("8"),verbose = TRUE, logfc.threshold = FALSE,min.pct = FALSE) DEK4351.20 <- FindMarkers(KO.cells, ident.1 = c("3","4","5","1" ), ident.2 = c("2","0"),verbose = TRUE, logfc.threshold = FALSE,min.pct = FALSE) DEK4.3 <- FindMarkers(KO.cells, ident.1 = "4", ident.2 = "3",verbose = TRUE, logfc.threshold = FALSE,min.pct = FALSE) DEK4.2 <- FindMarkers(KO.cells, ident.1 = "4", ident.2 = "2",verbose = TRUE, logfc.threshold = FALSE,min.pct = FALSE) DEK4.0 <- FindMarkers(KO.cells, ident.1 = "4", ident.2 = "0",verbose = TRUE, logfc.threshold = FALSE,min.pct = FALSE) DEK4.1 <- FindMarkers(KO.cells, ident.1 = "4", ident.2 = "1",verbose = TRUE, logfc.threshold = FALSE,min.pct = FALSE) DEK4.5 <- FindMarkers(KO.cells, ident.1 = "4", ident.2 = "5",verbose = TRUE, logfc.threshold = FALSE,min.pct = FALSE) ##### comaprison between G171B vs G171C for KO.cell clusters of 0+1+12 ######33 DEK4_C.B <- FindMarkers(KO.cells, ident.1 = "4_G171C", ident.2 = "4_G171B",verbose = TRUE, logfc.threshold = FALSE,min.pct = FALSE) DEK3_C.B <- FindMarkers(KO.cells, ident.1 = "3_G171C", ident.2 = "3_G171B",verbose = TRUE, logfc.threshold = FALSE,min.pct = FALSE) DEK1_C.B <- FindMarkers(KO.cells, ident.1 = "1_G171C", ident.2 = "1_G171B",verbose = TRUE, logfc.threshold = FALSE,min.pct = FALSE) DEK2_C.B <- FindMarkers(KO.cells, ident.1 = "2_G171C", ident.2 = "2_G171B",verbose = TRUE, logfc.threshold = FALSE,min.pct = FALSE) DEK0_C.B <- FindMarkers(KO.cells, ident.1 = "0_G171C", ident.2 = "0_G171B",verbose = TRUE, logfc.threshold = FALSE,min.pct = FALSE) DEK5_C.B <- FindMarkers(KO.cells, ident.1 = "5_G171C", ident.2 = "5_G171B",verbose = TRUE, logfc.threshold = FALSE,min.pct = FALSE) ################################### Write the results ######################### write.csv(DE0112.258, "CountResult/Markers/DE0112.258") write.csv(DE0112.All, "CountResult/Markers/DE0112.All") write.csv(DE1.2, "CountResult/Markers/DE1.2") write.csv(DE1.5, "CountResult/Markers/DE1.5") write.csv(DE1.8, "CountResult/Markers/DE1.8") write.csv(DE2.5, "CountResult/Markers/DE2.5") write.csv(DE2.8, "CountResult/Markers/DE2.8") write.csv(DE5.8, "CountResult/Markers/DE5.8") write.csv(DEK4351.20, "CountResult/Markers/DEK4351.20") write.csv(DEK4.3, "CountResult/Markers/DEK4.3") write.csv(DEK4.2, "CountResult/Markers/DEK4.2") write.csv(DEK4.0, "CountResult/Markers/DEK4.0") write.csv(DEK4.1, "CountResult/Markers/DEK4.1") write.csv(DEK4.5, "CountResult/Markers/DEK4.5") write.csv(DEK4_C.B, "CountResult/Markers/DEK4_C.B") write.csv(DEK3_C.B, "CountResult/Markers/DEK3_C.B") write.csv(DEK1_C.B, "CountResult/Markers/DEK1_C.B") write.csv(DEK2_C.B, "CountResult/Markers/DEK2_C.B") write.csv(DEK0_C.B, "CountResult/Markers/DEK0_C.B") write.csv(DEK5_C.B, "CountResult/Markers/DEK5_C.B") combined$celltype.stim <- paste(Idents(combined), combined$stim, sep = "_") combined$celltype <- Idents(combined) Idents(combined) <- "celltype.stim" test.combined$celltype.stim <- paste(Idents(test.combined), test.combined$stim, sep = "_") test.combined$celltype <- Idents(test.combined) Idents(test.combined) <- "celltype.stim" Combined_G171B_vs_G171C <- FindMarkers(combined, ident.1 = c("0_G171B","1_G171B","2_G171B","3_G171B","4_G171B","5_G171B","6_G171B","7_G171B","8_G171B","9_G171B","10_G171B","11_G171B","12_G171B" ), ident.2 = c("0_G171C","1_G171C","2_G171C","3_G171C","4_G171C","5_G171C","6_G171C","7_G171C","8_G171C","9_G171C","10_G171C","11_G171C","12_G171C" ), verbose = TRUE, logfc.threshold = FALSE,min.pct = FALSE) Combined_G171C_vs_G171B_Clust1 <- FindMarkers(combined, ident.1 = c("1_G171C"), ident.2 = c("1_G171B"), verbose = TRUE, logfc.threshold = FALSE,min.pct = FALSE) Combined_G171C_vs_G171B_Clust0 <- FindMarkers(combined, ident.1 = c("0_G171C"), ident.2 = c("0_G171B"), verbose = TRUE, logfc.threshold = FALSE,min.pct = FALSE) Combined_G171C_vs_G171B_Clust12 <- FindMarkers(combined, ident.1 = c("12_G171C"), ident.2 = c("12_G171B"), verbose = TRUE, logfc.threshold = FALSE,min.pct = FALSE) Combined_G171C_vs_G171B_Clust2 <- FindMarkers(combined, ident.1 = c("2_G171C"), ident.2 = c("2_G171B"), verbose = TRUE, logfc.threshold = FALSE,min.pct = FALSE) Combined_G171C_vs_G171B_Clust5 <- FindMarkers(combined, ident.1 = c("5_G171C"), ident.2 = c("5_G171B"), verbose = TRUE, logfc.threshold = FALSE,min.pct = FALSE) Combined_G171C_vs_G171B_Clust8 <- FindMarkers(combined, ident.1 = c("8_G171C"), ident.2 = c("8_G171B"), verbose = TRUE, logfc.threshold = FALSE,min.pct = FALSE) Combined_G171C_vs_G171B_Clust3 <- FindMarkers(combined, ident.1 = c("3_G171C"), ident.2 = c("3_G171B"), verbose = TRUE, logfc.threshold = FALSE,min.pct = FALSE) Combined_G171C_vs_G171B_Clust4 <- FindMarkers(combined, ident.1 = c("4_G171C"), ident.2 = c("4_G171B"), verbose = TRUE, logfc.threshold = FALSE,min.pct = FALSE) Combined_G171C_vs_G171B_Clust6 <- FindMarkers(combined, ident.1 = c("6_G171C"), ident.2 = c("6_G171B"), verbose = TRUE, logfc.threshold = FALSE,min.pct = FALSE) Combined_G171C_vs_G171B_Clust7 <- FindMarkers(combined, ident.1 = c("7_G171C"), ident.2 = c("7_G171B"), verbose = TRUE, logfc.threshold = FALSE,min.pct = FALSE) Combined_G171C_vs_G171B_Clust9 <- FindMarkers(combined, ident.1 = c("9_G171C"), ident.2 = c("9_G171B"), verbose = TRUE, logfc.threshold = FALSE,min.pct = FALSE) Combined_G171C_vs_G171B_Clust10 <- FindMarkers(combined, ident.1 = c("10_G171C"), ident.2 = c("10_G171B"), verbose = TRUE, logfc.threshold = FALSE,min.pct = FALSE) Combined_G171C_vs_G171B_Clust11 <- FindMarkers(combined, ident.1 = c("11_G171C"), ident.2 = c("11_G171B"), verbose = TRUE, logfc.threshold = FALSE,min.pct = FALSE) write.csv(Combined_G171C_vs_G171B_Clust1, "CountResult/Markers/Combined_G171C_vs_G171B_Clust1") write.csv(Combined_G171C_vs_G171B_Clust0, "CountResult/Markers/Combined_G171C_vs_G171B_Clust0") write.csv(Combined_G171C_vs_G171B_Clust12, "CountResult/Markers/Combined_G171C_vs_G171B_Clust12") write.csv(Combined_G171C_vs_G171B_Clust2, "CountResult/Markers/Combined_G171C_vs_G171B_Clust2") write.csv(Combined_G171C_vs_G171B_Clust5, "CountResult/Markers/Combined_G171C_vs_G171B_Clust5") write.csv(Combined_G171C_vs_G171B_Clust8, "CountResult/Markers/Combined_G171C_vs_G171B_Clust8") write.csv(Combined_G171C_vs_G171B_Clust3, "CountResult/Markers/Combined_G171C_vs_G171B_Clust3") write.csv(Combined_G171C_vs_G171B_Clust4, "CountResult/Markers/Combined_G171C_vs_G171B_Clust4") write.csv(Combined_G171C_vs_G171B_Clust6, "CountResult/Markers/Combined_G171C_vs_G171B_Clust6") write.csv(Combined_G171C_vs_G171B_Clust7, "CountResult/Markers/Combined_G171C_vs_G171B_Clust7") write.csv(Combined_G171C_vs_G171B_Clust9, "CountResult/Markers/Combined_G171C_vs_G171B_Clust9") write.csv(Combined_G171C_vs_G171B_Clust10, "CountResult/Markers/Combined_G171C_vs_G171B_Clust10") write.csv(Combined_G171C_vs_G171B_Clust11, "CountResult/Markers/Combined_G171C_vs_G171B_Clust11") PC_vs_PP_G171B <- FindMarkers(KO.cells, ident.1 = c("4_G171B","3_G171B"), ident.2 = c("0_G171B","2_G171B"), verbose = TRUE, test.use = "MAST", logfc.threshold = FALSE,min.pct = FALSE) head(PC_vs_PP_G171B, n = 15) PC_vs_PP_G171C <- FindMarkers(KO.cells, ident.1 = c("4_G171C","3_G171C"), ident.2 = c("0_G171C","2_G171C"), verbose = TRUE, test.use = "MAST", logfc.threshold = FALSE,min.pct = FALSE) head(PC_vs_PP_G171C, n = 15) PC_vs_PP_G171BC <- FindMarkers(KO.cells, ident.1 = c("4_G171B","3_G171B","4_G171C","3_G171C"), ident.2 = c("0_G171B","2_G171B","0_G171C","2_G171C"), verbose = TRUE, test.use = "MAST", logfc.threshold = FALSE,min.pct = FALSE) head(PC_vs_PP_G171C, n = 15) lnc5998_KO_DE_1 <- FindMarkers(KO.cells, ident.1 = c("4_G171B","3_G171B","5_G171B"), ident.2 = c("4_G171C","3_G171C","5_G171C"), verbose = TRUE, test.use = "MAST", logfc.threshold = FALSE,min.pct = FALSE) head(lnc5998_KO_DE, n = 15) lnc5998_KO_DE_2 <- FindMarkers(KO.cells, ident.1 = c("4_G171C","3_G171C","5_G171C"), ident.2 =c("4_G171B","3_G171B","5_G171B") , verbose = TRUE, test.use = "MAST", logfc.threshold = FALSE,min.pct = FALSE) head(lnc5998_KO_DE, n = 15) cell.type.genes <- (PC_vs_PP_G171BC[1]) # Takes all the unique cell type specific genes GOterms = topGOterms(fg.genes = cell.type.genes, bg.genes = rownames(KO.cells@assays$RNA@dataKO.cells@assays$RNA@data), organism = "Mouse") cell.type.genes <- (PC_vs_PP_G171BC[1]) # Takes all the unique cell type specific genes GOterms = topGOterms(fg.genes = rownames(cell.type.genes), bg.genes = rownames(KO.cells@assays$RNA@data), organism = "Mouse") ``` ```{r} AvergeExpression2 <- function (object, assays = NULL, features = NULL, return.seurat = FALSE, add.ident = NULL, slot = "data", use.scale = FALSE, use.counts = FALSE, verbose = TRUE, ...) { fxn.average <- switch(EXPR = slot, data = function(x) { return(mean(x = x)) }, mean) object.assays <- FilterObjects(object = object, classes.keep = "Assay") assays <- assays %||% object.assays ident.orig <- Idents(object = object) orig.levels <- levels(x = Idents(object = object)) ident.new <- c() if (!all(assays %in% object.assays)) { assays <- assays[assays %in% object.assays] if (length(assays) == 0) { stop("None of the requested assays are present in the object") } else { warning("Requested assays that do not exist in object. Proceeding with existing assays only.") } } if (!is.null(x = add.ident)) { new.data <- FetchData(object = object, vars = add.ident) new.ident <- paste(Idents(object)[rownames(x = new.data)], new.data[, 1], sep = "_") Idents(object, cells = rownames(new.data)) <- new.ident } data.return <- list() for (i in 1:length(x = assays)) { data.use <- GetAssayData(object = object, assay = assays[i], slot = slot) features.assay <- features if (length(x = intersect(x = features, y = rownames(x = data.use))) < 1) { features.assay <- rownames(x = data.use) } data.all <- data.frame(row.names = features.assay) for (j in levels(x = Idents(object))) { temp.cells <- WhichCells(object = object, idents = j) features.assay <- unique(x = intersect(x = features.assay, y = rownames(x = data.use))) if (length(x = temp.cells) == 1) { data.temp <- (data.use[features.assay, temp.cells]) if (slot == "data") { data.temp <- data.temp } } if (length(x = temp.cells) > 1) { data.temp <- apply(X = data.use[features.assay, temp.cells, drop = FALSE], MARGIN = 1, FUN = fxn.average) } data.all <- cbind(data.all, data.temp) colnames(x = data.all)[ncol(x = data.all)] <- j if (verbose) { message(paste("Finished averaging", assays[i], "for cluster", j)) } if (i == 1) { ident.new <- c(ident.new, as.character(x = ident.orig[temp.cells[1]])) } } names(x = ident.new) <- levels(x = Idents(object)) data.return[[i]] <- data.all names(x = data.return)[i] <- assays[[i]] } if (return.seurat) { toRet <- CreateSeuratObject(counts = data.return[[1]], project = "Average", assay = names(x = data.return)[1], ...) if (length(x = data.return) > 1) { for (i in 2:length(x = data.return)) { toRet[[names(x = data.return)[i]]] <- CreateAssayObject(counts = data.return[[i]]) } } if (DefaultAssay(object = object) %in% names(x = data.return)) { DefaultAssay(object = toRet) <- DefaultAssay(object = object) } Idents(toRet, cells = colnames(x = toRet)) <- ident.new[colnames(x = toRet)] Idents(object = toRet) <- factor(x = Idents(object = toRet), levels = as.character(x = orig.levels), ordered = TRUE) toRet <- NormalizeData(object = toRet, verbose = verbose) toRet <- ScaleData(object = toRet, verbose = verbose) return(toRet) } else { return(data.return) } } ``` Find differential expression markers ```{r} combined.markers <- FindAllMarkers(object = combined, only.pos = TRUE, min.pct = 0.25, thresh.use = 0.25) combined.markers %>% group_by(cluster) %>% top_n(2, avg_logFC) KO.markers <- FindAllMarkers(object = KO.cells, only.pos = TRUE, min.pct = 0.25, thresh.use = 0.25) response3 <- FindMarkers(combined, ident.1 = c("1_G171B","0_G171B","2_G171B"), ident.2 = c("1_G171C", "0_G171C","2_G171C"), verbose = TRUE, test.use = "MAST", logfc.threshold = FALSE,min.pct = FALSE) head(response3, n = 15) ``` ```{r} genes_hep_main =c('Alb', 'Ttr', 'Apoa1', 'Serpina1c') genes_hep = c('Alb', 'Ttr', 'Apoa1', 'Serpina1c', 'Cyp2e1', 'Glul', 'Oat', 'Gulo', 'Ass1', 'Hamp', 'Gstp1', 'Ubb', 'Cyp2f2', 'Pck1', 'Hal', 'Cdh1') genes_endo = c('Pecam1', 'Nrp1', 'Kdr','Oit3') genes_kuppfer = c( 'Clec4f', 'Cd68') genes_nk = c('Il2rb', 'Nkg7', 'Cxcr6', 'Gzma') genes_b = c('Cd79a', 'Cd79b') genes_bec = c('Epcam', 'Krt19', 'Krt7') genes_immune = 'Ptprc' HSC = c("Acta2", "Gfap","Des","Ngfr","Col2a1","Vim","Lama1","Nes") all_genes = c(genes_hep, genes_endo, genes_kuppfer, genes_nk, genes_b, genes_bec, genes_immune, HSC) genes_bec_b_immune = c(genes_bec,genes_b,genes_immune) genes_zones = c('Cyp2e1', 'Glul', 'Oat', 'Gulo', 'Ass1', 'Hamp', 'Gstp1', 'Ubb', 'Cyp2f2', 'Pck1', 'Hal', 'Cdh1') receptor_KO <- c("ncRNA-inter-chr7-5998","Cyp2b10","Nr1i2","Nr1i3","Ppara","Pparg","Ppargc1b","Ppard") ``` Visualize top genes in principal components ```{r, echo=FALSE, fig.height=4, fig.width=8} PCHeatmap(object = tiss1, pc.use = 1:3, cells.use = 500, do.balanced = TRUE, label.columns = FALSE, num.genes = 8) ``` Later on (in FindClusters and TSNE) you will pick a number of principal components to use. This has the effect of keeping the major directions of variation in the data and, ideally, supressing noise. There is no correct answer to the number to use, but a decent rule of thumb is to go until the plot plateaus. ```{r} PCElbowPlot(object = tiss1) ``` Choose the number of principal components to use. ```{r} # Set number of principal components. n.pcs = 10 ``` The clustering is performed based on a nearest neighbors graph. Cells that have similar expression will be joined together. The Louvain algorithm looks for groups of cells with high modularity--more connections within the group than between groups. The resolution parameter determines the scale. Higher resolution will give more clusters, lower resolution will give fewer. For the top-level clustering, aim to under-cluster instead of over-cluster. It will be easy to subset groups and further analyze them below. ```{r} # Set resolution res.used <- 4 tiss1 <- FindClusters(object = tiss1, reduction.type = "pca", dims.use = 1:n.pcs, resolution = res.used, print.output = 0, save.SNN = TRUE, force.recalc = TRUE) ``` We use TSNE solely to visualize the data. ```{r} # If cells are too spread out, you can raise the perplexity. If you have few cells, try a lower perplexity (but never less than 10). tiss1 <- RunTSNE(object = tiss1, dims.use = 1:n.pcs, seed.use = 10, perplexity=30) ``` ```{r} TSNEPlot(object = tiss1, do.label = T, pt.size = 1.2, label.size = 4) ``` ## Compare to previous annotations ```{r} previous_annotation = read.csv("/Users/kkarri/Documents/Lab/Single_cell_project/dropseq/Liver_droplet_annotation.csv", stringsAsFactors = FALSE) cols = c('free_annotation', 'cell_ontology_class') for (col in cols){ previous_col = paste0('previous_', col) tiss1@meta.data[, previous_col] <- "NA" tiss1@meta.data[as.character(previous_annotation$X), previous_col] <- previous_annotation[, col] print(table(tiss1@meta.data[, previous_col])) print(table(tiss1@meta.data[, previous_col], tiss@ident)) } tiss1 = compare_previous_annotation(tiss1, tissue_of_interest, "droplet") TSNEPlot(object = tiss1, do.return = TRUE, group.by = "previous_cell_ontology_class") table(tiss1@meta.data[, "previous_cell_ontology_class"], tiss@ident) ``` ```{r} tiss1 = compare_previous_annotation(tiss1, tissue_of_interest, "droplet") TSNEPlot(object = tiss1, do.return = TRUE, group.by = "previous_cell_ontology_class") table(tiss1@meta.data[, "previous_cell_ontology_class"], tiss1@ident) ``` ```{r} TSNEPlot(tiss1, group.by="mouse.sex") TSNEPlot(tiss1, group.by="mouse.id") ``` Significant genes: hepatocyte: Alb, Ttr, Apoa1, and Serpina1c pericentral: Cyp2e1, Glul, Oat, Gulo midlobular: Ass1, Hamp, Gstp1, Ubb periportal: Cyp2f2, Pck1, Hal, Cdh1 endothelial cells: Pecam1, Nrp1, Kdr+ and Oit3+ Kuppfer cells: Emr1, Clec4f, Cd68, Irf7 NK/NKT cells: Zap70, Il2rb, Nkg7, Cxcr6, Klr1c, Gzma B cells: Cd79a, Cd79b, Cd74 and Cd19 Immune cells: Ptprc ```{r, echo=FALSE, fig.height=16, fig.width=12} # Hepatic marker FeaturePlot(tiss1, c(genes_hep), pt.size = 1, nCol = 4, cols.use = c("grey", "red")) # Endothelial markers FeaturePlot(tiss1, c(genes_endo), pt.size = 1, nCol = 4, cols.use = c("grey", "red")) # Kupffer cells FeaturePlot(tiss1, c(genes_kuppfer), pt.size = 1, nCol = 4, cols.use = c("grey", "red")) # genes_nk FeaturePlot(tiss1, c(genes_nk), pt.size = 1, nCol = 4, cols.use = c("grey", "red")) # genes_b FeaturePlot(tiss1, c(genes_b), pt.size = 1, nCol = 4, cols.use = c("grey", "red")) # genes bile duct endo cells FeaturePlot(tiss1, c(genes_bec), pt.size = 1, nCol = 4, cols.use = c("grey", "red")) # genes immune FeaturePlot(tiss1, c(genes_immune), pt.size = 1, nCol = 4, cols.use = c("grey", "red")) ``` Dotplots let you see the intensity of exppression and the fraction of cells expressing for each of your genes of interest. The radius shows you the percent of cells in that cluster with at least one read sequenced from that gene. The color level indicates the average Z-score of gene expression for cells in that cluster, where the scaling is done over taken over all cells in the sample. #We have various immune cell types in the last cluster ```{r, echo=FALSE, fig.height=4, fig.width=10} DotPlot(tiss1, c(genes_kuppfer, genes_nk, genes_b, "Ptprc"), plot.legend = T, col.max = 2.5, do.return = T) + coord_flip() ``` ```{r, echo=FALSE, fig.height=8, fig.width=10} DotPlot(tiss1, c(genes_hep_main, genes_endo, genes_nk, genes_kuppfer, genes_bec_b_immune), plot.legend = T, col.max = 2.5, do.return = T) + coord_flip() ``` Using the markers above, we can confidentaly label many of the clusters: 19: endothelial cells 20: bile duct epithelial cells 21: immune cells rest are hepatocytes We will add those cell_ontology_classes to the dataset. ```{r} tiss1 <- StashIdent(object = tiss1, save.name = "cluster.ids") cluster.ids <- c(0,1,2,3,4,5,6,7,8,9,10,11,12,13,14,15,16,17,18,19,20) free_annotation <- c( NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, "bile duct epithelial cells", "endothelial cell of hepatic sinusoid", NA ) cell_ontology_class <- c( "hepatocyte", "hepatocyte", "hepatocyte", "hepatocyte", "hepatocyte", "hepatocyte", "hepatocyte", "hepatocyte", "hepatocyte", "hepatocyte", "hepatocyte", "hepatocyte", "hepatocyte", "hepatocyte", "hepatocyte", "hepatocyte", "hepatocyte", "hepatocyte", "duct epithelial cell", "endothelial cell of hepatic sinusoid", "hepatocyte") tiss1 = stash_annotations(tiss1, cluster.ids, free_annotation, cell_ontology_class) ``` ## Checking for batch effects Color by metadata, like plate barcode, to check for batch effects. ```{r} TSNEPlot(object = tiss1, do.return = TRUE, group.by = "channel") TSNEPlot(object = tiss1, do.return = TRUE, group.by = "free_annotation") ``` ## Subcluster Let's drill down on the hepatocytes. ```{r} subtiss1 = SubsetData(tiss1, ident.use = c(0,1,2,3,4,5,6,7,8,9,10,11,12,13,14,15,16,17,20)) ``` ```{r} subtiss1 <- subtiss1 %>% ScaleData() %>% FindVariableGenes(do.plot = FALSE, x.high.cutoff = Inf, y.cutoff = 0.5) %>% RunPCA(do.print = FALSE) ``` ```{r} PCHeatmap(object = subtiss1, pc.use = 1:3, cells.use = 20, do.balanced = TRUE, label.columns = FALSE, num.genes = 8) PCElbowPlot(subtiss1) ``` ```{r} sub.n.pcs = 8 sub.res.use = 0.5 subtiss1 <- subtiss1 %>% FindClusters(reduction.type = "pca", dims.use = 1:sub.n.pcs, resolution = sub.res.use, print.output = 0, save.SNN = TRUE, force.recalc = TRUE) %>% RunTSNE(dims.use = 1:sub.n.pcs, seed.use = 10, perplexity=8) TSNEPlot(object = subtiss1, do.label = T, pt.size = 1, label.size = 4) ``` ```{r, echo=FALSE, fig.height=25, fig.width=25} FeaturePlot(subtiss1, genes_hep,cols.use = c("grey", "red"), pt.size = 4, nCol = 4) ``` ```{r, echo=FALSE, fig.height=8, fig.width=10} DotPlot(subtiss1, all_genes, col.max = 2.5, plot.legend = T, do.return = T) + coord_flip() ``` ```{r} BuildClusterTree(subtiss1) ``` ```{r, echo=FALSE, fig.height=10, fig.width=8} #female genes have lower expression in cluster 6 relative to other female clusters, especally Xist FeaturePlot(subtiss1,c('Mup20', 'Mup1','Mup12', 'Mup21', 'Cyp2d9', 'Xist', 'A1bg', 'Cyp2c69'),cols.use = c("grey", "red"), pt.size = 3, nCol = 2) DotPlot(tiss1,c('Mup20', 'Mup1','Mup12', 'Mup21', 'Cyp2d9', 'Xist', 'A1bg', 'Cyp2c69'), plot.legend = T, col.max = 2.5, do.return = T) + coord_flip() ``` From these genes, it appears that the clusters represent: 0: midlobular male 1: pericentral female 2: periportal female 3: periportal male 4: midlobular male 5: pericentral male 6: midlobular female 7: midlobular female The multitude of clusters of each type correspond mostly to individual animals/sexes. ```{r} table(FetchData(subtiss1, c('mouse.sex','ident')) %>% droplevels()) ``` ```{r} sub.cluster.ids <- c(0, 1, 2, 3, 4, 5, 6, 7) sub.free_annotation <- c("periportal female", "midlobular male", "pericentral female", "periportal male", "midlobular male", "pericentral male", "midlobular female", "midlobular female") sub.cell_ontology_class <- c("hepatocyte", "hepatocyte", "hepatocyte", "hepatocyte", "hepatocyte", "hepatocyte", "hepatocyte", "hepatocyte") subtiss1 = stash_annotations(subtiss1, sub.cluster.ids, sub.free_annotation, sub.cell_ontology_class) tiss1 = stash_subtiss_in_tiss(tiss1, subtiss1) ``` Liver zonation markers ```{r} genes_zones = c('Cyp2e1', 'Glul', 'Oat', 'Gulo', 'Ass1', 'Hamp', 'Gstp1', 'Ubb', 'Cyp2f2', 'Pck1', 'Hal', 'Cdh1') FeaturePlot(subtiss1,c(genes_zones),cols.use = c("grey", "red"), pt.size = 1, nCol = 4) DotPlot(subtiss1,c(genes_zones), plot.legend = T, col.max = 2.5, do.return = T) + coord_flip() TSNEPlot(object = subtiss1, do.label = T, pt.size = 1, label.size = 4, group.by="free_annotation") TSNEPlot(object = tiss1, do.label = T, pt.size = 1, label.size = 4, group.by="free_annotation") ``` ########## Find cluster markers for lncRNAs ```{r} MIN_LOGFOLD_CHANGE = 1 # set to minimum required average log fold change in gene expression. MIN_PCT_CELLS_EXPR_GENE = 0.1 all.markers = FindAllMarkers(tiss1, min.pct = MIN_PCT_CELLS_EXPR_GENE, logfc.threshold = MIN_LOGFOLD_CHANGE, only.pos = TRUE, test.use="bimod") # likelihood ratio test lnc_all_markers <- grep(pattern = "^ncRNA", x= rownames(all.markers), value = TRUE) lnc_all_markers #[1] "ncRNA_inter_chr10_92081" "ncRNA_intra_chr16_13383" "ncRNA_inter_chr17_13605" "ncRNA_inter_chr14_11815" #[5] "ncRNA_inter_chr18_14344" FeaturePlot(subtiss1,c(lnc_all_markers),cols.use = c("grey", "red"), pt.size = 1, nCol = 4) ######################### lncRNA markers- CELL TYPE MARKER ############ markers.hep <- FindMarkers(object = tiss1, ident.1 = c(0,1,2,3,4,5,6,7,8,9,10,11,12,13,14,15,16,17,20), ident.2 = c(18,19),only.pos = TRUE, min.pct = 0.25, thresh.use = 0.25) lnc_markers_hep <- grep(pattern = "^ncRNA", x= rownames(markers.hep), value = TRUE) lnc_markers_hep FeaturePlot(tiss1,c(lnc_markers_hep),cols.use = c("grey", "red"), pt.size = 1, nCol = 4) DotPlot(tiss1,lnc_markers_hep, plot.legend = T, col.max = 2.5, do.return = T) + coord_flip() #[1] "ncRNA_as_chr11_9423" "ncRNA_as_chr7_6166" "ncRNA_inter_chr4_3295" "ncRNA_inter_chr17_14026" #[5] "ncRNA_inter_chr3_2915" "ncRNA_inter_chr5_4547" "ncRNA_inter_chr15_12684" markers.hep.MAST <- FindMarkers(object = tiss1, ident.1 = c(0,1,2,3,4,5,6,7,8,9,10,11,12,13,14,15,16,17,20), ident.2 = c(18,19),only.pos = TRUE, test.use = "MAST") lnc_markers_hep_MAST_TABLE <- subset(markers.hep.MAST, grepl("^ncRNA", rownames(markers.hep.MAST))) lnc_markers_hep_MAST <- grep(pattern = "^ncRNA", x= rownames(markers.hep.MAST), value = TRUE) lnc_markers_hep_MAST markers.endo <- FindMarkers(object = tiss1, ident.1 = c(18,19), only.pos = TRUE, min.pct = 0.25, thresh.use = 0.5) lnc_markers_endo <- grep(pattern = "^ncRNA", x= rownames(markers.endo), value = TRUE) lnc_markers_endo FeaturePlot(tiss1,c(lnc_markers_endo),cols.use = c("grey", "red"), pt.size = 1, nCol = 4) DotPlot(tiss1,lnc_markers_endo, plot.legend = T, col.max = 2.5, do.return = T) + coord_flip() #"ncRNA_inter_chr15_12770", "ncRNA_inter_chr12_10817", "ncRNA_as_chr13_11451", markers.endo.MAST <- FindMarkers(object = tiss1, ident.1 = 19, test.use = "MAST" ,only.pos = TRUE) lnc_markers_endo_MAST_TABLE <- subset(markers.endo.MAST, grepl("^ncRNA", rownames(markers.endo.MAST))) lnc_markers_endo_MAST <- grep(pattern = "^ncRNA", x= rownames(markers.endo.MAST), value = TRUE) lnc_markers_endo_MAST ################## lncRNA expression ########################3 # "ncRNA_inter_chr17_13605" , "ncRNA_intra_chr16_13383" ########## Periporal markers- zonation markers ############ markers.pc <- FindMarkers(object = subtiss1, ident.1 = c(2,5), only.pos = FALSE, min.pct = 0.001, thresh.use = 0.001, test.use = "bimod" ) markers.pc.MAST <- FindMarkers(object = subtiss1, ident.1 = c(2,5), ident.2 = c(0,3), test.use = "MAST" ,only.pos = TRUE) lnc_markers_pc <- subset(markers.pc, grepl("^ncRNA", rownames(markers.pc))) lnc_markers_pc <- grep(pattern = "^ncRNA", x= rownames(markers.pc), value = TRUE) lnc_markers_pc markers.pc.MAST <- FindMarkers(object = subtiss1, ident.1 = c(2,5), ident.2 = c(0,3), test.use = "MAST" ,only.pos = TRUE) lnc_markers_pc_MAST <- subset(markers.pc.MAST, grepl("^ncRNA", rownames(markers.pc.MAST))) lnc_markers_pc_MAST <- grep(pattern = "^ncRNA", x= rownames(markers.pc.MAST), value = TRUE) lnc_markers_pc_MAST DotPlot(tiss1, lnc_markers_pc, plot.legend = T, col.max = 2.5, do.return = T, group.by="free_annotation") + coord_flip() FeaturePlot(subtiss1,c(lnc_markers_pc),cols.use = c("grey", "red"), pt.size = 1, nCol = 4) ############################### midlobular genes ############# markers.mid <- FindMarkers(object = subtiss1, ident.1 = c(1,4,6,7), only.pos = FALSE, min.pct = 0.001, thresh.use = 0.05) lnc_markers_mid <- subset(markers.mid, grepl("^ncRNA", rownames(markers.mid))) lnc_markers_mid <- grep(pattern = "^ncRNA", x= rownames(markers.mid), value = TRUE) lnc_markers_mid DotPlot(tiss1, lnc_markers_mid, plot.legend = T, col.max = 2.5, do.return = T, group.by="free_annotation") + coord_flip() FeaturePlot(subtiss1,c(lnc_markers_mid),cols.use = c("grey", "red"), pt.size = 1, nCol = 4) markers.mid.MAST <- FindMarkers(object = subtiss1, ident.1 = c(1,4,6,7),test.use = "MAST",only.pos = TRUE ) lnc_markers_mid_MAST_TABLE <- subset(markers.mid.MAST, grepl("^ncRNA", rownames(markers.mid.MAST))) lnc_markers_mid_MAST <- grep(pattern = "^ncRNA", x= rownames(markers.mid.MAST), value = TRUE) lnc_markers_mid_MAST #####3 periportalmarker genes############3 markers.pp <- FindMarkers(object = subtiss1, ident.1 = c(0,3), only.pos = FALSE, min.pct = 0.001, thresh.use = 0.05) lnc_markers_pp <- subset(markers.pp, grepl("^ncRNA", rownames(markers.pp))) lnc_markers_pp <- grep(pattern = "^ncRNA", x= rownames(markers.pp), value = TRUE) lnc_markers_pp markers.pp.MAST <- FindMarkers(object = subtiss1, ident.1 = c(0,3), ident.2 = c(2,5),test.use = "MAST",only.pos = TRUE ) lnc_markers_pp_MAST_TABLE <- subset(markers.pp.MAST, grepl("^ncRNA", rownames(markers.pp.MAST))) lnc_markers_pp_MAST <- grep(pattern = "^ncRNA", x= rownames(markers.pp.MAST), value = TRUE) lnc_markers_pp_MAST FeaturePlot(subtiss1,c(lnc_markers_pp),cols.use = c("grey", "red"), pt.size = 1, nCol = 4, max.cutoff = 1) DotPlot(tiss1, c(lnc_markers_pp,"Cyp2e1","Cyp2f2"), plot.legend = T, col.max = 2.5, do.return = T, group.by= "free_annotation") + coord_flip() ################## amle and female specific ############################ markers.female <- FindMarkers(object = subtiss1, ident.1 = c(0,2,6,7), only.pos = TRUE, min.pct = 0.1, logfc.threshold = 1) lnc_markers_female <- subset(markers.female, grepl("^ncRNA", rownames(markers.female))) lnc_markers_female <- grep(pattern = "^ncRNA", x= rownames(markers.female), value = TRUE) lnc_markers_female FeaturePlot(subtiss1,c(lnc_markers_female),cols.use = c("grey", "red"), pt.size = 1, nCol = 4, max.cutoff = 1) DotPlot(tiss1, c(lnc_markers_female,"Cyp2e1","Cyp2f2"), plot.legend = T, col.max = 2.5, do.return = T, group.by= "free_annotation") + coord_flip() markers.male <- FindMarkers(object = subtiss1, ident.1 = c(1,3,4,5), only.pos = TRUE, min.pct = 0.001, thresh.use = 0.05) lnc_markers_male <- subset(markers.male, grepl("^ncRNA", rownames(markers.male))) lnc_markers_male <- grep(pattern = "^ncRNA", x= rownames(markers.male), value = TRUE) lnc_markers_male FeaturePlot(subtiss1,c(lnc_markers_male),cols.use = c("grey", "red"), pt.size = 1, nCol = 4, max.cutoff = 1) DotPlot(tiss1, c(lnc_markers_male,"Cyp2e1","Cyp2f2"), plot.legend = T, col.max = 2.5, do.return = T, group.by= "free_annotation") + coord_flip() ############################ Female zonate specific genes ################################### markers.pericentral.female <- FindMarkers(object = tiss1, ident.1 = c(6,11,14,20), test.use = "MAST", only.pos = TRUE, min.pct = 0.1, ident.2 = c(2,3,15,12,13,8,5,16), logfc.threshold = 1) markers.periportal.female <- FindMarkers(object = tiss1, ident.1 = c(2,3,15), only.pos = TRUE, min.pct = 0.1, ident.2 = c(6,11,14,20,12,13,8,5,16), logfc.threshold = 1) markers.pericentral.male <- FindMarkers(object = tiss1, ident.1 = c(13,12), test.use = "MAST", only.pos = TRUE, min.pct = 0.1, ident.2 = c(2,3,15,8,5,16,6,11,14,20), logfc.threshold = 1) markers.periportal.male <- FindMarkers(object = tiss1, ident.1 = c(8,5,16), test.use = "MAST", only.pos = TRUE, min.pct = 0.1, ident.2 = c(2,3,15,13,12,6,11,14,20), logfc.threshold = 1) ############### xeno-lncs CAR?RXR ################## FeaturePlot(tiss1,c("ncRNA_inter_chr15_12684","ncRNA_inter_chr8_7430","ncRNA_inter_chr7_6222"),cols.use = c("grey", "red"), pt.size = 1, nCol = 4, max.cutoff = 1) ####################################################################### markers.endo.2 <- FindMarkers(object = seurat_drop, logfc.threshold = 2,ident.1 = "Endothelial", only.pos = TRUE, min.pct = 0.25, thresh.use = 0.25) lnc.endo.2 <- grep(pattern = "^ncRNA", x= rownames(markers.endo.2), value = TRUE) lnc.endo.2 ``` Zonated lncRNAs ```{r} pp_zontaed <- c('ncRNA_inter_chr14_12016','ncRNA_as_chr19_15090','ncRNA_inter_chr10_9351','ncRNA_inter_chr16_13170', 'ncRNA_inter_chr3_2697','ncRNA_inter_chr1_274','ncRNA_as_chr6_5518','ncRNA_inter_chr14_12066','ncRNA_intra_chr12_10871', 'ncRNA_inter_chr16_13510','ncRNA_inter_chr3_2314','ncRNA_inter_chr10_9264,'ncRNA_inter_chr9_8122') ``` ## Checking for batch effects Color by metadata, like plate barcode, to check for batch effects. ```{r} TSNEPlot(object = subtiss1, do.return = TRUE, group.by = "mouse.sex") ``` # Final coloring Color by cell ontology class on the original TSNE. ```{r} TSNEPlot(object = tiss1, do.return = TRUE, group.by = "cell_ontology_class") ``` # Save the Robject for later ```{r} filename = here('00_data_ingest', '04_tiss1ue_robj_generated', paste0("droplet_", tiss1ue_of_interest, "refinedcells_seurat_tiss1.Robj")) print(filename) save(tiss1, file=filename) ``` ```{r} # To reload a saved object filename = here('00_data_ingest', '04_tiss1ue_robj_generated', paste0("droplet_", tissue_of_interest, "seurat_smartdrop-integrated-8272019.Robj")) load(file=filename) ``` # Export the final metadata ```{r} save_annotation_csv(tiss1, tiss1ue_of_interest, "droplet") ```