---
title: "Rat scRNAseq"
format: html
editor: visual
---

```{r}
# INSTALL
#remotes::install_github("satijalab/seurat", "seurat5", quiet = TRUE)
#install.packages("paletteer")
#BiocManager::install("TOAST")
#devtools::install_github('xuranw/MuSiC')
#BiocManager::install("SingleCellExperiment")
#install.packages("FactoMineR")
#install.packages("factoextra")
#BiocManager::install("glmGamPoi")
#install.packages('devtools')
#devtools::install_github('immunogenomics/presto')

library(Seurat)
library(tidyverse)
library(readxl)
library(paletteer)
library(pheatmap)
library(MuSiC)
library(SingleCellExperiment)
library("FactoMineR")
library("factoextra")
library(ggforce)
library(glmGamPoi)
library(cowplot)
 library(ggpubr)


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

folder_project="/home/rstudio/mydatalocal/"
dir.create(paste0(folder_project,"/results"))
setwd(paste0(folder_project,"/results"))
dir.create(paste0(folder_project,"/data"))


```

## Rat data

Load count

```{r}


# get metadata
 system("wget 'https://www.ncbi.nlm.nih.gov/geo/download/?acc=GSE209821&format=file&file=GSE209821%5FEXPORT%5FGEO%5Fmeta%2Edata%2Ecsv%2Egz' -O ~/mydatalocal/data/GSE209821_EXPORT_GEO_meta.data.csv")

 #get counts / or direct upload from local machine
#system("wget 'https://www.ncbi.nlm.nih.gov/geo/download/?acc=GSE209821&format=file&file=GSE209821%5FEXPORT%5FGEO%5Fcounts%5Fpostfilter%2Erds%2Egz' -O ~/mydatalocal/data/GSE209821_EXPORT_GEO_counts_postfilter.Rds")

scRNA.df <- readRDS("../data/GSE209821_EXPORT_GEO_counts_postfilter.rds")
rownames(scRNA.df) <- toupper(rownames(scRNA.df))

scRNA.matrix <- data.matrix(scRNA.df)

phenotable <- read.table(file = "../data/GSE209821_EXPORT_GEO_meta.data.csv",sep=",",h=T,row.names=1)

scRNA.matrix<-scRNA.matrix[,phenotable$exp.cond=="Lean"]
phenotable<-phenotable[phenotable$exp.cond=="Lean",]



rat<-CreateSeuratObject(counts=scRNA.matrix, meta.data=phenotable )


```

Process rat scRNAseq data

```{r}

rat <- subset(rat,
              subset = nCount_RNA > 200 &
                nCount_RNA <  median(rat$nCount_RNA) * 10 ) %>%
  Seurat::SCTransform(
    vst.flavor = "v2",
    verbose = FALSE,
    min_cell = 5,
    vars.to.regress = NULL
  ) %>%
  Seurat::RunPCA(npcs = 50, verbose = FALSE) %>%
  Seurat::RunUMAP(reduction = "pca",
                  dims = 1:50,
                  verbose = FALSE) %>%
  Seurat::FindNeighbors(reduction = "pca",
                        dims = 1:50,
                        verbose = FALSE) %>%
  Seurat::FindClusters(resolution = 0.8, verbose = FALSE)



gc()


rat[["percent.mt"]] <- PercentageFeatureSet(rat, pattern = "^MT-")

Idents(rat)=rat$orig.ident
virat <- VlnPlot(rat, features = c("nFeature_RNA", "nCount_RNA", "percent.mt"), ncol = 3, alpha = 0.05) & theme(legend.position = 'none',axis.title.x=element_blank(),axis.text.x=element_blank()) 
ggsave(virat, file="rat_violin.png",width=15,height=10,units="cm")
ggsave(virat, file="rat_violin.pdf",width=15,height=10,units="cm")

# nb of genes per cell
summary(rat$nFeature_RNA)
# nb of cells per gene

save(rat, file = "rat.RData")


rat1 <- DimPlot(object = rat,group.by = "orig.ident", split.by="orig.ident",ncol=3)+   ggtitle(NULL)  +
  coord_fixed() + NoLegend()
ggsave(rat1, file="rat_umap1.png")
ggsave(rat1, file="rat_umap1.pdf")



fig_rat_qc=plot_grid(virat, rat1, labels = c('A', 'B'), ncol = 1, rel_heights =  c(0.35, 1), rel_widths  =  c(0.35, 1))
ggsave(fig_rat_qc, file="fig_rat_qc.png",width=20,height=15,units="cm")
ggsave(fig_rat_qc, file="fig_rat_qc.pdf",width=20,height=15,units="cm")



Idents(rat)=rat$clusters2
rat1 =DimPlot(object = rat, label = TRUE) +
  coord_fixed()
ggsave(rat1, file="rat_umap_ident2.png")
ggsave(rat1, file="rat_umap_ident2.pdf")



Idents(rat)=rat$clusters3
rat1 =DimPlot(object = rat, label = TRUE) +
  coord_fixed()
ggsave(rat1, file="rat_umap_ident3.png")
ggsave(rat1, file="rat_umap_ident3.pdf")


rat1 =DimPlot(object = rat, label = FALSE) +
  coord_fixed()+ NoAxes() + NoLegend()
ggsave(rat1, file="rat_umap_ident3_wolabel.png")
ggsave(rat1, file="rat_umap_ident3_wolabel.pdf")


# get marker genes for these cell types
rat.markers <- FindAllMarkers(rat)
rat.markers_filter <- rat.markers %>%
    group_by(cluster) %>%
    dplyr::filter(avg_log2FC > 1)


rat.markers_filter %>%
    group_by(cluster) %>%
    dplyr::filter(avg_log2FC > 1) %>%
    slice_head(n = 10) %>%
    ungroup() -> top10
heatmk <- DoHeatmap(rat, features = top10$gene, size=5)+guides(colour=FALSE)+ scale_y_discrete(labels=NULL)

ggsave(heatmk, file="rat_markers_top10.png",width=35,height=25,units="cm")
ggsave(heatmk, file="rat_markers_top10.pdf",width=35,height=25,units="cm")





rat.markers_filter %>%
    group_by(cluster) %>%
    dplyr::filter(avg_log2FC > 1) %>%
    slice_head(n = 50) %>%
    ungroup() -> top50
heatmk <- DoHeatmap(rat, features = top50$gene)+guides(colour=FALSE)+ scale_y_discrete(labels=NULL)

ggsave(heatmk, file="rat_markers_top50.png",width=35,height=25,units="cm")
ggsave(heatmk, file="rat_markers_top50.pdf",width=35,height=25,units="cm")




rat.markers_filter %>%
    group_by(cluster) %>%
    dplyr::filter(avg_log2FC > 1) %>%
    slice_head(n = 100) %>%
    ungroup() -> top100
heatmk <- DoHeatmap(rat, features = top100$gene)+guides(colour=FALSE)+ scale_y_discrete(labels=NULL) 

ggsave(heatmk, file="rat_markers_top100.png",width=35,height=25,units="cm")
#ggsave(heatmk, file="rat_markers_top100.pdf",width=35,height=25,units="cm")



# get average expression for the whole dataset
Idents(rat)=rat$exp.cond
rat_avg_tot <- AggregateExpression(object = rat)$RNA

# get average expression individual per individual
Idents(rat)=rat$orig.ident
rat_avg_ind <- AggregateExpression(object = rat)$RNA

rat_avg_ind <- as.data.frame(rat_avg_ind)
rat_avg_ind$tot=rat_avg_tot[,1]

write.table(file="rat_pseudobulk.txt",rat_avg_ind,row.names=T,sep="\t",quote=F)

# get average expression individual per cell type
Idents(rat)=rat$clusters3
rat_avg <- AggregateExpression(object = rat)$RNA


```

Represent expression levels for markers genes on rat clusters

