--- title: "Liver FACS Notebook" output: pdf_document: default html_notebook: default --- Enter the name of the tissue you want to analyze. ```{r} tissue_of_interest = "Liver" library(here) source(here("00_data_ingest", "02_tissue_analysis_rmd", "boilerplate.R")) tiss = load_tissue_facs(tissue_of_interest) library(scater) library(MAST) ``` ```{r} plate_metadata <- read.csv("/restricted/projectnb/waxmanlab/kkarri/scRNAseq_data_integration/Liver_facs_annotation.csv", sep=",", header = TRUE) colnames(plate_metadata)[1] <- "plate.barcode" raw.data = read.csv("/restricted/projectnb/waxmanlab/kkarri/scRNAseq_data_integration/liver_facs_scrna_data.csv", sep=",", row.names=1) colnames(plate_metadata)[1] <- "plate.barcode" plate.barcodes = lapply(colnames(raw.data), function(x) strsplit(strsplit(x, "_")[[1]][1], '.', fixed=TRUE)[[1]][2]) barcode.df = t.data.frame(as.data.frame(plate.barcodes)) rownames(barcode.df) = colnames(raw.data) barcode.df= cbind(barcode.df, colnames(raw.data)) colnames(barcode.df) = c('plate.barcode1', 'plate.barcode') # old annotations #rnames = row.names(barcode.df) #meta.data <- barcode.df #row.names(meta.data) <- rnames rnames = row.names(barcode.df) meta.data <- merge(barcode.df, plate_metadata, by='plate.barcode', sort = F) row.names(meta.data) <- rnames # Sort cells by cell name meta.data = meta.data[order(rownames(meta.data)), ] raw.data = raw.data[,rownames(meta.data)] # 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 smartseq <- CreateSeuratObject(raw.data) smartseq@meta.data$tech <- "smartseq2" tiss <- AddMetaData(object = tiss, meta.data) tiss <- AddMetaData(object = tiss, percent.ercc, col.name = "percent.ercc") # Change default name for sums of counts from nUMI to nReads colnames(tiss@meta.data)[colnames(tiss@meta.data) == 'nUMI'] <- 'nReads' # Create metadata columns for cell_ontology_class tiss@meta.data[,'free_annotation'] <- NA tiss@meta.data[,'cell_ontology_class'] <- NA #Calculate percent ribosomal genes. ribo.genes <- grep(pattern = "^Rp[sl][[:digit:]]", x = rownames(x = tiss@data), value = TRUE) percent.ribo <- Matrix::colSums(tiss@raw.data[ribo.genes, ])/Matrix::colSums(tiss@raw.data) tiss <- AddMetaData(object = tiss, metadata = percent.ribo, col.name = "percent.ribo") tiss <- FilterCells(object = tiss, subset.names = c("nGene", "nReads"), low.thresholds = c(300, 50000)) tiss <- process_tissue(tiss, 1e6) ``` We can visualize top genes in each principal component. ```{r, echo=FALSE} PCHeatmap(object = tiss, pc.use = 1:3, cells.use = 500, do.balanced = TRUE, label.columns = FALSE, num.genes = 8) ``` We then project onto just the top principal components. This has the effect of keeping the major directions of variation in the data and, ideally, supressing noise. A decent rule of thumb is to pick the elbow in the plot below. ```{r} PCElbowPlot(object = tiss) ``` Choose the number of principal components to use. ```{r} n.pcs = 11 ``` ## Cluster ```{r} # Set resolution res.used <- 1 tiss <- FindClusters(object = tiss, reduction.type = "pca", dims.use = 1:n.pcs, resolution = res.used, print.output = 0, save.SNN = TRUE) ``` We use tSNE solely to visualize the data. ```{r} tiss <- RunTSNE(object = tiss, dims.use = 1:n.pcs, seed.use = 10, perplexity=30) ``` ```{r} TSNEPlot(object = tiss, do.label = T, pt.size = 1.2, label.size = 4) TSNEPlot(object = tiss, do.label = T, pt.size = 1.2, label.size = 4, group.by="mouse.sex") ``` ## Label clusters using marker genes Check expression of genes useful for indicating cell type. ```{r} genes_hep_main = c('Alb', 'Ttr', 'Apoa1', 'Serpina1c') #hepatocyte genes_endo = c('Pecam1', 'Nrp1', 'Kdr','Oit3') # endothelial genes_kuppfer = c('Emr1', 'Clec4f', 'Cd68', 'Irf7') # Kuppfer cells genes_nk = c('Zap70', 'Il2rb', 'Nkg7', 'Cxcr6') # Natural Killer cells genes_b = c('Cd79a', 'Cd79b', 'Cd74', 'Cd19') # B Cells genes_all = c(genes_hep_main, genes_endo, genes_kuppfer, genes_nk, genes_b) ######## zonation markerks ##################### genes_pericentral = c('Cyp2e1', 'Glul', 'Oat', 'Gulo') genes_mid = c('Ass1', 'Hamp', 'Gstp1', 'Ubb') genes_periporatl= c('Cyp2f2', 'Pck1', 'Hal', 'Cdh1') ``` In the tSNE plots below, the intensity of each point represents the log-normalized gene expression $N_{ij}$. ```{r, echo=FALSE, fig.height=20, fig.width=16} FeaturePlot(tiss, genes_all, pt.size = 3, nCol = 4, cols.use = c("lightgrey", "blue"), no.legend = F) ``` Dotplots show, for each cluster and gene, the fraction of cells with at least one read for the gene (circle size) and the average scaled expression for that gene among the cells expressing it (circle color). ```{r, echo=FALSE} DotPlot(tiss, genes_pericentral, plot.legend = T, col.max = 2.5, x.lab.rot = T) ``` The low but nonzero levels of Albumin present in all clusters is consistent with a small amount of leakage, either through physical contamination or index hopping. Nevertheless, the absolute levels of expression confirm a sharp difference between the hepatocyte clusters and the others. ```{r, echo=FALSE, fig.height=3, fig.width=6} VlnPlot(tiss, 'Alb', use.raw = T, do.return = T) ``` To confirm the identity of a cluster, you can inspect the genes differentially expressed in that cluster compared to the others. Cell type markers - Hepatocytes, endothelial, kupffer, NK cell, B-cell ```{r} clust.markers.hep <- FindMarkers(object = tiss, ident.1 = c(1,2,3,4,6), ident.2 = c(0,5,8,7), only.pos = TRUE, min.pct = 0.25, thresh.use=0.25, logfc.threshold = 1.5) lnc_markers_hep <- grep(pattern = "^ncRNA", x= rownames(clust.markers.hep), value = TRUE) lnc_markers_hep lnc_markers_hep_table <- clust.markers.hep [lnc_markers_hep,] FeaturePlot(tiss,c(lnc_markers_hep),cols.use = c("grey", "red"), pt.size = 1, nCol = 4) DotPlot(tiss,lnc_markers_hep, plot.legend = T, col.max = 2.5, do.return = T) + coord_flip() top4_hep <- c("ncRNA_inter_chr1_291","ncRNA_inter_chr15_12684","ncRNA_inter_chr4_3295","ncRNA_as_chr7_6166") clust.markers.hep.MAST <- FindMarkers(object = tiss, ident.1 = c(1,2,3,4,6), ident.2 = c(0,5,8,7), only.pos = TRUE, test.use = 'MAST') lnc_markers_hep_MAST <- grep(pattern = "^ncRNA", x= rownames(clust.markers.hep.MAST), value = TRUE) lnc_markers_hep_MAST lnc_markers_hep_table <- clust.markers.hep.MAST [lnc_markers_hep_MAST,] ############################## endothelial #############################3 clust.markers.endo <- FindMarkers(object = tiss, ident.1 = 0, ident.2 = c(1,2,3,4,5,6,7,8), only.pos = TRUE, min.pct = 0.25, thresh.use=0.25, logfc.threshold = 1.5) clust.markers.endo <- FindMarkers(object = tiss, ident.1 = 0, ident.2 = c(1,2,3,4,5,6,7,8), only.pos = TRUE, min.pct = 0.25, thresh.use=0.25, logfc.threshold = 1.5, ) lnc_LM_Endo <- c('ncRNA_as_chr1_482','ncRNA_as_chr19_14777','ncRNA_as_chr11_9534','ncRNA_inter_chr13_11551','ncRNA_inter_chr1_129','ncRNA_inter_chr9_7979','ncRNA_as_chr2_1739','ncRNA_inter_chr17_14189','ncRNA_intra_chr17_13672','ncRNA_as_chr15_12584','ncRNA_inter_chr2_1432','ncRNA_as_chr2_1123','ncRNA_as_chr5_4263','ncRNA_as_chr15_12466') lnc_markers_endo <- grep(pattern = "^ncRNA", x= rownames(clust.markers.endo), value = TRUE) lnc_markers_endo lnc_markers_endo_table <- clust.markers.endo [lnc_markers_endo,] DotPlot(tiss,lnc_markers_endo, plot.legend = T, col.max = 2.5, do.return = T) + coord_flip() top4_endo <- c("ncRNA_as_chr4_3500","ncRNA_inter_chr10_9294","ncRNA_as_chr6_5266","ncRNA_inter_chr8_7126") FeaturePlot(tiss,top4_endo,cols.use = c("grey", "blue"), pt.size = 1, nCol = 4) ########################## Kupffer cell ################################### clust.markers.kupffer <- FindMarkers(object = tiss, ident.1 = 5, ident.2 = c(1,2,3,4,0,6,7,8), only.pos = TRUE, min.pct = 0.25, thresh.use=0.25, logfc.threshold = 1.5) lnc_markers_kupffer <- grep(pattern = "^ncRNA", x= rownames(clust.markers.kupffer), value = TRUE) lnc_markers_kupffer lnc_markers_kupffer_table <- clust.markers.kupffer [lnc_markers_kupffer,] lnc_markers_kupffer_table DotPlot(tiss,lnc_markers_kupffer, plot.legend = T, col.max = 2.5, do.return = T) + coord_flip() top2_kupffer <- c("ncRNA_inter_chr4_3805") FeaturePlot(tiss,top2_kupffer,cols.use = c("grey", "blue"), pt.size = 1, nCol = 4) ########################### NK cell markers clust.markers.nk <- FindMarkers(object = tiss, ident.1 = 8, ident.2 = c(1,2,3,4,0,6,7,5), only.pos = TRUE, min.pct = 0.25, thresh.use=0.25, logfc.threshold = 1.5) lnc_markers_nk <- grep(pattern = "^ncRNA", x= rownames(clust.markers.nk), value = TRUE) lnc_markers_nk lnc_markers_nk_table <- clust.markers.nk [lnc_markers_nk,] lnc_markers_nk_table DotPlot(tiss,lnc_markers_nk, plot.legend = T, col.max = 2.5, do.return = T) + coord_flip() top4_nk <- c("ncRNA_inter_chr17_13999","ncRNA_inter_chr15_12575","ncRNA_inter_chr12_10439","ncRNA_inter_chr10_8747") FeaturePlot(tiss,top4_nk,cols.use = c("grey", "blue"), pt.size = 1, nCol = 4) ##########################B cell markers ########################## clust.markers.b <- FindMarkers(object = tiss, ident.1 = 7, ident.2 = c(1,2,3,4,0,6,8,5), only.pos = TRUE, min.pct = 0.25, thresh.use=0.25, logfc.threshold = 1.5) lnc_markers_b <- grep(pattern = "^ncRNA", x= rownames(clust.markers.b), value = TRUE) lnc_markers_b lnc_markers_b_table <- clust.markers.b [lnc_markers_b,] lnc_markers_b_table DotPlot(tiss,lnc_markers_b, plot.legend = T, col.max = 2.5, do.return = T) + coord_flip() top4_b <- c("ncRNA_inter_chr6_5241","ncRNA_inter_chr16_13113","ncRNA_inter_chr16_13114","ncRNA_inter_chr5_4788") FeaturePlot(tiss,top4_b,cols.use = c("grey", "blue"), pt.size = 1, nCol = 4) ############### xeno responsive ############################ xeno_mCAR <- c ("ncRNA_inter_chr15_12684", "ncRNA_inter_chr2_1830", "ncRNA_as_chr16_13146", "ncRNA_as_chr10_8791", "ncRNA_inter_chr18_14478", "ncRNA_inter_chr10_9193", "ncRNA_as_chr7_5921", "ncRNA_inter_chr8_7430", "ncRNA_as_chr1_762", "ncRNA_as_chr19_14977", "ncRNA_inter_chr16_13177", "ncRNA_inter_chr7_6222", "ncRNA_as_chr1_979", "ncRNA_as_chr6_5635", "ncRNA_inter_chr10_9418", "ncRNA_as_chr15_12874", "ncRNA_as_chr4_3496", "ncRNA_as_chr16_13153", "ncRNA_as_chr7_6048", "ncRNA_inter_chr10_9302", "ncRNA_as_chr4_3298", "ncRNA_as_chr10_8962", "ncRNA_inter_chr8_6938", "ncRNA_as_chr1_1040", "ncRNA_inter_chr17_13999", "ncRNA_inter_chr15_12575", "ncRNA_as_chr17_13888", "ncRNA_as_chr10_9385" ) xenocar_12 <- c( "ncRNA_as_chr10_9385","ncRNA_as_chr19_14977","ncRNA_as_chr6_5635","ncRNA_as_chr7_5921","ncRNA_inter_chr10_9193","ncRNA_inter_chr10_9418","ncRNA_inter_chr15_12684","ncRNA_inter_chr16_13177","ncRNA_inter_chr7_6222","ncRNA_inter_chr8_7430") features.plot = c("ncRNA_inter_chr2_1830","ncRNA_as_chr16_13146","ncRNA_as_chr1_1040","ncRNA_inter_chr17_13999") ``` Using the markers, we can confidentaly label the clusters. We provide both a free annotation (where anything name can be used) and a cell ontology class. The latter uses a controlled vocabulary for easy comparison between studies and different levels of the taxonomy. ```{r} tiss <- StashIdent(object = tiss, save.name = "cluster.ids") cluster.ids <- c(0, 1, 2, 3, 4, 5, 6, 7, 8) free_annotation <- c( "endothelial cell of hepatic sinusoid", NA, NA, NA, NA, "Kupffer cell", "hepatocyte", "B cell", "NK/NKT cells") cell_ontology_class <-c( "endothelial cell of hepatic sinusoid", "hepatocyte", "hepatocyte", "hepatocyte", "hepatocyte", "Kupffer cell", "hepatocyte", "B cell", "natural killer cell") tiss = stash_annotations(tiss, cluster.ids, free_annotation, cell_ontology_class) ``` ## Subcluster Let's drill down on the hepatocytes. ```{r} subtiss = SubsetData(tiss, ident.use = c(1,2,3,4,6)) subtiss_endo = SubsetData(tiss, ident.use = 0) ``` ```{r} subtiss <- subtiss %>% ScaleData() %>% FindVariableGenes(do.plot = FALSE, x.high.cutoff = Inf, y.cutoff = 0.5) %>% RunPCA(do.print = FALSE) subtiss_endo <- subtiss_endo %>% ScaleData() %>% FindVariableGenes(do.plot = FALSE, x.high.cutoff = Inf, y.cutoff = 0.5) %>% RunPCA(do.print = FALSE,pc.genes = c(Endo_zonated_lncs,lnc_LM_Endo)) ``` ```{r} PCHeatmap(object = subtiss, pc.use = 1:3, cells.use = 20, do.balanced = TRUE, label.columns = FALSE, num.genes = 8) PCElbowPlot(subtiss) ``` zonated endothelial cells ```{r} zonated_endo <- c('Serinc3', 'Igfbp7', 'Ptprb', 'Dnase1l3', 'MALAT1', 'Rpl39', 'Dbp', 'Ehd3', 'Eif3e', 'Jak1', 'Tmem64', 'Pls3', 'Fabp4', 'Clec4g', 'Pfkfb3', 'Il6st', 'Rpl22l1', 'Rap1b', 'Ppap2a', 'Bmp2', 'Mrc1', 'Stt3a', 'Slc43a3', 'Gja4', 'Gnaq', 'Oit3', 'Cd36', 'Gpr116', 'Fcgr2b', 'Med13', 'Rbbp9', 'Snx2', 'Sparc', 'Bgn', 'Tmem2', 'Tgfbr2', 'Rasa1', 'Maf', 'Plxnc1', 'Kdr', 'Gpihbp1', 'Kit', 'Stab2', 'Sgms1', 'Med21', 'Rnd3', 'Crim1', 'Flt1', 'Cyp4b1', 'Jam2', 'Tspan7', 'Xpo1', 'Cdc27', 'Flt4', 'Ccnd1', 'Thbd', 'Prelp', 'Wnt2', 'Clec14a', 'Mmrn2', 'C1qtnf1', 'Ptgs1', 'Cd55', 'Tmx3', 'Marcks', '4931406P16Rik', 'Stab1', 'Mylip', 'Cpd', 'Ccnd2', 'Cd38', 'Pfdn4', 'Fcho2', 'Ets1', 'Lyve1', 'Tjp1', 'Rasgrp3', 'Klf6', 'Gimap4', 'Sema6a', 'Ly6a', 'Lrrc32', 'Gimap6', 'Acer2', 'Cnn2', 'Olfm1', 'Gpr182', 'Il13ra1', 'Notch1', 'Cotl1', 'Slfn5', 'Sgk1', 'Amotl1', 'Esam', 'Ccdc80', 'Cdh13', 'Plxnd1', 'Msrb3', 'Nid1', 'Kcnb1', 'Ecm1', 'Emilin1', 'Tpm4', 'Akap2', 'Tmem204', 'Msn', 'Ifi44', 'Cmtm3', 'Tm4sf1', 'Klf2', 'Tns3', 'Ier3', 'Mospd1', 'Nkd1', 'Nrp2', 'Tmem88', 'Sox18', 'Atp2b4', 'Mndal', 'Plk2', 'Gatm', 'ncRNA_as_chr1_482', 'Pea15a', 'Hspg2', 'Rasip1', 'Arrb2', 'Pcdh17', 'AI607873', 'Ppic', 'Tmem106a', 'Epb4.1l2', 'Rab3b', 'Fam189a2', 'Mef2c', 'Rin2', 'Tmtc3', 'Chst15', 'Tgfb1i1', 'ncRNA_as_chr19_14777', 'Fam174b', 'Sema3f', 'Colec11', 'Slc43a2', 'Csnk1e', 'Hecw2', 'Robo1', 'Sipa1', 'P2ry1', 'Serpina3g', 'Trim47', 'Fmnl2', 'Ltbp4', 'Itga9', 'Fgd5', 'Ppp1r16b', 'Lfng', 'Gimap1', 'Lrat', 'Wnt9b', 'Pear1', 'Jag1', 'Ctdspl', 'Dkk3', 'Clec1a', 'Dpy19l4', 'Cldn5', 'Cav1', 'Abcc4', 'Lmcd1', 'Ckap4', 'Adam23', 'Trove2', 'Plscr4', 'Rgl1', 'Cyyr1', 'Vipr1', 'Dusp7', 'Zfp715', 'Zeb2', 'Trim30a', 'Ifi204', 'Plekho1', 'Rnf144a', 'Osmr', 'Lphn1', '9430020K01Rik', 'Fgfr1', 'Mfge8', 'ncRNA_as_chr11_9534', 'Cc2d2a', 'Sh3rf1', 'ncRNA_inter_chr13_11551', 'Hip1', 'Ldb2', 'Tbc1d19', 'Serpina3f', 'Hlx', 'Fmn1', 'Pld1', 'Fndc5', 'Ifit2', 'Syt12', 'Gpr56', 'Lama4', 'Col14a1', 'ncRNA_inter_chr1_129', 'Kank3', 'Gnai1', 'ncRNA_inter_chr9_7979', 'Dchs1', 'Impdh1', 'Anxa3', 'Igfbp3', 'ncRNA_as_chr2_1739', 'Slc7a7', 'Exoc3l', 'Plekhg1', 'Chst2', 'B4galt4', 'ncRNA_inter_chr17_14189', 'Clasp2', 'Rab31', 'Efnb2', 'ncRNA_intra_chr17_13672', 'Ntf3', 'Maml3', 'Selp', 'Ppp1r9a', 'Paqr4', 'ncRNA_as_chr15_12584', 'Col1a2', 'Dip2a', 'ncRNA_inter_chr2_1432', 'Chst7', 'ncRNA_as_chr2_1123', 'Gata2', 'ncRNA_as_chr5_4263', 'Klf4', 'Dll4', 'ncRNA_as_chr15_12466', 'Sox17', 'Flrt1', 'Me2', 'Tmem44', 'Rasd1', 'Slc41a1', 'Enpp6', 'Nova2', 'Nrarp', 'Twist1', 'Pla2r1', 'Rbfox3', 'Txndc16', 'Itgb3', 'Des', 'Msr1', 'Cadm3', 'Ddx26b', 'Ntm', 'Mecom', 'Map3k8', 'Atpaf1', 'Bicc1', 'Cd34') ``` ```{r} sub.n.pcs = 8 sub.res.use =3.5 subtiss <- subtiss %>% FindClusters(reduction.type = "pca", dims.use = 1:sub.n.pcs, resolution = 1, print.output = 0, save.SNN = TRUE, force.recalc = TRUE) %>% RunTSNE(dims.use = 1:sub.n.pcs, seed.use = 10, perplexity=30) TSNEPlot(object = subtiss, do.label = T, pt.size = 1, label.size = 4) sub.n.pcs = 8 sub.res.use =3.5 subtiss_endo <- subtiss_endo %>% FindClusters(reduction.type = "pca", dims.use = 1:sub.n.pcs, resolution = 1, print.output = 0, save.SNN = TRUE, force.recalc = TRUE,genes.use = c(Endo_zonated_lncs,lnc_LM_Endo)) %>% RunTSNE(dims.use = 1:sub.n.pcs, seed.use = 10, perplexity=30) TSNEPlot(object = subtiss_endo, do.label = T, pt.size = 1, label.size = 4) ``` ```{r} Endo_zonated_lncs <- c('ncRNA_as_chr10_8920', 'ncRNA_as_chr10_8922', 'ncRNA_as_chr10_8929', 'ncRNA_as_chr10_9411', 'ncRNA_as_chr11_10121', 'ncRNA_as_chr11_9897', 'ncRNA_as_chr12_10617', 'ncRNA_as_chr13_11477', 'ncRNA_as_chr15_12584', 'ncRNA_as_chr18_14332', 'ncRNA_as_chr1_38', 'ncRNA_as_chr1_733', 'ncRNA_as_chr2_1188', 'ncRNA_as_chr2_1893', 'ncRNA_as_chr4_3220', 'ncRNA_as_chr4_3316', 'ncRNA_as_chr5_4263', 'ncRNA_as_chr6_5147', 'ncRNA_as_chr7_6440', 'ncRNA_as_chr8_6954', 'ncRNA_as_chrX_15320', 'ncRNA_inter_chr10_8466', 'ncRNA_inter_chr10_8771', 'ncRNA_inter_chr10_9037', 'ncRNA_inter_chr10_9093', 'ncRNA_inter_chr11_10099', 'ncRNA_inter_chr11_9575', 'ncRNA_inter_chr11_9966', 'ncRNA_inter_chr12_10529', 'ncRNA_inter_chr12_10657', 'ncRNA_inter_chr13_11551', 'ncRNA_inter_chr15_12347', 'ncRNA_inter_chr15_12517', 'ncRNA_inter_chr16_13214', 'ncRNA_inter_chr16_13449', 'ncRNA_inter_chr17_14034', 'ncRNA_inter_chr19_14875', 'ncRNA_inter_chr1_602', 'ncRNA_inter_chr1_812', 'ncRNA_inter_chr1_900', 'ncRNA_inter_chr1_937', 'ncRNA_inter_chr2_1432', 'ncRNA_inter_chr2_1906', 'ncRNA_inter_chr2_2001', 'ncRNA_inter_chr3_2237', 'ncRNA_inter_chr3_2486', 'ncRNA_inter_chr3_2508', 'ncRNA_inter_chr3_2697', 'ncRNA_inter_chr4_3058', 'ncRNA_inter_chr4_3350', 'ncRNA_inter_chr4_3379', 'ncRNA_inter_chr4_3499', 'ncRNA_inter_chr5_4065', 'ncRNA_inter_chr5_4347', 'ncRNA_inter_chr5_4705', 'ncRNA_inter_chr5_4791', 'ncRNA_inter_chr6_4914', 'ncRNA_inter_chr6_4966', 'ncRNA_inter_chr6_5217', 'ncRNA_inter_chr6_5612', 'ncRNA_inter_chr7_5879', 'ncRNA_inter_chr7_5931', 'ncRNA_inter_chr7_6006', 'ncRNA_inter_chr7_6061', 'ncRNA_inter_chr8_6826', 'ncRNA_inter_chr8_7357', 'ncRNA_inter_chr9_8323', 'ncRNA_inter_chr9_8339', 'ncRNA_intra_chr19_14833') sub.cluster.ids <- c(0, 2, 1) sub.free_annotation <- c("periportal","pericentral","periportal") sub.cell_ontology_class <- c("endothelial","endothelial","endothelial") subtiss_endo = stash_annotations(subtiss_endo, sub.cluster.ids, sub.free_annotation, sub.cell_ontology_class) tiss = stash_subtiss_in_tiss(tiss, subtiss) ``` ```{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(subtiss,c('Mup20', 'Mup1','Mup12', 'Mup21', 'Cyp2d9', 'Xist', 'A1bg', 'Cyp2c69'),cols.use = c("grey", "red"), pt.size = 3, nCol = 2) DotPlot(subtiss,c('Mup20', 'Mup1','Mup12', 'Mup21', 'Cyp2d9', 'Xist', 'A1bg', 'Cyp2c69'), plot.legend = T, col.max = 2.5, do.return = T) + coord_flip() ``` Liver zonation markers ```{r} genes_zones = c('Cyp2e1', 'Glul', 'Oat', 'Gulo', 'Ass1', 'Hamp', 'Gstp1', 'Ubb', 'Cyp2f2', 'Pck1', 'Hal', 'Cdh1') FeaturePlot(subtiss,c(genes_zones),cols.use = c("grey", "red"), pt.size = 1, nCol = 4) DotPlot(subtiss,c(genes_zones), plot.legend = T, col.max = 2.5, do.return = T) + coord_flip() TSNEPlot(object = subtiss, do.label = T, pt.size = 1, label.size = 4, group.by="free_annotation") TSNEPlot(object = tiss, do.label = T, pt.size = 1, label.size = 4, group.by="free_annotation") ``` ```{r} sub.cluster.ids <- c(0, 1, 2, 3, 4) sub.free_annotation <- c("midlobular hepatocyte -1(F)", "midlobular hepatocyte (M)", "pericentral hepatocyte (M)", "periportal hepatocyte (M)", "midlobular hepatocyte-2 (F)") sub.cell_ontology_class <- c("hepatocyte", "hepatocyte", "hepatocyte", "hepatocyte", "hepatocyte") subtiss = stash_annotations(subtiss, sub.cluster.ids, sub.free_annotation, sub.cell_ontology_class) tiss = stash_subtiss_in_tiss(tiss, subtiss) ``` The multitude of clusters of each type correspond mostly to individual animals/sexes. ```{r} table(FetchData(subtiss1, c('mouse.sex','ident')) %>% droplevels()) ``` ## Checking for batch effects Color by metadata, like plate barcode, to check for batch effects. Here we see that the clusters are segregated by sex. ```{r} TSNEPlot(object = tiss, do.return = TRUE, group.by = "mouse.id") ``` Nevertheless, every cluster contains cells from multiple mice. ```{r} table(FetchData(tiss, c('mouse.id','ident')) %>% droplevels()) ``` # Final coloring Color by cell ontology class on the original tSNE. ```{r} TSNEPlot(object = tiss, group.by = "cell_ontology_class") ```