pow_baudin.Rmd

---
title: "Title: Barley (Baudin) RNAseq analysis during powdery mildew infection"
author: "Author: Paula Moolhuijzen"
date: "Date: `r format(Sys.time(), '%d %B, %Y')`"
output: html_document
---

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

### DESeq2 analysis of Baudin isogenic lines

```{r, echo=FALSE, message=F, warning=F}
getwd()
setwd(getwd())
print(getwd())

library(DESeq2)
library(ggplot2)
library(genefilter)
library(gplots)
library(affy)
library(pheatmap)
library(dplyr)
library(pca3d)
library(ggfortify)
library("Hmisc")
library(corrplot)
require(scatterplot3d)
library(reshape2)
library(ggpubr)
library(dendextend)
#library(plyr)
```

Raw data set = H201SC19093309
```{r, echo=FALSE, message=F, warning=F}

meta <- read.table("meta_baudin", sep="\t", header = TRUE, row.names = 1)
mat <- read.table("matrix_baudin", sep="\t", header = TRUE, row.names=1)
mat <- mat[rowSums(mat) > 90, ]

data <- DESeqDataSetFromMatrix(countData=mat, colData=meta, design=~line + treatment)

rld <- rlogTransformation(data)
vst <- vst(data, blind = FALSE)
dds <- DESeq(data)
vsd <- vst(dds, blind = FALSE)
```



#### PCA plot log transformed data set

Variance data set based on line and treatment.

```{r , echo=FALSE}
rv <- rowVars(assay(vst))
# select the ntop genes by variance
select <- head(order(rv, decreasing=TRUE), 1000)
plotPCA(vsd, intgroup = "line_treatment")
```

PCA 3-Dimensional plot
```{r, echo=FALSE}

pca <- prcomp(t(assay(vsd)[select,]))



```

PCA 2-Dimensional
```{r, echo=FALSE}

pca2d(pca, group=meta$line_treatment, legend="topright", show.ellipses=TRUE, ellipse.ci=0.75,show.plane=FALSE)

#dev.off
```



```{r, echo=FALSE}
res <- cor(mat[select,])
res <- round(res, 3)

#install.packages("Hmisc")

ass <- assay(vsd)[select,]
res2 <-rcorr(res, type = c("pearson","spearman"))

```


Flattened correlation coefficients and correlation p-values
```{r, echo=FALSE}

# ++++++++++++++++++++++++++++
# flattenCorrMatrix
# ++++++++++++++++++++++++++++
# cormat : matrix of the correlation coefficients
# pmat : matrix of the correlation p-values
flattenCorrMatrix <- function(cormat, pmat) {
  ut <- upper.tri(cormat)
  data.frame(
    row = rownames(cormat)[row(cormat)[ut]],
    column = rownames(cormat)[col(cormat)[ut]],
    cor  =(cormat)[ut],
    p = pmat[ut]
    )
}

```


Correlation plot with significance

```{r, echo=FALSE}

pdf("baudin-corr-plot.pdf")
corrplot(res2$r, type = "upper", order = "hclust", 
         tl.col = "black", sig.level=0.01, tl.srt = 45, insig="blank",  cl.lim=c(min(res2$r), max(res2$r)), is.corr=FALSE)

dev.off()


#?corrplot
```

Plot of the correlation coefficent 

```{r, echo=FALSE}
pdf(file="baudin-heatmap.pdf")

library(RColorBrewer)
heatmap.2(res2$r, scale = "row", key=TRUE, col= colorRampPalette(brewer.pal(8, "Oranges"))(10))
dev.off()


```



Bar plot of Principal components

```{r, echo=FALSE}
#Determine the proportion of variance of each component
#Proportion of variance equals (PC stdev^2) / (sum all PCs stdev^2)
pca.proportionvariances <- ((pca$sdev^2) / (sum(pca$sdev^2)))*100
barplot(pca.proportionvariances, cex.names=1, xlab=paste("Principal component (PC), 1-", length(pca$sdev)), ylab="Proportion of variation (%)", main="Screen plot", ylim=c(0,100))

```


Scatterplot of PCA samples