```{r}

markers <- as.data.frame(read_xlsx(path = "../data/cell-types-markers_associated-biblio.xlsx",sheet ="All-markers"))
markers$MGI<-toupper(markers$`Genes MGI`)
#non-unique markers: 
nnu <- names(table(markers$MGI)[table(markers$MGI)==2])
markers<-markers[!markers$MGI%in%nnu,]
markers<-markers[!is.na(markers$Segment),]
markers$Segment[markers$Segment=="-"]="Immune"
# Markers annotated with "-" correspond to Immune cells. We remove markers annotated with NA as they are not identified

## DE genes in murine
de <- read.table("../data/Table_DE_results_cond_season_murinae_deseq.tsv",sep="\t",h=T)

markers$Genes=toupper(markers$`Genes MGI`)

rat_avg_mk <- merge(rat_avg,markers,by.x=0,by.y="Genes")
rownames(rat_avg_mk)<-rat_avg_mk$Row.names
rat_avg_mk = with(rat_avg_mk,  rat_avg_mk[order(rat_avg_mk$Segment) , ])

de$X2<-toupper(row.names(de))
datade <- merge(rat_avg_mk,de,by.x="Row.names",by.y=0,all.x=T)
datade$xm<-ifelse(datade$log2FoldChange>0,"xeric","mesic")
datade$xm[is.na(datade$log2FoldChange)]="NA"
datade = datade[order(paste0(datade$Segments,datade$xm)),]



gene_functions_df <- data.frame(Segment = rat_avg_mk$Segment)
row.names(gene_functions_df) <- rownames(rat_avg_mk)

newCols <- paletteer_d("ggsci::springfield_simpsons")[1:length(unique(rat_avg_mk$Segment))]
annoRow <- newCols

ann_colors = list(Segments = annoRow)


heat_plot <- pheatmap(rat_avg_mk[,2:9],
scale="row",
cluster_rows = FALSE,
cluster_cols = TRUE, # set to FALSE if you want to remove the dendograms
fontsize_row = 10, # row label font size
fontsize_col = 7, # column label font size
angle_col = 45, # sample names at an angle
show_colnames = T, show_rownames = F, # displaying column and row names,
annotation_row = gene_functions_df,
#annotation_colors = ann_colors,
main = "Marker gene expression in Rat scRNAseq data") # a title for our heatmap
pdf("Rat_markers_per_celltype.pdf", height = 20, width = 10)
heat_plot
dev.off()

pdf("Rat_markers_per_celltype.pdf", height = 7, width = 5)
heat_plot
dev.off()
  
  

## with DE genes
ann_df <- data.frame(CT = colnames(datade)[2:9])
row.names(ann_df) <- colnames(datade)[2:9]
head(ann_df)

gene_functions_df <- data.frame(Segments = datade$Segments,  Condition=datade$xm)
row.names(gene_functions_df) <- rownames(datade)
# Markers annotated with "-" correspond to Immune cells. We remove markers annotated with NA as they are not identified
gene_functions_df$Segments[gene_functions_df$Segments=="-"]="Immune"
gene_functions_df<-gene_functions_df[!is.na(gene_functions_df$Segment),]
datade<-datade[!is.na(datade$Segment),]
datade1<-datade[datade$Segments!="Immune"&datade$Segments!="-",]
datade1<-datade1[,c("Row.names","Endo","Prox tub","IC","LOH","DCT/CNT/PC","Stroma","Immune","Podo","Segments")]

newCols <- paletteer_d("ggsci::springfield_simpsons")[1:length(unique(datade1$Segments))]
annoRow <- newCols
names(annoRow) <- unique(datade1$Segment)
newCols <- paletteer_d("ggthemes::stata_s2color")[1:8] 
annoCol <- newCols
names(annoCol) <- names(datade1)[2:9]

ann_colors = list(
  Segments = annoRow, CT = annoCol, Condition = c(xeric="orange", mesic="green","NA"="white"))


heat_plot <- pheatmap(datade1[,2:9], 
                      scale="row",
          annotation_colors = ann_colors,
                      cluster_rows = FALSE, 
                      cluster_cols = FALSE, # set to FALSE if you want to remove the dendograms
                      annotation_row = gene_functions_df, # row (gene) annotations
                      annotation_col = ann_df, # column (sample) annotations
                   
                      annotation_names_row = F, 
                      annotation_names_col = F,
                      fontsize_row = 10,          # row label font size
                      fontsize_col = 7,          # column label font size 
                      angle_col = 45, # sample names at an angle
                      #legend_breaks = c(-2, 0, 2), # legend customisation
                      #legend_labels = c("Low", "Medium", "High"), # legend customisation
                      show_colnames = T, show_rownames = F, # displaying column and row names
                      main = "Marker gene expression in rat data") # a title for our heatmap
  
pdf("Marker_rat_deconvolution_plot_short.pdf", height = 10, width = 8)
heat_plot
dev.off()


ordered=c("IC","DCT/CNT/PC","Endo","Stroma","LOH","Podo","Prox tub","Immune")    

ann_df <- data.frame(CT = ordered)
row.names(ann_df) <- ordered
head(ann_df)

heat_plot <- pheatmap(datade1[,ordered], 
                      scale="row",
          annotation_colors = ann_colors,
                      cluster_rows = FALSE, 
                      cluster_cols = FALSE, # set to FALSE if you want to remove the dendograms
                      annotation_row = gene_functions_df, # row (gene) annotations
                      annotation_col = ann_df, # column (sample) annotations
                   
                      annotation_names_row = F, 
                      annotation_names_col = F,
                      fontsize_row = 10,          # row label font size
                      fontsize_col = 7,          # column label font size 
                      angle_col = 45, # sample names at an angle
                      #legend_breaks = c(-2, 0, 2), # legend customisation
                      #legend_labels = c("Low", "Medium", "High"), # legend customisation
                      show_colnames = T, show_rownames = F, # displaying column and row names
                      main = "Marker gene expression in rat data") # a title for our heatmap
  
pdf("Marker_rat_deconvolution_plot_short_reorder.pdf", height = 10, width = 8)
heat_plot
dev.off()



```

### Check quality of deconvolution on pseudobulks: for rat data

```{r}
# we use 3 individuals separately or the total


sample.ids <- phenotable$orig.ident
# individual.ids and cell.types should be in the same order as in sample.ids
individual.labels <- row.names(phenotable)
cell.type.labels <- phenotable$clusters3
  

sc.pheno <- data.frame(check.names=F, check.rows=F,
                       stringsAsFactors=F,
                       row.names=individual.labels,
                       SubjectName=sample.ids,
                       cellType=cell.type.labels)

sc.meta <- data.frame(labelDescription=c("SubjectName",
                                         "cellType"),
                      row.names=c("SubjectName",
                                  "cellType"))
sc.pdata <- new("AnnotatedDataFrame",
                data=sc.pheno,
                varMetadata=sc.meta)

sc.eset <- ExpressionSet(assayData=scRNA.matrix,
                                  phenoData=sc.pdata)


sc.eset2 <- SingleCellExperiment(list(counts=scRNA.matrix),
    colData=DataFrame(cellTypeID=sc.pheno$cellType,cellType=sc.pheno$cellType,sampleID=sc.pheno$SubjectName,SubjectName=sc.pheno$SubjectName))



pseudobulk.matrix <- data.matrix(cbind(round(rat_avg_ind,0)))
pseudobulk.matrix[is.na(pseudobulk.matrix)] <- 0 
pseudobulk.eset <- ExpressionSet(assayData = pseudobulk.matrix)
bulk.mtx = exprs(pseudobulk.eset)

## Deconvolve with selected CTs for pseudobulk
prop_music_pseudorat=music_prop(bulk.mtx = bulk.mtx, sc.sce = sc.eset2, 
                      clusters = 'cellType', samples = 'sampleID',
                      select.ct = c("DCT/CNT/PC","Endo","IC","LOH","Podo","Prox tub","Stroma"), 
                      verbose = F)$Est.prop.weighted
write.table(file="prop_music_pseudobulk_rat.txt",prop_music_pseudorat,quote=F)

prop_music_pseudorat=as.data.frame(prop_music_pseudorat)

trueprop<-table(rat$clusters3[rat$clusters3!="Immune"],rat$exp.cond[rat$clusters3!="Immune"])
pind<-table(rat$clusters3[rat$clusters3!="Immune"],rat$orig.ident[rat$clusters3!="Immune"])
trueprop_pseudo_rat<-data.frame(pind/apply(pind,2,sum))
trueprop_pseudo_rat_tot<-data.frame(trueprop/sum(trueprop))
trueprop_pseudo_rat<-rbind(trueprop_pseudo_rat,trueprop_pseudo_rat_tot)

prop_music_pseudorat$sample=row.names(prop_music_pseudorat)
prop_music_pseudorat=prop_music_pseudorat%>%as_tibble%>%pivot_longer(-sample)

prop_music_pseudorat$sample=gsub("-","_",prop_music_pseudorat$sample)
prop_music_pseudorat$sample[prop_music_pseudorat$sample=="tot"]="Lean"


compareprop<-trueprop_pseudo_rat%>%left_join(prop_music_pseudorat,by=c("Var2"="sample","Var1"="name"))%>%
  rename("Var1"="Cell_Type")%>%
  mutate("Pseudobulk" = case_when(str_detect(Var2,"Lean") ~ "Total",str_detect(Var2,"RK") ~ "Individual"))

  


g<-ggplot(compareprop,aes(Freq,value,col=Cell_Type,shape=Pseudobulk))+geom_point()+xlab("true proportion")+ylab("estimated proportion")+ theme_bw()
ggsave(g,file="proportions_frompseudorat_trueest.png")
ggsave(g,file="proportions_frompseudorat_trueest.pdf")

```

### Check quality of deconvolution on pseudobulks: for mouse data

