---
title: "script_pairwise"
format: html
editor: visual
---

## Pairwise

Script to compute the distance between pairs of samples of different aridity condition

## Running Code

You can add options to executable code like this

```{r}

library(tidyverse)
library(readxl)
library(paletteer)
library(pheatmap)
library(MuSiC)
library(SingleCellExperiment)
library("FactoMineR")
library("factoextra")
library(wrMisc)
library(DESeq2)

options(future.globals.maxSize = 15 * 1e9)  # allow R to use 256 GB of RAM

#knitr::opts_chunk$set(root.dir = normalizePath("/home/rstudio/ifbdata/mydatalocal/results/pairwise"))  
```

```{r}



# read expression data and metadata
# count_table
Alltable_onl = read.csv('../data/datalocalandonline_cds_counts_2025_final.tsv', h=T, sep = "\t", row.names=1) #15796   125
colnames(Alltable_onl) = gsub("loc_cds_SHPC43_Mastomys_kollmanspergeri", "loc_cds_SHPC43_Mastomys_kollmannspergeri", colnames(Alltable_onl))
colnames(Alltable_onl) = gsub("loc_cds_SHPC74_Mastomys_kollmanspergeri", "loc_cds_SHPC74_Mastomys_kollmannspergeri", colnames(Alltable_onl))
colnames(Alltable_onl) = gsub("loc_cds_SHPC75_Mastomys_kollmanspergeri", "loc_cds_SHPC75_Mastomys_kollmannspergeri", colnames(Alltable_onl))
colnames(Alltable_onl) = gsub("loc_cds_SHPC76_Mastomys_kollmanspergeri", "loc_cds_SHPC76_Mastomys_kollmannspergeri", colnames(Alltable_onl))
colnames(Alltable_onl) = gsub("SHPC46_Arvicanthis_niloticus_lmo", "SHPC46_Mastomys_natalensis", colnames(Alltable_onl)) #change Ani_lmo to Mna even if not using in expression)
colnames(Alltable_onl) = gsub("onl5_cds_SRR17216323_Peromyscus_maniculatus", "onl_cds_SRR17216323_Peromyscus_maniculatus", colnames(Alltable_onl))
colnames(Alltable_onl) = gsub("onl5_cds_SRR17216333_Peromyscus_maniculatus", "onl_cds_SRR17216333_Peromyscus_maniculatus", colnames(Alltable_onl))
colnames(Alltable_onl) = gsub("onl5_cds_SRR17216352_Peromyscus_maniculatus", "onl_cds_SRR17216352_Peromyscus_maniculatus", colnames(Alltable_onl))
# coldata


coldata <- read.csv("../data/coldata_July2025_family.tsv", sep = "\t")
coldata_tot <- coldata[which(coldata$total_set == "yes"),]
coldata_mur <- coldata[which(coldata$muridae_set == "yes"),]
coldata_rec <- coldata[which(coldata$recent_set == "yes"),]
coldata_anc <- coldata[which(coldata$ancient_set == "yes"),]

# DE genes

Demur = read.csv('../DESeq/Table_DE_results_cond_season_murinae_deseq.tsv', h=T, sep = "\t", row.names=1) #15796   125
Demur$gene=row.names(Demur)
Demur = Demur[Demur$sig=="FDR<0.1",]

Detot = read.csv('../DESeq/Table_DE_results_cond_season_total_deseq.tsv', h=T, sep = "\t", row.names=1) #15796   125
Detot$gene=row.names(Detot)
Detot = Detot[Detot$sig=="FDR<0.1",]


Deanc = read.csv('../DESeq/Table_DE_results_cond_season_ancient_deseq.tsv', h=T, sep = "\t", row.names=1) #15796   125
Deanc$gene=row.names(Deanc)
Deanc = Deanc[Deanc$sig=="FDR<0.1",]



Derec = read.csv('../DESeq/Table_DE_results_cond_season_recent_deseq.tsv', h=T, sep = "\t", row.names=1) #15796   125
Derec$gene=row.names(Derec)
Derec = Derec[Derec$sig=="FDR<0.1",]

# results from deconvolution

prop= read.table("../Deconvolutions/results/ref_prop_music_total_rat.txt",sep="\t",h=T,row.names=1)
prop$IDshort = prop$ID
prop$ID = row.names(prop)

```

