--- title: "G171_Downsampled" author: "KK" date: "January 28, 2020" output: html_document --- ```{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) ``` ```{r} 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.down0.5 <- Read10X("/net/waxman-server/mnt/data/waxmanlabvm_home/kkarri/G171/Analysis/not_bycol/G171B_down0.5/") #raw.data.down0.5_bycol <- Read10X("/net/waxman-server/mnt/data/waxmanlabvm_home/kkarri/G171/Analysis/Transcript_downsampled_0.5_bycol/Liver-10X_G171B") #raw.data.down0.9_bycol <- Read10X("/net/waxman-server/mnt/data/waxmanlabvm_home/kkarri/G171/Analysis/Transcript_downsampled_0.9_bycol/Liver-10X_G171B") raw.data.down0.25 <- Read10X("/net/waxman-server/mnt/data/waxmanlabvm_home/kkarri/G171/Analysis/not_bycol/G171B_down0.25/") raw.data.down1<- Read10X("/net/waxman-server/mnt/data/waxmanlabvm_home/kkarri/G171/Analysis/not_bycol/G171B_down1/") joint.bcs <- intersect(colnames(raw.data), colnames(raw.data.down0.5)) # Subset RNA and HTO counts by joint cell barcodes raw.data <- raw.data[, joint.bcs] raw.data.down0.5 <- raw.data.down0.5[, joint.bcs] raw.data.down0.5_bycol <- raw.data.down0.5_bycol[, joint.bcs] raw.data.down0.9_bycol <- raw.data.down0.9_bycol[, joint.bcs] raw.data.down0.25 <- raw.data.down0.25[, joint.bcs] raw.data.down1<- raw.data.down1[, joint.bcs] 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,] colnames(raw.data.down0.5) <- lapply(colnames(raw.data.down0.5), function(x) paste0(tissue_metadata_G171B$channel[1],'_',x)) meta.data0.5 = data.frame(row.names = colnames(raw.data.down0.5)) meta.data0.5['channel'] = tissue_metadata_G171B$channel[1] rnames = row.names(meta.data0.5) meta.data0.5 <- merge(meta.data0.5, tissue_metadata_G171B, sort = F) row.names(meta.data0.5) <- rnames # Order the cells alphabetically to ensure consistency. ordered_cell_names = order(colnames(raw.data.down0.5)) raw.data.down0.5 = raw.data.down0.5[,ordered_cell_names] meta.data0.5 = meta.data0.5[ordered_cell_names,] erccs <- grep(pattern = "^ERCC-", x = rownames(x = raw.data.down0.5), value = TRUE) percent.ercc <- Matrix::colSums(raw.data.down0.5[erccs, ])/Matrix::colSums(raw.data.down0.5) ercc.index <- grep(pattern = "^ERCC-", x = rownames(x = raw.data.down0.5), value = FALSE) raw.data.down0.5 <- raw.data.down0.5[-ercc.index,] ncRNA0.5 <- grep(pattern = "^ncRNA", x = rownames(x = raw.data.down0.5), value = TRUE) #ncountsum(raw.data.down0.5[ncRNA0.5, ])>0 lncRNA0.5 <- Matrix::colSums(raw.data.down0.5[ncRNA0.5, ]>0) mt.index.0.5 <- grep(pattern = "^mt-", x = rownames(x = raw.data.down0.5), value = FALSE) raw.data.down0.5 <- raw.data.down0.5[-mt.index.0.5,] #colnames(raw.data.down1) <- lapply(colnames(raw.data.down1), function(x) paste0(tissue_metadata_G171B$channel[1],'_',x)) meta.data1 = data.frame(row.names = colnames(raw.data.down1)) 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.down1)) raw.data.down1 = raw.data.down1[,ordered_cell_names] meta.data1 = meta.data1[ordered_cell_names,] erccs <- grep(pattern = "^ERCC-", x = rownames(x = raw.data.down1), value = TRUE) percent.ercc <- Matrix::colSums(raw.data.down1[erccs, ])/Matrix::colSums(raw.data.down1) ercc.index <- grep(pattern = "^ERCC-", x = rownames(x = raw.data.down1), value = FALSE) raw.data.down1 <- raw.data.down1[-ercc.index,] ncRNA1 <- grep(pattern = "^ncRNA", x = rownames(x = raw.data.down1), value = TRUE) lncRNA1 <- Matrix::colSums(raw.data.down1[ncRNA1, ]>0) mt.index.1 <- grep(pattern = "^mt-", x = rownames(x = raw.data.down1), value = FALSE) raw.data.down1 <- raw.data.down1[-mt.index.1,] colnames(raw.data.down0.25) <- lapply(colnames(raw.data.down0.25), function(x) paste0(tissue_metadata_G171B$channel[1],'_',x)) meta.data0.25 = data.frame(row.names = colnames(raw.data.down0.25)) meta.data0.25['channel'] = tissue_metadata_G171B$channel[1] rnames = row.names(meta.data0.25) meta.data0.25 <- merge(meta.data0.25, tissue_metadata_G171B, sort = F) row.names(meta.data0.25) <- rnames # Order the cells alphabetically to ensure consistency. ordered_cell_names = order(colnames(raw.data.down0.25)) raw.data.down0.25 = raw.data.down0.25[,ordered_cell_names] meta.data0.25 = meta.data0.25[ordered_cell_names,] erccs <- grep(pattern = "^ERCC-", x = rownames(x = raw.data.down0.25), value = TRUE) percent.ercc <- Matrix::colSums(raw.data.down0.25[erccs, ])/Matrix::colSums(raw.data.down0.25) ercc.index <- grep(pattern = "^ERCC-", x = rownames(x = raw.data.down0.25), value = FALSE) raw.data.down0.25 <- raw.data.down0.25[-ercc.index,] ncRNA0.25 <- grep(pattern = "^ncRNA", x = rownames(x = raw.data.down0.25), value = TRUE) lncRNA0.25 <- Matrix::colSums(raw.data.down0.25[ncRNA0.25, ]>0) mt.index.0.25 <- grep(pattern = "^mt-", x = rownames(x = raw.data.down0.25), value = FALSE) raw.data.down0.25 <- raw.data.down0.25[-mt.index.0.25,] droplet <- CreateSeuratObject(raw.data) # dropseq #droplet <- AddMetaData(object = droplet, meta.data1) droplet@meta.data$tech <- "G171B" #droplet <- subset(droplet, subset = nCount_RNA > 500 & nCount_RNA <50000) droplet <- NormalizeData(droplet, verbose = FALSE) droplet <- FindVariableFeatures(droplet, selection.method = "vst", nfeatures = 2000) droplet$stim <- "G171B" droplet <- SCTransform(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:25) droplet <- FindClusters(droplet, resolution = 0.4 ) droplet <- RunTSNE(droplet, reduction = "pca", dims = 1:25) droplet_down0.5 <- CreateSeuratObject(raw.data.down0.5) # dropseq droplet_down0.5 <- AddMetaData(object = droplet_down0.5, metadata = lncRNA0.5, col.name = "nlncRNA") droplet_down0.5 <- AddMetaData(object = droplet_down0.5, meta.data0.5) droplet_down0.5@meta.data$tech <- "G171B" #droplet_down0.5 <- subset(droplet_down0.5, subset = nCount_RNA > 500 & nCount_RNA <100000) #droplet_down0.5 <- subset(droplet_down0.5, subset = nFeature_RNA > 200 & nCount_RNA > 500) droplet_down0.5 <- NormalizeData(droplet_down0.5, verbose = FALSE) droplet_down0.5 <- FindVariableFeatures(droplet_down0.5, selection.method = "vst", nfeatures = 2000) droplet_down0.5$stim <- "G171B" droplet_down0.5 <- SCTransform(droplet_down0.5) droplet_down0.5 <- ScaleData(droplet_down0.5, verbose = FALSE) droplet_down0.5 <- RunPCA(droplet_down0.5, npcs = 30, verbose = FALSE) droplet_down0.5 <- RunUMAP(droplet_down0.5, reduction = "pca", dims = 1:25) droplet_down0.5 <- FindNeighbors(droplet_down0.5, reduction = "pca", dims = 1:25) droplet_down0.5 <- FindClusters(droplet_down0.5, resolution = 1 ) droplet_down0.5 <- RunTSNE(droplet_down0.5, reduction = "pca", dims = 1:25) droplet_down0.25 <- CreateSeuratObject(raw.data.down0.25) # dropseq droplet_down0.25 <- AddMetaData(object = droplet_down0.25, metadata = lncRNA0.25, col.name = "nlncRNA") droplet_down0.25 <- AddMetaData(object = droplet_down0.25, meta.data0.25) droplet_down0.25@meta.data$tech <- "G171B" #droplet_down0.25 <- subset(droplet_down0.25, subset = nCount_RNA > 500 & nCount_RNA <100000) #droplet_down0.25 <- subset(droplet_down0.25, subset = nFeature_RNA > 200 & nCount_RNA > 500) droplet_down0.25 <- NormalizeData(droplet_down0.25, verbose = FALSE) droplet_down0.25 <- FindVariableFeatures(droplet_down0.25, selection.method = "vst", nfeatures = 2000) droplet_down0.25$stim <- "G171B" droplet_down0.25 <- SCTransform(droplet_down0.25) droplet_down0.25 <- ScaleData(droplet_down0.25, verbose = FALSE) droplet_down0.25 <- RunPCA(droplet_down0.25, npcs = 30, verbose = FALSE) droplet_down0.25 <- RunUMAP(droplet_down0.25, reduction = "pca", dims = 1:25) droplet_down0.25 <- FindNeighbors(droplet_down0.25, reduction = "pca", dims = 1:25) droplet_down0.25 <- FindClusters(droplet_down0.25, resolution = 1 ) droplet_down0.25 <- RunTSNE(droplet_down0.25, reduction = "pca", dims = 1:25) droplet_down1 <- CreateSeuratObject(raw.data.down1) # dropseq droplet_down1 <- AddMetaData(object = droplet_down1, metadata = lncRNA1, col.name = "nlncRNA") droplet_down1 <- AddMetaData(object = droplet_down1, meta.data1) droplet_down1@meta.data$tech <- "G171B" #droplet_down1 <- subset(droplet_down1, subset = nCount_RNA > 500 & nCount_RNA <100000) #droplet_down1 <- subset(droplet_down1, subset = nFeature_RNA > 200 & nCount_RNA > 500) droplet_down1 <- NormalizeData(droplet_down1, verbose = FALSE) droplet_down1 <- FindVariableFeatures(droplet_down1, selection.method = "vst", nfeatures = 2000) droplet_down1$stim <- "G171B" droplet_down1 <- SCTransform(droplet_down1) droplet_down1 <- ScaleData(droplet_down1, verbose = FALSE) droplet_down1 <- RunPCA(droplet_down1, npcs = 30, verbose = FALSE) droplet_down1 <- RunUMAP(droplet_down1, reduction = "pca", dims = 1:25) droplet_down1 <- FindNeighbors(droplet_down1, reduction = "pca", dims = 1:25) droplet_down1 <- FindClusters(droplet_down1, resolution = 1) droplet_down1 <- RunTSNE(droplet_down1, reduction = "pca", dims = 1:25) ``` ```{r} pd1 <- UMAPPlot(droplet, reduction = "umap", label=TRUE, label.size=5,pt.size=0.5) DefaultAssay(droplet) <- "RNA" droplet <- NormalizeData(droplet) droplet <- ScaleData(droplet) d1 <- DotPlot(droplet, features = all_genes)+RotatedAxis() pd2.1 <- UMAPPlot(droplet_down1, reduction = "umap", label=TRUE, label.size=5, pt.size=0.5) DefaultAssay(droplet_down1) <- "RNA" droplet_down1 <- NormalizeData(droplet_down1) d2.1 <- DotPlot(droplet_down1, features = all_genes,label.size=2)+RotatedAxis() pd2.0.5 <- UMAPPlot(droplet_down0.5, reduction = "umap", label=TRUE, label.size=5, pt.size=0.5) DefaultAssay(droplet_down0.5) <- "RNA" droplet_down0.5 <- NormalizeData(droplet_down0.5) d2.0.5 <- DotPlot(droplet_down0.5, features = all_genes,label.size=2)+RotatedAxis() pd2.0.25 <- UMAPPlot(droplet_down0.25, reduction = "umap", label=TRUE, label.size=5, pt.size=0.5) DefaultAssay(droplet_down0.25) <- "RNA" droplet_down0.25 <- NormalizeData(droplet_down0.25) d2.0.25 <- DotPlot(droplet_down0.25, features = all_genes,label.size=2)+RotatedAxis() plot_grid(pd2.1,pd2.0.5 ,pd2.0.25, d2.1,d2.0.5,d2.0.25) vplot1 <- VlnPlot(droplet_down0.5, features = c("nFeature_RNA", "nlncRNA"), ncol = 2,pt.size = 0) plot1 <- FeatureScatter(droplet_down0.5, feature1 = "nCount_RNA", feature2 = "nFeature_RNA") plot1.1 <- FeatureScatter(droplet_down0.5, feature1 = "nCount_RNA", feature2 = "nlncRNA", pt.size = 0.1) vplot1 plot1 vplot2 <- VlnPlot(droplet_down0.25, features = c("nFeature_RNA", "nlncRNA"), ncol = 2, pt.size = 0) plot2 <- FeatureScatter(droplet_down0.25, feature1 = "nCount_RNA", feature2 = "nFeature_RNA") plot2.1 <- FeatureScatter(droplet_down0.25, feature1 = "nCount_RNA", feature2 = "nlncRNA", pt.size = 0.1) vplot2 plot2 vplot3 <- VlnPlot(droplet_down1, features = c("nFeature_RNA", "nlncRNA"), ncol = 2, pt.size = 0) plot3 <- FeatureScatter(droplet_down1, feature1 = "nCount_RNA", feature2 = "nFeature_RNA") plot3.1 <- FeatureScatter(droplet_down1, feature1 = "nCount_RNA", feature2 = "nlncRNA",pt.size = 0.1) vplot3 plot3 ``` ```{r} hep <- subset(droplet_down1, idents = c(0,1,3,6)) endo <- subset(droplet_down1, idents = 3) IK <- subset(droplet_down1, idents = 5) HSC <- subset(droplet_down1, idents = 6) unknown <- subset(droplet_down1, idents = 8) raw.data.hep <- as.matrix(GetAssayData(droplet_down1, slot = c("counts"))[, WhichCells(droplet_down1, ident = c(0,1,3,6))]) hep.sum <- as.data.frame(rowSums(raw.data.hep)) hep.mean<- as.data.frame(rowMeans(raw.data.hep)) hep <- merge(hep.sum, hep.mean) ``` ```{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','Igfbp7','Aqp1') 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("Dcn","Lama1","Nes") Dividing = "Top2a" Bplasma= "Jchain" Mac= "Csf1r" Chol="Sox9" Y_genes <- c("Uty","Ddx3y","Kdm5d","Eif2s3y", "Gm47283") sex <- c("Cyp2d9", "ncRNA-inter-chrX-15394","Cyp2c69", 'Mup20', 'Mup1','Mup12', 'Mup21', 'Cyp2d9') F_sex <- c('Sult3a1', 'A1bg', 'Fmo3', 'Cyp2b9', 'Sult2a1','Cyp2b13') all_genes = c(genes_hep, genes_endo, genes_kuppfer, Mac,genes_nk, genes_b, genes_bec, genes_immune, HSC, Dividing) ```