--- title: "pseudotime analysis" output: trajectories --- Constructing single-cell trajectories Order cells in pseudotime along a trajectory ```{r} library(monocle3) library(ggplot2) library(dplyr) library(Seurat) library(cowplot) library(patchwork) library(tidyverse) library(ggpubr) library(pheatmap) library(RColorBrewer) library(viridis) library(scales) library(remotes) library(SeuratWrappers) library(shiny) ``` ```{r} setwd('Z:/16-Our_Papers/In Preparation/MetabolicModeling-Miro1_Claudia&Alise/Figures/Figure 2/Originals') #saveRDS(combined.all2, file="combined.all_integrated_annotated.rds") combined.all2<-readRDS("combined.all_integrated_annotated.rds") ``` #Double check the object ```{r} combined.all2$CELL_ID <-Idents(combined.all2) Idents(combined.all2) = combined.all2$CELL_ID # Visualization p <- DimPlot(combined.all2, reduction = "umap", label = TRUE) p1 <- DimPlot(combined.all2, reduction = "umap", split.by="Genotype",label =TRUE, pt.size = 1.5) p2 <- DimPlot(combined.all2, reduction = "umap", label = TRUE, pt.size = 1.5 ) p3 <- DimPlot(combined.all2, reduction = "umap", group.by = "Day",label = TRUE) p4 <- DimPlot(combined.all2, reduction = "umap", group.by = "Genotype",pt.size = 1.5) p5 <- DimPlot(combined.all2, reduction = "umap", split.by = "Genotype",pt.size = 1.5) p6 <- DimPlot(combined.all2, reduction = "umap", split.by = "Day",pt.size = 1.5) print(p) print(p1) print(p2) print(p3) print(p4) print(p5) print(p6) ``` #Subset the object removing GC samples ```{r} #subset WT and PD Idents(combined.all2) = combined.all2$Genotype WTPDobj <- subset(combined.all2, idents = c("WT_D30","WT_D60","PD_D30","PD_D60")) Idents(WTPDobj) = WTPDobj$CELL_ID ``` #Subset object to analyse only dopaminergic neuron clusters ```{r} levels(WTPDobj) DNobj <- subset(WTPDobj, idents = c("TH-high Dopaminergic neurons","Dopaminergic neurons")) levels(DNobj) ``` #create object B for the whole object or only DNs!!! The same script applies ```{r} my.counts = WTPDobj@assays$RNA@counts # <- change object here WTPDobj or DNobj my.counts= as.matrix(my.counts) cell_meta = data.frame("Cell" = colnames(my.counts)) rownames(cell_meta) = colnames(my.counts) gene_meta = data.frame("gene_short_name" = rownames(my.counts)) rownames(gene_meta) = rownames(my.counts) B = new_cell_data_set(my.counts, cell_metadata = cell_meta, gene_metadata = gene_meta) B <- as.cell_data_set(WTPDobj) # <- change object here WTPDobj or DNobj ``` ```{r} head(colData(B)) fData(B) rownames(fData(B))[1:10] head(counts(B)) ``` #keep metadata from seurat object for the pseudotime ```{r} recreate.partitions <- c(rep(1, length(B@colData@rownames))) names(recreate.partitions) <- B@colData@rownames recreate.partitions <- as.factor(recreate.partitions) recreate.partitions B@clusters@listData[["UMAP"]][["partitions"]] <- recreate.partitions list.cluster <- WTPDobj@active.ident # <- change object here WTPDobj or DNobj B@clusters@listData[["UMAP"]][["clusters"]] <- list.cluster B@int_colData@listData[["reducedDims"]]@listData[["UMAP"]] <- WTPDobj@reductions$umap@cell.embeddings # <- change object here WTPDobj or DNobj ``` ```{r} cluster.before.traj <-plot_cells(B, color_cells_by = "cluster", label_groups_by_cluster = F, group_label_size = 5) + theme(legend.position = "right") cluster.before.traj ``` ```{r} B <- learn_graph(B, use_partition = F) ``` ```{r} plot_cells(B, color_cells_by = "cluster",label_groups_by_cluster = T, label_branch_points = T, label_roots = T, label_leaves = T, group_label_size = 5) ``` #Order cells along the pseudotime and subset for DEG analysis ```{r} B <- order_cells(B, reduction_method = "UMAP") ps1<-plot_cells(B, color_cells_by = "pseudotime", label_groups_by_cluster = T, label_branch_points = F, label_roots = T, label_leaves = F, graph_label_size = 5) + theme( axis.line = element_line(colour = 'black', linewidth = 0.5),axis.title.y = element_text(size = 18),axis.text.y = element_text(size = 16), axis.title.x = element_text(size = 18),axis.text.x = element_text(size = 16),legend.text = element_text(size=16),legend.title = element_text(size=18) ) ``` #Visualize pseudotime for cell clusters in each Genotype ```{r} B$monocle3_pseudotime <- pseudotime(B) head(B$cluster) data.pseudo <- as.data.frame(colData(B)) dataWT<-subset(data.pseudo, data.pseudo$Genotype %in% c("WT_D30","WT_D60")) dataPD<-subset(data.pseudo, data.pseudo$Genotype %in% c("PD_D30","PD_D60")) #For all cell types psWT<-ggplot(dataWT, aes(monocle3_pseudotime, CELL_ID, fill = CELL_ID)) + scale_fill_manual(values=c("#CC79A7", "#66CC99", "lightskyblue","#9999CC", "#E69F00","coral3","black")) + geom_boxplot() + theme_bw() + theme(legend.position = "none", axis.title.y = element_text(size = 16),axis.text.y = element_text(size = 14), axis.title.x = element_text(size = 16),axis.text.x = element_text(size = 14),title = element_text(size=16) )+ labs (title = "Pseudotime_WT") psPD<-ggplot(dataPD, aes(monocle3_pseudotime, CELL_ID, fill = CELL_ID)) + scale_fill_manual(values=c("#CC79A7", "#66CC99", "lightskyblue","#9999CC", "#E69F00","coral3","black")) + geom_boxplot() + theme_bw() + theme(legend.position = "none", axis.title.y = element_text(size = 16),axis.text.y = element_text(size = 14), axis.title.x = element_text(size = 16),axis.text.x = element_text(size = 14),title = element_text(size=16) )+ labs (title = "Pseudotime_PD") print(psWT) print(psPD) ``` ```{r} #Subset dopaminergic neurons DNs<-subset(data.pseudo, data.pseudo$CELL_ID %in% c("Dopaminergic neurons","TH-high Dopaminergic neurons")) DNs$Genotype<-factor(DNs$Genotype, levels = c("PD_D60","PD_D30","WT_D60","WT_D30")) dopa<-ggplot(DNs, aes(monocle3_pseudotime, Genotype, fill = Genotype)) + geom_boxplot()+ geom_point() + scale_fill_manual(values=c("darkred","rosybrown","cornflowerblue", "lightskyblue")) + theme_bw() + theme(legend.position = "none", axis.title.y = element_text(size = 20),axis.text.y = element_text(size = 18), axis.title.x = element_text(size = 18),axis.text.x = element_text(size = 16)) print(dopa) ``` #DEGs over pseudotime #For dopaminergic neurons, repeat the script with DNobj ```{r} B1<-subset(B, B@colData$Genotype %in% c("WT_D30","WT_D60")) B2<-subset(B, B@colData$Genotype %in% c("PD_D30","PD_D60")) #Find DEGs over pseudotime in WT samples degWT <- graph_test(B1, neighbor_graph = "principal_graph") degWTfiltr<-subset(degWT,degWT$status == "OK") degWTfiltrTOP<-subset(degWTfiltr, degWTfiltr$q_value < 0.05) #write.csv(degWTfiltrTOP,"pseudotime_genes_WT30D60_DNs.csv") #Find DEGs over pseudotime in PD samples degPD <- graph_test(B2, neighbor_graph = "principal_graph") degPDfiltr<-subset(degPD,degPD$status == "OK") degPDfiltrTOP<-subset(degPDfiltr, degPDfiltr$q_value < 0.05) #write.csv(degPDfiltrTOP,"pseudotime_genes_PD30D60_DNs.csv") ``` #Plot WT significant pseudotime genes across all samples in DNobj ```{r} dataWT<-read.csv("Z:/16-Our_Papers/In Preparation/MetabolicModeling-Miro1_Claudia&Alise/Figures/Figure 3/pseudotime_genes_WT30D60_DNs.csv") dataPD<-read.csv("Z:/16-Our_Papers/In Preparation/MetabolicModeling-Miro1_Claudia&Alise/Figures/Figure 3/pseudotime_genes_PD30D60_DNs.csv") #combine lists dataA<-full_join(dataWT,dataPD, by="X") Idents(DNobj) = DNobj$Genotype dataAll<-subset(DNobj, features= dataA$X) genes_exprA <- AverageExpression(object=dataAll)$RNA Z_Normalize <- function(x){return((x-mean(x))/(sd(x)))} NormalisedData_Zscore <- t(apply(as.matrix(genes_exprA),1,Z_Normalize)) phA<-pheatmap(NormalisedData_Zscore, colorRampPalette(rev(brewer.pal(n = 5, name = "RdYlBu")))(100), show_rownames = F, angle_col = 90, fontsize = 10, cellwidth = 20, cellheight = 0.08) print(phA) ``` #Subset and plot the top by q value for WT ```{r} #subset the top by q value for WT data_orderedWT <- dataWT[order(dataWT$q_value),] dataTopWT<-data_orderedWT[1:100,] #write.csv(dataTopWT,"TOP100_DNs_WT.csv") data2WT<-subset(DNobj, features= dataTopWT$X) genes_exprWT <- AverageExpression(object=data2WT)$RNA Z_Normalize <- function(x){return((x-mean(x))/(sd(x)))} NormalisedData_Zscore <- t(apply(as.matrix(genes_exprWT),1,Z_Normalize)) phWT<-pheatmap(NormalisedData_Zscore, colorRampPalette(rev(brewer.pal(n = 5, name = "RdYlBu")))(100), show_rownames = T, angle_col = 90, cutree_rows = 2, fontsize_row = 5, fontsize_col=12, cellwidth = 20, cellheight = 5.5) print(phWT) ``` #Subset and plot the top by q value for PD ```{r} #subset the top by q value for PD data_orderedPD <- dataPD[order(dataPD$q_value),] dataTopPD<-data_orderedPD[1:100,] #write.csv(dataTopPD,"TOP100_DNs_PD.csv") data2PD<-subset(DNobj, features= dataTopPD$X) genes_exprPD <- AverageExpression(object=data2PD)$RNA Z_Normalize <- function(x){return((x-mean(x))/(sd(x)))} NormalisedData_Zscore <- t(apply(as.matrix(genes_exprPD),1,Z_Normalize)) phPD<-pheatmap(NormalisedData_Zscore, colorRampPalette(rev(brewer.pal(n = 5, name = "RdYlBu")))(100), show_rownames = T, angle_col = 90, cutree_rows = 2, fontsize_row = 5, fontsize_col=12, cellwidth = 20, cellheight = 5.5) print(phPD) ``` #Find specific DEGs over pseudotime for WT and PD ```{r} WTspecific<- dplyr::setdiff(dataWT$X,dataPD$X) WTsp_subs<-subset(dataWT, dataWT$X %in% WTspecific) #write.csv(WTsp_subs,"specific_pseudoTgenes_DNs_WT.csv") PDspecific<- dplyr::setdiff(dataPD$X,dataWT$X) PDsp_subs<-subset(dataPD, dataPD$X %in% PDspecific) #write.csv(PDspecific,"specific_pseudoTgenes_DNs_PD.csv") ``` #Visulaise enriched molecular functions and reactome pathways of the TOP 100 genes of WT and PD pseudotime ```{r} MF_WT<-read.csv("Z:/16-Our_Papers/In Preparation/MetabolicModeling-Miro1_Claudia&Alise/Figures/Figure 3/Molec_function_WT.csv") MF_PD<-read.csv("Z:/16-Our_Papers/In Preparation/MetabolicModeling-Miro1_Claudia&Alise/Figures/Figure 3/Molec_function_PD.csv") RC_WT<-read.csv("Z:/16-Our_Papers/In Preparation/MetabolicModeling-Miro1_Claudia&Alise/Figures/Figure 3/Reactome_WT.csv") RC_PD<-read.csv("Z:/16-Our_Papers/In Preparation/MetabolicModeling-Miro1_Claudia&Alise/Figures/Figure 3/Reactome_PD.csv") ``` ```{r} Plot <- ggplot(data=MF_WT[1:5,], aes(x=FDR, y= reorder(term.description,FDR),fill = term.description))+ geom_bar(width=0.7,stat="identity") + scale_fill_manual(values=c("khaki","steelblue3","royalblue3","goldenrod2","#414487FF"))+ #coord_flip()+ theme_classic() + # scale_fill_viridis(discrete=TRUE, option = "cividis") + labs(y = " ", x="MF terms enrichment -Log10(FDR)")+ theme(axis.text.x = element_text(size=14), axis.title.x=element_text(size=16), axis.title.y=element_blank(), axis.text.y=element_text(size=16, colour = "black"), legend.position = "none")+ scale_y_discrete(labels = label_wrap(50)) print(Plot) ``` ```{r} Plot1 <- ggplot(data=RC_WT[1:5,], aes(x=FDR, y= reorder(term.description,FDR),fill = term.description))+ geom_bar(width=0.7,stat="identity") + scale_fill_manual(values=c("#414487FF","goldenrod2","steelblue2","royalblue3","khaki"))+ #coord_flip()+ theme_classic() + # scale_fill_viridis(discrete=TRUE, option = "cividis") + labs(y = " ", x="Reactome terms enrichment -Log10(FDR)")+ theme(axis.text.x = element_text(size=14), axis.title.x=element_text(size=14), axis.title.y=element_blank(), axis.text.y=element_text(size=16, colour = "black"), legend.position = "none")+ scale_y_discrete(labels = label_wrap(62)) print(Plot1) ``` ```{r} Plot2 <- ggplot(data=MF_PD[1:5,], aes(x=FDR, y= reorder(term.description,FDR),fill = term.description))+ geom_bar(width=0.7,stat="identity") + scale_fill_manual(values=c("#7e4e90ff","#eb8055ff","coral3", "#a65c85ff","#f9a242ff"))+ #coord_flip()+ "#7e4e90ff","#eb8055ff","coral3", "#a65c85ff","#f9a242ff" theme_classic() + # scale_fill_viridis(discrete=TRUE, option = "cividis") + labs(y = " ", x="MF terms enrichment -Log10(FDR)")+ theme(axis.text.x = element_text(size=14), axis.title.x=element_text(size=16), axis.title.y=element_blank(), axis.text.y=element_text(size=16, colour = "black"), legend.position = "none")+ scale_y_discrete(labels = label_wrap(50)) print(Plot2) ``` ```{r} Plot3 <- ggplot(data=RC_PD[1:5,], aes(x=FDR, y= reorder(term.description,FDR),fill = term.description))+ geom_bar(width=0.7,stat="identity") + scale_fill_manual(values=c("#414487FF","#7AD151FF","#22A884FF","#2A788EFF","#FDE725FF"))+ #coord_flip()+ theme_classic() + # scale_fill_viridis(discrete=TRUE, option = "cividis") + labs(y = " ", x="Reactome terms enrichment -Log10(FDR)")+ theme(axis.text.x = element_text(size=14), axis.title.x=element_text(size=14), axis.title.y=element_blank(), axis.text.y=element_text(size=16, colour = "black"), legend.position = "none")+ scale_y_discrete(labels = label_wrap(62)) print(Plot3) ``` #Plot genes of interest over time ```{r} Idents(DNobj) = (DNobj)$Genotype my_levels = c("WT_D30","WT_D60", "PD_D30","PD_D60") levels(DNobj)=my_levels c1=DNobj@meta.data %>% mutate(Models= factor(Genotype, levels=my_levels)) %>% mutate(expression = DNobj@assays$RNA@data['HK2',]) %>% ggplot(aes(x = Models, y = expression)) + geom_jitter(aes(color = expression )) + scale_color_viridis(discrete=FALSE, option = "viridis") + geom_smooth(method="loess",se=FALSE, aes(group=1), linewidth=2, color="darkred") + #stat_summary(fun = mean, geom = "line",aes(group = Genotype), position = position_dodge(width = 0.5))+ theme_classic()+ #facet_grid(~CELL_ID, labeller=label_wrap_gen(width = 2, multi_line = TRUE))+ stat_compare_means(comparisons = list(c("WT_D30", "WT_D60"),c("WT_D30", "PD_D30"),c("WT_D60", "PD_D60"),c("PD_D30", "PD_D60")),method = "wilcox.test", label = "p.signif") + # scale_fill_manual(values= c( "darkgreen", "red", "skyblue"), name = "Condition") + labs(title = "HK2", y="RNA expression level")+ theme(plot.title=element_text(size=22, face="bold", hjust=0, color = "black"), axis.text.x=element_text(size=18,color = "black", angle=90), axis.title.x=element_blank(), axis.title.y=element_text(size=20,color = "black"), axis.text.y=element_text(size=16, colour = "black"), legend.text = element_text(size=14, color = "black"), legend.title = element_text(size=15)) c1 c2=DNobj@meta.data %>% mutate(Models= factor(Genotype, levels=my_levels)) %>% mutate(expression = DNobj@assays$RNA@data['PKM',]) %>% ggplot(aes(x = Models, y = expression)) + geom_jitter(aes(color = expression )) + scale_color_viridis(discrete=FALSE, option = "viridis") + geom_smooth(method="loess",se=FALSE, aes(group=1), linewidth=2, color="darkred") + #stat_summary(fun = mean, geom = "line",aes(group = Genotype), position = position_dodge(width = 0.5))+ theme_classic()+ #facet_grid(~CELL_ID, labeller=label_wrap_gen(width = 2, multi_line = TRUE))+ stat_compare_means(comparisons = list(c("WT_D30", "WT_D60"),c("WT_D30", "PD_D30"),c("WT_D60", "PD_D60"),c("PD_D30", "PD_D60")),method = "wilcox.test", label = "p.signif") + # scale_fill_manual(values= c( "darkgreen", "red", "skyblue"), name = "Condition") + labs(title = "PKM", y="RNA expression level")+ theme(plot.title=element_text(size=22, face="bold", hjust=0, color = "black"), axis.text.x=element_text(size=18,color = "black", angle=90), axis.title.x=element_blank(), axis.title.y=element_text(size=20,color = "black"), axis.text.y=element_text(size=16, colour = "black"), legend.text = element_text(size=14, color = "black"), legend.title = element_text(size=15)) c2 c3=DNobj@meta.data %>% mutate(Models= factor(Genotype, levels=my_levels)) %>% mutate(expression = DNobj@assays$RNA@data['NDUFA1',]) %>% ggplot(aes(x = Models, y = expression)) + geom_jitter(aes(color = expression )) + scale_color_viridis(discrete=FALSE, option = "viridis") + geom_smooth(method="loess",se=FALSE, aes(group=1), linewidth=2, color="darkred") + #stat_summary(fun = mean, geom = "line",aes(group = Genotype), position = position_dodge(width = 0.5))+ theme_classic()+ #facet_grid(~CELL_ID, labeller=label_wrap_gen(width = 2, multi_line = TRUE))+ stat_compare_means(comparisons = list(c("WT_D30", "WT_D60"),c("WT_D30", "PD_D30"),c("WT_D60", "PD_D60"),c("PD_D30", "PD_D60")),method = "wilcox.test", label = "p.signif") + # scale_fill_manual(values= c( "darkgreen", "red", "skyblue"), name = "Condition") + labs(title = "NDUFA1", y="RNA expression level")+ theme(plot.title=element_text(size=22, face="bold", hjust=0, color = "black"), axis.text.x=element_text(size=18,color = "black", angle=90), axis.title.x=element_blank(), axis.title.y=element_text(size=20,color = "black"), axis.text.y=element_text(size=16, colour = "black"), legend.text = element_text(size=14, color = "black"), legend.title = element_text(size=15)) c3 c4=DNobj@meta.data %>% mutate(Models= factor(Genotype, levels=my_levels)) %>% mutate(expression = DNobj@assays$RNA@data['ATP5F1A',]) %>% ggplot(aes(x = Models, y = expression)) + geom_jitter(aes(color = expression )) + scale_color_viridis(discrete=FALSE, option = "viridis") + geom_smooth(method="loess",se=FALSE, aes(group=1), linewidth=2, color="darkred") + #stat_summary(fun = mean, geom = "line",aes(group = Genotype), position = position_dodge(width = 0.5))+ theme_classic()+ #facet_grid(~CELL_ID, labeller=label_wrap_gen(width = 2, multi_line = TRUE))+ stat_compare_means(comparisons = list(c("WT_D30", "WT_D60"),c("WT_D30", "PD_D30"),c("WT_D60", "PD_D60"),c("PD_D30", "PD_D60")),method = "wilcox.test", label = "p.signif") + # scale_fill_manual(values= c( "darkgreen", "red", "skyblue"), name = "Condition") + labs(title = "ATP5F1A", y="RNA expression level")+ theme(plot.title=element_text(size=22, face="bold", hjust=0, color = "black"), axis.text.x=element_text(size=18,color = "black", angle=90), axis.title.x=element_blank(), axis.title.y=element_text(size=20,color = "black"), axis.text.y=element_text(size=16, colour = "black"), legend.text = element_text(size=14, color = "black"), legend.title = element_text(size=15)) c4 c5=DNobj@meta.data %>% mutate(Models= factor(Genotype, levels=my_levels)) %>% mutate(expression = DNobj@assays$RNA@data['LDHA',]) %>% ggplot(aes(x = Models, y = expression)) + geom_jitter(aes(color = expression )) + scale_color_viridis(discrete=FALSE, option = "viridis") + geom_smooth(method="loess",se=FALSE, aes(group=1), linewidth=2, color="darkred") + #stat_summary(fun = mean, geom = "line",aes(group = Genotype), position = position_dodge(width = 0.5))+ theme_classic()+ #facet_grid(~CELL_ID, labeller=label_wrap_gen(width = 2, multi_line = TRUE))+ stat_compare_means(comparisons = list(c("WT_D30", "WT_D60"),c("WT_D30", "PD_D30"),c("WT_D60", "PD_D60"),c("PD_D30", "PD_D60")),method = "wilcox.test", label = "p.signif") + # scale_fill_manual(values= c( "darkgreen", "red", "skyblue"), name = "Condition") + labs(title = "LDHA", y="RNA expression level")+ theme(plot.title=element_text(size=22, face="bold", hjust=0, color = "black"), axis.text.x=element_text(size=18,color = "black", angle=90), axis.title.x=element_blank(), axis.title.y=element_text(size=20,color = "black"), axis.text.y=element_text(size=16, colour = "black"), legend.text = element_text(size=14, color = "black"), legend.title = element_text(size=15)) c5 ``` #Plot mitochondrial genes ```{r} mt.genes <- rownames(DNobj)[grep("^MT-",rownames(DNobj))] Idents(DNobj) = DNobj$Genotype my_levels = c("WT_D30","WT_D60","PD_D30","PD_D60") Idents(DNobj) <- factor(Idents(DNobj), levels= my_levels) mt<-DotPlot(DNobj, features = mt.genes, cols = c("lightgrey", "coral3"),dot.scale = 12) + RotatedAxis()+theme(axis.text.x = element_text(size=20, color="black", angle = 90),axis.text.y = element_text(size=13, color="black"), axis.title.x = element_blank(),axis.title.y = element_blank())+coord_flip() mt ``` #Perform differential expression analysis between dopaminergic neuron clusters WTD60vsWTD30, PDD60vs PDD30 and WTD60vsPDD60 ```{r} #DEGs TH-high WT D30 vs WT D60 Idents(DNobj) = DNobj$CELL_ID TH_high<-subset(DNobj,idents = "TH-high Dopaminergic neurons" ) Idents(TH_high) = TH_high$Genotype DEGS.THhigh.WTD60vsWTD30 <- FindMarkers(TH_high, ident.1 = "WT_D60", ident.2 = "WT_D30", only.pos = F,test.use = "DESeq2") #write.csv(DEGS.THhigh.WTD60vsWTD30, file = "Z:/16-Our_Papers/In #Preparation/MetabolicModeling-Miro1_Claudia&Alise/Figures/Figure 3/DEGS.THhigh.WTD60vsWTD30.csv") Idents(DNobj) = DNobj$CELL_ID DNs<-subset(DNobj,idents = "Dopaminergic neurons" ) Idents(DNs) = DNs$Genotype #DEGs DNs WT D30 vs WT D60 DEGS.DNs.WTD60vsWTD30 <- FindMarkers(DNs, ident.1 = "WT_D60", ident.2 = "WT_D30", only.pos = F,test.use = "DESeq2") #write.csv(DEGS.DNs.WTD60vsWTD30, file = "Z:/16-Our_Papers/In #Preparation/MetabolicModeling-Miro1_Claudia&Alise/Figures/Figure 3/DEGS.DNs.WTD60vsWTD30.csv") #DEGs DNs PD D30 vs PD D60 DEGS.DNs.PDD60vsPDD30 <- FindMarkers(DNs, ident.1 = "PD_D60", ident.2 = "PD_D30", only.pos = F,test.use = "DESeq2") #write.csv(DEGS.DNs.PDD60vsPDD30, file = "Z:/16-Our_Papers/In #Preparation/MetabolicModeling-Miro1_Claudia&Alise/Figures/Figure 3/DEGS.DNs.PDD60vsPDD30.csv") #compare PD D60 vs WT D60 Idents(DNobj) = DNobj$Genotype DEGS.DNS_all.PDD60vsWTD60 <- FindMarkers(DNobj, ident.1 = "PD_D60", ident.2 = "WT_D60", only.pos = F,test.use = "DESeq2") #write.csv(DEGS.DNS_all.PDD60vsWTD60, file = "Z:/16-Our_Papers/In #Preparation/MetabolicModeling-Miro1_Claudia&Alise/Figures/Figure 3/DEGS.DNS_all.PDD60vsWTD60.csv") ``` ```{r} #library(EnhancedVolcano) df <- as.data.frame(DEGS.DNS_all.PDD60vsWTD60) df_sig <- subset(df, df$p_val<0.05) df1<-as.data.frame(df_sig$avg_log2FC) df2 <- as.data.frame(df_sig$p_val_adj) df3<-as.data.frame(df_sig$p_val) genes <- rownames(df_sig) dataVolc <- cbind.data.frame(genes,df1,df2,df3) colnames(dataVolc) <- c("genes","log2FoldChange", "pvalue_adj", "pvalue") DEGS.DNS_all.PDD60vsWTD60 <- EnhancedVolcano(dataVolc, lab = dataVolc$genes, x = 'log2FoldChange', y = 'pvalue_adj' , xlim = c(-2, 2), ylim = c(0, 150), #selectLab = "SNCA", title = ' ', subtitle = NULL, FCcutoff = 1, pCutoff = 0.05, pointSize =2, colAlpha = 1, legendLabels = c("NS","log2 FC > 1","p.adj. < 0.05","p.adj. < 0.05 & log2 FC > 1"), legendLabSize = 10, drawConnectors = TRUE, widthConnectors = 0.3, lengthConnectors = unit(0.01,'npc'), legendPosition = "none", labSize = 3) plot(DEGS.DNS_all.PDD60vsWTD60) ``` #Top 10 biological processes from DEGs between WT D60 and PD D60 ```{r} BP_PDvsWT<-read.csv("Z:/16-Our_Papers/In Preparation/MetabolicModeling-Miro1_Claudia&Alise/Figures/Figure 3/PDvsWT_D60_DNS.csv") library(paletteer) Plot5 <- ggplot(data=BP_PDvsWT[1:10,], aes(x=FDR, y= reorder(term.description,FDR),fill = FDR))+ geom_bar(width=0.7,stat="identity") + scale_fill_gradient2( low = "mistyrose", mid = "rosybrown3", high = "indianred", midpoint = 7.79 )+ #coord_flip()+ theme_classic() + # scale_fill_viridis(discrete=TRUE, option = "cividis") + labs(y = " ", x="BP terms enrichment -Log10(FDR)")+ theme(axis.text.x = element_text(size=14), axis.title.x=element_text(size=15), axis.title.y=element_blank(), axis.text.y=element_text(size=15, colour = "black"), legend.position = "none")+ scale_y_discrete(labels = label_wrap(40)) print(Plot5) ```