```{r}
## with mouse data pseudobulk
+

pseudobulk.matrix <- data.matrix(cbind(round(mouse_avg_ind,0)))
pseudobulk.matrix[is.na(pseudobulk.matrix)] <- 0 
pseudobulk.eset <- ExpressionSet(assayData = pseudobulk.matrix)
bulk.mtx = exprs(pseudobulk.eset)

## Deconvolve with selected CTs for pseudobulk
prop_music_pseudomus=music_prop(bulk.mtx = bulk.mtx, sc.sce = sc.eset2, 
                      clusters = 'cellType', samples = 'sampleID',
                      select.ct = c("DCT/CNT/PC","Endo","IC","LOH","Podo","Prox tub","Stroma"), 
                      verbose = F)$Est.prop.weighted
write.table(file="prop_music_pseudobulk_mouse.txt",prop_music_pseudomus,quote=F)

prop_music_pseudomus=as.data.frame(prop_music_pseudomus)

trueprop<-table(mouse$V2,mouse$orig.ident)
pind<-table(mouse$V2,mouse$sample)
trueprop_pseudo_mouse<-data.frame(pind/apply(pind,2,sum))
trueprop_pseudo_mouse_tot<-data.frame(trueprop/sum(trueprop))
trueprop_pseudo_mouse<-rbind(trueprop_pseudo_mouse,trueprop_pseudo_mouse_tot)

prop_music_pseudomus$sample=paste0("g",row.names(prop_music_pseudomus))
prop_music_pseudomus=prop_music_pseudomus%>%as_tibble%>%pivot_longer(-sample)

prop_music_pseudomus$name[prop_music_pseudomus$name=="Prox tub"]="PT"
prop_music_pseudomus$name[prop_music_pseudomus$name=="DCT"]="DCT/CNT/PC"
trueprop_pseudo_mouse$Var2=paste0("g",trueprop_pseudo_mouse$Var2)
trueprop_pseudo_mouse$Var2=gsub("SeuratProject","tot",trueprop_pseudo_mouse$Var2)



compareprop<-trueprop_pseudo_mouse%>%left_join(prop_music_pseudomus,by=c("Var2"="sample","Var1"="name"))


g<-ggplot(compareprop,aes(Freq,value,col=Var1,shape=Var2=="Lean"))+geom_point()+xlab("true proportion")+ylab("estimated proportion")+ theme_minimal()
ggsave(g,file="proportions_frompseudomus_trueest.png")
ggsave(g,file="proportions_frompseudomus_trueest.pdf")


```

# Deconvolutions with rat reference data

```{r}


group.colors <- c(arid =  "#FFA500", mesic =  "#5EA336")


# Prepare input bulk data for deconvolution 
bulk.df <- read.table("../data/datalocalandonline_cds_counts_2025_final.tsv",h=T,row.names=1)
rownames(bulk.df) <- toupper(rownames(bulk.df))

# select Murinae species
coldata <- read.table("../data/coldata_July2025_family.tsv",h=T,sep="\t")
coldataMurinae <- coldata[coldata$muridae_set=="yes",]
coldataMurinae$ID <- gsub("kollmannspergeri","kollmanspergeri",coldataMurinae$ID)
bulk.df_Murinae <- bulk.df[,coldataMurinae$ID]
#coldataMurinae <- coldataMurinae[-grep("kollmannspergeri",coldataMurinae$ID),]
bulk.df_Murinae_sort <- bulk.df_Murinae[,coldataMurinae$ID]


bulk.matrix <- data.matrix(round(bulk.df_Murinae_sort,0))
bulk.matrix[is.na(bulk.matrix)] <- 0 
bulk.eset <- ExpressionSet(assayData = bulk.matrix)

row.names(coldataMurinae)=coldataMurinae$ID
metadata <- data.frame(labelDescription= names(coldataMurinae)[1:20], row.names=names(coldataMurinae)[1:20])

bulk.eset <- ExpressionSet(assayData = data.matrix(bulk.matrix),
                           phenoData = new("AnnotatedDataFrame", data = coldataMurinae[,1:20], 
                                           varMetadata = metadata) )



# individual.ids and cell.types should be in the same order as in sample.ids
sample.ids <- phenotable$orig.ident
individual.labels <- row.names(phenotable)
cell.type.labels <- phenotable$clusters3

sc.pheno <- data.frame(check.names=F, check.rows=F,
                       stringsAsFactors=F,
                       row.names=individual.labels,
                       SubjectName=sample.ids,
                       cellType=cell.type.labels)

sc.meta <- data.frame(labelDescription=c("SubjectName",
                                         "cellType"),
                      row.names=c("SubjectName",
                                  "cellType"))
sc.pdata <- new("AnnotatedDataFrame",
                data=sc.pheno,
                varMetadata=sc.meta)

sc.eset <- ExpressionSet(assayData=scRNA.matrix,
                                  phenoData=sc.pdata)


sc.eset2 <- SingleCellExperiment(list(counts=scRNA.matrix),
    colData=DataFrame(cellTypeID=sc.pheno$cellType,cellType=sc.pheno$cellType,sampleID=sc.pheno$SubjectName,SubjectName=sc.pheno$SubjectName))

bulk.mtx = exprs(bulk.eset)

## With selected CTs

prop_music=music_prop(bulk.mtx = bulk.mtx, sc.sce = sc.eset2, 
                      clusters = 'cellType', samples = 'sampleID',
                      select.ct = c("DCT/CNT/PC","Endo","IC","LOH","Podo","Prox tub","Stroma"), 
                      verbose = F)$Est.prop.weighted
colnames(prop_music)=c("DCT/CNT/PC","Endo","IC","LOH","Podo","Prox_tub","Stroma")
write.table(file="prop_music_murinae_rat.txt",prop_music,quote=F)
## With all CTs


# prop_music


## PCA with selected cell types
ref.MurinaeM_Music <- merge(prop_music,coldataMurinae,by.x=0,by.y="ID")
#ref.MurinaeM_Music$DCT=ref.MurinaeM_Music$DCT+ref.MurinaeM_Music$LOH
table_ACPM = ref.MurinaeM_Music
table_ACPM$ID=paste(table_ACPM$sp_short_name,row.names(table_ACPM))
rownames(table_ACPM)=table_ACPM$ID
res_ACPM = table_ACPM %>%  PCA(quali.sup = c(1,9:38), graph = FALSE)
var <- get_pca_var(res_ACPM)

pdf("Murinae_Biplot_per_cond_season_PCA_Music_selectedCellTypes_ind.pdf")
fviz_pca_biplot (res_ACPM,col.var = "black", col.ind = as.factor(table_ACPM$cond_season), 
             addEllipses = FALSE,  alpha.ind=0.8, 
             palette = group.colors, legend.title = "Groups",
            ggtheme 
            = theme_bw())
dev.off()


### represent percentage of cell types

p<- table_ACPM%>%
  select(c("Row.names","DCT/CNT/PC","Endo","IC","LOH","Podo","Prox_tub", "Stroma","species","cond_season","sp_short_name","ID"))%>%
  pivot_longer(-c("species","sp_short_name","cond_season","Row.names","ID"),names_to="CT",values_to="prop")%>%
  ggplot(aes(x=cond_season,y=prop,fill=cond_season))+
  geom_boxplot(outlier.shape = NA)+
  geom_jitter(alpha=0.2,height=0,width=0.1)+
  scale_fill_manual(values=c("arid"="#FFA500","mesic"="#5EA336"))+
  facet_wrap(~CT, scales = "free")+theme_bw()+
  ylab("Proportion")+
  ylab("Condition")+
  theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1),strip.text = element_text(
    size = 6, angle = 90))
   
ggsave(p,file="boxplot_celltypes_Music_selectedCellTypes_Murinae.pdf",width = 10,
  height = 20,units="cm")



p<- table_ACPM%>%
  select(c("Row.names","DCT/CNT/PC","Endo","IC","LOH","Podo","Prox_tub", "Stroma","species","cond_season","sp_short_name","ID"))%>%
  pivot_longer(-c("species","sp_short_name","cond_season","Row.names","ID"),names_to="CT",values_to="prop")%>%
  ggplot(aes(x=cond_season,y=prop,fill=cond_season))+
  geom_boxplot(outlier.shape = NA)+
  geom_jitter(alpha=0.2,height=0,width=0.1)+
  scale_fill_manual(values=c("arid"="#FFA500","mesic"="#5EA336"))+
  facet_row(~CT, scales = "free")+theme_bw()+
  ylab("Proportion")+
  ylab("Condition")+
  theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1),strip.text = element_text(
    size = 6, angle = 90))
   
ggsave(p,file="boxplot_celltypes_Music_selectedCellTypes_Murinae.pdf",width = 30,
  height = 10,units="cm")



# Compare cell types inferred for mouse and rat to pseudobulk counts and inferences on pseudobulks


bulkprop_mouserat <- table_ACPM%>%filter(species%in%c("Mus_musculus","Rattus_norvegicus"))
compareprop<-trueprop_pseudo_rat%>%left_join(prop_music_pseudorat,by=c("Var2"="sample","Var1"="name"))

compareprop_rat_long<-compareprop%>%
  pivot_longer(-c(Var1,Var2)) %>%
  rename("Var1"="Celltype") %>%
  rename("Var2" = "Sample") %>%
  rename("name" = "Type") %>%
  relocate("Celltype", .after = "Sample") %>%
  mutate(Type=fct_recode(Type, sc_count = "Freq", sc_deconvolution = "value"))%>%
  mutate(Celltype=fct_recode(Celltype, Prox_tub = "Prox tub"))

  
bulkprop_mouserat_long<-bulkprop_mouserat%>%
  select("DCT/CNT/PC","Endo","IC","LOH","Podo","Prox_tub", "Stroma","ID")%>%
  pivot_longer(-ID) %>%
  rename("name"="Celltype") %>%
  rename("ID" = "Sample") %>%
  mutate("Type" = case_when(str_detect(Sample,"Mmu") ~ "bulk deconvolution Mus",str_detect(Sample,"Rno") ~ "bulk deconvolution Rat"))


compareprop_rat_bulk_sc <- bulkprop_mouserat_long%>%bind_rows(compareprop_rat_long)



p<- compareprop_rat_bulk_sc%>%
    ggplot(aes(x=Celltype,y=value,color=Type,group_by=Sample))+
  geom_jitter(alpha=0.4,height=0,width=0.1) +
  theme_bw()+
  ylab("Cell-type proportion")
ggsave(p,file="CellProp_comparisons_musrat.pdf",width = 20,
  height = 15,units="cm")
ggsave(p,file="CellProp_comparisons_musrat.png",width = 20,
  height = 15,units="cm")


compareprop_rat_bulk_sc_selec <- compareprop_rat_bulk_sc%>%
  filter(Celltype%in%c("Endo","LOH","Prox_tub", "Stroma"))%>%
  group_by(Sample,Type)%>%
  mutate(sumprop=sum(value))%>%
   mutate(value2=value/sumprop)
  
p2<- compareprop_rat_bulk_sc_selec%>%
    ggplot(aes(x=Celltype,y=value2,color=Type,group_by=Sample))+
  geom_jitter(alpha=0.4,height=0,width=0.1) +
  theme_bw()+
  ylab("Cell-type proportion")
ggsave(p2,file="CellProp_comparisons_musrat_selec.pdf",width = 20,
  height = 15,units="cm")
ggsave(p2,file="CellProp_comparisons_musrat_selec.png",width = 20,
  height = 15,units="cm")



#fig_deconvo_qc=plot_grid(p + theme(legend.position="none"), p2, labels = c('A', 'B'),rel_widths = c(1,1.1))
#ggsave(fig_deconvo_qc, file="fig_deconvo_qc.png",width=25,height=15,units="cm")
#ggsave(fig_deconvo_qc, file="fig_deconvo_qc.pdf",width=25,height=15,units="cm")
  
  
   

```

