---
title: "Figure 4"
author: "Sang Liu"
date: '2023-10-14'
output: html_document
---

```{r setup, include=FALSE}
knitr::opts_chunk$set(echo = TRUE)
```

## Load packages

```{r}
library(ggplot2)
library(Seurat)
library(fastcluster)
library(ggpubr)
```

## Figure 4C

```{r}
load("~/arbor_features")
arbornames<-rownames(allenarbors_dataframe)[which(allenarbors_dataframe$AxonSubtype %in% c(45:64)&allenarbors_dataframe$annot %in% calannot("TH"))]
statistics<-allenarbors_dataframe[arbornames,c("branchpoints","cable.length","Volume","Density","Lx","Ly","Lz","LxvsLy","LzvsLxy","Mean_Segment_Length","maxbranchorder",
                                               "terminal_branch_length","tortuosity","Mean_PartitionAsymmetry_Degree","Ratio_AsymmetryBranch","Mean_BifurAngle","node_root_maxEucDist","node_root_meanEucDist","Fractal_Dim")]
delete<-which(!colnames(allenarbors_dataframe) %in% colnames(statistics))
statistics_labels<-allenarbors_dataframe[arbornames,delete]
statistics$Mean_PartitionAsymmetry_Degree[which(is.na(statistics$Mean_PartitionAsymmetry_Degree))]<-0
statistics$Mean_BifurAngle[which(is.na(statistics$Mean_BifurAngle))]<-0
statistics_labels<-statistics_labels[apply(statistics,1,function(x)all(!is.na(x)&x!=Inf)),]
statistics<-statistics[apply(statistics,1,function(x)all(!is.na(x)&x!=Inf)),]

targets<-c("RT","SPA","PP",calannot("TH",3))
for(i in targets){
  statistics_labels$annot[which(statistics_labels$annot %in% calannot(i))]<-i
}
statistics_labels$Type[which(statistics_labels$AxonSubtype %in% 45:52)]<-"CT"
statistics_labels$Type[which(statistics_labels$AxonSubtype %in% 53:64)]<-"PT"

#remove highly correlated features
cormatx=matrix(0,ncol(statistics),ncol(statistics))
for(i in 1:(ncol(statistics)-1)){
  for(j in (i+1):ncol(statistics)){
    cormatx[i,j]=cor(statistics[,i],statistics[,j])
  }
}
statistics[,unique(which(abs(cormatx)>0.9,arr.ind = T)[,1])]<-NULL

samples<-1:nrow(statistics)
seurat<-CreateSeuratObject(t(statistics[samples,]))
seurat<-ScaleData(seurat)
seurat<-FindVariableFeatures(seurat)
seurat <- RunPCA(seurat)
seurat <- RunUMAP(seurat,dims = 1:(nrow(seurat)-1))
seurat <- RunTSNE(seurat,dims = 1:(nrow(seurat)-1))
seurat@meta.data$Type<-ifelse(statistics_labels$AxonSubtype %in% 45:52,"CT",ifelse(statistics_labels$AxonSubtype %in% 53:64,"PT","TC"))
seurat@meta.data$neuronsubtype<-statistics_labels$AxonSubtype[samples]
seurat@meta.data$annots<-statistics_labels$annot[samples]
seurat@meta.data$annots<-factor(seurat@meta.data$annots,levels=targets[which(targets %in% unique(seurat@meta.data$annots))])
seurat@meta.data$somaposition<-statistics_labels$somaposition[samples]

out.dist<-dist(seurat@reductions$umap@cell.embeddings[,1:2]) 
out.hclust<-fastcluster::hclust(out.dist,method="ward.D")
seurat@meta.data$Subclass<-cutree(out.hclust,k=4,order_clusters_as_data = F)[order(out.hclust$order)]
```

```{r,fig.cap="Figure 4C (right)",fig.align='center',fig.width=3,fig.height=3}
UMAPPlot(seurat,group.by = "Subclass",pt.size=1,raster=T)+
  theme_classic()+
  theme(plot.title = element_blank(),
        text=element_text(size=5),
        axis.text.x = element_blank(),
        axis.text.y =element_blank(),
        axis.title.x =element_blank(),
        axis.title.y =element_blank(),
        axis.ticks = element_blank(),
        axis.line = element_blank(),
        strip.background = element_blank(),
        strip.text.x = element_text(angle=0,size=5),
        legend.title = element_text(hjust = 0.5,size=7),
        legend.text = element_text(hjust = 0.5,size=5),
        legend.key.width =unit(0.1,"mm"),
        legend.key.height =unit(0.1,"mm"),
        legend.text.align=0)

ggsave(paste0(FIGUREDIR,"Figure 4C.pdf"),device = "pdf",width = 60,height = 50,units = "mm")
```

