Liver Droplet- Raw data Notebook

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

# 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)

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"

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)

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)

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

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

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"

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)

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"

```

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)

 

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

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")
 

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

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)

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

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.

PCElbowPlot(object = tiss1)

Choose the number of principal components to use.

# 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.

# 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.

# 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)

TSNEPlot(object = tiss1, do.label = T, pt.size = 1.2, label.size = 4)

Compare to previous annotations

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)

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)

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

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

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.

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.

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.

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))

subtiss1 <- subtiss1 %>% ScaleData() %>%
  FindVariableGenes(do.plot = FALSE, x.high.cutoff = Inf, y.cutoff = 0.5) %>%
  RunPCA(do.print = FALSE)

PCHeatmap(object = subtiss1, pc.use = 1:3, cells.use = 20, do.balanced = TRUE, label.columns = FALSE, num.genes = 8)
PCElbowPlot(subtiss1)

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)

BuildClusterTree(subtiss1)

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.

table(FetchData(subtiss1, c('mouse.sex','ident')) %>% droplevels())

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

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

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

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.

TSNEPlot(object = subtiss1, do.return = TRUE, group.by = "mouse.sex")

Final coloring

Color by cell ontology class on the original TSNE.

TSNEPlot(object = tiss1, do.return = TRUE, group.by = "cell_ontology_class")

Save the Robject for later

filename = here('00_data_ingest', '04_tiss1ue_robj_generated', 
                     paste0("droplet_", tiss1ue_of_interest, "refinedcells_seurat_tiss1.Robj"))
print(filename)
save(tiss1, file=filename)

# 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

save_annotation_csv(tiss1, tiss1ue_of_interest, "droplet")