All pairs, union of DE genes

```{r}

unionDE = unique(c(Demur$gene, Deanc$gene, Derec$gene, Detot$gene))


############################ DATASET TOT (80 samples)
# Prepare dataset for coldata_tot
coldata=coldata_tot
All_table = Alltable_onl[,colnames(Alltable_onl) %in% coldata$ID]
Alltableok1 <- round(All_table[complete.cases(All_table),])
Alltableok1[is.na(Alltableok1)]=0
n=apply(Alltableok1,1,function(x){sum(x==0)})
Alltableok1=Alltableok1[n==0,]
# Prepare coldata
coldata <- coldata[ which(coldata$ID %in% colnames(Alltableok1)), ]
coldata <- coldata[ order(match(coldata$ID, colnames(Alltableok1))), ]#Reorder table such coldata and matrix has same order
coldata$Batch_number <- as.factor(coldata$Batch_number)
coldata$cond_season = factor(coldata$cond_season, level = c("arid","mesic"))

combExprMat = as.matrix(Alltableok1)
colnames(combExprMat) = coldata$species
speciesLabels= coldata$species
row.names(coldata_tot)=coldata_tot$ID
coldata_tot=coldata_tot[names(Alltableok1),]

# normalize the data
ddsInput <- DESeqDataSetFromMatrix(countData = as.matrix(round(Alltableok1)),
                                   colData = coldata_tot,
                                   design = ~ 1)
dds <- DESeq(ddsInput)

rlg <- vst(dds, blind=TRUE)
ntd <- normTransform(dds)

#se <- SummarizedExperiment(log2(counts(dds, normalized=TRUE) + 1),
#                           colData=colData(dds))
resBasemean <- counts(dds, norm=T)
names(resBasemean) = coldata_tot$species

counts_sp_norm <- as.data.frame(counts(dds, norm=T))

countsDE = counts_sp_norm %>% 
  subset(select=coldata_tot$ID)  %>% 
  filter(rownames(counts_sp_norm) %in% unionDE) %>% 
  rowGrpMeans(grp = coldata_tot$species, na.rm = TRUE)
  

coldata_exp_sp <- coldata_tot%>%
  select("species","family","cond_season")%>%
  distinct()

distPairs <- as.data.frame(cor(countsDE, method="pearson"))

distPairs$Species1 <- rownames(distPairs)
distPairs_long_XM <- distPairs%>%
  pivot_longer(-Species1,names_to="Species2",values_to="cor")%>%
  inner_join(coldata_exp_sp, join_by(Species1 == species))%>%
  inner_join(coldata_exp_sp, join_by(Species2 == species))%>%
  filter(Species1!=Species2)%>%
  filter(cond_season.x!=cond_season.y)%>%
  mutate(same_family = family.x==family.y)

# select members of pair that belong to same family / different family 


distPairs_long_XM %>% 
  ggplot(aes(x=same_family, y=cor)) + geom_boxplot() + geom_jitter(shape=16, position=position_jitter(0.2))



distPairs_long_all <- distPairs%>%
  pivot_longer(-Species1,names_to="Species2",values_to="cor")%>%
  inner_join(coldata_exp_sp, join_by(Species1 == species))%>%
  inner_join(coldata_exp_sp, join_by(Species2 == species))%>%
  filter(Species1!=Species2)%>%
  mutate(same_cond = paste(cond_season.x,cond_season.y))%>%
  mutate(same_family = family.x==family.y)%>%
  mutate(intra = Species1 < Species2)%>%
  filter(intra==TRUE)%>%
  mutate(condition = case_when(cond_season.x=="mesic"&cond_season.y=="arid"~"mesic arid",
                   cond_season.x=="arid"&cond_season.y=="mesic"~"mesic arid",
                   cond_season.x=="arid"&cond_season.y=="arid"~"arid arid",
                   cond_season.x=="mesic"&cond_season.y=="mesic"~"mesic mesic"))%>%
  separate(Species1,c("genera1","sp1"), sep = "_",remove=FALSE)%>%
  separate(Species2,c("genera2","sp2"), sep = "_",remove=FALSE)%>%
  mutate(same_genus = genera1==genera2)

# select members of pair that belong to same family / different family 


p <- distPairs_long_all %>% 
  ggplot(aes(x=same_family, y=cor,fill=condition)) + geom_boxplot(outlier.shape = NA)+ 
geom_point(position = position_jitterdodge(jitter.width = 0.1), alpha=0.3)+ 
  scale_fill_manual(values=c("orange", "grey","darkgreen" )) + theme_bw()

ggsave("total_pair_by_family_distpearson.pdf",p,width=20,height=20,units="cm")




p <- distPairs_long_all %>% 
  ggplot(aes(x=same_genus, y=cor,fill=condition)) + geom_boxplot(outlier.shape = NA)+ 
geom_point(position = position_jitterdodge(jitter.width = 0.1), alpha=0.3)+ 
  scale_fill_manual(values=c("orange", "grey","darkgreen")) + theme_bw()

ggsave("total_pair_by_genus_distpearson.pdf",p,width=20,height=20,units="cm")





p <- distPairs_long_all %>% 
  ggplot(aes(x=paste(same_genus,same_family), y=cor,fill=condition)) + geom_boxplot(outlier.shape = NA)+ 
geom_point(position = position_jitterdodge(jitter.width = 0.1), alpha=0.3)+ 
  scale_fill_manual(values=c("orange", "grey","darkgreen")) + theme_bw()


ggsave("total_pair_by_genus_by_family_distpearson.pdf",p,width=20,height=20,units="cm")



```

