M0 <- read.csv("../ETEToolkit/ETE3_results_genetree_M0.csv", stringsAsFactors = F, header = F)
colnames(M0) <- c("gene","dNdS")
M12 <- read.csv("../ETEToolkit/ETE3_genetree_12_results.csv", stringsAsFactors = F, header = F)
colnames(M12) <- c("gene","lnL_M1","lnL_M2","p_value_12","w_M1","w_M2")
M12[4] <- lapply(M12[4], gsub, pattern="\\*", replacement="")
M12$p_value_12 <- as.numeric(M12$p_value_12)
M78 <- read.csv("../ETEToolkit/ETE3_results_genetree_87.csv", stringsAsFactors = F, header = F)
colnames(M78) <- c("gene","lnL_M7","lnL_M8","p_value_78","w_M7","w_M8")
M78[4] <- lapply(M78[4], gsub, pattern="\\*", replacement="")
M78$p_value_78 <- as.numeric(M78$p_value_78)

M12$'M1-vs-M2' <- M12$lnL_M2 - M12$lnL_M1
M12.s <- M12[c(1,7,4)]
M12.s$fdr_12 <- p.adjust(M12$p_value_12, method = "fdr")
M78$'M7-vs-M8' <- M78$lnL_M8 - M78$lnL_M7
M78.s <- M78[c(1,7,4)]
M78.s$fdr_78 <- p.adjust(M78$p_value_78, method = "fdr")

library(dplyr)
gene.summary.complete <- dplyr::left_join(gene.d.summary, M12.s, by = "gene")
gene.summary.complete <- dplyr::left_join(gene.summary.complete, M78.s, by = "gene")
gene.summary.complete <- dplyr::left_join(gene.summary.complete, M0, by = "gene")
gene.summary.complete <- gene.summary.complete[c(1,3,4,9,10,11,12,13,14,15,5,6,7,8)]

library(knitr)
low.p <- gene.summary.complete[order(gene.summary.complete$fdr_12),]
write.csv(low.p,"full.gene.diversity.csv")
signif.78 <- gene.summary.complete[complete.cases(gene.summary.complete$`M7-vs-M8`),]
signif.78 <- signif.78[signif.78$fdr_78 <= 0.05,]
signif.78 <- signif.78[order(signif.78$fdr_78),]
write.csv(signif.78, "positively.selected.genes.csv")

writeFasta<-function(data, filename){
  fastaLines = c()
  for (rowNum in 1:nrow(data)){
    fastaLines = c(fastaLines, as.character(paste(">", data[rowNum,"name"], sep = "")))
    fastaLines = c(fastaLines,as.character(data[rowNum,"seq"]))
  }
  fileConn<-file(filename)
  writeLines(fastaLines, fileConn)
  close(fileConn)
}

signif.78.seq <- cbind(signif.78$gene,signif.78$seq)
signif.78.seq <- unique(signif.78.seq)
colnames(signif.78.seq) <- c("name","seq")
writeFasta(signif.78.seq,"positively.selected.genes.faa")
kable(low.p[c(1:20),c(1:10)])

library(ggplot2)
ggplot(gene.summary.complete, aes(x = fdr_12)) + geom_histogram(bins = 100) + labs(title="Distribution of M1-vs-M2 fdr-corrected pvalues")+ theme_classic()

ggplot(gene.summary.complete, aes(x = fdr_78)) + geom_histogram(bins = 100) + labs(title="Distribution of M7-vs-M8 fdr-corrected pvalues")+ theme_classic()

gene.summary.M78.logical <- gene.summary.complete[complete.cases(gene.summary.complete$`M7-vs-M8`),]
gene.summary.M78.logical$M78.logical <- gene.summary.M78.logical$fdr_78 <= 0.05
dn.ds.box <- ggplot(gene.summary.M78.logical, aes(x=M78.logical,y=dNdS)) + theme_classic() +
  geom_jitter(size = 01,shape=21,aes(x=M78.logical,y=dNdS, color=M78.logical)) +
  geom_boxplot(outlier.alpha = 0, alpha = 0) + ylim(0,2) +
  labs(y="dN/dS", 
       x="M7-vs-M8 test fdr-corrected p-value <0.05", title = "Boxplots of dN/dS values") + guides(color = guide_legend(override.aes = list(size=5)))
plot(dn.ds.box)

