
#.libPaths("/Users/eli/Documents/Tosches_Lab/Seurat5/Seurat5/renv/")

#setwd("~/Documents/Tosches_Lab/Tosches_R/")

library(Seurat)




eggnog_salamander <- read.table("pleurodeles_eggnog_pruned.txt" , header=T)
eggnog_mouse <- read.table("mus_eggnog_pruned.txt" , header=T)



CR_pw <- readRDS("TP73_all_clean_broad_ident_MT_120624.rds")


DimPlot(CR_pw, label = T)

#FeaturePlot(CR_pw, features = c("NFIX","PCP4","ASIC4","EOMES"), label = T)

#CR_pw@meta.data$integrate_ident<-paste0(("pw-"),CR_pw@meta.data$seurat_clusters)

CR_pw@meta.data$integrate_ident<-Idents(CR_pw)



#DimPlot(CR_mouse_allen, label = T, group.by = "subclass") 
#FeaturePlot(CR_mouse_allen, features = c("Trp73","Ebf3","Eomes","Stk32a"))

CR_mouse_allen <- readRDS("mouse_allen_TP73_MT_merge_111124.rds")
DimPlot(CR_mouse_allen, label = T)





CR_mouse_allen@meta.data$integrate_ident<-Idents(CR_mouse_allen)




DimPlot(CR_pw, label = T, group.by = "integrate_ident")

DimPlot(CR_mouse_allen, label = T, group.by = "integrate_ident")

DimPlot(CR_mouse_allen, label = T, group.by = "broad_ident")
pleuro.0<-UpdateSeuratObject(CR_pw)

allen.0<-UpdateSeuratObject(CR_mouse_allen)


Idents(pleuro.0)<-"broad_ident"
Idents(allen.0)<-"broad_ident"


pleuro.0<-subset(pleuro.0, idents ="salamander MT", invert = T)

allen.0<-subset(allen.0, idents ="mouse MT", invert = T)




shared_genes<-intersect(eggnog_mouse$eggnog, eggnog_salamander$eggnog)

eggnog_mouse<-subset(eggnog_mouse, eggnog_mouse$eggnog %in% shared_genes)
eggnog_salamander<-subset(eggnog_salamander, eggnog_salamander$eggnog %in% shared_genes)

eggnog_allen<-eggnog_mouse


s_data<-pleuro.0[["RNA"]]$counts
s_data <- s_data[rownames(s_data) %in% eggnog_salamander[,1],]

allen_data<-allen.0[["RNA"]]$counts
## this used to say y_data[]... but not sure where y_wold come from. changed to allen_data
allen_data <- allen_data[rownames(allen_data) %in% eggnog_mouse[,1],]





eggnog_s <- eggnog_salamander[eggnog_salamander[,1] %in% rownames(s_data),]
eggnog_allen <- eggnog_mouse[eggnog_mouse[,1] %in% rownames(allen_data),]

s_data <- s_data[order(match(rownames(s_data), eggnog_s[,1])),]
allen_data <- allen_data[order(match(rownames(allen_data), eggnog_allen[,1])),]



all.equal(rownames(s_data),as.character(eggnog_s[,1]))
all.equal(rownames(allen_data),as.character(eggnog_allen[,1]))






rownames(s_data) <- eggnog_s$eggnog
rownames(allen_data) <- eggnog_allen$eggnog 






pleuro.meta<-pleuro.0@meta.data
allen.meta<-allen.0@meta.data




allen.meta$species<-"mouse"

pleuro.meta$species<-"pleuro"



allen.meta$source<-"allen"

pleuro.meta$source<-"pleuro"






pleuro.1<-CreateSeuratObject(counts = s_data, meta.data = pleuro.meta)


allen.1<-CreateSeuratObject(counts = allen_data, meta.data = allen.meta)



pleuro.1 <- SCTransform(pleuro.1, variable.features.n = 8000 )
pleuro.1 <- RunPCA(pleuro.1, features = VariableFeatures(object = pleuro.1))
pleuro.1 <- FindNeighbors(pleuro.1, dims = 1:35)
pleuro.1 <- FindClusters(pleuro.1, resolution = 1)
pleuro.1 <- RunUMAP(pleuro.1, dims = 1:35)



allen.1 <- SCTransform(allen.1, variable.features.n = 8000 )
allen.1 <- RunPCA(allen.1, features = VariableFeatures(object = allen.1))
allen.1 <- FindNeighbors(allen.1, dims = 1:35)
allen.1 <- FindClusters(allen.1, resolution = 1)
allen.1 <- RunUMAP(allen.1, dims = 1:35)




options(future.globals.maxSize = 60000 * 1024^2)

object.list <- list(pleuro.1, allen.1)




Idents(pleuro.1)<-"orig.ident"

pleuro.1.dp<-DotPlot(pleuro.1, features = VariableFeatures(pleuro.1))

pleuro.1.dp.clean<-subset(pleuro.1.dp$data, pleuro.1.dp$data$pct.exp>.1)


var.feat.pleuro.1<-pleuro.1.dp.clean$features.plot



allen.1@meta.data$conform<-"mouse"
Idents(allen.1)<-"conform"

allen.1.dp<-DotPlot(allen.1, features = VariableFeatures(allen.1))

allen.1.dp.clean<-subset(allen.1.dp$data, allen.1.dp$data$pct.exp>.1)


var.feat.allen.1<-allen.1.dp.clean$features.plot


shared_variable_genes<-intersect(var.feat.allen.1, var.feat.pleuro.1)

shared_variable_genes


int.features<-SelectIntegrationFeatures(object.list = object.list, nfeatures = 5000)
int.features<-subset(int.features, int.features%in%shared_variable_genes)
object.list <-  PrepSCTIntegration(object.list = object.list, anchor.features = int.features)
anchors <- FindIntegrationAnchors(object.list = object.list, normalization.method = "SCT", reduction="cca", anchor.features = int.features)
integrated_v4 <- IntegrateData(anchorset = anchors, normalization.method = "SCT", dims=1:50)





integrated_v4 <- RunPCA(integrated_v4, npcs = 100 )


integrated_v4 <- FindNeighbors(integrated_v4, dims = 1:50)
integrated_v4 <- FindClusters(integrated_v4, resolution = .1)

integrated_v4 <- RunUMAP(integrated_v4, dims = 1:50)



DimPlot(integrated_v4,   group.by = "species", raster = F)


DimPlot(integrated_v4,   group.by = "broad_ident", raster = F)

DimPlot(integrated_v4, group.by = "integrate_ident", raster = F, label = T)


DimPlot(integrated_v4, group.by = "seurat_clusters", raster = F, label = T)


saveRDS(integrated_v4,file="seurat5_clean_Allen_Woych_TP73_if5000_8000SVG_SCT_OG_d50_040125")