Annotate pairs, split DE genes

```{r}

# annotate non-convergent pairs :

toremove = list( 
  c("Fukomys_damarensis","Fukomys_mechowii"),
  c("Fukomys_damarensis","Fukomys_micklemi"),
  c("Fukomys_mechowii","Fukomys_micklemi"),
  c("Nannomys_cot","Nannomys_indutus"),
  c("Mus_musculus","Mus_musculus_ddo"),
  c("Mus_musculus","Mus_musculus_wla"),
  c("Mus_musculus_ddo","Mus_musculus_wla"),
  c("Mastomys_coucha","Mastomys_erythroleucus"),
  c("Mastomys_coucha","Mastomys_kollmannspergeri"),
  c("Mastomys_erythroleucus","Mastomys_kollmannspergeri"),
  c("Abrothrix_longipilis","Abrothrix_olivaceus")
)
toremove <- data.frame(t(data.frame(toremove)))
toremove$pair <- paste(toremove$X1,toremove$X2)


# annotate pairs from within/between datasets: 

age = list( 
  c("Peromyscus_eremicus","within"),
  c("Apodemus_mystacinus","within"),
  c("Nannomys_cot","within"),
  c("Nannomys_indutus","within"),
  c("Mus_macedonicus","within"),
  c("Chaetodipus_baileyi","between"),
  c("Micaelamys_namaquensis","between"),
  c("Meriones_unguiculatus","between"),
  c("Chinchilla_lanigera","between")
)
age <- data.frame(t(data.frame(age)))
names(age)=c("Species","age")


compute_cor=function(distaname="pearson",genesetname="unionDE",geneset=unionDE)
{
  countsDE = counts_sp_norm %>% 
  subset(select=coldata_tot$ID)  %>% 
  filter(rownames(counts_sp_norm) %in% geneset) %>% 
  rowGrpMeans(grp = coldata_tot$species, na.rm = TRUE)
  
distPairs <- as.data.frame(cor(countsDE, method=distaname))
distPairs$Species1 <- rownames(distPairs)


distPairs_long_all <- distPairs%>%
  pivot_longer(-Species1,names_to="Species2",values_to="cor")%>%
  inner_join(coldata_exp_sp, join_by(Species1 == species))%>%
  inner_join(coldata_exp_sp, join_by(Species2 == species))%>%
  filter(Species1!=Species2)%>%
  mutate(pair = paste(Species1,Species2))%>%
  anti_join(toremove, join_by(pair==pair))%>%
  mutate(same_cond = paste(cond_season.x,cond_season.y))%>%
  mutate(same_family = family.x==family.y)%>%
  mutate(intra = Species1 < Species2)%>%
  filter(intra==TRUE)%>%
  mutate(condition = case_when(cond_season.x=="mesic"&cond_season.y=="arid"~"mesic arid",
                   cond_season.x=="arid"&cond_season.y=="mesic"~"mesic arid",
                   cond_season.x=="arid"&cond_season.y=="arid"~"arid arid",
                   cond_season.x=="mesic"&cond_season.y=="mesic"~"mesic mesic"))%>%
  separate(Species1,c("genera1","sp1"), sep = "_",remove=FALSE)%>%
  separate(Species2,c("genera2","sp2"), sep = "_",remove=FALSE)%>%
  mutate(same_genus = genera1==genera2)

# select members of pair that belong to same family / different family 


p1 <- distPairs_long_all %>% 
  ggplot(aes(x=same_family, y=cor,fill=condition)) + geom_boxplot(outlier.shape = NA)+ 
geom_point(position = position_jitterdodge(jitter.width = 0.1), alpha=0.3)+ 
  scale_fill_manual(values=c("orange", "grey","darkgreen" )) + theme_bw()#+ylim(0,1)

name_graph <- paste0(paste("distance",distaname,"genes",genesetname,"allpairs_family",sep="_"),".pdf")
ggsave(name_graph,p,width=20,height=20,units="cm")

name_graph <- paste0(paste("distance",distaname,"genes",genesetname,"allpairs_family",sep="_"),".png")
ggsave(name_graph,p,width=20,height=20,units="cm")




p2 <- distPairs_long_all %>% 
  ggplot(aes(x=same_genus, y=cor,fill=condition)) + geom_boxplot(outlier.shape = NA)+ 
geom_point(position = position_jitterdodge(jitter.width = 0.1), alpha=0.3)+ 
  scale_fill_manual(values=c("orange", "grey","darkgreen" )) + theme_bw()#+ylim(0,1)

name_graph <- paste0(paste("distance",distaname,"genes",genesetname,"allpairs_genus",sep="_"),".pdf")
ggsave(name_graph,p,width=20,height=20,units="cm")

name_graph <- paste0(paste("distance",distaname,"genes",genesetname,"allpairs_genus",sep="_"),".png")
ggsave(name_graph,p,width=20,height=20,units="cm")




distPairs_long_subset <- distPairs_long_all%>%
   inner_join(age, join_by(Species1 == Species))%>%
   inner_join(age, join_by(Species2 == Species))%>%
  mutate(age_category = case_when(age.x=="between"&age.y=="between"~"between between",
                   age.x=="between"&age.y=="within"~"between within",
                   age.x=="within"&age.y=="between"~"between within",
                   age.x=="within"&age.y=="within"~"within within"))


p3 <- distPairs_long_subset %>% 
  ggplot(aes(x=age_category, y=cor,fill=same_family)) + geom_boxplot(outlier.shape = NA)+ 
geom_point(position = position_jitterdodge(jitter.width = 0.1), alpha=0.3)+ theme_bw()#+ylim(0,1)

name_graph <- paste0(paste("distance",distaname,"genes",genesetname,"agecategory",sep="_"),".pdf")
ggsave(name_graph,p,width=20,height=20,units="cm")

name_graph <- paste0(paste("distance",distaname,"genes",genesetname,"agecategory",sep="_"),".png")
ggsave(name_graph,p,width=20,height=20,units="cm")


return(list(p1,p2,p3))
}
  
  


unionDE = unique(c(Demur$gene, Deanc$gene, Derec$gene, Detot$gene))
murinaeexcluDE = Demur$gene[!Demur$gene%in%unique(c(Deanc$gene, Derec$gene, Detot$gene))]
ancientexcluDE = Deanc$gene[!Deanc$gene%in%unique(c(Demur$gene, Derec$gene, Detot$Gene))]


plots1<-compute_cor(distaname="pearson",genesetname="unionDE",geneset=unionDE)
plots2<-compute_cor(distaname="spearman",genesetname="unionDE",geneset=unionDE)
# compute_cor(distaname="spearman",genesetname="murinae",geneset=deMurinae$Gene)
# compute_cor(distaname="spearman",genesetname="ancient",geneset=deAncient$Gene)
# compute_cor(distaname="spearman",genesetname="recent",geneset=deRecent$Gene)
# compute_cor(distaname="spearman",genesetname="total",geneset=deTotal$Gene)
# compute_cor(distaname="spearman",genesetname="murinaeexcluDE",geneset=murinaeexcluDE)
# compute_cor(distaname="spearman",genesetname="ancientexcluDE",geneset=ancientexcluDE)
plots3<-compute_cor(distaname="spearman",genesetname="all",geneset=rownames(counts_sp_norm))

# compute_cor(distaname="pearson",genesetname="murinae",geneset=deMurinae$Gene)
# compute_cor(distaname="pearson",genesetname="ancient",geneset=deAncient$Gene)
# compute_cor(distaname="pearson",genesetname="recent",geneset=deRecent$Gene)
# compute_cor(distaname="pearson",genesetname="total",geneset=deTotal$Gene)
# compute_cor(distaname="pearson",genesetname="murinaeexcluDE",geneset=murinaeexcluDE)
# compute_cor(distaname="pearson",genesetname="ancientexcluDE",geneset=ancientexcluDE)
plots4<-compute_cor(distaname="pearson",genesetname="all",geneset=rownames(counts_sp_norm))




```