```{r, echo=FALSE}

#install.packages("scatterplot3d")

par(mar=c(4,4,4,4), cex=1.0, cex.main=0.8, cex.axis=0.8)

colors <- c("blue", "cornflowerblue", "aquamarine4", "aquamarine1", "darkorchid4", "darkorchid3")
colors <- colors[as.numeric(meta$line_treatment)]
s3d<-scatterplot3d(pca$x[,1:3],  main="", color=colors, pch=17, xlab=paste("PC1, ", round(pca.proportionvariances[1], 2), "%"), ylab=paste("PC2, ", round(pca.proportionvariances[2], 2), "%"), zlab=paste("PC3, ", round(pca.proportionvariances[3], 2), "%"), grid=TRUE, box=FALSE)

```


#### Plots for raw read count per sample
```{r, echo=FALSE}
bpt <- barplot(colSums(mat)/1000000, 
               main="Total no. reads per sample (million)",
               #col=phenoTable$color, 
               las=1,  horiz=TRUE,
               ylab="Samples", cex.names=0.5,
               xlab="Million counts")

```

#### Log2 gene read count
```{r, echo=FALSE}
epsilon <- 1 # pseudo-count to avoid problems with log(0)
hg <- hist(as.matrix(log2(mat + epsilon)), breaks=100, col="blue", border="white",
           main="Log2 counts per gene", xlab="log2(counts+1)", ylab="Number of genes", 
           las=1, cex.axis=0.7)
```

#### Log2 read count per sample
```{r mat, echo=FALSE}
bx <- boxplot(log2(mat + epsilon), pch=".", 
              horizontal=TRUE, cex.axis=0.5, main="Log2 counts per sample",
              las=1, ylab="Samples", xlab="log2(counts +1)")
```

#### Normalised read counts samples

```{r ,echo=FALSE }
# plot 6
epsilon <- 1 
bx <- boxplot(log2(counts(dds, normalized=TRUE)+epsilon),  col=meta$Group, cex.axis=0.7, 
              las=1, xlab="log2(normalized counts)", horizontal=TRUE, main="Normalized counts")
```

#### Density plot 
```{r,echo=FALSE }
ds <- plotDensity(log2(counts(dds, normalized=TRUE)+epsilon),  
                  xlab="log2(normalized counts)", cex.lab=0.7, panel.first=grid()) 
```

#### Mean and variance relationship - Computing mean and variance
```{r, echo=FALSE }
norm.counts <- counts(dds, normalized=TRUE)
mean.counts <- rowMeans(norm.counts)
variance.counts <- apply(norm.counts, 1, var)

mean.var.col <- densCols(x=log2(mean.counts), y=log2(variance.counts))
plot(x=log2(mean.counts), y=log2(variance.counts), pch=16, cex=0.5, 
     col=mean.var.col, main="Mean-variance relationship",
     xlab="Mean log2(normalized counts) per gene",
     ylab="Variance of log2(normalized counts)",
     panel.first = grid())
abline(a=0, b=1, col="brown")

```

#### MA-plot

```{r , echo=FALSE}

de <- results(dds)
mostsig <- subset(de, (padj <= 0.05 & !is.na(padj)) & abs(log2FoldChange) >= 1 )
ggmaplot(de, fdr = 0.05, fc = 2, size = 0.4, palette = c("#B31B21", "#1465AC", "darkgray"),
         genenames = as.vector(mostsig$name), legend = "top", top = 5, font.label = c("bold", 11),
         font.legend = "bold", font.main = "bold",
         ggtheme = ggplot2::theme_minimal())
```


#### Plot dispersion estimates

```{r, echo=FALSE}
eds <- plotDispEsts(dds)
```

#### Plot expression of top 20 genes with largest variance (vst)
```{r, echo=FALSE }
topVarGenes <- head(order(rowVars(assay(vsd)), decreasing = TRUE), 100)
top  <- assay(vst)[ topVarGenes, ]
#top  <- top - rowMeans(top)

anno <- as.data.frame(colData(dds)[, c("line","treatment")])

write.table(top, file="Baudin-topVarGenes.txt")
pdf("Baudin-top100-expression.pdf")

pheatmap(top, annotation_col = anno, show_rownames = F)

dev.off()
```

#### Plot expression of top 20 genes with largest variance (Var-MeanVar)

```{r,echo=FALSE }
topVarGenes <- head(order(rowVars(assay(vst)), decreasing = TRUE),100)
top  <- assay(vst)[ topVarGenes, ]
top  <- top - rowMeans(top)

anno <- as.data.frame(colData(vst)[, c("line","treatment")])

pdf("Baudin-top100-express-mean.pdf")

pheatmap(top, annotation_col = anno, show_rownames = F)

dev.off()
```