```{r}
# select all species
coldata <- read.table("../data/coldata_July2025_family.tsv",h=T,sep="\t")
coldataTotal <- coldata[coldata$total_set=="yes",]
coldataTotal$ID <- gsub("kollmannspergeri","kollmanspergeri",coldataTotal$ID)
bulk.df_Total <- bulk.df[,coldataTotal$ID]
#coldataMurinae <- coldataMurinae[-grep("kollmannspergeri",coldataMurinae$ID),]
bulk.df_Total_sort <- bulk.df_Total[,coldataTotal$ID]


bulk.matrix <- data.matrix(round(bulk.df_Total_sort,0))
bulk.matrix[is.na(bulk.matrix)] <- 0 
bulk.eset <- ExpressionSet(assayData = bulk.matrix)

row.names(coldataTotal)=coldataTotal$ID
metadata <- data.frame(labelDescription= names(coldataTotal)[1:20], row.names=names(coldataTotal)[1:20])

bulk.eset <- ExpressionSet(assayData = data.matrix(bulk.matrix),
                           phenoData = new("AnnotatedDataFrame", data = coldataTotal[,1:20], 
                                           varMetadata = metadata) )



sample.ids <- phenotable$orig.ident

# individual.ids and cell.types should be in the same order as in sample.ids
individual.labels <- row.names(phenotable)

  
cell.type.labels <- phenotable$clusters3

sc.pheno <- data.frame(check.names=F, check.rows=F,
                       stringsAsFactors=F,
                       row.names=individual.labels,
                       SubjectName=sample.ids,
                       cellType=cell.type.labels)

sc.meta <- data.frame(labelDescription=c("SubjectName",
                                         "cellType"),
                      row.names=c("SubjectName",
                                  "cellType"))
sc.pdata <- new("AnnotatedDataFrame",
                data=sc.pheno,
                varMetadata=sc.meta)

sc.eset <- ExpressionSet(assayData=scRNA.matrix,
                                  phenoData=sc.pdata)


sc.eset2 <- SingleCellExperiment(list(counts=scRNA.matrix),
    colData=DataFrame(cellTypeID=sc.pheno$cellType,cellType=sc.pheno$cellType,sampleID=sc.pheno$SubjectName,SubjectName=sc.pheno$SubjectName))

bulk.mtx = exprs(bulk.eset)

## With selected CTs

prop_music=music_prop(bulk.mtx = bulk.mtx, sc.sce = sc.eset2, 
                      clusters = 'cellType', samples = 'sampleID',
                      select.ct = c("DCT/CNT/PC","Endo","IC","LOH","Podo","Prox tub","Stroma"), 
                      verbose = F)$Est.prop.weighted
colnames(prop_music)=c("DCT/CNT/PC","Endo","IC","LOH","Podo","Prox_tub","Stroma")
write.table(file="prop_music_total_rat.txt",prop_music,quote=F,row.names=F,sep="\t")


## With all CTs
#prop_music=read.table("prop_music.txt",h=T,row.names=1)

prop_music2=music_prop(bulk.mtx = bulk.mtx, sc.sce = sc.eset2, 
                      clusters = 'cellType', samples = 'sampleID',
                      verbose = F)$Est.prop.weighted
## With all CTs + marker genes only

#prop_music3_markers=music_prop(bulk.mtx = bulk.mtx, sc.sce = sc.eset2, markers=markers$`Genes MGI`,
#                      clusters = 'cellType', samples = 'sampleID',
#                      verbose = F)$Est.prop.weighted

#prop_music3=as.data.frame(prop_music3_markers)





## PCA with selected cell types
ref.Total_Music <- merge(prop_music,coldataTotal,by.x=0,by.y="ID")
#ref.MurinaeM_Music$DCT=ref.MurinaeM_Music$DCT+ref.MurinaeM_Music$LOH
table_ACPM = ref.Total_Music
table_ACPM$ID=paste(table_ACPM$sp_short_name,row.names(table_ACPM))
rownames(table_ACPM)=table_ACPM$ID
res_ACPM = table_ACPM %>%  PCA(quali.sup = c(1,9:38), graph = FALSE)
var <- get_pca_var(res_ACPM)

write.table(file="ref_prop_music_total_rat.txt",table_ACPM,quote=F,row.names=F,sep="\t")



pdf("Total_Biplot_per_cond_season_PCA_Music_selectedCellTypes_ind.pdf")
fviz_pca_biplot (res_ACPM,col.var = "black", col.ind = as.factor(table_ACPM$cond_season), 
             addEllipses = FALSE,  alpha.ind=0.8, 
             palette = group.colors, legend.title = "Groups",
            ggtheme 
            = theme_bw())
dev.off()


#https://www.cfholbert.com/blog/cluster-pca/
pca1 <- dudi.pca(table_ACPM[,2:8], scan = FALSE, nf = 3)
rand1 <- rtest(discrimin(pca1, as.factor(table_ACPM$cond_season), scan = FALSE), 1000)
rand1
pdf("Total_effect_season_mctest.pdf")
plot(rand1, main = "Monte-Carlo test")
dev.off()

rand2 <- rtest(discrimin(pca1, as.factor(table_ACPM$family), scan = FALSE), 1000)
rand2
pdf("Total_effect_family_mctest.pdf")
plot(rand2, main = "Monte-Carlo test")
dev.off()


### represent percentage of cell types

p<- table_ACPM%>%
  select(c("Row.names","DCT/CNT/PC","Endo","IC","LOH","Podo","Prox_tub", "Stroma","species","cond_season","sp_short_name","ID"))%>%
  pivot_longer(-c("species","sp_short_name","cond_season","Row.names","ID"),names_to="CT",values_to="prop")%>%
  ggplot(aes(x=cond_season,y=prop,fill=cond_season))+
  geom_boxplot(outlier.shape = NA)+
  geom_jitter(alpha=0.2,height=0,width=0.1)+
  scale_fill_manual(values=c("arid"="#FFA500","mesic"="#5EA336"))+
  facet_wrap(~CT, scales = "free")+theme_bw()+
  ylab("Proportion")+
  ylab("Condition")+
  theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1),strip.text = element_text(
    size = 6, angle = 90))
   

ggsave(p,file="boxplot_celltypes_Music_selectedCellTypes_Total.pdf",width = 10,
  height = 20,units="cm")




p<- table_ACPM%>%
  select(c("Row.names","DCT/CNT/PC","Endo","IC","LOH","Podo","Prox_tub", "Stroma","species","cond_season","sp_short_name","ID"))%>%
  pivot_longer(-c("species","sp_short_name","cond_season","Row.names","ID"),names_to="CT",values_to="prop")%>%
  ggplot(aes(x=cond_season,y=prop,fill=cond_season))+
  geom_boxplot(outlier.shape = NA)+
  geom_jitter(alpha=0.2,height=0,width=0.1)+
  scale_fill_manual(values=c("arid"="#FFA500","mesic"="#5EA336"))+
  facet_row(~CT, scales = "free")+theme_bw()+
  ylab("Proportion")+
  ylab("Condition")+
  theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1),strip.text = element_text(
    size = 6, angle = 90))
   
ggsave(p,file="boxplot_celltypes_Music_selectedCellTypes_Total.pdf",width = 30,
  height = 10,units="cm")




p<- table_ACPM%>%
  select(c("Row.names","DCT/CNT/PC","Endo","IC","LOH","Prox_tub", "Stroma","species","cond_season","family","sp_short_name","ID"))%>%
  pivot_longer(-c("species","sp_short_name","cond_season","Row.names","ID","family"),names_to="CT",values_to="prop")%>%
  ggplot(aes(x=cond_season,y=prop,fill=cond_season))+
  geom_boxplot(outlier.shape = NA)+
  geom_jitter(alpha=0.2,height=0,width=0.1)+
  scale_fill_manual(values=c("arid"="#FFA500","mesic"="#5EA336"))+
  facet_row(~CT, scales = "free")+theme_bw()+
  ylab("Proportion")+
  xlab("Condition")+theme_bw()
  theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1),strip.text = element_text(
    size = 6, angle = 90))
   
ggsave(p,file="boxplot_celltypes_Music_selectedCellTypes_Total.pdf",width = 20,
  height = 7,units="cm")
ggsave(p,file="boxplot_celltypes_Music_selectedCellTypes_Total.png",width = 20,
  height = 7,units="cm")



p<- table_ACPM%>%
  select(c("Row.names","DCT/CNT/PC","Endo","IC","LOH","Podo","Prox_tub", "Stroma","species","cond_season","family","sp_short_name","ID"))%>%
  pivot_longer(-c("species","sp_short_name","cond_season","Row.names","ID","family"),names_to="CT",values_to="prop")%>%
  ggplot(aes(x=cond_season,y=prop,fill=cond_season))+
  geom_boxplot(outlier.shape = NA)+
  geom_jitter(alpha=0.2,height=0,width=0.1)+
  scale_fill_manual(values=c("arid"="#FFA500","mesic"="#5EA336"))+
  facet_grid(family~CT, scales = "free")+theme_bw()+
  ylab("Proportion")+
  ylab("Condition")+
  theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1),strip.text = element_text(
    size = 6, angle = 90))
   
ggsave(p,file="boxplot_celltypes_Music_selectedCellTypes_Total_family.pdf",width = 20,
  height = 30,units="cm")


ggsave(p,file="boxplot_celltypes_Music_selectedCellTypes_Total_family.png",width = 20,
  height = 30,units="cm")






p<- table_ACPM%>%
  select(c("Row.names","DCT/CNT/PC","Endo","IC","LOH","Podo","Prox_tub", "Stroma","species","cond_season","family","sp_short_name","ID"))%>%
  pivot_longer(-c("species","sp_short_name","cond_season","Row.names","ID","family"),names_to="CT",values_to="prop")%>%
  filter(family%in%c("Cricetidae","Heteromyidae","Muridae","Nesomyidae"))%>%
  ggplot(aes(x=cond_season,y=prop,fill=cond_season))+
  geom_boxplot(outlier.shape = NA)+
  geom_jitter(alpha=0.2,height=0,width=0.1)+
  scale_fill_manual(values=c("arid"="#FFA500","mesic"="#5EA336"))+
  facet_grid(family~CT, scales = "free")+theme_bw()+
  ylab("Proportion")+
  ylab("Condition")+
  theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1),strip.text = element_text(
    size = 6, angle = 90))
   
ggsave(p,file="boxplot_celltypes_Music_selectedCellTypes_Total_family2cond.pdf",width = 20,
  height = 30,units="cm")


ggsave(p,file="boxplot_celltypes_Music_selectedCellTypes_Total_family2cond.png",width = 20,
  height = 20,units="cm")

 theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1),strip.text = element_text(
    size = 6, angle = 90))

theme_minimal

ggsave(p,file="boxplot_celltypes_Music_selectedCellTypes_Total_family2cond.png",width = 20,
  height = 20,units="cm")


p<- table_ACPM%>%
  select(c("Row.names","DCT/CNT/PC","Endo","IC","LOH","Podo","Prox_tub", "Stroma","species","cond_season","family","sp_short_name","ID"))%>%
  pivot_longer(-c("species","sp_short_name","cond_season","Row.names","ID","family"),names_to="CT",values_to="prop")%>%
  filter(family%in%c("Cricetidae","Heteromyidae","Muridae","Nesomyidae"))%>%
  ggplot(aes(x=cond_season,y=prop,fill=cond_season))+
  geom_boxplot(outlier.shape = NA)+
  geom_jitter(alpha=0.2,height=0,width=0.1)+
  scale_fill_manual(values=c("arid"="#FFA500","mesic"="#5EA336"))+
  facet_grid(family~CT, scales = "free")+theme_bw()+
  ylab("Proportion")+
  ylab("Condition")+
  theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1),strip.text = element_text(
    size = 6, angle = 90))+
  stat_compare_means(method = "wilcox.test",aes(label = paste0("p = ", after_stat(p.format))))

ggsave(p,file="boxplot_celltypes_Music_selectedCellTypes_Total_family2cond.pdf",width = 20,
  height = 30,units="cm")


ggsave(p,file="boxplot_celltypes_Music_selectedCellTypes_Total_family2cond.png",width = 20,
  height = 20,units="cm")



table_ACPM=read.table("ref_prop_music_total_rat.txt",sep="\t",h=T)

sapply(c("DCT.CNT.PC","Endo","IC","LOH","Podo","Prox_tub","Stroma"),function(x){p=anova(lm(table_ACPM[,x]~table_ACPM$family+table_ACPM$cond_season))["Pr(>F)"][2,];signif(p,digits=2)})
#DCT/CNT/PC       Endo         IC        LOH       Podo   Prox tub     Stroma 
# OLD   0.0300     0.0064     0.7900     0.9400        NaN     0.0071     0.0030 
# NEW    0.0420     0.0530     0.9700     0.7900     0.4500     0.0530     0.0029 


sapply(c("DCT.CNT.PC","Endo","IC","LOH","Podo","Prox_tub","Stroma"),function(x){p=anova(lm(table_ACPM[,x]~table_ACPM$cond_season))["Pr(>F)"][1,];signif(p,digits=2)})
#DCT.CNT.PC       Endo         IC        LOH       Podo   Prox_tub     Stroma 
# OLD   0.00560    0.12000    0.41000    0.59000        NaN    0.08800    0.00054 
# NEW   0.01100    0.41000    0.66000    0.42000    0.25000    0.27000    0.00038 


## to test on the 4 families with species in both habitats

 table_ACPM_fam=table_ACPM[table_ACPM$family%in%c("Cricetidae","Heteromyidae","Muridae","Nesomyidae"),]

sapply(c("DCT.CNT.PC","Endo","IC","LOH","Podo","Prox_tub","Stroma"),function(x){p=anova(lm(table_ACPM_fam[,x]~table_ACPM_fam$family+table_ACPM_fam$cond_season))["Pr(>F)"][2,];signif(p,digits=2)})

#DCT.CNT.PC       Endo         IC        LOH       Podo   Prox_tub     Stroma 
#    0.0430     0.0340     0.9600     0.7600     0.5000     0.0630     0.0017 


```