Pairwise cell proportions

```{r}



compute_cor_prop=function(distaname="pearson",listcell=c("Endo","IC","LOH","Prox_tub","Stroma"))
{
  prop_sp = prop %>% 
  select(all_of(listcell),species,cond_season,family) %>%
    group_by(species,cond_season,family)%>%
    summarize_all(mean)
  
  prop_sp=as.data.frame(prop_sp)
  row.names(prop_sp)=prop_sp$species

distPairs <- as.data.frame(cor(t(prop_sp[,listcell]), method=distaname))
distPairs$Species1 <- prop_sp$species

coldata_exp_sp <- coldata_tot%>%
  select("species","family","cond_season")%>%
  distinct()

distPairs_long_all <- distPairs%>%
  pivot_longer(-Species1,names_to="Species2",values_to="cor")%>%
  inner_join(coldata_exp_sp, join_by(Species1 == species))%>%
  inner_join(coldata_exp_sp, join_by(Species2 == species))%>%
  filter(Species1!=Species2)%>%
  mutate(pair = paste(Species1,Species2))%>%
  anti_join(toremove, join_by(pair==pair))%>%
  mutate(same_cond = paste(cond_season.x,cond_season.y))%>%
  mutate(same_family = family.x==family.y)%>%
  mutate(intra = Species1 < Species2)%>%
  filter(intra==TRUE)%>%
  mutate(condition = case_when(cond_season.x=="mesic"&cond_season.y=="arid"~"mesic arid",
                   cond_season.x=="arid"&cond_season.y=="mesic"~"mesic arid",
                   cond_season.x=="arid"&cond_season.y=="arid"~"arid arid",
                   cond_season.x=="mesic"&cond_season.y=="mesic"~"mesic mesic"))%>%
  separate(Species1,c("genera1","sp1"), sep = "_",remove=FALSE)%>%
  separate(Species2,c("genera2","sp2"), sep = "_",remove=FALSE)%>%
  mutate(same_genus = genera1==genera2)

# select members of pair that belong to same family / different family 

write.table(distPairs_long_all,file=paste0("distPairs_long_all",distaname,".txt"),sep="\t",quote=F)


p1 <- distPairs_long_all %>% 
  ggplot(aes(x=same_family, y=cor,fill=condition)) + geom_boxplot(outlier.shape = NA)+ 
geom_point(position = position_jitterdodge(jitter.width = 0.1), alpha=0.3)+ 
  scale_fill_manual(values=c("orange", "grey","darkgreen" )) + theme_bw()+ylim(0,1)

name_graph <- paste0(paste("distance",distaname,"prop_allpairs_family",sep="_"),".pdf")
ggsave(name_graph,p,width=20,height=20,units="cm")

name_graph <- paste0(paste("distance",distaname,"prop_allpairs_family",sep="_"),".png")
ggsave(name_graph,p,width=20,height=20,units="cm")




p2 <- distPairs_long_all %>% 
  ggplot(aes(x=same_genus, y=cor,fill=condition)) + geom_boxplot(outlier.shape = NA)+ 
geom_point(position = position_jitterdodge(jitter.width = 0.1), alpha=0.3)+ 
  scale_fill_manual(values=c("orange", "grey","darkgreen" )) + 
  theme_bw()+ylim(0,1)

name_graph <- paste0(paste("distance",distaname,"prop_allpairs_genus",sep="_"),".pdf")
ggsave(name_graph,p,width=20,height=20,units="cm")

name_graph <- paste0(paste("distance",distaname,"prop_allpairs_genus",sep="_"),".png")
ggsave(name_graph,p,width=20,height=20,units="cm")


distPairs_long_subset <- distPairs_long_all%>%
   inner_join(age, join_by(Species1 == Species))%>%
   inner_join(age, join_by(Species2 == Species))%>%
  mutate(age_category = case_when(age.x=="between"&age.y=="between"~"between between",
                   age.x=="between"&age.y=="within"~"between within",
                   age.x=="within"&age.y=="between"~"between within",
                   age.x=="within"&age.y=="within"~"within within"))

write.table(distPairs_long_subset,file=paste0("distPairs_long_subset",distaname,".txt"),sep="\t",quote=F)


p3 <- distPairs_long_subset %>% 
  ggplot(aes(x=age_category, y=cor,fill=same_family)) + 
  geom_boxplot(outlier.shape = NA)+ 
geom_point(position = position_jitterdodge(jitter.width = 0.1), alpha=0.3)+ 
  theme_bw()+ylim(0,1)

name_graph <- paste0(paste("distance",distaname,"prop_agecategory",sep="_"),".pdf")
ggsave(name_graph,p,width=20,height=20,units="cm")

name_graph <- paste0(paste("distance",distaname,"prop_agecategory",sep="_"),".png")
ggsave(name_graph,p,width=20,height=20,units="cm")


return(list(p1,p2,p3))
}
  
  
resprop<-compute_cor_prop(distaname="pearson")

#distPairs_long_subset[distPairs_long_subset$cor<0.8,]

 library(ggpubr)
#spearman
p<-ggarrange(plots3[[1]],resprop[[1]],plots3[[3]],plots2[[3]], ncol = 2, nrow = 2,  common.legend = TRUE,labels=c("all","prop","all","DE"))
ggsave("cor_spearman.pdf",p,width=20,height=20,units="cm")
ggsave("cor_spearman.png",p,width=20,height=20,units="cm")



#pearson
p<-ggarrange(plots4[[1]],resprop[[1]],plots4[[3]], plots1[[3]], ncol = 2, nrow = 2,  common.legend = TRUE,labels=c("all","prop","all","DE"))
ggsave("cor_pearson.pdf",p,width=20,height=20,units="cm")
ggsave("cor_pearson.png",p,width=20,height=20,units="cm")

```

