
#.libPaths("/burg/tosches/users/eg3092/Seurat_5_2024/")
#setwd("/burg/tosches/users/eg3092_v2")


library(Seurat)
 
pleuro.0<- readRDS("clean_annotate_v2_sal.dev_mat_neur_TP73_iso_060825_v3.rds")

mouse.0   <- readRDS("FMH_Neuron_mouse_develop_clean_2000-_1CLR_072624_annotate_V2.rds")


Idents(mouse.0)<-"r_anno"
mouse.0<-subset(mouse.0, idents= c("Mouse Midbrain"), invert = T)



Idents(pleuro.0)<-"r_anno"
pleuro.0<-subset(pleuro.0, idents= c("Pleuro Midbrain"), invert = T)




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


Ss = "P"

Sm = "M"


pleuro.0<-UpdateSeuratObject(pleuro.0)

mouse.0<-UpdateSeuratObject(mouse.0)


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)

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

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


eggnog_s <- eggnog_salamander[eggnog_salamander[,1] %in% rownames(s_data),]
eggnog_m <- eggnog_mouse[eggnog_mouse[,1] %in% rownames(m_data),]


s_data <- s_data[order(match(rownames(s_data), eggnog_s[,1])),]
m_data <- m_data[order(match(rownames(m_data), eggnog_m[,1])),]

all.equal(rownames(s_data),as.character(eggnog_s[,1]))
all.equal(rownames(m_data),as.character(eggnog_m[,1]))






rownames(s_data) <- eggnog_s$eggnog
rownames(m_data) <- eggnog_m$eggnog 




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



mouse.meta$species<-"mouse"
pleuro.meta$species<-"pleuro"

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


mouse.1<-CreateSeuratObject(counts = m_data, meta.data = mouse.meta)


nDims=90


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


pleuro.1 <- SCTransform(pleuro.1, variable.features.n = 4000, vars.to.regress = c("nCount_RNA","Stage","percent.mt","G2M.Score","S.Score") )
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)

print("checkpoint_1")

mouse.1 <- SCTransform(mouse.1, variable.features.n = 4000, vars.to.regress =c("percent.mt", "CellCycle", "Age") )
mouse.1 <- RunPCA(mouse.1, features = VariableFeatures(object = mouse.1))
mouse.1 <- FindNeighbors(mouse.1, dims = 1:35)
mouse.1 <- FindClusters(mouse.1, resolution = 1)
mouse.1 <- RunUMAP(mouse.1, dims = 1:35)

print("checkpoint_2")

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

object.list <- list(pleuro.1, mouse.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>.25)


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



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

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

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


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


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

#shared.variable.genes
shared_variable_genes



int.features<-SelectIntegrationFeatures(object.list = object.list, nfeatures = 2000)
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:90)

print("checkpoint_3")

integrated_v4 <- RunPCA(integrated_v4,npcs = 100 )
integrated_v4 <- RunUMAP(integrated_v4, dims = 1:90)




#int.features
length(int.features)
int.features



gc()
DefaultAssay(integrated_v4)<-"integrated"

integrated_v4 <- FindNeighbors(integrated_v4, dims = 1:90)
integrated_v4 <- FindClusters(integrated_v4, resolution = 1)
integrated_v4 <- RunUMAP(integrated_v4, dims = 1:90)

saveRDS(integrated_v4, file = "seurat5_clean_Lamanno_Woych_tp73_iso_FORE_mat_V2R_ivf4000_f2000_SCT_OG_SVG_filter_d90_060925_clr1.rds")

#pdf(paste0("plots_seurat5_clean_Lamanno_Woych_tp73_iso_FORE_mat_V2R_ivf4000_f2000_SCT_OG_SVG_filter_d90_060925_clr1.pdf"), width=20, height=20)

DimPlot(integrated_v4, label = T, raster = F)+NoAxes()

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

DimPlot(integrated_v4, group.by = "r_anno", raster = F)+NoAxes()
DimPlot(integrated_v4, group.by = "TP73_type_allen", raster = F)+NoAxes()
DimPlot(integrated_v4, group.by = "integrate_ident", raster = F)+NoAxes()
DefaultAssay(integrated_v4)<-"SCT"
FeaturePlot(integrated_v4, features = c("TP73","EBF3","EBF2","STK32A","EMX2","EOMES","LHX1"))

#FeaturePlot(integrated_v4, features = c("TP73","EBF3","EBF2","STK32A","EMX2","EOMES","LHX1"))
dev.off()




