---
title: "JSLAB19.Rmd"
author: "MAKINO Manon"
date: '2022-05-31'
output: 
  html_document:
    toc: TRUE
    toc_float: TRUE
---

```{r, include=FALSE}
library(TCC)
library(MBCluster.Seq)
library(openxlsx)
library(ggplot2)
library(ggsci)
library(reshape2)
```

# サンプルのクラスタリング

```{r, fig.width=5, fig.height=4}
data <- read.xlsx("JSLAB18.xlsx", rowNames=TRUE)
out <- clusterSample(data, dist.method="spearman",
            hclust.method="average", unique.pattern=TRUE)
plot(out)
```

# MBCdeg2の実行

## 1. ファイルの読込とsubsetting
```{r}
data_all <- read.xlsx("JSLAB18.xlsx", rowNames=TRUE)
data <- data_all[, 4:9]
```

## 2. 群ラベル情報の作成とTCC正規化
```{r}
data.cl <- c(1,1,1,2,2,2)
tcc <- new('TCC',data,data.cl)
tcc <- calcNormFactors(tcc,norm.method="tmm",test.method="edger", 
                     iteration=3, FDR=0.1,floorPDEG=0.05)
norm.factors <- tcc$norm.factors
ef.libsizes <- colSums(data)*norm.factors
size.factors <- ef.libsizes/mean(ef.libsizes)
```

## 3. MBCluster.Seq実行の共通部分
```{r}
set.seed(2022)
hoge <- RNASeq.Data(Count = data, Treatment = data.cl,
                  Normalizer= log(size.factors))
```

## 4. MBCluster.Seq本番(K=3)
```{r}
K <- 3                    #クラスター数を指定
out_f <- "MBCdeg_K3.xlsx" #出力ファイル名を指定
c0 <- KmeansPlus.RNASeq(data=hoge, nK=K,
           model="nbinom", print.steps=F)
cls <- Cluster.RNASeq(data=hoge, model="nbinom",
           centers=c0$centers, method="EM")
PP <- cls$probability
colnames(PP) <- paste("c", 1:K, sep="")
tmp <- cbind(rownames(data), data, PP, cls$cluster)
write.xlsx(x=tmp, file=out_f)
cls$centers
table(cls$cluster)
cls_k3 <- cls
```

## 5. MBCluster.Seq本番(K=4)
```{r}
K <- 4                    #クラスター数を指定
out_f <- "MBCdeg_K4.xlsx" #出力ファイル名を指定
c0 <- KmeansPlus.RNASeq(data=hoge, nK=K,
           model="nbinom", print.steps=F)
cls <- Cluster.RNASeq(data=hoge, model="nbinom",
           centers=c0$centers, method="EM")
PP <- cls$probability
colnames(PP) <- paste("c", 1:K, sep="")
tmp <- cbind(rownames(data), data, PP, cls$cluster)
write.xlsx(x=tmp, file=out_f)
cls$centers
table(cls$cluster)
cls_k4 <- cls
```

## 6. MBCluster.Seq本番(K=5)
```{r}
K <- 5                    #クラスター数を指定
out_f <- "MBCdeg_K5.xlsx" #出力ファイル名を指定
c0 <- KmeansPlus.RNASeq(data=hoge, nK=K,
           model="nbinom", print.steps=F)
cls <- Cluster.RNASeq(data=hoge, model="nbinom",
           centers=c0$centers, method="EM")
PP <- cls$probability
colnames(PP) <- paste("c", 1:K, sep="")
tmp <- cbind(rownames(data), data, PP, cls$cluster)
write.xlsx(x=tmp, file=out_f)
cls$centers
table(cls$cluster)
cls_k5 <- cls
```

## 7. 作図(元データの作成)
```{r}
dat <- data.frame(
  c1 = c(cls_k3$centers[1,],cls_k4$centers[1,],cls_k5$centers[1,]),
  c2 = c(cls_k3$centers[2,],cls_k4$centers[2,],cls_k5$centers[2,]),
  c3 = c(cls_k3$centers[3,],cls_k4$centers[3,],cls_k5$centers[3,]),
  c4 = c(rep(NA,2),cls_k4$centers[4,],cls_k5$centers[4,]),
  c5 = c(rep(NA,2),rep(NA,2),cls_k5$centers[5,])
)
melt(dat)
matome <- data.frame(melt(dat),k=rep(c(rep('K=3',2),rep('K=4',2),rep('K=5',2)),5),u=rep(1:2,15))
```

## 8. 作図（本番)
```{r, fig.width=6, fig.height=4}
g <- ggplot(data = matome)
g <- g + geom_line(aes(x = u, y = value),size = 1)
g <- g + theme_bw()
g <- g + theme(panel.grid.major =element_line(colour = "gray" ),
               aspect.ratio = 1)
g <- g + labs(x='Group name',y='log(pH7 / pH4.5)')
g <- g + scale_x_continuous(breaks=1:2,labels = c('pH4.5','pH7'))
g <- g + facet_grid(k ~ variable)
plot(g)
```