Annotate age of transitions for all total dataset

```{r}


# annotate non-convergent pairs :

toremove = list( 
  c("Fukomys_damarensis","Fukomys_mechowii"),
  c("Fukomys_damarensis","Fukomys_micklemi"),
  c("Fukomys_mechowii","Fukomys_micklemi"),
  c("Nannomys_cot","Nannomys_indutus"),
  c("Mus_musculus","Mus_musculus_ddo"),
  c("Mus_musculus","Mus_musculus_wla"),
  c("Mus_musculus_ddo","Mus_musculus_wla"),
  c("Mastomys_coucha","Mastomys_erythroleucus"),
  c("Mastomys_coucha","Mastomys_kollmannspergeri"),
  c("Mastomys_erythroleucus","Mastomys_kollmannspergeri"),
  c("Abrothrix_longipilis","Abrothrix_olivaceus")
)
toremove <- data.frame(t(data.frame(toremove)))
toremove$pair <- paste(toremove$X1,toremove$X2)


# annotate pairs from within/between datasets: 

age = list(c("Acomys_dimidiatus","old"),
c("Apodemus_mystacinus","young"),
c("Arvicanthis_niloticus","young"),    
c("Chaetodipus_baileyi","old"),
c("Chinchilla_lanigera","old"),
c("Dipodomys_spectabilis2","old"),   
c("Eospalax_fontanierii","old"),
c("Fukomys_damarensis","old"),
c("Fukomys_mechowii","old"),     
c("Fukomys_micklemi","old"),
c("Heterocephalus_glaber","old"),
c("Mastomys_coucha","young"),          
c("Mastomys_erythroleucus","young"),
c("Mastomys_kollmannspergeri","young"),
c("Meriones_unguiculatus","old"),    
c("Mesocricetus_auratus","old"),
c("Micaelamys_namaquensis","old"),
c("Mus_macedonicus","young"),          
c("Nannomys_cot","young"),
c("Nannomys_indutus","young"),
c("Peromyscus_eremicus","young")      
# ("Saccostomus_campestris","NA") REMOVED because confidence interval is too large
)
age <- data.frame(t(data.frame(age)))
names(age)=c("Species","age")



compute_cor=function(distaname="pearson",genesetname="unionDE",geneset=unionDE)
{
  countsDE = counts_sp_norm %>% 
  subset(select=coldata_tot$ID)  %>% 
  filter(rownames(counts_sp_norm) %in% geneset) %>% 
  rowGrpMeans(grp = coldata_tot$species, na.rm = TRUE)
  
distPairs <- as.data.frame(cor(countsDE, method=distaname))
distPairs$Species1 <- rownames(distPairs)


distPairs_long_all <- distPairs%>%
  pivot_longer(-Species1,names_to="Species2",values_to="cor")%>%
  inner_join(coldata_exp_sp, join_by(Species1 == species))%>%
  inner_join(coldata_exp_sp, join_by(Species2 == species))%>%
  filter(Species1!=Species2)%>%
  mutate(pair = paste(Species1,Species2))%>%
  anti_join(toremove, join_by(pair==pair))%>%
  mutate(same_cond = paste(cond_season.x,cond_season.y))%>%
  mutate(same_family = family.x==family.y)%>%
  mutate(intra = Species1 < Species2)%>%
  filter(intra==TRUE)%>%
  mutate(condition = case_when(cond_season.x=="mesic"&cond_season.y=="xeric"~"mesic xeric",
                   cond_season.x=="xeric"&cond_season.y=="mesic"~"mesic xeric",
                   cond_season.x=="xeric"&cond_season.y=="xeric"~"xeric xeric",
                   cond_season.x=="mesic"&cond_season.y=="mesic"~"mesic mesic"))%>%
  separate(Species1,c("genera1","sp1"), sep = "_",remove=FALSE)%>%
  separate(Species2,c("genera2","sp2"), sep = "_",remove=FALSE)%>%
  mutate(same_genus = genera1==genera2)

# select members of pair that belong to same family / different family 


p1 <- distPairs_long_all %>% 
  ggplot(aes(x=same_family, y=cor,fill=condition)) + geom_boxplot(outlier.shape = NA)+ 
geom_point(position = position_jitterdodge(jitter.width = 0.1), alpha=0.3)+ 
  scale_fill_manual(values=c("darkgreen", "grey","orange" )) + theme_bw()#+ylim(0,1)

name_graph <- paste0(paste("distance",distaname,"genes",genesetname,"allpairs_family",sep="_"),".pdf")
ggsave(name_graph,p,width=20,height=20,units="cm")

name_graph <- paste0(paste("distance",distaname,"genes",genesetname,"allpairs_family",sep="_"),".png")
ggsave(name_graph,p,width=20,height=20,units="cm")




p2 <- distPairs_long_all %>% 
  ggplot(aes(x=same_genus, y=cor,fill=condition)) + geom_boxplot(outlier.shape = NA)+ 
geom_point(position = position_jitterdodge(jitter.width = 0.1), alpha=0.3)+ 
  scale_fill_manual(values=c("darkgreen", "grey","orange" )) + theme_bw()#+ylim(0,1)

name_graph <- paste0(paste("distance",distaname,"genes",genesetname,"allpairs_genus",sep="_"),".pdf")
ggsave(name_graph,p,width=20,height=20,units="cm")

name_graph <- paste0(paste("distance",distaname,"genes",genesetname,"allpairs_genus",sep="_"),".png")
ggsave(name_graph,p,width=20,height=20,units="cm")




distPairs_long_subset <- distPairs_long_all%>%
   inner_join(age, join_by(Species1 == Species))%>%
   inner_join(age, join_by(Species2 == Species))%>%
  mutate(age_category = case_when(age.x=="old"&age.y=="old"~"old old",
                   age.x=="old"&age.y=="young"~"old young",
                   age.x=="young"&age.y=="old"~"old young",
                   age.x=="young"&age.y=="young"~"young young"))


p3 <- distPairs_long_subset %>% 
  ggplot(aes(x=age_category, y=cor,fill=same_family)) + geom_boxplot(outlier.shape = NA)+ 
geom_point(position = position_jitterdodge(jitter.width = 0.1), alpha=0.3)+ theme_bw()#+ylim(0,1)

name_graph <- paste0(paste("distance",distaname,"genes",genesetname,"agecategory_allpairs",sep="_"),".pdf")
ggsave(name_graph,p,width=20,height=20,units="cm")

name_graph <- paste0(paste("distance",distaname,"genes",genesetname,"agecategory_allpairs",sep="_"),".png")
ggsave(name_graph,p,width=20,height=20,units="cm")


return(list(p1,p2,p3))
}
  
  


plots5<-compute_cor(distaname="pearson",genesetname="unionDE",geneset=unionDE)
plots6<-compute_cor(distaname="spearman",genesetname="unionDE",geneset=unionDE)
# compute_cor(distaname="spearman",genesetname="murinae",geneset=deMurinae$Gene)
# compute_cor(distaname="spearman",genesetname="ancient",geneset=deAncient$Gene)
# compute_cor(distaname="spearman",genesetname="recent",geneset=deRecent$Gene)
# compute_cor(distaname="spearman",genesetname="total",geneset=deTotal$Gene)
# compute_cor(distaname="spearman",genesetname="murinaeexcluDE",geneset=murinaeexcluDE)
# compute_cor(distaname="spearman",genesetname="ancientexcluDE",geneset=ancientexcluDE)
plots7<-compute_cor(distaname="spearman",genesetname="all",geneset=rownames(counts_sp_norm))

# compute_cor(distaname="pearson",genesetname="murinae",geneset=deMurinae$Gene)
# compute_cor(distaname="pearson",genesetname="ancient",geneset=deAncient$Gene)
# compute_cor(distaname="pearson",genesetname="recent",geneset=deRecent$Gene)
# compute_cor(distaname="pearson",genesetname="total",geneset=deTotal$Gene)
# compute_cor(distaname="pearson",genesetname="murinaeexcluDE",geneset=murinaeexcluDE)
# compute_cor(distaname="pearson",genesetname="ancientexcluDE",geneset=ancientexcluDE)
plots8<-compute_cor(distaname="pearson",genesetname="all",geneset=rownames(counts_sp_norm))



 library(ggpubr)
#spearman
p<-ggarrange(plots3[[1]],resprop[[1]],plots7[[3]],plots6[[3]], ncol = 2, nrow = 2,  common.legend = TRUE,labels=c("all","prop","all","DE"))
ggsave("cor_spearman_allpairs.pdf",p,width=20,height=20,units="cm")
ggsave("cor_spearman_allpairs.png",p,width=20,height=20,units="cm")



#pearson
p<-ggarrange(plots4[[1]],resprop[[1]],plots8[[3]], plots5[[3]], ncol = 2, nrow = 2,  common.legend = TRUE,labels=c("all","prop","all","DE"))
ggsave("cor_pearson_allpairs.pdf",p,width=20,height=20,units="cm")
ggsave("cor_pearson_allpairs.png",p,width=20,height=20,units="cm")


```