### Check quality of deconvolution on pseudobulks: for mouse data

```{r}


scRNA.df <- readRDS("../../data/sou/Park_scRNA.rds")
rownames(scRNA.df) <- scRNA.df$MGI
scRNA.df <- scRNA.df[, -which(names (scRNA.df) %in% c("MGI"))]

#non-unique markers: 
nn <- toupper(rownames(scRNA.df))
nnu <- names(table(nn)[table(nn)>1])
scRNA.df<-scRNA.df[!toupper(row.names(scRNA.df))%in%nnu,]

row.names(scRNA.df)=toupper(row.names(scRNA.df))


scRNA.matrix <- data.matrix(scRNA.df)
sample.ids <- colnames(scRNA.matrix)

phenotable <- read.table(file = "../../data/sou/Park_12_CT_phenoData.txt")

mouse<-CreateSeuratObject(counts=scRNA.matrix, meta.data=phenotable )

mouse <- mouse %>%
  Seurat::SCTransform(
    vst.flavor = "v2",
    verbose = FALSE,
    min_cell = 5,
    vars.to.regress = NULL
  ) %>%
  Seurat::RunPCA(npcs = 50, verbose = FALSE) %>%
  Seurat::RunUMAP(reduction = "pca",
                  dims = 1:50,
                  verbose = FALSE)

Idents(mouse)=mouse$V2

mou <-DimPlot(object = mouse, label = TRUE) +
  coord_fixed()
ggsave(mou, file="mou_umap.png")
ggsave(mou, file="mou_umap.pdf")


# to plot colors corresponding to rat map
ggplotColours <- function(n = 6, h = c(0, 360) + 15){
     if ((diff(h) %% 360) < 1) h[2] <- h[2] - 360/n
     hcl(h = (seq(h[1], h[2], length = n)), c = 100, l = 65)
 }
 
 color_list <- ggplotColours(n=8)
 color_list
names(color_list)=c("Endo","Prox tub","IC","LOH","DCT","Stroma","Immune","Podo")
 
color_list_order=color_list[c("Immune","IC","DCT","Endo","Stroma","LOH","Immune","Immune","Immune","grey","Podo","Prox tub")]
color_list_order[is.na(color_list_order)]="grey"
names(color_list_order)=levels(as.factor(mouse$V2))

mou_col <-DimPlot(object = mouse, label = TRUE,group.by = "V2",cols = color_list_order)+
  coord_fixed()
ggsave(mou_col, file="mouse_umap_identV2colorat.png",width=20,height=20,units="cm")
ggsave(mou_col, file="mouse_umap_identV2colorat.pdf",width=20,height=20,units="cm")


mouse_avg <- AggregateExpression(object = mouse)$RNA


Idents(mouse)=mouse$orig.ident
mouse_avg_tot <- AggregateExpression(object = mouse)$RNA

ide<-data.frame(t(data.frame(strsplit(phenotable$V1,"\\."))))
sample.ids <- paste0("sample",ide$X2)

# individual.ids and cell.types should be in the same order as in sample.ids
individual.labels <- row.names(phenotable)

  
cell.type.labels <- phenotable$V2

sc.pheno <- data.frame(check.names=F, check.rows=F,
                       stringsAsFactors=F,
                       row.names=individual.labels,
                       SubjectName=sample.ids,
                       cellType=cell.type.labels)

sc.meta <- data.frame(labelDescription=c("SubjectName",
                                         "cellType"),
                      row.names=c("SubjectName",
                                  "cellType"))
sc.pdata <- new("AnnotatedDataFrame",
                data=sc.pheno,
                varMetadata=sc.meta)


sc.eset2 <- SingleCellExperiment(list(counts=scRNA.matrix),
    colData=DataFrame(cellTypeID=sc.pheno$cellType,cellType=sc.pheno$cellType,sampleID=sc.pheno$SubjectName,SubjectName=sc.pheno$SubjectName))

bulk.mtx
pseudobulk.matrix <- data.matrix(cbind(round(mouse_avg_tot,0),round(mouse_avg_tot,0)))
pseudobulk.matrix[is.na(pseudobulk.matrix)] <- 0 
pseudobulk.eset <- ExpressionSet(assayData = pseudobulk.matrix)
bulk.mtx = exprs(pseudobulk.eset)

## With selected CTs for pseudobulk
prop_music_pseudomouse=music_prop(bulk.mtx = bulk.mtx, sc.sce = sc.eset2, 
                      clusters = 'cellType', samples = 'sampleID',
                      select.ct = c("CD","DCT","Endo","Fib","LOH","Novel1","Podo","PT"), 
                      verbose = F)$Est.prop.weighted
write.table(file="prop_music_pseudobulk_mouse.txt",prop_music_pseudomouse,quote=F)


trueprop<-table(mouse$V2[mouse$V2%in%c("CD","DCT","Endo","Fib","LOH","Novel1","Podo","PT")])
trueprop<-trueprop/sum(trueprop)
typecell=names(trueprop)

mousebulkprop<-data.frame(trueprop=trueprop)
mousebulkprop$estprop=prop_music_pseudomouse[1,]
mousebulkprop$`cell type`=typecell
g<-ggplot(mousebulkprop,aes(trueprop,estprop,col=`cell type`))+geom_point()+xlab("true proportion")+ylab("estimated proportion") + theme_minimal()
ggsave(g,file="proportions_frompseudomouse_trueest.png")
ggsave(g,file="proportions_frompseudomouse_trueest.pdf")

cor.test(mousebulkprop$trueprop.Freq,mousebulkprop$estprop,meth="spearman")
#0.9047619
cor.test(ratbulkprop$trueprop.Freq,ratbulkprop$estprop,meth="spearman")
#0.9642857 



mouseexp<-bulk.df_Total_sort[,coldataTotal$species=="Mus_musculus"]
mouseexp$averagemouse<-apply(mouseexp,1,mean)
mouseexp_dat<-merge(mouseexp,mouse_avg_tot,by.x=0,by.y=0)
g<-ggplot(mouseexp_dat,aes(log(averagemouse),log(V1),alpha=0.3))+geom_point()+xlab("log(bulk expression)")+ylab("log(pseudo bulk expression)")
ggsave(g,file="expression_bulk_pseudo_mouse.png")
ggsave(g,file="expression_bulk_pseudo_mouse.pdf")

cor.test(mouseexp_dat$averagemouse,mouseexp_dat$V1,meth="spearman")
#0.6877439 




ratexp<-bulk.df_Total_sort[,coldataTotal$species=="Rattus_norvegicus"]
ratexp$averagerat<-apply(ratexp,1,mean)
ratexp_dat<-merge(ratexp,rat_avg_tot,by.x=0,by.y=0)
g<-ggplot(ratexp_dat,aes(log(averagerat),log(V1),alpha=0.3))+geom_point()+xlab("log(bulk expression)")+ylab("log(pseudo bulk expression)")
ggsave(g,file="expression_bulk_pseudo_rat.png")
ggsave(g,file="expression_bulk_pseudo_rat.pdf")

cor.test(ratexp_dat$averagerat,ratexp_dat$V1,meth="spearman")
#0.5698054 


```