```{r, echo=FALSE}
resultsNames(dds)
```

### Identify significant differentially expressed (DE) genes

```{r , echo=FALSE, message=F, warning=F }
meta <- read.table("meta_baudin", sep="\t", header = TRUE, row.names = 1)
#str(meta)
mat <- read.table("matrix_baudin", sep="\t", header = TRUE, row.names=1)
#mat <- mat[rowSums(mat) > 80, ]

data <- DESeqDataSetFromMatrix(countData=mat, colData=meta, design=~line_treatment)
dds <- DESeq(data)
resultsNames(dds)
```

Printing most significant at (padj <= 0.05 & !is.na(padj)) & abs(log2FoldChange) >= 1 )
```{r , echo=FALSE, message=F, warning=F }
res <- results(dds)
mostsig <- subset(res, (padj <= 0.05 & !is.na(padj)) & abs(log2FoldChange) >= 1 )
dim(mostsig)
write.table(mostsig, file="Baudin-most-significant.txt", sep="\t", col.names=TRUE, row.names=TRUE)
```

```{r, echo=FALSE}

ggmaplot(res, fdr = 0.05, fc = 2, size = 0.4, palette = c("#B31B21", "#1465AC", "darkgray"),
         genenames = as.vector(mostsig$name), legend = "top", top = 5, font.label = c("bold", 11),
         font.legend = "bold", font.main = "bold",
         ggtheme = ggplot2::theme_minimal())
```


#### Most signifcant gene
```{r}
topGene <- rownames(res)[which.min(res$padj)]
topGene
```

```{r}
library("ggbeeswarm")
geneCounts <- plotCounts(dds, gene = topGene, intgroup = c("treatment","line"),
                         returnData = TRUE)
ggplot(geneCounts, aes(x = treatment, y = count, color = line)) +
  scale_y_log10() +  geom_beeswarm(cex = 3)

```

#### Z83834.1 H.vulgare mRNA for Mlo protein HORVU4Hr1G082710 counts
```{r, echo=FALSE}
#topGene <- rownames(res)[which.max(res$log2FoldChange)]

Mlo <- "HORVU4Hr1G082710" # Z83834.1 H.vulgare mRNA for Mlo protein
```



```{r}
library("ggbeeswarm")
geneCounts <- plotCounts(dds, gene = Mlo, intgroup = c("treatment","line"),
                         returnData = TRUE)
ggplot(geneCounts, aes(x = treatment, y = count, title = Mlo, color = line)) +
  scale_y_log10() +  geom_beeswarm(cex = 3)

```


#### Gene with largest fold change

```{r, echo=FALSE}
LargestFC <- rownames(res)[which.max(res$log2FoldChange)]
LargestFC
```



```{r}
library("ggbeeswarm")
geneCounts <- plotCounts(dds, gene = LargestFC, intgroup = c("treatment","line"),
                         returnData = TRUE)
ggplot(geneCounts, aes(x = treatment, y = count, title = LargestFC , color = line)) +
  scale_y_log10() +  geom_beeswarm(cex = 3)

```


```{r, echo=FALSE}
vst <- vst(dds, blind = FALSE)
#head(assay(vst), 3)
```

### Simple contrast analysis of uninfected (control) Baudin isogenic lines 

##### 1. Contrast isogenic lines S12 (12 copy point mutation mlo 0%) versus S0 (Wt mlo 100%) controls
```{r, echo=FALSE }

de_S12_S0_c <- results(dds, contrast=c("line_treatment","S12_Control","S0_Control"), lfcThreshold=1, altHypothesis="greaterAbs", pAdjustMethod = "BH")

mostsig <- subset(de_S12_S0_c, (padj <= 0.05 & !is.na(padj)) & abs(log2FoldChange) >= 1 )
ggmaplot(de_S12_S0_c, fdr = 0.05, fc = 2, size = 0.4, palette = c("#B31B21", "#1465AC", "darkgray"),
         genenames = as.vector(mostsig$name), legend = "top", top = 5, font.label = c("bold", 11),
         font.legend = "bold", font.main = "bold",
         ggtheme = ggplot2::theme_minimal())
```


