library(Tempora) library(Seurat) library(ggplot2) library(cowplot) library(stringr) library(dplyr) library(tidyverse) library(pbapply) library(clusterProfiler) library(enrichplot) library(data.table) library(writexl) library(EnhancedVolcano) library(ggsunburst) library(ggrepel) library(RColorBrewer) library(svglite) library(ggpubr) source("functions.R") save_plot = function(plotobj,outname, fig.width, fig.height) { print(paste(outname, fig.width, fig.height)) fname=paste(outname, "png", sep=".") print(paste("Saving to file", fname)) png(filename=fname, width = fig.width, height = fig.height, units = 'in', res = 300)#width = fig.width*100, height=fig.height*100) plot(plotobj) dev.off() fname=paste(outname, "pdf", sep=".") print(paste("Saving to file", fname)) pdf(file=fname, width = fig.width, height=fig.height) plot(plotobj) dev.off() fname=paste(outname, "svg", sep=".") print(paste("Saving to file", fname)) svglite::svglite(file = fname, width = fig.width, height = fig.height) plot(plotobj) dev.off() return(plotobj) } files <- Sys.glob("../21068/calls_full/*_M*/outs/raw*/*.mtx.gz") allfiles.raw = list() allABs.raw = list() for (file in files) { samplename = str_split(dirname(file), "/")[[1]][4] foldername = dirname(file) print(paste(samplename, foldername)) h5file = Read10X(foldername,unique.features = TRUE) allfiles.raw[[samplename]] = h5file$`Gene Expression` # the antibody capture from cellranger is incorrect! allABs.raw[[samplename]] = h5file$`Antibody Capture` } statDF = data.frame(matrix(ncol=2,nrow=0, dimnames=list(NULL, c("RawAll", "RawFiltered")))) toCellPrefix = function(x) { datasetprefix = x datasetprefix = str_replace_all(datasetprefix, "[./]", "") datasetprefix = str_split(datasetprefix, "_")[[1]][1] return(paste(datasetprefix, sep="")) } makesum = function(a, suffix) { out = {} out["sum"] = sum(a) return(out) } makeSeuratObj = function(matrix, proj) { obj = CreateSeuratObject(matrix, project=proj) print("Renaming Cells") obj <- RenameCells(obj, add.cell.id=proj) print(paste("Seurat obj project", obj@project.name)) mtPattern = "^mt-" rplPattern = "^Rpl" rpsPattern = "^Rps" rpPattern = "^Rps|^Rpl" selGenes = rownames(obj)[grepl(rownames(obj), pattern=mtPattern)] print(paste("Got a total of mt-Genes:", length(selGenes), paste(head(selGenes), collapse=", "))) selGenes = rownames(obj)[grepl(rownames(obj), pattern=rplPattern)] print(paste("Got a total of Rpl-Genes:", length(selGenes), paste(head(selGenes), collapse=", "))) selGenes = rownames(obj)[grepl(rownames(obj), pattern=rpsPattern)] print(paste("Got a total of Rps-Genes:", length(selGenes), paste(head(selGenes), collapse=", "))) selGenes = rownames(obj)[grepl(rownames(obj), pattern=rpPattern)] print(paste("Got a total of Rp-Genes:", length(selGenes), paste(head(selGenes), collapse =", "))) obj[["percent.mt"]] <- PercentageFeatureSet(obj, pattern = mtPattern) obj[["percent.rpl"]] <- PercentageFeatureSet(obj, pattern = rplPattern) obj[["percent.rps"]] <- PercentageFeatureSet(obj, pattern = rpsPattern) obj[["percent.rp"]] <- PercentageFeatureSet(obj, pattern = rpPattern) return(obj) } objlist = list() for (x in names(allfiles.raw)) { matrix = allfiles.raw[[x]] origSize = ncol(matrix) filteredAll = dim(matrix)[2] filteredObj = makeSeuratObj(matrix, x) rawFilteredSize = ncol(filteredObj) filteredObj <- NormalizeData(filteredObj, verbose = FALSE) filteredObj <- FindVariableFeatures(filteredObj, verbose = FALSE) objlist[[x]] = filteredObj statDF <- rbind(statDF, data.frame("RawAll"=origSize, "RawFiltered"=rawFilteredSize)) } for (name in names(objlist)) { print(name) #p=VlnPlot(objlist[[name]], features = c("nFeature_RNA", "nCount_RNA", "percent.mt"), ncol = 3, pt.size = 0) #save_plot(p, paste(name, "violins_qc", sep="_"), fig.width=10, fig.height=6) plot1 <- FeatureScatter(objlist[[name]], feature1 = "nCount_RNA", feature2 = "percent.mt") plot2 <- FeatureScatter(objlist[[name]], feature1 = "nCount_RNA", feature2 = "nFeature_RNA") save_plot(plot1 + plot2, paste("plots_qc/", paste(name, "scatter_ncount_mt", sep="_"), sep=""), fig.width=10, fig.height=6) plot1 <- FeatureScatter(objlist[[name]], feature1 = "nCount_RNA", feature2 = "percent.rp") plot2 <- FeatureScatter(objlist[[name]], feature1 = "nCount_RNA", feature2 = "nFeature_RNA") save_plot(plot1 + plot2, paste("plots_qc/", paste(name, "scatter_ncount_rp", sep="_"), sep=""), fig.width=10, fig.height=6) } objlist <- lapply(X = objlist, FUN = function(obj) { # mt content: https://www.ncbi.nlm.nih.gov/pmc/articles/PMC6072887/ obj <- subset(obj, subset = nFeature_RNA > 100 & nFeature_RNA < 6000 & nCount_RNA > 500) obj <- subset(obj, subset = percent.mt < 15) print(obj) return(obj) }) for (name in names(objlist)) { p=VlnPlot(objlist[[name]], features = c("nFeature_RNA", "nCount_RNA", "percent.mt"), ncol = 3, pt.size = 0) save_plot(p, paste("plots_qc/", paste(name, "filtered_violins_qc", sep="_"), sep=""), fig.width=10, fig.height=6) plot1 <- FeatureScatter(objlist[[name]], feature1 = "nCount_RNA", feature2 = "percent.mt") plot2 <- FeatureScatter(objlist[[name]], feature1 = "nCount_RNA", feature2 = "nFeature_RNA") save_plot(plot1 + plot2, paste("plots_qc/", paste(name, "filtered_scatter_ncount_mt", sep="_"), sep=""), fig.width=10, fig.height=6) plot1 <- FeatureScatter(objlist[[name]], feature1 = "nCount_RNA", feature2 = "percent.rp") plot2 <- FeatureScatter(objlist[[name]], feature1 = "nCount_RNA", feature2 = "nFeature_RNA") save_plot(plot1 + plot2, paste("plots_qc/", paste(name, "filtered_scatter_ncount_rp", sep="_"), sep=""), fig.width=10, fig.height=6) } htoObjList = list() for (name in names(objlist)) { cellsGEX = length(colnames(objlist[[name]])) cellsAB = length(colnames(allABs.raw[[name]])) gexnames = substring(colnames(objlist[[name]]), 11) htosKeep = setdiff(rownames(allABs.raw[[name]]), c("Mtag3", "Mtag10")) htoCountMatrix = allABs.raw[[name]] htoCountMatrix = htoCountMatrix[htosKeep, ] htoCountMatrix = htoCountMatrix[, colSums(htoCountMatrix) >= 1] cellsJoint = intersect(gexnames, colnames(htoCountMatrix)) htoCountMatrix = htoCountMatrix[, cellsJoint] colnames(htoCountMatrix) = paste(objlist[[name]]@project.name, colnames(htoCountMatrix), sep="_") xobj = subset(objlist[[name]], cells=colnames(htoCountMatrix)) # Normalize RNA data with log normalization xobj <- NormalizeData(xobj) # Find and scale variable features xobj <- FindVariableFeatures(xobj, selection.method = "mean.var.plot") xobj <- ScaleData(xobj, features = VariableFeatures(xobj)) xobj[["HTO"]] <- CreateAssayObject(counts = htoCountMatrix) xobj <- NormalizeData(xobj, assay = "HTO", normalization.method = "CLR") xobj <- HTODemux(xobj, assay = "HTO", positive.quantile = 0.99) print(table(xobj$HTO_classification.global)) htoObjList[[name]] = xobj } for (name in names(objlist)) { print(table(htoObjList[[name]]$HTO_classification.global)) } for (name in names(htoObjList)) { p=VlnPlot(htoObjList[[name]], features = c("nFeature_RNA", "nCount_RNA", "percent.mt"), ncol = 3, pt.size = 0, group.by = "HTO_classification.global") save_plot(p, paste("plots_hto_qc/", name, "_hto_violins_qc", sep=""), fig.width=10, fig.height=6) r=RidgePlot(htoObjList[[name]], assay = "HTO", features = rownames(htoObjList[[name]][["HTO"]]), ncol = 2) save_plot(r, paste("plots_hto_qc/", name, "_hto_ridge_qc", sep=""), fig.width=10, fig.height=14) #f = FeatureScatter(htoObjList[[name]], feature1 = "hto_Blood", feature2 = "hto_Thrombus") #save_plot(f, paste("plots_hto_qc/", name, "_hto_scatter_qc", sep=""), fig.width=10, fig.height=6) p=VlnPlot(htoObjList[[name]], features = "nCount_RNA", pt.size = 0, log = TRUE, group.by= "HTO_classification.global") save_plot(p, paste("plots_hto_qc/", name, "_hto_vln_classification_qc", sep=""), fig.width=10, fig.height=6) print(name) print(table(htoObjList[[name]]$HTO_classification)) } htoFiltered <- lapply(X = htoObjList, FUN = function(obj) { # mt content: https://www.ncbi.nlm.nih.gov/pmc/articles/PMC6072887/ print(obj) print(table(obj$HTO_classification)) obj <- subset(obj, subset = HTO_classification.global == "Singlet") print(obj) print(table(obj$HTO_classification)) return(obj) }) print("cells per experiment") print(mapply(sum, lapply(htoFiltered, function(x) {dim(x)[2]}))) print("total cells") print(sum(mapply(sum, lapply(htoFiltered, function(x) {dim(x)[2]})))) save.image("htodemux.v4.RData", compress = FALSE) sample2info = list( "sample_M1"= list(tags_lung=c("Mtag1", "Mtag2.1", "Mtag4", "Mtag5"), tags_kidney=c("Mtag6", "Mtag7", "Mtag8", "Mtag9"), animals=c(346, 345, 344, 343), conditions=c("cre+", "cre-", "cre+", "cre-")), "sample_M2" = list(tags_lung=c("Mtag1", "Mtag2.1", "Mtag4", "Mtag5"), tags_kidney=c("Mtag6", "Mtag7", "Mtag8", "Mtag9"), animals=c(342, 341, 340, 338), conditions=c("cre+", "cre-", "cre+", "cre-")) ) lungSampleList = list() kidneySampleList = list() for (name in names(htoFiltered)) { print(name) sampleTags = sample2info[[name]]$tags_lung sampleLung = subset(htoFiltered[[name]], subset=HTO_classification %in% sampleTags) toAnimal = sample2info[[name]]$animals toCondition = sample2info[[name]]$conditions names(toAnimal) = sampleTags names(toCondition) = sampleTags animalFeature = toAnimal[sampleLung$HTO_classification] conditionFeature = toCondition[sampleLung$HTO_classification] names(animalFeature) = names(sampleLung$HTO_classification) names(conditionFeature) = names(sampleLung$HTO_classification) sampleLung$animalFeature = animalFeature sampleLung$conditionFeature = conditionFeature lungSampleList[[name]] = sampleLung sampleTags = sample2info[[name]]$tags_kidney sampleKidney = subset(htoFiltered[[name]], subset=HTO_classification %in% sampleTags) names(toAnimal) = sampleTags names(toCondition) = sampleTags animalFeature = toAnimal[sampleKidney$HTO_classification] conditionFeature = toCondition[sampleKidney$HTO_classification] names(animalFeature) = names(sampleKidney$HTO_classification) names(conditionFeature) = names(sampleKidney$HTO_classification) sampleKidney$animalFeature = animalFeature sampleKidney$conditionFeature = conditionFeature kidneySampleList[[name]] = sampleKidney } for (name in names(lungSampleList)) { sampleElem = lungSampleList[[name]] print(sampleElem) print(unique(paste(name, sampleElem$animalFeature, sampleElem$conditionFeature, sampleElem$HTO_classification))) } for (name in names(kidneySampleList)) { sampleElem = kidneySampleList[[name]] print(sampleElem) print(unique(paste(name, sampleElem$animalFeature, sampleElem$conditionFeature, sampleElem$HTO_classification))) } conditions.lung = c() animals.lung = c() for (name in names(lungSampleList)) { conditions.lung = c(conditions.lung, lungSampleList[[name]]$conditionFeature) animals.lung = c(animals.lung, lungSampleList[[name]]$animalFeature) } conditions.kidney = c() animals.kidney = c() for (name in names(kidneySampleList)) { conditions.kidney = c(conditions.kidney, kidneySampleList[[name]]$conditionFeature) animals.kidney = c(animals.kidney, kidneySampleList[[name]]$animalFeature) } # Basic function to convert human to mouse gene names #https://github.com/satijalab/seurat/issues/2493 s.genes <- convertHumanGeneList(cc.genes$s.genes) g2m.genes <- convertHumanGeneList(cc.genes$g2m.genes) #### # # # DO LUNG ANALYSIS # # #### print(s.genes) print(g2m.genes) objlist.lung <- lapply(X = lungSampleList, FUN = function(x) { DefaultAssay(x) = "RNA" x }) objlist.lung <- lapply(X = objlist.lung, FUN = function(x) { DefaultAssay(x) = "RNA" print(x) x <- CellCycleScoring( x, g2m.features = g2m.genes, s.features = s.genes) #"nCount_RNA","S.Score", "G2M.Score" suppressWarnings(x <- SCTransform(x,vars.to.regress = c('percent.mt', 'percent.rp'), verbose = T)) }) features <- SelectIntegrationFeatures(object.list = objlist.lung, nfeatures = 2000) print(length(features)) objlist.lung <- PrepSCTIntegration(object.list = objlist.lung, anchor.features = features) objlist.lung.anchors <- FindIntegrationAnchors(object.list = objlist.lung, normalization.method = "SCT", anchor.features = features) obj.lung.integrated <- IntegrateData(anchorset = objlist.lung.anchors, normalization.method = "SCT",verbose=T) print(obj.lung.integrated) #write_rds(obj.lung.integrated, "lung.integrated.rds", compress="none") #save.image(file="lung.integration.v4.final.RData", compress = F) DefaultAssay(obj.lung.integrated) = "integrated" #obj.lung.integrated <- ScaleData(obj.lung.integrated, verbose = FALSE) obj.lung.integrated <- RunPCA(obj.lung.integrated, verbose = FALSE) obj.lung.integrated <- FindNeighbors(obj.lung.integrated, dims = 1:50) obj.lung.integrated <- FindClusters(obj.lung.integrated, resolution = 0.2) obj.lung.integrated$idents = Idents(obj.lung.integrated) obj.lung.integrated <- RunUMAP(obj.lung.integrated, dims = 1:50, n.neighbors = 50, return.model = TRUE) p=DimPlot(obj.lung.integrated, pt.size = 0.001, label=T) save_plot(p, "plots_lung/dimplot_umap", fig.width=12, fig.height=8) p=VlnPlot(obj.lung.integrated, c("percent.mt"), pt.size=0, group.by="idents") save_plot(p, "plots_lung/vln_percentmt", fig.width=12, fig.height=5) p=FeaturePlot(obj.lung.integrated, c("percent.mt"), pt.size=0, label=T) save_plot(p, "plots_lung/fplt_percentmt", fig.width=12, fig.height=8) # as suggested here: https://github.com/satijalab/seurat/issues/2180#issuecomment-832064768 DefaultAssay(obj.lung.integrated) <- "RNA" obj.lung.integrated <- NormalizeData(obj.lung.integrated, verbose = FALSE) p=FeaturePlot(obj.lung.integrated, c("Ccr2"), pt.size=0, label=T) save_plot(p, "plots_lung/fplt_ccr2", fig.width=12, fig.height=8) obj.lung.integrated$conditionFeature = conditions.lung obj.lung.integrated$animalFeature = animals.lung print(obj.lung.integrated) p=DimPlot(obj.lung.integrated, pt.size=0.001, label=T, split.by="conditionFeature") save_plot(p, "plots_lung/dimplot_umap_by_cond", fig.width=20, fig.height=8) p=DimPlot(obj.lung.integrated, pt.size=0.001, label=T, split.by="animalFeature") save_plot(p, "plots_lung/dimplot_umap_by_animal", fig.width=40, fig.height=8) p=DimPlot(obj.lung.integrated, pt.size=0.001, label=T, split.by="HTO_classification") save_plot(p, "plots_lung/dimplot_umap_by_hto", fig.width=40, fig.height=8) countByOIAbs = getCellCountDF(obj.lung.integrated, prefix="", select_by = "idents", group_by="orig.ident", relative=F, show.percent=F, outname="plots_lung/ia_count_by_origident.tsv") countByAFAbs = getCellCountDF(obj.lung.integrated, prefix="", select_by = "idents", group_by="animalFeature", relative=F, show.percent=F, outname="plots_lung/ia_count_by_animal.tsv") countByCFAbs = getCellCountDF(obj.lung.integrated, prefix="", select_by = "idents", group_by="conditionFeature", relative=F, show.percent=F, outname="plots_lung/ia_count_by_condition.tsv") DefaultAssay(obj.lung.integrated) = "RNA" deResTT = makeDEResults(obj.lung.integrated, assay="RNA", test="t") exprdfTT = getDEXpressionDF(obj.lung.integrated, deResTT, assay="RNA") write.table(exprdfTT, "plots_lung/expr_test_t.tsv", sep="\t", row.names=F, quote = F) write_xlsx(exprdfTT, "plots_lung/expr_test_t.xlsx") markers.use.tt=subset(exprdfTT ,avg_log2FC>0&p_val_adj<0.01&!startsWith(gene, "mt-")&!startsWith(gene, "Rp")) finalMarkers.use.tt = markers.use.tt %>% arrange(p_val_adj, desc(abs(pct.1)*abs(avg_log2FC),)) %>% group_by(clusterID) %>% dplyr::slice(1:20) finalMarkers.use.tt finalMarkers.all.tt = markers.use.tt %>% arrange(p_val_adj, desc(abs(pct.1)*abs(avg_log2FC),)) %>% group_by(clusterID) finalMarkers.all.tt = finalMarkers.all.tt[order(finalMarkers.all.tt$clusterID),] finalMarkers.all.tt p_dp_genes_idents = DotPlot(obj.lung.integrated, features = unique(finalMarkers.use.tt$gene), assay="RNA")+ theme(axis.text.x = element_text(angle = 90, hjust = 1.0, vjust = 0.5)) save_plot(p_dp_genes_idents, "plots_lung/dotplot_idents_genes", 45, 8) cluster_mean_expr = get_population_expression_data(obj.lung.integrated, group="idents", outname="plots_lung/cluster_mean_expr_lung") cells.lung.cre_pos = names(obj.lung.integrated$conditionFeature[obj.lung.integrated$conditionFeature == "cre+"]) cells.lung.cre_neg = names(obj.lung.integrated$conditionFeature[obj.lung.integrated$conditionFeature == "cre-"]) print(paste(length(cells.lung.cre_pos), length(cells.lung.cre_neg), length(obj.lung.integrated$conditionFeature))) de_lung_cre = compareCellsByCluster(obj.lung.integrated, cells.lung.cre_pos, cells.lung.cre_neg, "lung_cre_pos", "lung_cre_neg", group.by="idents", test="MAST", outfolder="de_lung", fcCutoff=0.1) de_lung_cre_t = compareCellsByCluster(obj.lung.integrated, cells.lung.cre_pos, cells.lung.cre_neg, "lung_cre_pos", "lung_cre_neg", group.by="idents", test="t", outfolder="de_lung_t", fcCutoff=0.1) de_lung_cre_all = compareClusters(obj.lung.integrated, cells.lung.cre_pos, cells.lung.cre_neg, "lung_cre_pos", "lung_cre_neg", test="MAST", outfolder="de_lung", fcCutoff=0.1) makeVolcanos(de_lung_cre, "Volcano", "lung_volcanos/volcanos", FCcutoff=0.1) makeVolcanos(de_lung_cre_t, "Volcano", "lung_volcanos_t/volcanos", FCcutoff=0.1) makeVolcanos(list("de_lung_cre_all"=de_lung_cre_all), "Volcano", "lung_volcanos/volcanos_all", FCcutoff=0.1) write_rds(obj.lung.integrated, "lung.integrated.rds", compress="none") save.image(file="lung.integration.v4.final.RData", compress = F) #### # # # DO KIDNEY ANALYSIS # # #### objlist.kidney <- lapply(X = kidneySampleList, FUN = function(x) { DefaultAssay(x) = "RNA" x }) objlist.kidney <- lapply(X = objlist.kidney, FUN = function(x) { DefaultAssay(x) = "RNA" print(x) x <- CellCycleScoring( x, g2m.features = g2m.genes, s.features = s.genes) #"nCount_RNA","S.Score", "G2M.Score" suppressWarnings(x <- SCTransform(x,vars.to.regress = c('percent.mt', 'percent.rp'), verbose = T)) }) features <- SelectIntegrationFeatures(object.list = objlist.kidney, nfeatures = 2000) print(length(features)) objlist.kidney <- PrepSCTIntegration(object.list = objlist.kidney, anchor.features = features) objlist.kidney.anchors <- FindIntegrationAnchors(object.list = objlist.kidney, normalization.method = "SCT", anchor.features = features) obj.kidney.integrated <- IntegrateData(anchorset = objlist.kidney.anchors, normalization.method = "SCT",verbose=T) print(obj.kidney.integrated) DefaultAssay(obj.kidney.integrated) = "integrated" #obj.lung.integrated <- ScaleData(obj.lung.integrated, verbose = FALSE) obj.kidney.integrated <- RunPCA(obj.kidney.integrated, verbose = FALSE) obj.kidney.integrated <- FindNeighbors(obj.kidney.integrated, dims = 1:50) obj.kidney.integrated <- FindClusters(obj.kidney.integrated, resolution = 0.2) obj.kidney.integrated$idents = Idents(obj.kidney.integrated) obj.kidney.integrated <- RunUMAP(obj.kidney.integrated, dims = 1:50, n.neighbors = 50, return.model = TRUE) DefaultAssay(obj.kidney.integrated) = "integrated" obj.kidney.integrated <- FindClusters(obj.kidney.integrated, resolution = 0.5) obj.kidney.integrated$idents = Idents(obj.kidney.integrated) p=DimPlot(obj.kidney.integrated, pt.size = 0.001, label=T) save_plot(p, "plots_kidney/dimplot_umap", fig.width=12, fig.height=8) p=VlnPlot(obj.kidney.integrated, c("percent.mt"), pt.size=0, group.by="idents") save_plot(p, "plots_kidney/vln_percentmt", fig.width=12, fig.height=5) p=FeaturePlot(obj.kidney.integrated, c("percent.mt"), pt.size=0, label=T) save_plot(p, "plots_kidney/fplt_percentmt", fig.width=12, fig.height=8) DefaultAssay(obj.kidney.integrated) <- "RNA" obj.kidney.integrated <- NormalizeData(obj.kidney.integrated, verbose = FALSE) p=FeaturePlot(obj.lung.integrated, c("Ccr2"), pt.size=0, label=T) save_plot(p, "plots_lung/fplt_ccr2", fig.width=12, fig.height=8) obj.kidney.integrated$conditionFeature = conditions.kidney obj.kidney.integrated$animalFeature = animals.kidney print(obj.kidney.integrated) p=DimPlot(obj.kidney.integrated, pt.size=0.001, label=T, split.by="conditionFeature") save_plot(p, "plots_kidney/dimplot_umap_by_cond", fig.width=20, fig.height=8) p=DimPlot(obj.kidney.integrated, pt.size=0.001, label=T, split.by="animalFeature") save_plot(p, "plots_kidney/dimplot_umap_by_animal", fig.width=40, fig.height=8) p=DimPlot(obj.kidney.integrated, pt.size=0.001, label=T, split.by="HTO_classification") save_plot(p, "plots_kidney/dimplot_umap_by_hto", fig.width=40, fig.height=8) write_rds(obj.kidney.integrated, "kidney.integrated.rds", compress="none") save.image(file="both.integration.v4.final.RData", compress = F) countByOIAbs = getCellCountDF(obj.kidney.integrated, prefix="", select_by = "idents", group_by="orig.ident", relative=F, show.percent=F, outname="plots_kidney/ia_count_by_origident.tsv") countByAFAbs = getCellCountDF(obj.kidney.integrated, prefix="", select_by = "idents", group_by="animalFeature", relative=F, show.percent=F, outname="plots_kidney/ia_count_by_animal.tsv") countByCFAbs = getCellCountDF(obj.kidney.integrated, prefix="", select_by = "idents", group_by="conditionFeature", relative=F, show.percent=F, outname="plots_kidney/ia_count_by_condition.tsv") DefaultAssay(obj.kidney.integrated) = "RNA" deResTT.kidney = makeDEResults(obj.kidney.integrated, assay="RNA", test="t") exprdfTT.kidney = getDEXpressionDF(obj.kidney.integrated, deResTT.kidney, assay="RNA") write.table(exprdfTT.kidney, "plots_kidney/expr_test_t.tsv", sep="\t", row.names=F, quote = F) write_xlsx(exprdfTT.kidney, "plots_kidney/expr_test_t.xlsx") markers.use.tt.kidney=subset(exprdfTT.kidney ,avg_log2FC>0&p_val_adj<0.01&!startsWith(gene, "mt-")&!startsWith(gene, "RP")) finalMarkers.use.tt.kidney = markers.use.tt.kidney %>% arrange(p_val_adj, desc(abs(pct.1)*abs(avg_log2FC),)) %>% group_by(clusterID) %>% dplyr::slice(1:20) finalMarkers.use.tt.kidney finalMarkers.all.tt.kidney = markers.use.tt.kidney %>% arrange(p_val_adj, desc(abs(pct.1)*abs(avg_log2FC),)) %>% group_by(clusterID) finalMarkers.all.tt.kidney = finalMarkers.all.tt.kidney[order(finalMarkers.all.tt$clusterID),] finalMarkers.all.tt.kidney p_dp_genes_idents = DotPlot(obj.kidney.integrated, features = unique(finalMarkers.use.tt.kidney$gene), assay="RNA")+ theme(axis.text.x = element_text(angle = 90, hjust = 1.0, vjust = 0.5)) save_plot(p_dp_genes_idents, "plots_kidney/dotplot_idents_genes", 45, 8) cluster_mean_expr.kidney = get_population_expression_data(obj.kidney.integrated, group="idents", outname="plots_kidney/cluster_mean_expr_kidney") cells.kidney.cre_pos = names(obj.kidney.integrated$conditionFeature[obj.kidney.integrated$conditionFeature == "cre+"]) cells.kidney.cre_neg = names(obj.kidney.integrated$conditionFeature[obj.kidney.integrated$conditionFeature == "cre-"]) print(paste(length(cells.kidney.cre_pos), length(cells.kidney.cre_neg), length(obj.kidney.integrated$conditionFeature))) de_kidney_cre = compareCellsByCluster(obj.kidney.integrated, cells.kidney.cre_pos, cells.kidney.cre_neg, "kidney_cre_pos", "kidney_cre_neg", group.by="idents", test="MAST", outfolder="de_kidney", fcCutoff=0.1) de_kidney_cre_t = compareCellsByCluster(obj.kidney.integrated, cells.kidney.cre_pos, cells.kidney.cre_neg, "kidney_cre_pos", "kidney_cre_neg", group.by="idents", test="t", outfolder="de_kidney_t", fcCutoff=0.1) de_kidney_cre_all = compareClusters(obj.kidney.integrated, cells.kidney.cre_pos, cells.kidney.cre_neg, "kidney_cre_pos", "kidney_cre_neg", test="MAST", outfolder="de_kidney", fcCutoff=0.1) makeVolcanos(de_kidney_cre, "Volcano", "kidney_volcanos/volcanos", FCcutoff=0.1) makeVolcanos(de_kidney_cre_t, "Volcano", "kidney_volcanos_t/volcanos", FCcutoff=0.1) makeVolcanos(list("de_kidney_cre_all"=de_kidney_cre_all), "Volcano", "kidney_volcanos/volcanos_all", FCcutoff=0.1) write_rds(obj.kidney.integrated, "kidney.integrated.rds", compress="none") save.image(file="both.integration.v4.final.RData", compress = F) #load("both.integration.v4.final.RData") source("functions.R") # # ## ### Additional Plots Lung ## # # obj.lung.integrated$conditionFeature2 = factor(obj.lung.integrated$conditionFeature , levels=c("cre-", "cre+") ) for (clusterID in unique(obj.lung.integrated$idents)) { print(clusterID) clusterCells = names(obj.lung.integrated$idents[obj.lung.integrated$idents == clusterID]) print(length(clusterCells)) plotElems = list( "Lung"=list(cells=clusterCells, label="Lung") ) titleText = paste("Scaled Average Expression of selected genes (Cluster ", clusterID, ")", sep="") plotGenes = c("Jun", "Cst3", "Fkbp5", "Cebpb","Irf8","Lst1","Cx3cr1","Ccr5") p=makeSideBySideDotPlot(obj.lung.integrated, plotElems, featureGenes = plotGenes, group.by = "conditionFeature2", col.min = -2, col.max = 2, title=titleText, scaled=T, featureWiseScale=F, cols = c("blue","yellow", "red")) save_plot(p, paste("dotplots_lung/dotplot_cre_updown", clusterID, sep="."), fig.width=16, fig.height=8) p=makeSideBySideDotPlot(obj.lung.integrated, plotElems, featureGenes = plotGenes, group.by = "conditionFeature2", col.min = 0, col.max = 30, title=titleText, scaled=F, featureWiseScale=F, cols = c("blue","yellow", "red")) save_plot(p, paste("dotplots_lung/dotplot_abs_cre_updown", clusterID, sep="."), fig.width=16, fig.height=8) plots=VlnPlot(subset(obj.lung.integrated, cells=clusterCells), features=plotGenes, ncol=2, group.by="idents", split.by="conditionFeature2", pt.size=0, combine = FALSE, assay="RNA") for(i in 1:length(plots)) { plots[[i]] <- plots[[i]] + geom_boxplot(color="grey", alpha=0.4, position =position_dodge(width = 0.9)) + stat_summary(fun=mean, geom="point", aes(group=split), position=position_dodge(.9), color="black", size=4) } plots=CombinePlots(plots) save_plot(plots, paste("dotplots_lung/vlnplot_cre_updown", clusterID, sep="."), fig.width=10, fig.height=12) } # # ## ### Additional Plots Kidney ## # # obj.kidney.integrated$conditionFeature2 = factor(obj.kidney.integrated$conditionFeature , levels=c("cre-", "cre+") ) for (clusterID in unique(obj.kidney.integrated$idents)) { print(clusterID) clusterCells = names(obj.kidney.integrated$idents[obj.kidney.integrated$idents == clusterID]) print(length(clusterCells)) plotElems = list( "Kidney"=list(cells=clusterCells, label="Kidney") ) titleText = paste("Scaled Average Expression of selected genes (Cluster ", clusterID, ")", sep="") plotGenes = c("Ddit4", "Ptprc", "Atf3", "Dusp1","Fos","Ccr5","Cx3cr1","Syk") p=makeSideBySideDotPlot(obj.kidney.integrated, plotElems, featureGenes = plotGenes, group.by = "conditionFeature2", col.min = -2, col.max = 2, title=titleText, scaled=T, featureWiseScale=F, cols = c("blue","yellow", "red")) save_plot(p, paste("dotplots_kidney/dotplot_cre_updown", clusterID, sep="."), fig.width=16, fig.height=8) p=makeSideBySideDotPlot(obj.kidney.integrated, plotElems, featureGenes = plotGenes, group.by = "conditionFeature2", col.min = 0, col.max = 30, title=titleText, scaled=F, featureWiseScale=F, cols = c("blue","yellow", "red")) save_plot(p, paste("dotplots_kidney/dotplot_abs_cre_updown", clusterID, sep="."), fig.width=16, fig.height=8) plots=VlnPlot(subset(obj.kidney.integrated, cells=clusterCells), features=plotGenes, ncol=2, group.by="idents", split.by="conditionFeature2", pt.size=0, combine = FALSE, assay="RNA") for(i in 1:length(plots)) { plots[[i]] <- plots[[i]] + geom_boxplot(color="grey", alpha=0.4, position =position_dodge(width = 0.9)) + stat_summary(fun=mean, geom="point", aes(group=split), position=position_dodge(.9), color="black", size=4) } plots=CombinePlots(plots) save_plot(plots, paste("dotplots_kidney/vlnplot_cre_updown", clusterID, sep="."), fig.width=10, fig.height=12) } # ## ### Reviewer plots ## # #Könntest Du einen Violin, einen Dot Plot und einen Feature Plot für folgende Marker auf alle Zellen in Kidney und Lunge machen? #Adgre1 (F4/80) #CX3CR1 #CD11c -> Itgax #CD16 -> Fcgr3, Fcgr4 #CD11b -> Itgam load("both.integration.v4.final.RData") source("functions.R") makeReviewerPlots = function(obj.in, outfolder) { DefaultAssay(obj.in) = "RNA" p=DimPlot(obj.in, group.by="idents") save_plot(p, paste(outfolder, "dimplot", sep="/"), fig.width=10, fig.height=8) gois = c("Cx3cr1", "Itgax", "Fcgr1", "Fcgr3", "Fcgr4", "Itgam", "Adgre1") objGenes = rownames(obj.in) print(setdiff(gois, objGenes)) p=VlnPlot(obj.in, gois, group.by="idents", ncol=1) save_plot(p, paste(outfolder, "vlnplot", sep="/"), fig.width=6, fig.height=14) p=DotPlot(obj.in, features=gois, group.by="idents") save_plot(p, paste(outfolder, "dotplot", sep="/"), fig.width=10, fig.height=4) for (gene in gois) { print(gene) p=FeaturePlot(obj.in, gene) save_plot(p, paste(outfolder, paste("featureplot", gene, sep="_"), sep="/"), fig.width=10, fig.height=8) } } makeReviewerPlots(obj.lung.integrated,"reviewer/lung") makeReviewerPlots(obj.kidney.integrated,"reviewer/kidney") lung.loom <- as.loom(obj.lung.integrated, filename = "lung.loom", verbose = FALSE) library(SeuratDisk) DefaultAssay(obj.lung.integrated) <- "RNA" SaveH5Seurat(obj.lung.integrated, "lung.h5seurat") Convert("lung.h5seurat", dest="h5ad") # or set assay="RNA". library(monocle3) library("SeuratWrappers") mc3_lung <- as.cell_data_set(obj.lung.integrated) mc3_lung <- preprocess_cds(mc3_lung, num_dim = 50) mc3_lung <- align_cds(mc3_lung, alignment_group = "animalFeature") mc3_lung <- reduce_dimension(mc3_lung, reduction_method="UMAP", preprocess_method="PCA") p=plot_cells(mc3_lung, label_groups_by_cluster=FALSE,color_cells_by = "idents") save_plot(p, "monocle/dim_reduce", fig.width=8, fig.height=6) mc3_lung <- cluster_cells(cds = mc3_lung, reduction_method = "UMAP") p=plot_cells(cds, color_cells_by = "partition") save_plot(p, "monocle/partitions", fig.width=8, fig.height=6) mc3_lung <- learn_graph(mc3_lung) p=plot_cells(mc3_lung, color_cells_by = "partition", label_groups_by_cluster=T, label_leaves=T, label_branch_points=T) save_plot(p, "monocle/partitions_graph", fig.width=8, fig.height=6) get_earliest_principal_node <- function(cds, time_bin="130-170"){ cell_ids <- which(colData(cds)[, "embryo.time.bin"] == time_bin) closest_vertex <- cds@principal_graph_aux[["UMAP"]]$pr_graph_cell_proj_closest_vertex closest_vertex <- as.matrix(closest_vertex[colnames(cds), ]) root_pr_nodes <- igraph::V(principal_graph(cds)[["UMAP"]])$name[as.numeric(names (which.max(table(closest_vertex[cell_ids,]))))] root_pr_nodes } mc3_lung <- order_cells(mc3_lung, root_pr_nodes=get_earliest_principal_node(mc3_lung)) p = plot_cells(mc3_lung, color_cells_by = "partition", group_cells_by = "partition", group_label_size = 4) save_plot(p, "monocle/clusters", fig.width=8, fig.height=6) mc3_lung <- learn_graph(mc3_lung, use_partition = TRUE) p= plot_cells(mc3_lung, color_cells_by = "pseudotime", label_groups_by_cluster = FALSE, group_label_size = 3.5, graph_label_size = 2) + scale_color_d3(palette = "category20b") save_plot(p, "monocle/clusters_cols", fig.width=8, fig.height=6) # a helper function to identify the root principal points: get_earliest_principal_node <- function(cds, time_bin="130-170"){ cell_ids <- which(colData(cds)[, "embryo.time.bin"] == time_bin) closest_vertex <- cds@principal_graph_aux[["UMAP"]]$pr_graph_cell_proj_closest_vertex closest_vertex <- as.matrix(closest_vertex[colnames(cds), ]) root_pr_nodes <- igraph::V(principal_graph(cds)[["UMAP"]])$name[as.numeric(names (which.max(table(closest_vertex[cell_ids,]))))] root_pr_nodes } mc3_lung <- order_cells(mc3_lung, root_pr_nodes=get_earliest_principal_node(mc3_lung)) plot_cells(mc3_lung, color_cells_by = "pseudotime", label_cell_groups=FALSE, label_leaves=FALSE, label_branch_points=FALSE, graph_label_size=1.5)