## Figure 4E

```{r,fig.cap="Figure 4E (left)",fig.align='center',fig.width=3,fig.height=3}
type<-c("CT","PT")
arborsfreq=rep(NA,length(unique(seurat@meta.data$Subclass))*length(type))
P<-rep(NA,length(unique(seurat@meta.data$Subclass))*length(type))
n<-1
for(i in unique(seurat@meta.data$Subclass)){
  for(j in type){
    arborsfreq[n]=length(which(seurat@meta.data$Subclass==i&seurat@meta.data$Type==j))/length(which(seurat@meta.data$Type==j))
    ma<-matrix(nrow=2,ncol=2)
    ma[1,1]<-length(which(seurat@meta.data$Subclass==i&seurat@meta.data$Type==j))
    ma[2,1]<-length(which(seurat@meta.data$Subclass!=i&seurat@meta.data$Type==j))
    ma[1,2]<-length(which(seurat@meta.data$Subclass==i&seurat@meta.data$Type!=j))
    ma[2,2]<-length(which(seurat@meta.data$Subclass!=i&seurat@meta.data$Type!=j))
    if(ma[1,1]/ma[1,2]>ma[2,1]/ma[2,2]){
      P[n]<-fisher.test(ma)$p.value
    }else{
      P[n]<-1
    }
    
    n<-n+1
  }
}
P<-p.adjust(P,"fdr")
marker<-rep("",length(P))
marker[which(P<=0.05)]<-"*"
marker[which(P<=0.01)]<-"**"
marker[which(P<=0.001)]<-"***"
marker[which(P<=0.0001)]<-"****"
ggplot(data.frame(Proportion=arborsfreq,Subclass=as.factor(rep(unique(seurat@meta.data$Subclass),each=length(type))),Type=as.factor(rep(type,times=length(unique(seurat@meta.data$Subclass))))), aes(x = Type, y = Proportion,fill=Subclass)) +
  geom_bar(stat = "identity", position = "stack")+
  geom_text(aes(x=Type,y=Proportion,label=marker),size=1.5,position= "stack",stat = "identity")+
  scale_y_continuous(n.breaks = 2)+
  theme_classic()+
  theme(plot.title = element_blank(),
        text=element_text(size=5),
        axis.text.x = element_text(size=7),
        axis.text.y =element_text(size=5),
        axis.title.x =element_text(size=7),
        axis.title.y =element_text(size=7),
        strip.background = element_blank(),
        strip.text.x = element_text(angle=0,size=5),
        legend.title = element_text(hjust = 0.5,size=7),
        legend.text = element_text(hjust = 0.5,size=5),
        legend.key.size = unit(1,"mm"),
        legend.text.align=0)
ggsave(paste0(FIGUREDIR,"Figure 4E"),width = 40,height = 30,units = "mm")
```