Number of DE genes at log2FoldChange >= 1 & padj <= 0.05 for S12C_control versus S0C_control
```{r , echo=FALSE}

de_S12_S0_c <- subset(de_S12_S0_c, abs(de_S12_S0_c$log2FoldChange) >= 1 & de_S12_S0_c$padj <= 0.05)
write.table(de_S12_S0_c, file="most-significant-S12C_control-S0C_control.txt", sep="\t", col.names=TRUE, row.names=TRUE)
dim(de_S12_S0_c)
```

```{r , echo=FALSE}
#head(de_S12_S0_c)
de_S12_S0_cl <- rownames(de_S12_S0_c)
deS12S0cl <- log2(assay(dds)[ de_S12_S0_cl, ] +1)

anno <- as.data.frame(colData(dds)[, c("line","treatment")])

pheatmap(deS12S0cl, annotation_col = anno, show_rownames = F, main = "S12S / S0S", cutree_rows = 4)

```

```{r}

cal_z_score <- function(deS12S0cl){
  (deS12S0cl - mean(deS12S0cl)) / sd(deS12S0cl)
}
 deS12S0cln <- t(apply(deS12S0cl, 1, cal_z_score))
anno <- as.data.frame(colData(dds)[, c("line","treatment")])

pheatmap(deS12S0cln, annotation_col = anno, show_rownames = F, main = "S12S / S0S", cutree_rows = 4)
?pheatmap
```

##### 2. Simple contrast isogenic lines S2 (two copy point mutation mlo 20%) and S12 (12 copy mutation mlo 0%) controls

Number of DE genes at abs(log2FoldChange) >= 1 & padj <= 0.05 for S2_control versus S12_control
```{r}
de_S2S12_c <-results(dds, contrast=c("line_treatment","S2_Control","S12_Control"), lfcThreshold=1, altHypothesis="greaterAbs", pAdjustMethod = "BH")

mostsig <- subset(de_S2S12_c, (padj <= 0.05 & !is.na(padj)) & abs(log2FoldChange) >= 1 )
ggmaplot(de_S2S12_c, fdr = 0.05, fc = 2, size = 0.4, palette = c("#B31B21", "#1465AC", "darkgray"),
         genenames = as.vector(mostsig$name), legend = "top", top = 5, font.label = c("bold", 11),
         font.legend = "bold", font.main = "bold",
         ggtheme = ggplot2::theme_minimal())
```

```{r}
de_S2S12_c <- subset(de_S2S12_c, !is.na(padj) & abs(log2FoldChange) >= 1 & padj <= 0.05)
write.table(de_S2S12_c, file="most-significant-S2_control-S12_control.txt", sep="\t", col.names=TRUE, row.names=TRUE)
dim(de_S2S12_c)
```


```{r}
de_S2_S12_cl <- rownames(de_S2S12_c)
deS2S12cl <- log2(assay(dds)[de_S2_S12_cl, ] +1)
#deS12S0cl  <- deS12S0cl - rowMeans(deS12S0cl)

anno <- as.data.frame(colData(dds)[, c("line","treatment")])

pheatmap(deS2S12cl, annotation_col = anno, show_rownames = F, main = "S2S / S12S", cutree_rows = 5)

```

```{r}

cal_z_score <- function(deS2S12cl){
  (deS2S12cl - mean(deS2S12cl)) / sd(deS2S12cl)
}
 deS2S12cln <- t(apply(deS2S12cl, 1, cal_z_score))
anno <- as.data.frame(colData(dds)[, c("line","treatment")])

pheatmap(deS2S12cln, annotation_col = anno, show_rownames = F, main = "S2S / S12S", cutree_rows = 5)
?pheatmap
```

##### 3. Simple contrast isogenic lines S2 (two copy point mutation mlo 20%) and S0 (Wt mlo 100%) controls

```{r, echo=FALSE }
# e.g. for P18 infection vs control
de_S2_S0_c <- results(dds, contrast=c("line_treatment","S2_Control","S0_Control"), lfcThreshold=1, altHypothesis="greaterAbs", pAdjustMethod = "BH")


mostsig <- subset(de_S2_S0_c, (padj <= 0.05 & !is.na(padj)) & abs(log2FoldChange) >= 1 )
ggmaplot(de_S2_S0_c, fdr = 0.05, fc = 2, size = 0.4, palette = c("#B31B21", "#1465AC", "darkgray"),
         genenames = as.vector(mostsig$name), legend = "top", top = 5, font.label = c("bold", 11),
         font.legend = "bold", font.main = "bold",
         ggtheme = ggplot2::theme_minimal())
```