### Project hamster data

```{r}
# toupper gene symbol 10x
feat = "../../data/ham/filtered/features.tsv.gz"
read.delim(feat, header = FALSE) %>%
      dplyr::mutate(dplyr::across(everything(), ~ toupper(.))) %>%
      readr::write_tsv(file = feat, col_names = FALSE)



mat<-Read10X("../../data/ham/filtered/")
ham<-CreateSeuratObject(counts=mat)


mt_gene <- mtham
mt_gene <- mt_gene[mt_gene %in% row.names(ham)]

rb_gene <- rps_gene_ham
rb_gene <- rb_gene[rb_gene %in% row.names(ham)]


ham[["percent.mt"]] <-
  Seurat::PercentageFeatureSet(
    object = ham,
    features = mt_gene,
    assay = "RNA",
    col.name = "percent"
  )$percent
ham[["percent.rb"]] <-
  Seurat::PercentageFeatureSet(
    object = ham,
    features = rb_gene,
    assay = "RNA",
    col.name = "percent"
  )$percent


s=apply(ham@assays$RNA$counts[mt_gene,],1,sum)
hist(s,nclass=100)
s=apply(ham@assays$RNA$counts[mt_gene,],2,sum)
hist(s,nclass=100)
n=apply(ham@assays$RNA$counts,2,sum)
hist(n,nclass=100)

pdf("prop_mito.pdf")
hist(s/n,nclass=100,main="mtRNA/ncounts")
dev.off()

png("prop_mito.png")
hist(s/n,nclass=100,main="mtRNA/ncounts")
dev.off()


ham <- subset(ham,
              subset = nCount_RNA > 200 &
                nCount_RNA <  median(ham$nCount_RNA) * 10 &
                percent.mt < 80) %>%
  Seurat::SCTransform(
    vst.flavor = "v2",
    verbose = FALSE,
    min_cell = 5,
    vars.to.regress = NULL
  ) %>%
  Seurat::RunPCA(npcs = 50, verbose = FALSE) %>%
  Seurat::RunUMAP(reduction = "pca",
                  dims = 1:50,
                  verbose = FALSE) %>%
  Seurat::FindNeighbors(reduction = "pca",
                        dims = 1:50,
                        verbose = FALSE) %>%
  Seurat::FindClusters(resolution = 0.8, verbose = FALSE)

ham <- CellCycleScoring(ham,
                        s.features = cc.genes.updated.2019$s.genes,
                        g2m.features = cc.genes.updated.2019$g2m.genes)

ham1 <- DimPlot(object = ham, label = TRUE) +
  coord_fixed()
ggsave(ham1, file="ham_umap1.png")
ggsave(ham1, file="ham_umap1.pdf")

ham2 <- DimPlot(object = ham,
        group.by = "Phase",
        shuffle = TRUE) +
  coord_fixed()
ggsave(ham2, file="ham_umap2.png")
ggsave(ham2, file="ham_umap2.pdf")


ham.query <- ham

 
map_query <- function(obj.query, obj.ref, meta_col = "V2") {
  anchors <- FindTransferAnchors(
    reference = obj.ref,
    query = obj.query,
    dims = 1:30,
    reference.reduction = "pca",
    normalization.method = "SCT",
    recompute.residuals = FALSE,
    features = intersect(row.names(obj.ref), row.names(obj.query))
  )
  
  obj.ref <- RunUMAP(
    obj.ref,
    dims = 1:30,
    reduction = "pca",
    return.model = TRUE
  )
  
  obj.query <- MapQuery(
    anchorset = anchors,
    reference = obj.ref,
    query = obj.query,
    refdata = list(celltype = "V2"),
    reference.reduction = "pca",
    reduction.model = "umap"
  )
  
  p1 <- DimPlot(
    obj.ref,
    reduction = "umap",
    group.by = meta_col,
    label = TRUE,
    label.size = 3,
    repel = TRUE
  ) + NoLegend() + ggtitle("Reference annotations")
  
  p2 <- DimPlot(
    obj.query,
    reduction = "ref.umap",
    group.by = "predicted.celltype",
    label = TRUE,
    label.size = 3,
    repel = TRUE
  ) + NoLegend() + ggtitle("Query transferred labels")
  
  p3 <-  FeaturePlot(obj.query, reduction = "ref.umap", features = "predicted.celltype.score")
  print(table(obj.query$predicted.celltype, obj.query$seurat_clusters))
  
  print(p1)
  print(p2)
  print(p3)
  
  ggsave(p1, file="mouse_ref.png")
ggsave(p1, file="mouse_ref.pdf")

ggsave(p2, file="ham_query.png")
ggsave(p2, file="ham_query.pdf")

ggsave(p3, file="ham_queryscore.png")
ggsave(p3, file="ham_queryscore.pdf")

  obj.query
  
}

ham_avg <- map_query(obj.query = ham.query, obj.ref = mouse.ref, meta_col = "V2")
Idents(ham_avg)=ham_avg$predicted.celltype


ham_avg_exp <- AggregateExpression(object = ham_avg)$RNA

mus_avg_exp <- AggregateExpression(object = mus_avg)$RNA


markers <- as.data.frame(read_xlsx(path = "../../data/markers/cell-types-markers_associated-biblio.xlsx",sheet ="All-markers"))
markers$MGI<-toupper(markers$`Genes MGI`)
#non-unique markers: 
nnu <- names(table(markers$MGI)[table(markers$MGI)==2])
markers<-markers[!markers$MGI%in%nnu,]
markers<-markers[!is.na(markers$Segment),]
markers$Segment[markers$Segment=="-"]="Immune"
# Markers annotated with "-" correspond to Immune cells. We remove markers annotated with NA as they are not identified

## DE genes in murine
de <- read.table("../../data/expr_table/env_seasonal_MvsX_Murinae_invert.tsv",sep="\t",h=T)



ham_avg_mk <- merge(ham_avg_exp,markers,by.x=0,by.y="MGI")
rownames(ham_avg_mk)<-ham_avg_mk$Row.names
ham_avg_mk = with(ham_avg_mk,  ham_avg_mk[order(ham_avg_mk$Segment) , ])

de$X2<-toupper(de$X)
datade <- merge(ham_avg_mk,de,by.x="Row.names",by.y="X2",all.x=T)
datade$xm<-ifelse(datade$log2FoldChange_OK>0,"xeric","mesic")
datade$xm[is.na(datade$log2FoldChange_OK)]="NA"
datade = datade[order(paste0(datade$Segments,datade$xm)),]


gene_functions_df <- data.frame(Segment = ham_avg_mk$Segment)
row.names(gene_functions_df) <- rownames(ham_avg_mk)

newCols <- paletteer_d("ggsci::springfield_simpsons")[1:length(unique(ham_avg_mk$Segment))]
annoRow <- newCols

ann_colors = list(Segments = annoRow)


heat_plot <- pheatmap(ham_avg_mk[,2:9],
scale="row",
cluster_rows = FALSE,
cluster_cols = TRUE, # set to FALSE if you want to remove the dendograms
fontsize_row = 10, # row label font size
fontsize_col = 7, # column label font size
angle_col = 45, # sample names at an angle
show_colnames = T, show_rownames = F, # displaying column and row names,
annotation_row = gene_functions_df,
#annotation_colors = ann_colors,
main = "Marker gene expression in Rat scRNAseq data") # a title for our heatmap
pdf("Ham_markers_per_celltype.pdf", height = 20, width = 10)
heat_plot
dev.off()

pdf("Ham_markers_per_celltype.pdf", height = 7, width = 5)
heat_plot
dev.off()
  
  

## with DE genes
ann_df <- data.frame(CT = colnames(datade)[2:11])
row.names(ann_df) <- colnames(datade)[2:11]
head(ann_df)

gene_functions_df <- data.frame(Segments = datade$Segments,  Condition=datade$xm)
row.names(gene_functions_df) <- rownames(datade)
# Markers annotated with "-" correspond to Immune cells. We remove markers annotated with NA as they are not identified
gene_functions_df$Segments[gene_functions_df$Segments=="-"]="Immune"
gene_functions_df<-gene_functions_df[!is.na(gene_functions_df$Segment),]
datade<-datade[!is.na(datade$Segment),]
datade1<-datade[datade$Segments!="Immune"&datade$Segments!="-",]
datade1<-datade1[,c("Row.names","CD","DCT","Endo","Fib","LOH","PT","Segments")]

newCols <- paletteer_d("ggsci::springfield_simpsons")[1:length(unique(datade1$Segments))]
annoRow <- newCols
names(annoRow) <- unique(datade1$Segment)
newCols <- paletteer_d("ggthemes::stata_s2color")[1:6] 
annoCol <- newCols
names(annoCol) <- names(datade1)[2:7]

ann_colors = list(
  Segments = annoRow, CT = annoCol, Condition = c(xeric="orange", mesic="green","NA"="white"))


heat_plot <- pheatmap(datade1[,2:7], 
                      scale="row",
          annotation_colors = ann_colors,
                      cluster_rows = FALSE, 
                      cluster_cols = FALSE, # set to FALSE if you want to remove the dendograms
                      annotation_row = gene_functions_df, # row (gene) annotations
                      annotation_col = ann_df, # column (sample) annotations
                   
                      annotation_names_row = F, 
                      annotation_names_col = F,
                      fontsize_row = 10,          # row label font size
                      fontsize_col = 7,          # column label font size 
                      angle_col = 45, # sample names at an angle
                      #legend_breaks = c(-2, 0, 2), # legend customisation
                      #legend_labels = c("Low", "Medium", "High"), # legend customisation
                      show_colnames = T, show_rownames = F, # displaying column and row names
                      main = "Marker gene expression in hamster data") # a title for our heatmap
  
pdf("Marker_ham_deconvolution_plot_short.pdf", height = 10, width = 8)
heat_plot
dev.off()




```

