# 3d R plot
#saveat <- "/mnt/NAS_coruscant_datashare/haasf/madland_RNA-seq_Schippers"
saveat <- "R:/haasf/madland_RNA-seq_Schippers"

#file.rpkm <- '/mnt/Teamdrive/projects/MAdLand/cooperations/Schippers_Lab_RNA-seq_analyses/Mp.p.sort.wotRNA.rpkm'
file.rpkm <- "S:/projects/MAdLand/cooperations/Schippers_Lab_RNA-seq_analyses/Mp.p.sort.wotRNA.rpkm"
# file.rpkm <- "S:/projects/MAdLand/cooperations/Schippers_Lab_RNA-seq_analyses/31315_all.p.sort.rpkm"

data.rpkm <- read.table(file.rpkm, header=T, sep="\t", row.names=1)


# sort by colnames
data.rpkm <- data.rpkm[,order(colnames(data.rpkm))]


librariesName <- list(
  control_6h = c("control_6h", "red"),
  control_24h = c("control_24h", "blue"),
  LCS_6h = c("LCS_6h", "green"),
  LCS_24h = c("LCS_24h", "yellow")
)

header.ori <- c("MF82", "MF81", "MF80", "MF76", "MF75", "MF74", "MF68", "MF67", "MF66", "MF62", "MF61", "MF60")
header.new <- c("LCS_24h.2", "LCS_24h.1", "LCS_24h", "control_24h.2", "control_24h.1", "control_24h", "control_6h.2", "LCS_6h.1", "LCS_6h", "LCS_6h.2", "control_6h.1", "control_6h")
# header.ori <- c("P97", "P96", "P95", "P91", "P90", "P89", "P111", "P110", "P109", "P105", "P104", "P103")
# header.new <- c("LCS_6h.2", "LCS_6h.1", "LCS_6h", "control_6h.2", "control_6h.1", "control_6h", "LCS_24h.2", "LCS_24h.1", "LCS_24h", "control_24h.2", "control_24h.1", "control_24h")


header.new <- header.new[order(header.ori)]
header.ori <- header.ori[order(header.ori)]

col.header <- header.new

colnames(data.rpkm) <- col.header

library("DESeq2")
library("ggplot2")
library("RColorBrewer")
library("pheatmap")
library("BiocGenerics")
library("rgl")
library("magick")
library("sjmisc")


################### running ######################
### PCA RPKM ###

set.seed(0)

data.inv <- t(data.rpkm)
data.dist <-dist(data.inv, method="euclidean") # "euclidean", "maximum", "manhattan", "canberra", "binary" or "minkowski"
data.dist.hc <- hclust(data.dist,method="ward.D2")
data.dist.pca <-princomp(data.dist,cor=T)

pc1 <- data.dist.pca$scores[,1]
pc2 <- data.dist.pca$scores[,2]
pc3 <- data.dist.pca$scores[,3]

# create data frame for pc1 pc2 pc3
data.dist.pca.frame = data.frame(pc1,pc2,pc3)
rownames(data.dist.pca.frame) <- names(data.dist.pca$scale)
colnames(data.dist.pca.frame) <- c("pc1","pc2","pc3")

condition.values <- c()
condition.values.color <- c()


for(a in colnames(data.rpkm)) {

  v <- substr(a, nchar(a), nchar(a))
  if(str_contains(v, c(1,2,3,4,5,6,7,8,9,0), logic = "or")) {
    if (substr(substr(a, 1, nchar(a)-2), nchar(substr(a, 1, nchar(a)-2)), nchar(substr(a, 1, nchar(a)-2))) == ".") {
      n <- substr(a, 1, nchar(a)-3)
    } else {
      n <- substr(a, 1, nchar(a)-2)
    }

  } else {
    n <- a

  }
  condition.values <- c(condition.values, librariesName[n][[1]][1])
  condition.values.color <- c(condition.values.color, librariesName[n][[1]][2])
}


data.dist.pca.frame["tissue"] <- condition.values
data.dist.pca.frame["color"] <- condition.values.color
data.dist.pca.frame["name"] <- names(data.dist.pca$scale)
attr(data.dist.pca.frame, "percentVar") <- (data.dist.pca$sdev)^2 / sum(data.dist.pca$sdev^2) # cumsum()

# simple plot
png(filename=paste0(saveat, "/HC_RPKM_normalized.png"))
plot(data.dist.hc) # hc plot
dev.off()
png(filename=paste0(saveat, "/PCA_variance_RPKM_normalized.png"))
plot(data.dist.pca) # variances; var(data.dist.pca$sdev[1:9])
dev.off()
# get the parcent variation
percentVar <- round(100 * attr(data.dist.pca.frame, "percentVar"))

# 3d plot
plot3d(pc1, pc2, pc3,
       type = "s", # p, s, l, h, n
       #pch = c(1:3),
       col = condition.values.color,
       size = 1,
       xlab = paste0("PC1: ",percentVar[1],"% variance"),
       ylab = paste0("PC2: ",percentVar[2],"% variance"),
       zlab = paste0("PC3: ",percentVar[3],"% variance"),
       cex = 2,
       main = "", # -> princomp",
)

# shift <- matrix(4, 4, 4, byrow = TRUE)
# text3d(shift,texts=1:3)
grid3d(c("x", "y", "z"))

## add legend
legend3d("right", unique(condition.values), pch = 19, col = unique(condition.values.color))

#### video #####
M <- par3d("userMatrix")
play3d( par3dinterp( userMatrix=list(M,
                                     rotate3d(M, pi/2, 1, 0, 0),
                                     rotate3d(M, pi/2, 0, 1, 0) ) ),
        duration=2 )

movie3d(spin3d(axis = c(1, 2, 1)), duration = 5,
        dir = saveat)

#### video end ####

# pc1, pc2
png(filename=paste0(saveat, "/PCA_RPKM_normalized.png"))
ggplot(
  data.dist.pca.frame,
  aes(
    pc1,
    pc2,
    color=tissue)
) +
  geom_point(size=2.5) +
  xlab(
    paste0("PC1: ",percentVar[1],"% variance")
  ) +
  ylab(
    paste0("PC2: ",percentVar[2],"% variance")
  ) +
  #theme() + #, face="bold"
  scale_colour_manual(
    values= c("blue", "red", "yellow", "green") # dodgerblue3
  ) +
  ggtitle("PCA of all samples (RPKM normalized)") +
  theme(
    plot.title = element_text(lineheight=.8),
    panel.background = element_rect(fill = "gray95")
  )
dev.off()





######## Some experiments ######
# for(i in colnames(data)) {
#     colnames(data)[colnames(data)==grep(i, names(data), value = TRUE)] <- librariesName[[i]]
# }
#
# for(i in colnames(data)) {
#     colnames(data)[colnames(data)==grep(i, names(data), value = TRUE)] <- names(librariesName[librariesName==grep(i, librariesName, value = TRUE)])
# }