```{r, echo=FALSE}
#Number of DE genes at log2FoldChange >= 1 & padj <= 0.05 for S2_control versus S0_control

de_S2_S0_c <- subset(de_S2_S0_c, abs(de_S2_S0_c$log2FoldChange) >= 1 & de_S2_S0_c$padj <= 0.05)
write.table(de_S2_S0_c, file="most-significant-S2_control-S0_control.txt", sep="\t", col.names=TRUE, row.names=TRUE)
dim(de_S2_S0_c)
```

```{r}
de_S2_S0_cl <- rownames(de_S2_S0_c)
deS2S0cl <- log2(assay(dds)[de_S2_S0_cl, ] +1)
#deS12S0cl  <- deS12S0cl - rowMeans(deS12S0cl)

anno <- as.data.frame(colData(dds)[, c("line","treatment")])

pheatmap(deS2S0cl, annotation_col = anno, show_rownames = F, main = "S2S / S0S", cutree_rows = 6)

```

```{r}

cal_z_score <- function(deS2S0cl){
  (deS2S0cl - mean(deS2S0cl)) / sd(deS2S0cl)
}
 deS2S0cln <- t(apply(deS2S0cl, 1, cal_z_score))
anno <- as.data.frame(colData(dds)[, c("line","treatment")])

pheatmap(deS2S0cln, annotation_col = anno, show_rownames = F, main = "S2S / S0S", cutree_rows = 6)

```

##### 4. Simple contrast isogenic lines S2 (two copy point mutation mlo 20%) and S0 (Wt mlo 100%) Inoculation

```{r, echo=FALSE }
# e.g. for P18 infection vs control
de_S2_S0_p <- results(dds, contrast=c("line_treatment","S2_Inoculation","S0_Inoculation"), lfcThreshold=1, altHypothesis="greaterAbs", pAdjustMethod = "BH")

mostsig <- subset(de_S2_S0_p, (padj <= 0.05 & !is.na(padj)) & abs(log2FoldChange) >= 1 )
ggmaplot(de_S2_S0_p, fdr = 0.05, fc = 2, size = 0.4, palette = c("#B31B21", "#1465AC", "darkgray"),
         genenames = as.vector(mostsig$name), legend = "top", top = 5, font.label = c("bold", 11),
         font.legend = "bold", font.main = "bold",
         ggtheme = ggplot2::theme_minimal())
```


```{r, echo=FALSE }
#Number of DE genes at log2FoldChange >= 1 & padj <= 0.05 for S2_control versus S0_control

de_S2_S0_p <- subset(de_S2_S0_p, abs(de_S2_S0_p$log2FoldChange) >= 1 & de_S2_S0_p$padj <= 0.05)
write.table(de_S2_S0_p, file="most-significant-S2_Inoculation-S0_Inoculation.txt", sep="\t", col.names=TRUE, row.names=TRUE)
dim(de_S2_S0_p)
```

```{r}
de_S2_S0_pl <- rownames(de_S2_S0_p)
deS2S0pl <- log2(assay(dds)[de_S2_S0_pl, ] +1)
#deS12S0cl  <- deS12S0cl - rowMeans(deS12S0cl)

anno <- as.data.frame(colData(dds)[, c("line","treatment")])

pheatmap(deS2S0pl, annotation_col = anno, show_rownames = F, main = "S2P / S0P", cutree_rows = 4)

```

```{r}

cal_z_score <- function(deS2S0pl){
  (deS2S0pl - mean(deS2S0pl)) / sd(deS2S0pl)
}
 deS2S0pln <- t(apply(deS2S0pl, 1, cal_z_score))
anno <- as.data.frame(colData(dds)[, c("line","treatment")])

pheatmap(deS2S0pln, annotation_col = anno, show_rownames = F, main = "S2P / S0P", cutree_rows = 5)
```

##### 5. Simple contrast isogenic lines S2 (two copy point mutation mlo 20%) and S12 (11 copy) infection

```{r, echo=FALSE }
# e.g. for P18 infection vs control
de_S2_S12_p <- results(dds, contrast=c("line_treatment","S2_Inoculation","S12_Inoculation"), lfcThreshold=1, altHypothesis="greaterAbs", pAdjustMethod = "BH")

mostsig <- subset(de_S2_S12_p, (padj <= 0.05 & !is.na(padj)) & abs(log2FoldChange) >= 1 )
ggmaplot(de_S2_S12_p, fdr = 0.05, fc = 2, size = 0.4, palette = c("#B31B21", "#1465AC", "darkgray"),
         genenames = as.vector(mostsig$name), legend = "top", top = 5, font.label = c("bold", 11),
         font.legend = "bold", font.main = "bold",
         ggtheme = ggplot2::theme_minimal())

```