### deconvolution with mouse scrnaseq data

```{r}

# Prepare input bulk data for deconvolution 
bulk.df <- read.table("../../data/expr_table/datalocalandonline_counts.tsv",h=T,row.names=1)
rownames(bulk.df) <- toupper(rownames(bulk.df))


# to select Murinae species
#coldata <- read.table("../../data/expr_table/datalocalandonline_coldata.tsv",h=T,sep="\t")
#coldataMurinae <- coldata[coldata$exp_murinae=="yes",]
#coldataMurinae$ID <- gsub("kollmannspergeri","kollmanspergeri",coldataMurinae$ID)
#bulk.df_Murinae <- bulk.df[,coldataMurinae$ID]
#coldataMurinae <- coldataMurinae[-grep("kollmannspergeri",coldataMurinae$ID),]
#bulk.df_Murinae_sort <- bulk.df_Murinae[,coldataMurinae$ID]


bulk.matrix <- data.matrix(round(bulk.df_Total_sort,0))
bulk.matrix[is.na(bulk.matrix)] <- 0 
bulk.eset <- ExpressionSet(assayData = bulk.matrix)

row.names(coldataTotal)=coldataTotal$ID
metadata <- data.frame(labelDescription= names(coldataTotal)[1:20], row.names=names(coldataTotal)[1:20])

bulk.eset <- ExpressionSet(assayData = data.matrix(bulk.matrix),
                           phenoData = new("AnnotatedDataFrame", data = coldataTotal[,1:20], 
                                           varMetadata = metadata) )

bulk.mtx<-exprs(bulk.eset )
## With selected CTs

prop_music_mouse=music_prop(bulk.mtx = bulk.mtx, sc.sce = sc.eset2, 
                      clusters = 'cellType', samples = 'sampleID',
                      select.ct = c("CD","DCT","Endo","Fib","LOH","Novel1","Podo","PT"), 
                      verbose = F)$Est.prop.weighted
write.table(file="prop_music_mouse.txt",prop_music,quote=F)
## With all CTs

prop_music2_mouse=music_prop(bulk.mtx = bulk.mtx, sc.sce = sc.eset2, 
                      clusters = 'cellType', samples = 'sampleID',
                      verbose = F)$Est.prop.weighted
## With all CTs + marker genes only

#prop_music3_markers=music_prop(bulk.mtx = bulk.mtx, sc.sce = sc.eset2, markers=markers$`Genes MGI`,
#                      clusters = 'cellType', samples = 'sampleID',
#                      verbose = F)$Est.prop.weighted

#prop_music3=as.data.frame(prop_music3_markers)


## PCA with selected cell types
ref.Total_Music <- merge(prop_music_mouse,coldataTotal,by.x=0,by.y="ID")
#ref.MurinaeM_Music$DCT=ref.MurinaeM_Music$DCT+ref.MurinaeM_Music$LOH
table_ACPM = ref.Total_Music
table_ACPM$ID=paste(table_ACPM$sp_short_name,row.names(table_ACPM))
rownames(table_ACPM)=table_ACPM$ID
res_ACPM = table_ACPM %>%  PCA(quali.sup = c(1,10:42), graph = FALSE)
var <- get_pca_var(res_ACPM)

pdf("Total_Biplot_per_cond_season_PCA_Music_selectedCellTypes_mouse_ind.pdf")
fviz_pca_biplot (res_ACPM,col.var = "black", col.ind = as.factor(table_ACPM$cond_season), 
             addEllipses = FALSE,  alpha.ind=0.8, 
             palette = group.colors, legend.title = "Groups",
            ggtheme 
            = theme_bw())
dev.off()


#https://www.cfholbert.com/blog/cluster-pca/
pca1 <- dudi.pca(table_ACPM[,2:9], scan = FALSE, nf = 3)
rand1 <- rtest(discrimin(pca1, as.factor(table_ACPM$cond_season), scan = FALSE), 1000)
rand1
pdf("Total_effect_season_mctest_mouse.pdf")
plot(rand1, main = "Monte-Carlo test")
dev.off()

rand1 <- rtest(discrimin(pca1, as.factor(table_ACPM$family), scan = FALSE), 1000)
rand1
pdf("Total_effect_family_mctest_mouse.pdf")
plot(rand1, main = "Monte-Carlo test")
dev.off()


### represent percentage of cell types

p<- table_ACPM%>%
  select(c("Row.names","CD","DCT","Endo","Fib","LOH","Novel1","Podo","PT","species","cond_season","family","sp_short_name","ID"))%>%
  pivot_longer(-c("species","sp_short_name","cond_season","Row.names","ID","family"),names_to="CT",values_to="prop")%>%
  ggplot(aes(x=cond_season,y=prop,fill=cond_season))+
  geom_boxplot(outlier.shape = NA)+
  geom_jitter(alpha=0.2,height=0,width=0.1)+
  scale_fill_manual(values=c("xeric"="#FFA500","mesic"="#5EA336"))+
  facet_wrap(~CT, scales = "free")+theme_bw()+
  ylab("Proportion")+
  ylab("Condition")+
  theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1),strip.text = element_text(
    size = 6, angle = 90))
   
ggsave(p,file="boxplot_celltypes_Music_selectedCellTypes_Total_mouse.pdf",width = 10,
  height = 20,units="cm")




p<- table_ACPM%>%
  select(c("Row.names","CD","DCT","Endo","Fib","LOH","Novel1","Podo","PT","species","cond_season","family","sp_short_name","ID"))%>%
  pivot_longer(-c("species","sp_short_name","cond_season","Row.names","ID","family"),names_to="CT",values_to="prop")%>%
  ggplot(aes(x=cond_season,y=prop,fill=cond_season))+
  geom_boxplot(outlier.shape = NA)+
  geom_jitter(alpha=0.2,height=0,width=0.1)+
  scale_fill_manual(values=c("xeric"="#FFA500","mesic"="#5EA336"))+
  facet_row(~CT, scales = "free")+theme_bw()+
  ylab("Proportion")+
  ylab("Condition")+
  theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1),strip.text = element_text(
    size = 6, angle = 90))
   
ggsave(p,file="boxplot_celltypes_Music_selectedCellTypes_Total_mouse.pdf",width = 30,
  height = 10,units="cm")





p<- table_ACPM%>%
  select(c("Row.names","CD","DCT","Endo","Fib","LOH","Novel1","Podo","PT","species","cond_season","family","sp_short_name","ID"))%>%
  pivot_longer(-c("species","sp_short_name","cond_season","Row.names","ID","family"),names_to="CT",values_to="prop")%>%
  ggplot(aes(x=cond_season,y=prop,fill=cond_season))+
  geom_boxplot(outlier.shape = NA)+
  geom_jitter(alpha=0.2,height=0,width=0.1)+
  scale_fill_manual(values=c("xeric"="#FFA500","mesic"="#5EA336"))+
  facet_grid(family~CT, scales = "free")+theme_bw()+
  ylab("Proportion")+
  ylab("Condition")+
  theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1),strip.text = element_text(
    size = 6, angle = 90))
   
ggsave(p,file="boxplot_celltypes_Music_selectedCellTypes_Total_family_mouse.pdf",width = 20,
  height = 30,units="cm")


ggsave(p,file="boxplot_celltypes_Music_selectedCellTypes_Total_family_mouse.png",width = 20,
  height = 30,units="cm")






p<- table_ACPM%>%
  select(c("Row.names","CD","DCT","Endo","Fib","LOH","Novel1","Podo","PT","species","cond_season","family","sp_short_name","ID"))%>%
  pivot_longer(-c("species","sp_short_name","cond_season","Row.names","ID","family"),names_to="CT",values_to="prop")%>%
  filter(family%in%c("Cricetidae","Heteromyidae","Muridae","Nesomyidae"))%>%
  ggplot(aes(x=cond_season,y=prop,fill=cond_season))+
  geom_boxplot(outlier.shape = NA)+
  geom_jitter(alpha=0.2,height=0,width=0.1)+
  scale_fill_manual(values=c("xeric"="#FFA500","mesic"="#5EA336"))+
  facet_grid(family~CT, scales = "free")+theme_bw()+
  ylab("Proportion")+
  ylab("Condition")+
  theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1),strip.text = element_text(
    size = 6, angle = 90))
   
ggsave(p,file="boxplot_celltypes_Music_selectedCellTypes_Total_family2cond_mouse.pdf",width = 20,
  height = 30,units="cm")


ggsave(p,file="boxplot_celltypes_Music_selectedCellTypes_Total_family2cond_mouse.png",width = 20,
  height = 20,units="cm")



anova(lm(table_ACPM$Endo~table_ACPM$family+table_ACPM$cond_season))
#0.0048787
sapply(c("CD","DCT","Endo","Fib","LOH","Novel1","Podo","PT"),function(x){p=anova(lm(table_ACPM[,x]~table_ACPM$family+table_ACPM$cond_season))["Pr(>F)"][2,];signif(p,digits=2)}
+ )

#     CD     DCT    Endo     Fib     LOH  Novel1    Podo      PT 
#3.7e-01 5.9e-01 4.9e-03 1.3e-02 4.2e-01 1.3e-01 1.2e-05 3.7e-03 






# Compare cell types inferred for mouse and rat to pseudobulk counts and inferences on pseudobulks


bulkprop_mouserat <- table_ACPM%>%filter(species%in%c("Mus_musculus","Rattus_norvegicus"))
compareprop<-trueprop_pseudo_rat%>%left_join(prop_music_pseudorat,by=c("Var2"="sample","Var1"="name"))

compareprop_rat_long<-compareprop%>%
  pivot_longer(-c(Var1,Var2)) %>%
  rename("Var1"="Celltype") %>%
  rename("Var2" = "Sample") %>%
  rename("name" = "Type") %>%
  relocate("Celltype", .after = "Sample") %>%
  mutate(Type=fct_recode(Type, sc_count = "Freq", sc_deconvolution = "value"))%>%
  mutate(Celltype=fct_recode(Celltype, Prox_tub = "Prox tub"))

  
bulkprop_mouserat_long<-bulkprop_mouserat%>%
  select("DCT/CNT/PC","Endo","IC","LOH","Podo","Prox_tub", "Stroma","ID")%>%
  pivot_longer(-ID) %>%
  rename("name"="Celltype") %>%
  rename("ID" = "Sample") %>%
  mutate("Type" = case_when(str_detect(Sample,"Mmu") ~ "bulk deconvolution Mus",str_detect(Sample,"Rno") ~ "bulk deconvolution Rat"))

  
compareprop_rat_bulk_sc <- bulkprop_mouserat_long%>%bind_rows(compareprop_rat_long)



p<- compareprop_rat_bulk_sc%>%
    ggplot(aes(x=Celltype,y=value,color=Type,group_by=Sample))+
  geom_jitter(alpha=0.4,height=0,width=0.1) +
  theme_bw()+
  ylab("Cell-type proportion")
ggsave(p,file="CellProp_comparisons_musrat.pdf",width = 20,
  height = 15,units="cm")
ggsave(p,file="CellProp_comparisons_musrat.png",width = 20,
  height = 15,units="cm")


compareprop_rat_bulk_sc_selec <- compareprop_rat_bulk_sc%>%
  filter(Celltype%in%c("Endo","LOH","Prox_tub", "Stroma"))%>%
  group_by(Sample,Type)%>%
  mutate(sumprop=sum(value))%>%
   mutate(value2=value/sumprop)
  
p2<- compareprop_rat_bulk_sc_selec%>%
    ggplot(aes(x=Celltype,y=value2,color=Type,group_by=Sample))+
  geom_jitter(alpha=0.4,height=0,width=0.1) +
  theme_bw()+
  ylab("Cell-type proportion")
ggsave(p2,file="CellProp_comparisons_musrat_selec.pdf",width = 20,
  height = 15,units="cm")
ggsave(p2,file="CellProp_comparisons_musrat_selec.png",width = 20,
  height = 15,units="cm")



fig_deconvo_qc=plot_grid(p + theme(legend.position="none"), p2, labels = c('A', 'B'),rel_widths = c(1,1.1))
ggsave(fig_deconvo_qc, file="fig_deconvo_qc.png",width=25,height=15,units="cm")
ggsave(fig_deconvo_qc, file="fig_deconvo_qc.pdf",width=25,height=15,units="cm")
  
  




```

\`\`\`
