---
title: "R Notebook"
output: html_notebook
---

This is an [R Markdown](http://rmarkdown.rstudio.com) Notebook. When you execute code within the notebook, the results appear beneath the code. 

Try executing this chunk by clicking the *Run* button within the chunk or by placing your cursor inside it and pressing *Cmd+Shift+Enter*. 

```{r}
library(readxl)
library(ggplot2)
spot_counting <- read_excel("~/Desktop/scRNA-seq_snRNA-seq_adipocyte_lineages/revisions/spot_counting.xlsx",col_names = T)
spot_counting=spot_counting[rowSums(is.na(spot_counting[,1:7])) == 0, ]
spot_counting=spot_counting[,1:7]
spot_counting$foci=as.integer(spot_counting$`Number of foci`)
spot_counting=spot_counting[rowSums(is.na(spot_counting[,1:8])) == 0, ]
hist(log(spot_counting$foci,base=2))
ggplot(spot_counting,aes(x=foci))+geom_histogram(fill="white",color="black",bins=30,aes(y=..density..))+theme_classic(base_size = 20)+xlab("Number of spots")+geom_vline(xintercept = 45.72973,color="blue",size=1,linetype="dotted")+ geom_vline(xintercept = 178.54167	,color="blue",size=1,linetype="dotted")+geom_vline(xintercept = 371.66667,color="blue",size=1,linetype="dotted")
fitNB(spot_counting$foci)
```

```{r}
x=spot_counting[spot_counting$Annotation=="Target-/Target+",]
ggplot(x, aes(x=Annotation_v2, y=foci)) + geom_boxplot()+xlab("Cell-Type")+ylab("# of spots")+theme_classic(base_size = 20)+ geom_dotplot(binaxis='y', stackdir='center', dotsize=0.75)
wilcox.test(foci ~ Annotation_v2, data = x)
library(dplyr)
x %>% group_by(Annotation_v2) %>% summarise(mean(foci))
```

```{r}
library(mclust)
spot.gmm = Mclust(log(spot_counting$foci,base = 2))
summary(spot.gmm)
spot.gmm$parameters

spot.gmm.1 = Mclust(spot_counting$foci, G=1)
logLik(spot.gmm.1)
# 'log Lik.' -1226.241 (df=2)
logLik(spot.gmm)-logLik(spot.gmm.1)
1-pchisq(14.87156, df=3)
```

```{r}
 library(mixtools)
mixmdl = normalmixEM(spot_counting$foci)
plot(mixmdl,which=2)
```

```{r}
x=as.matrix(spot_counting[,"foci"])
y=dist(x)
z=hclust(y)
plot(z, cex = 0.6, hang = -1)
sub_grp <- cutree(z, k = 4)
table(sub_grp)

spot_counting$grp=sub_grp
spot_counting[spot_counting$grp==4,"grp"]=3 # Merging group 3 and group 4
spot_counting$grp=as.factor(spot_counting$grp)
ggplot(spot_counting,aes(x=grp,y=foci))+geom_boxplot()+xlab("Cluster")+ylab("Number of spots")+theme_classic(base_size = 20)+ geom_dotplot(binaxis='y', stackdir='center', dotsize=0.5)
library(dplyr)
spot_counting %>% group_by(grp) %>% summarise(mean(foci))
table(spot_counting$grp,spot_counting$Annotation_v2)
x=spot_counting[spot_counting$Annotation=="Target-/Target+",]
table(x$grp,x$Annotation_v2)
wilcox.test(x = spot_counting$foci[spot_counting$grp==2],y=spot_counting$foci[spot_counting$grp==3])
```