```{r, echo=FALSE }
#Number of DE genes at log2FoldChange >= 1 & padj <= 0.05 for S2_control versus S0_control

de_S2_S12_p <- subset(de_S2_S12_p, abs(de_S2_S12_p$log2FoldChange) >= 1 & de_S2_S12_p$padj <= 0.05)
write.table(de_S2_S12_p, file="most-significant-S2_Inoculation-S12_Inoculation.txt", sep="\t", col.names=TRUE, row.names=TRUE)
dim(de_S2_S12_p)
```

```{r}
de_S2_S12_pl <- rownames(de_S2_S12_p)
deS2S12pl <- log2(assay(dds)[de_S2_S12_pl, ] +1)
#deS12S0cl  <- deS12S0cl - rowMeans(deS12S0cl)

anno <- as.data.frame(colData(dds)[, c("line","treatment")])

pheatmap(deS2S12pl, annotation_col = anno, show_rownames = F, main = "S2P / S12P", cutree_rows = 7)

```

```{r}

cal_z_score <- function(deS2S12pl){
  (deS2S12pl - mean(deS2S12pl)) / sd(deS2S12pl)
}
 deS2S12pln <- t(apply(deS2S12pl, 1, cal_z_score))
anno <- as.data.frame(colData(dds)[, c("line","treatment")])

pheatmap(deS2S12pln, annotation_col = anno, show_rownames = F, main = "S2P / S12P", cutree_rows = 4)
```


##### 6. Simple contrast isogenic lines S12 (11  mutation mlo ) and S0 (Wild type) infection

```{r, echo=FALSE }
# e.g. for P18 infection vs control
de_S12_S0_p <- results(dds, contrast=c("line_treatment","S12_Inoculation","S0_Inoculation"), lfcThreshold=1, altHypothesis="greaterAbs", pAdjustMethod = "BH")

mostsig <- subset(de_S12_S0_p, (padj <= 0.05 & !is.na(padj)) & abs(log2FoldChange) >= 1 )
ggmaplot(de_S12_S0_p, fdr = 0.05, fc = 2, size = 0.4, palette = c("#B31B21", "#1465AC", "darkgray"),
         genenames = as.vector(mostsig$name), legend = "top", top = 5, font.label = c("bold", 11),
         font.legend = "bold", font.main = "bold",
         ggtheme = ggplot2::theme_minimal())

```

```{r, echo=FALSE }
#Number of DE genes at log2FoldChange >= 1 & padj <= 0.05 for S2_control versus S0_control

de_S12_S0_p <- subset(de_S12_S0_p, abs(de_S12_S0_p$log2FoldChange) >= 1 & de_S12_S0_p$padj <= 0.05)
write.table(de_S12_S0_p, file="most-significant-S12_Inoculation-S0_Inoculation.txt", sep="\t", col.names=TRUE, row.names=TRUE)
dim(de_S12_S0_p)
```

```{r}
de_S12_S0_pl <- rownames(de_S12_S0_p)
deS12S0pl <- log2(assay(dds)[de_S12_S0_pl, ] +1)
#deS12S0cl  <- deS12S0cl - rowMeans(deS12S0cl)

anno <- as.data.frame(colData(dds)[, c("line","treatment")])

pheatmap(deS12S0pl, annotation_col = anno, show_rownames = F, main = "S12P / S0P", cutree_rows = 5)

```

```{r}
de_S12_S0_pl <- rownames(de_S12_S0_p)
deS12S0pl <- log2(assay(dds)[de_S12_S0_pl, ] +1)
#deS12S0cl  <- deS12S0cl - rowMeans(deS12S0cl)
cal_z_score <- function(deS12S0pl){
  (deS12S0pl - mean(deS12S0pl)) / sd(deS12S0pl)
}
 deS12S0pln <- t(apply(deS12S0pl, 1, cal_z_score))
anno <- as.data.frame(colData(dds)[, c("line","treatment")])

pheatmap(deS12S0pln, annotation_col = anno, show_rownames = F, main = "S12P / S0P", cutree_rows = 2)
```