M1278.scatter <- ggplot(gene.summary.M78.logical, aes(x=fdr_12,y=fdr_78)) +
  geom_point(size=1, shape = 21, alpha = 0.4) +
  geom_hline(yintercept= 0.05, linetype = "dashed", color="red") +
  geom_vline(xintercept= 0.05, linetype = "dashed", color="red") +
  theme_bw()
plot(M1278.scatter)

library(dplyr)
geneclusters <- read.csv("../geosmithia_annotations/antismash/data/geneclusters.txt", sep="\t",header = F,stringsAsFactors = F) #hand-modify to extract scaffold name
colnames(geneclusters) <- c("AS_scaffold","scaffold","type","AS_genes")
geneclusters <- geneclusters[c(1,2)]
geneclusters <- unique(geneclusters)
AS_to_cluster <- read.csv("../geosmithia_annotations/antismash/AS_to_cluster.csv", header = T,stringsAsFactors = F) #custom build from index.html file
cluster.coords <- dplyr::left_join(AS_to_cluster,geneclusters,by="AS_scaffold")
cluster.coords <- cluster.coords[c(1,8,4,5,2,3,6,7)]

kable(cluster.coords) %>%
  kable_styling() %>%
  scroll_box(width = "100%", height = "400px")

library(stringr)
gff <- read.csv("../g.morbida_CDS.gff3", sep = "\t", row.names = NULL, skip = 2, header = FALSE)
gff <- gff[c(1,4,5,9)]
colnames(gff) <- c("scaffold","start","stop","gene")
gene.split <- str_split_fixed(gff$gene, ";",2)
gene.split <- str_split_fixed(gene.split[,1], "=",2)
gene.split <- str_split_fixed(gene.split[,2], ":",2)
gff$gene <- as.character(gene.split[,1])
gff$scaffold <- as.character(gff$scaffold)
kable(head(gff,n=50)) %>%
  kable_styling() %>%
  scroll_box(width = "100%", height = "400px")

cluster.coords$gene <- "NA"

gene.cluster.map <- subset(gff, scaffold == cluster.coords[2,"scaffold"] & start >= cluster.coords[2,"From"] & stop <=cluster.coords[2,"To"])
gene.cluster.map$Cluster <- cluster.coords[2,"Cluster"]

for (x in c(3:nrow(cluster.coords))) {
  x.sub1 <- subset(gff, scaffold == cluster.coords[x,"scaffold"] & start >= cluster.coords[x,"From"] & stop <=cluster.coords[x,"To"])
  x.sub1$Cluster <- cluster.coords[x,"Cluster"]
  gene.cluster.map <- rbind(gene.cluster.map,x.sub1)
}

kable(gene.cluster.map) %>%
  kable_styling() %>%
  scroll_box(width = "100%", height = "400px")

gene.cluster.map.concise <- unique(gene.cluster.map[c(4,5)])
cluster.type <- cluster.coords[c(5,6,7)]
colnames(cluster.type)[2] <- "Cluster_type"
gene.cluster.map.concise <- dplyr::left_join(gene.cluster.map.concise, cluster.type, by="Cluster")

kable(gene.cluster.map.concise) %>%
  kable_styling() %>%
  scroll_box(width = "100%", height = "400px")

low.p <- dplyr::left_join(low.p, gene.cluster.map.concise, by="gene")
low.p <- low.p[c(1:13,15,16,17,14)]
colnames(low.p)[c(14,15,16)] <- c("Cluster","Cluster_type","Most.similar.known.cluster")
write.csv(low.p,"full.gene.diversity.csv")
kable(head(low.p)) %>%
  kable_styling() %>%
  scroll_box(width = "100%", height = "400px")

signif.78 <- low.p[complete.cases(low.p$`M7-vs-M8`),]
signif.78 <- signif.78[signif.78$fdr_78 <= 0.05,]
signif.78 <- signif.78[order(signif.78$fdr_78),]
write.csv(signif.78, "positively.selected.genes.csv")

kable(signif.78) %>%
  kable_styling() %>%
  scroll_box(width = "100%", height = "400px")

SM.PS <- signif.78[complete.cases(signif.78$Cluster),]
SM.PS <- SM.PS[c(1,2,4),]
kable(SM.PS) %>%
  kable_styling() %>%
  scroll_box(width = "100%", height = "400px")

write.csv(SM.PS, "positively.selected.SM.cluster.csv")

signif.78 <- dplyr::left_join(signif.78, cazy.results, by="gene")

kable(signif.78) %>%
  kable_styling() %>%
  scroll_box(width = "100%", height = "400px")