```{r,fig.cap="Figure 4E (right)",fig.align='center',fig.width=3,fig.height=3}
Neurontype<-c("PT","CT")
data<-as.data.frame(t(as.matrix(seurat@assays$RNA@data)))
data$Type<-seurat@meta.data$Type
data<-reshape2::melt(data[,],id.vars = c("Type"))
features<-c("branchpoints","terminal-branch-length","Fractal-Dim")
data<-subset(data,variable %in% features)
colnames(data)[2]<-"features"
data$value<-as.numeric(data$value)
data$value[which(data$features %in% c("cable.length","terminal-branch-length"))]<-
  data$value[which(data$features %in% c("cable.length","terminal-branch-length"))]/1000
data$value[which(data$features %in% c("Volume"))]<-data$value[which(data$features %in% c("Volume"))]/1000000000

subdata<-data
ggviolin(
  subdata,
  x="Type",
  y = "value",
  add = "mean",
  xlab = "morphological features",
  size = 0.3,
  width = 0.9,
  add.params = list(size = 0.01),
  trim = T,
  scale = "width",
  color = "Type"
) +
  scale_y_log10(n.breaks = 4)+
  facet_wrap(
    "features",
    scales = "free",
    ncol = length(unique(data$features)),
    strip.position = "top",
) +
  stat_compare_means(
    comparisons = list(c("CT","PT")),
    label="p.signif",
    size = 2,
    method = "wilcox.test"
  ) +
  scale_color_manual(values=c("blue","red"))+
  labs(x = "Type") +
  theme(
    strip.background = element_blank(),
    strip.text = element_text(size = 7),
    axis.text.x = element_text(size = 5),
    axis.text.y = element_text(size = 5),
    axis.title.x = element_text(size = 0),
    axis.title.y = element_blank(),
    panel.spacing.y = unit(0, "mm"),
    legend.position = "none",
    legend.title = element_text(hjust = 0.5,size=5),
    legend.text = element_text(hjust = 0.5,size=5),
    legend.key.height=unit(0.1,"mm"),
    legend.key.width = unit(1,"mm"),legend.text.align=0,
    legend.box.spacing = unit(0.1,"mm")
)
```

## Figure 4F

```{r}
targets<-c("RT","SPA","PP",calannot("TH",3))
arborsfreq<-rep(NA,length(unique(seurat@meta.data$Subclass))*length(targets))
P<-rep(NA,length(unique(seurat@meta.data$Subclass))*length(targets))
arbornum=rep(NA,length(unique(seurat@meta.data$Subclass))*length(targets))
n<-1
for(i in unique(seurat@meta.data$Subclass)){
  for(j in targets){
    arborsfreq[n]=length(which(seurat@meta.data$Subclass==i&seurat@meta.data$annots==j))/length(which(seurat@meta.data$annots==j))
    arbornum[n]=length(which(seurat@meta.data$annots==j))
    ma<-matrix(nrow=2,ncol=2)
    ma[1,1]<-length(which(seurat@meta.data$Subclass==i&seurat@meta.data$annots==j))
    ma[2,1]<-length(which(seurat@meta.data$Subclass!=i&seurat@meta.data$annots==j))
    ma[1,2]<-length(which(seurat@meta.data$Subclass==i&seurat@meta.data$annots!=j))
    ma[2,2]<-length(which(seurat@meta.data$Subclass!=i&seurat@meta.data$annots!=j))
    P[n]<-fisher.test(ma,alternative = "greater")$p.value
    n<-n+1
  }
}
tmp<-data.frame(Proportion=arborsfreq,Num=arbornum,Subclass=as.factor(rep(unique(seurat@meta.data$Subclass),each=length(targets))),Region=factor(rep(targets,times=length(unique(seurat@meta.data$Subclass)))),P=P)
tmp<-subset(tmp,!is.na(Proportion)&Num>0)
tmp$Region<-factor(tmp$Region,levels=targets[which(targets %in% unique(tmp$Region))])
tmp$P<-p.adjust(tmp$P,"fdr")
tmp$marker<-rep("",length(tmp$P))
tmp$marker[which(tmp$P<=0.05)]<-"*"
tmp$marker[which(tmp$P<=0.01)]<-"**"
tmp$marker[which(tmp$P<=0.001)]<-"***"
tmp$marker[which(tmp$P<=0.0001)]<-"****"

ggplot(tmp, aes(x = Region, y = Proportion,fill=Subclass)) +
  geom_bar(stat = "identity", position = "stack")+
  geom_text(aes(x=Region,y=Proportion,label=marker),size=1.5,position= "stack",stat = "identity")+
  scale_y_continuous(n.breaks = 2)+
  theme_classic()+
  theme(plot.title = element_blank(),
        text=element_text(size=5),
        axis.text.x = element_text(size=5,angle=60,hjust=1),
        axis.text.y =element_text(size=5),
        axis.title.x =element_text(size=7),
        axis.title.y =element_text(size=7),
        strip.background = element_blank(),
        strip.text.x = element_text(angle=0,size=5),
        legend.title = element_text(hjust = 0.5,size=7),
        legend.text = element_text(hjust = 0.5,size=5),
        legend.key.size = unit(1,"mm"),
        legend.text.align=0)
ggsave(paste0(FIGUREDIR,"Figure 4F.pdf"),width = 100,height = 30,units = "mm")
```
