


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


library(Seurat)


# Take zebrafish object developent only (in Zebrafish  nomenclature), convert to salamander using zebrafish eggnog.


pleuro.0<- readRDS("clean_annotate_v2_sal.dev_mat_neur_FORE_TP73_iso_061625_v3.rds")
PD.0<- readRDS("ZF_pandley_anno_MAT_ident_dev_021015.rds")





eggnog_salamander <- read.table("pleurodeles_eggnog_pruned.txt" , header=T)
eggnog_ZF<-readRDS("ZF_eggnog_fix_020525.rds")

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


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

eggnog_ZF<-subset(eggnog_ZF, eggnog_ZF$eggnog %in% shared_genes)

eggnog_PD<-eggnog_ZF

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




PD_data<-PD.0[["RNA"]]$counts
PD_data <- PD_data[rownames(PD_data) %in% eggnog_PD[,1],]





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


eggnog_PD <- eggnog_PD[eggnog_PD[,1] %in% rownames(PD_data),]

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


PD_data <- PD_data[order(match(rownames(PD_data), eggnog_PD[,1])),]




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


all.equal(rownames(PD_data),as.character(eggnog_PD[,1]))




rownames(s_data) <- eggnog_s$eggnog


rownames(PD_data) <- eggnog_PD$eggnog 




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


PD.meta<-PD.0@meta.data





PD.meta$species<-"zebrafish"

pleuro.meta$species<-"pleuro"



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


PD_mtx<-as.matrix(PD_data)

#PD_df<-as.data.frame(sapply(PD_mtx, as.numeric))



PD_mtx_agg<-aggregate(PD_mtx, list(row.names(PD_mtx)), sum)

#PD_mtx_agg

rownames(PD_mtx_agg)<-PD_mtx_agg$Group.1


PD_mtx_agg<-PD_mtx_agg[, colnames(PD_mtx_agg) != "Group.1"]


#PD_mtx_agg_clean <- PD_mtx_agg[ , unlist(lapply(PD_mtx_agg, is.numeric))]

PD_mtx_agg_clean<-PD_mtx_agg[colSums(PD_mtx_agg) != 0]

PD_data.2<-as.sparse(PD_mtx_agg_clean)

PD.1<-CreateSeuratObject(counts = PD_data.2, meta.data = PD.meta)



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)


PD.1 <- SCTransform(PD.1, variable.features.n = 4000, vars.to.regress ="Stage" )
PD.1 <- RunPCA(PD.1, features = VariableFeatures(object = PD.1))
PD.1 <- FindNeighbors(PD.1, dims = 1:35)
PD.1 <- FindClusters(PD.1, resolution = 1)
PD.1 <- RunUMAP(PD.1, dims = 1:35)


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

object.list <- list(pleuro.1, PD.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



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

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

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


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


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







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



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


gc()

#integrated_v4<-readRDS("seurat5_Pandley_Woych_int_eggnog_MAT_f6000_SCT_OG_exact_d90_020525.rds")

#DefaultAssay(integrated_v4)<-"integrated"

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


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



pdf(paste0("plots_seurat5_Pandley_Woych_FORE_int_eggnog_MAT_V2R_vf4000_if2000_SCT_OG_SVG_exact_d80_061625.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 = "integrate_ident", raster = F)+NoAxes()

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

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


# 
# 
# integrated_v4_subset<-subset(integrated_v4, idents=c("49","53","20","41"))
# 
# saveRDS(integrated_v4_subset, file = "subset_seurat5_Pandley_Woych_int_eggnog_MAT_f6000_SCT_OG_exact_d90_020525_clr_1.5.rds")
# 
