library(ggplot2)
library(stringr)
#####
#This part of scripts produces the plots of alignment data in a directory "../report/images_assembly_zvl/" with format ".png". Only genome-aligned reverse reads were taken here. 
#For every sample was created a set of reads of equal capacity (min between all sample reads, estimated manually).
#####

#Uploading data in R
table_set <- data.frame(matrix(ncol = 24, nrow = 0))
colnames(table_set) <- unlist(strsplit("ID;Read;Count_SNPs;Cycles;GC;Mean_length_del;Mean_length_insert;Read_length;Align_length;Dels;Inserts;A>C;A>G;A>T;C>A;C>G;C>T;G>A;G>C;G>T;T>A;T>C;T>G",";"))

temp_name<-read.table("names_of_custom_chromosome.txt")
samples <- c()
for (i in temp_name[,c("V1")]) {
  samples<-c(samples, gsub("_genome_mapq20_table.txt", "", gsub("custom_", "", i)))
}

for (i in samples) {
  filename<-paste("custom_", i, "_genome_mapq20_table.txt", sep='')
  temp<-read.csv(filename, sep=";", header=TRUE)
  temp<-subset(temp, Read==2)
  temp<-temp[sample(nrow(temp), 200000),]
  temp<-cbind(temp, "Sample"=i)
  table_set<-rbind(table_set, temp)
}
rm(temp)
table_set$Sample<-as.factor(table_set$Sample)
dir.create("../report/images_assembly_zvl")

#GC-content

ggplot(data = table_set, aes(x = GC)) +
  geom_density(stat = "density", fill = "pink", alpha = 0.5) +
  facet_wrap(~Sample, nrow = 3) +
  geom_vline(xintercept = 50, linetype="dotted") + 
  theme(plot.title = element_text(hjust = 0.5))
ggsave("../report/images_assembly_zvl/gc.png", width = 2250, height = 1500, units = "px")

#Mismatches per cycle of sequence

df2<-data.frame("Cycles"=0, "Sample"=1)
for (i in samples) {
  cyc <-subset(table_set, Sample == i)[,4]
  cyc <- as.numeric(unlist(strsplit(cyc,",")))
  df2<-rbind(df2, cbind("Cycles" = cyc, "Sample" =i))
}
df2<-df2[-1,]
df2$Sample<-as.factor(df2$Sample)
df2$Cycles<-as.numeric(df2$Cycles)
rm (cyc, i)

ggplot(df2, aes(x = Cycles)) +
  geom_density(stat = "count", fill = "blue", alpha = 0.5) +
  #geom_line(stat = "count", fill = "blue", alpha = 0.5) +
  facet_wrap(~Sample, nrow = 3) +
  ggtitle("Mismatches per cycle") +
  theme(plot.title = element_text(hjust = 0.5))
ggsave("../report/images_assembly_zvl/mpc.png", width = 2250, height = 1500, units = "px")
ggplot(df2, aes(x = Cycles)) +
  geom_density(stat = "density", fill = "blue", alpha = 0.5) +
  facet_wrap(~Sample, nrow = 3) +
  ggtitle("Mismatches per cycle") +
  theme(plot.title = element_text(hjust = 0.5))
ggsave("../report/images_assembly_zvl/mpc_density.png", width = 2250, height = 1500, units = "px")
rm (df2)

#Mismatches depending on GC-content

mgc2<-cbind(table_set[, c("Count_SNPs", "Sample")], "GC" = round(table_set$GC))

ggplot(mgc2, aes(x=GC, y=Count_SNPs)) +
  stat_summary(fun = sum, geom="line", linewidth = 0.5) +
  facet_wrap(~Sample, nrow = 3) +
  ggtitle("Mismatches/GC") +
  geom_vline(xintercept = 50, linetype="dotted") + 
  theme_bw()+
  theme(plot.title = element_text(hjust = 0.5))
ggsave("../report/images_assembly_zvl/mgc.png", width = 2250, height = 1500, units = "px")
ggplot(mgc2, aes(x=GC, y=Count_SNPs)) +
  stat_summary(fun = mean, geom="line", size = 0.5) +
  facet_wrap(~Sample, nrow = 3) +
  ggtitle("Mismatches/GC") +
  geom_vline(xintercept = 50, linetype="dotted") + 
  theme_bw()+
  theme(plot.title = element_text(hjust = 0.5))
ggsave("../report/images_assembly_zvl/mgc_mean.png", width = 2250, height = 1500, units = "px")

rm (mgc2)

#Frequency of different types of mismatches

mtype2 <- data.frame(matrix(ncol = 3, nrow = 0))
colnames(mtype2)<- c("Sum", "type", "Sample")

for (i in samples) {
  temp <- apply(table_set[, c(12:23)][table_set$Sample==i,], 2, sum)
  temp <- data.frame("Sum" = temp)
  temp$type<-as.factor(colnames(table_set[, c(12:23)]))
  temp<-cbind(temp, "Sample" = i)
  mtype2<-rbind(mtype2, temp)
}
rm(temp)

ggplot(mtype2) +
  geom_line(aes(x = type, y = Sum, group = Sample)) +
  geom_point(aes(x = type, y = Sum)) +
  #facet_wrap(~Sample, nrow = 3, scales = "free_y") +
  facet_wrap(~Sample, nrow = 3, ) +
  ggtitle("Mismatch type") +
  theme(plot.title = element_text(hjust = 0.5),axis.text.x = element_text(angle = 90, hjust = 1))
ggsave("../report/images_assembly_zvl/mt.png", width = 2250, height = 1500, units = "px")

mtype2_2 <- data.frame(matrix(ncol = 3, nrow = 0))
colnames(mtype2_2)<- c("mean", "type", "Sample")
for (i in samples) {
  temp <- apply(table_set[table_set$Sample==i, c(12:23)], 2, mean)
  temp <- data.frame("mean" = temp)
  temp$type<-as.factor(colnames(table_set[, c(12:23)]))
  temp<-cbind(temp, "Sample" = i)
  mtype2_2<-rbind(mtype2_2, temp)
}
mtype2_2$Sample<-as.factor(mtype2_2$Sample)

ggplot(mtype2_2) +
  geom_line(aes(x = type, y = mean, group = Sample)) +
  geom_point(aes(x = type, y = mean)) +
  #facet_wrap(~Sample, nrow = 3, scales = "free_y") +
  facet_wrap(~Sample, nrow = 3) +
  ggtitle("Mismatch type") +
  theme(plot.title = element_text(hjust = 0.5),axis.text.x = element_text(angle = 90, hjust = 1))
ggsave("../report/images_assembly_zvl/mt_mean.png", width = 2250, height = 1500, units = "px")

rm (mtype2, mtype2_2)

#Number of deletions per sample

delets2 <- data.frame(matrix(ncol = 2, nrow = 0))
colnames(delets2)<- c("Nucleotide", "Sample")
for (i in samples) {
  temp <- table_set[table_set$Sample==i,][,10]
  temp <- unlist(strsplit(temp,""))
  temp<- temp[temp %in% c(letters, LETTERS)]
  temp<-cbind("Nucleotide"=temp, "Sample" = i)
  temp<-data.frame(temp)
  temp$Nucleotide <- as.factor(temp$Nucleotide)
  delets2<-rbind(temp, delets2)
}
delets2$Sample<-as.factor(delets2$Sample)
rm (temp)

ggplot(delets2) +
  geom_bar(aes(x = Nucleotide), width = 0.4, fill = "blue", alpha = 0.5, col = "black") +
  facet_wrap(~Sample, nrow = 3) +
  #ylim(0, 1500) +
  ggtitle ("Deletions") +
  theme(plot.title = element_text(hjust = 0.5))
ggsave("../report/images_assembly_zvl/delets.png", width = 2250, height = 1500, units = "px")

rm (delets2)

#Number of insertions per sample

inserts2 <- data.frame(matrix(ncol = 2, nrow = 0))
colnames(inserts2)<- c("Nucleotide", "Sample")
for (i in samples) {
  temp <- table_set[table_set$Sample==i,][,11]
  temp <- unlist(strsplit(temp,""))
  temp<- temp[temp %in% c(letters, LETTERS)]
  temp<-cbind("Nucleotide"=temp, "Sample" = i)
  temp<-data.frame(temp)
  temp$Nucleotide <- as.factor(temp$Nucleotide)
  inserts2<-rbind(temp, inserts2)
}
inserts2$Sample<-as.factor(inserts2$Sample)

ggplot(inserts2) +
  ggtitle("Insertions") +
  geom_bar(aes(x = Nucleotide), width = 0.4, fill = "blue", alpha = 0.5, col = "black") +
  facet_wrap(~Sample, nrow = 2) +
  #ylim(0, 2000) +
  theme(plot.title = element_text(hjust = 0.5))
ggsave("../report/images_assembly_zvl/inserts.png", width = 2250, height = 1500, units = "px")

rm (inserts2)

#Mean, min and max values of insertions per sample

insert_mean2 <- data.frame(matrix(ncol = 7, nrow = 0))
colnames(insert_mean2)<- c("Sample", "Min.", "1st Qu.", "Median", "Mean", "3rd Qu.", "Max.")
for (i in samples) {
  temp<-subset(table_set, Mean_length_insert!=0 & Sample == i)[,7]
  temp<-cbind("Sample"=i, t(summary(temp)))
  insert_mean2<-rbind(insert_mean2, temp)
}
insert_mean2[, 2]<-as.numeric(insert_mean2[, 2])
insert_mean2[, 4]<-as.numeric(insert_mean2[, 4])
insert_mean2[, 5]<-as.numeric(insert_mean2[, 5])
insert_mean2[, 7]<-as.numeric(insert_mean2[, 7])
ggplot(data.frame(insert_mean2), aes(x = factor(Sample), y = Mean), size = 0.5) +
  geom_pointrange(aes(ymin = 0, ymax = Max.)) +
  #scale_y_continuous(breaks = c(0:12), limits = c(0, 13)) +
  geom_point(aes(x = factor(Sample), y = 0), shape = 1) +
  geom_point(aes(x = factor(Sample), y = Max.), shape = 1) +
  ggtitle("Insertion mean length") +
  scale_x_discrete(name = "Sample") +
  #scale_y_continuous(name = "Length", breaks = seq(0, 18, 1)) +
  theme(plot.title = element_text(hjust = 0.5))
ggsave("../report/images_assembly_zvl/iml.png", width = 2250, height = 1500, units = "px")

rm (insert_mean2)

#Mean, min and max values of deletions per sample

delete_mean2 <- data.frame(matrix(ncol = 7, nrow = 0))
colnames(delete_mean2)<- c("Sample", "Min.", "1st Qu.", "Median", "Mean", "3rd Qu.", "Max.")
for (i in samples) {
  temp<-subset(table_set, Mean_length_del!=0 & Sample == i)[,6]
  temp<-cbind("Sample"=i, t(summary(temp)))
  delete_mean2<-rbind(delete_mean2, temp)
}
delete_mean2[, 2]<-as.numeric(delete_mean2[, 2])
delete_mean2[, 4]<-as.numeric(delete_mean2[, 4])
delete_mean2[, 5]<-as.numeric(delete_mean2[, 5])
delete_mean2[, 7]<-as.numeric(delete_mean2[, 7])
ggplot(data.frame(delete_mean2), aes(x = factor(Sample), y = Mean), size = 0.5) +
  geom_pointrange(aes(ymin = 0, ymax = Max.)) +
  #scale_y_continuous(breaks = c(0:22), limits = c(0, 22)) +
  geom_point(aes(x = factor(Sample), y = 0), shape = 1) +
  geom_point(aes(x = factor(Sample), y = Max.), shape = 1) +
  ggtitle("Deletion mean length") +
  scale_x_discrete(name = "Sample") +
  #scale_y_continuous(name = "Length", breaks = seq(0, 12, 1)) +
  theme(plot.title = element_text(hjust = 0.5))
ggsave("../report/images_assembly_zvl/dml.png", width = 2250, height = 1500, units = "px")

rm (delete_mean2)

#Number of mismatches per sample

ggplot(table_set, aes(x=Sample, y=Count_SNPs, group = 1)) +
  stat_summary(fun = sum, geom="bar", width = 0.7, fill = "white", color = "black", alpha = 0.7) +
  ggtitle("Mismatches per sample") +
  #geom_vline(xintercept = 50, linetype="dotted") + 
  theme(plot.title = element_text(hjust = 0.5))
ggsave("../report/images_assembly_zvl/mps_mean.png", width = 2250, height = 1500, units = "px")


ggplot(table_set, aes(x=Sample, y=Count_SNPs, group = 1)) +
  stat_summary(fun = mean, geom="point") +
  stat_summary(fun = mean, geom="line") +
  ggtitle("Mismatches per sample") +
  #geom_vline(xintercept = 50, linetype="dotted") + 
  theme(plot.title = element_text(hjust = 0.5))
ggsave("../report/images_assembly_zvl/mps.png", width = 2250, height = 1500, units = "px")

#Number of deletions per cycle of sequence

numbers_only <- function(x) !grepl("\\D", x)

dels2 <- data.frame(matrix(nrow = 0, ncol = 2))
colnames(dels2)<-c("Num", "Sample")
for (i in samples) {
  temp <- table_set[table_set$Sample==i,][,10]
  temp <- unlist(strsplit(temp, split="|",fixed = TRUE))
  temp <- unlist(strsplit(temp,","))
  temp<- temp[numbers_only(temp)]
  temp<-cbind("Num"=temp, "Sample" = i)
  temp<-data.frame(temp)
  dels2<-rbind(temp, dels2)
}
dels2$Num<-as.numeric(dels2$Num)
ggplot(dels2, aes(x = Num)) +
  #geom_density(stat = "count", alpha = 0.5) +
  stat_count(geom="line", position="identity") +
  stat_count(geom="point", position="identity", size = 1) +
  facet_wrap(~Sample, nrow = 3) +
  #ylim(0, 1500) +
  ggtitle ("Deletions per cycle") +
  theme_bw()+
  theme(plot.title = element_text(hjust = 0.5))
ggsave("../report/images_assembly_zvl/dpc.png", width = 2250, height = 1500, units = "px")

rm (dels2, temp)

#Number of insertions per cycle of sequence

ins2 <- data.frame(matrix(nrow = 0, ncol = 2))
colnames(ins2)<-c("Num", "Sample")
for (i in samples) {
  temp <- table_set[table_set$Sample==i,][,11]
  temp <- unlist(strsplit(temp, split="|",fixed = TRUE))
  temp <- unlist(strsplit(temp,","))
  temp<- temp[numbers_only(temp)]
  temp<-cbind("Num"=temp, "Sample" = i)
  temp<-data.frame(temp)
  ins2<-rbind(temp, ins2)
}
ins2$Num<-as.numeric(ins2$Num)
ggplot(ins2, aes(x = Num)) +
  stat_count(geom="line", position="identity") +
  stat_count(geom="point", position="identity", size = 1) +
  facet_wrap(~Sample, nrow = 3) +
  #ylim(0, 1500) +
  ggtitle ("Insertions per cycle") +
  scale_x_continuous(breaks=seq(0, 150, 25))+
  theme_bw()+
  theme(plot.title = element_text(hjust = 0.5)) 
ggsave("../report/images_assembly_zvl/ipc.png", width = 2250, height = 1500, units = "px")
rm (ins2, temp)

#Number of deletions per sample

parse_delets2 <- data.frame(matrix(nrow = 0, ncol = 2))
colnames(parse_delets2)<-c("Num", "Sample")
for (i in samples) {
  temp <- table_set[table_set$Sample==i,][,10]
  temp <- unlist(strsplit(temp,split="|",fixed = TRUE))
  temp<-cbind("Num"=temp, "Sample" = i)
  temp<-data.frame(temp)
  parse_delets2<-rbind(temp, parse_delets2)
}
ggplot(parse_delets2, aes(x=Sample, group = 1)) +
  stat_count(geom="bar", width = 0.7, fill = "white", color = "black", alpha = 0.7) +
  ggtitle("Deletions per sample") +
  #geom_vline(xintercept = 50, linetype="dotted") + 
  theme(plot.title = element_text(hjust = 0.5))
ggsave("../report/images_assembly_zvl/dps.png", width = 2250, height = 1500, units = "px")
rm (parse_delets2, temp)

#Number of insertions per sample

parse_inserts2 <- data.frame(matrix(nrow = 0, ncol = 2))
colnames(parse_inserts2)<-c("Num", "Sample")
for (i in samples) {
  temp <- table_set[table_set$Sample==i,][,11]
  temp <- unlist(strsplit(temp,split="|",fixed = TRUE))
  temp<-cbind("Num"=temp, "Sample" = i)
  temp<-data.frame(temp)
  parse_inserts2<-rbind(temp, parse_inserts2)
}
ggplot(parse_inserts2, aes(x=Sample, group = 1)) +
  stat_count(geom="bar", width = 0.7, fill = "white", color = "black", alpha = 0.7) +
  ggtitle("Insertions per sample") +
  #geom_vline(xintercept = 50, linetype="dotted") + 
  theme(plot.title = element_text(hjust = 0.5))
ggsave("../report/images_assembly_zvl/ips.png", width = 2250, height = 1500, units = "px")
rm (parse_inserts2, temp)

#Frequency of insertion motifs per sample

insert_motifs2 <- data.frame(matrix(nrow = 0, ncol = 2))
colnames(insert_motifs2)<-c("Num", "Sample")
for (i in samples) {
  temp <- table_set[table_set$Sample==i,][,11]
  temp <- unlist(strsplit(temp, split="|",fixed = TRUE))
  temp <- unlist(strsplit(temp,","))
  temp <- temp[!numbers_only(temp)]
  temp <- cbind("Num"=temp, "Sample" = i)
  temp <- data.frame(temp)
  insert_motifs2 <- rbind(temp, insert_motifs2)
}
ggplot(insert_motifs2, aes(x = Num)) +
  stat_count(geom='bar', width = 0.7, fill = "white", color = "black", alpha = 0.7) +
  facet_wrap(~Sample, nrow=4, scales = "free_y") +
  ggtitle("Insertion motifs per sample") +
  theme(plot.title = element_text(hjust = 0.5),axis.text.x = element_text(angle = 90, hjust = 1))
ggsave("../report/images_assembly_zvl/imps.png", width = 8250, height = 2500, units = "px")
rm (temp, insert_motifs2)

#Frequency of deletion motifs per sample

delete_motifs2 <- data.frame(matrix(nrow = 0, ncol = 2))
colnames(delete_motifs2)<-c("Num", "Sample")
for (i in samples) {
  temp <- table_set[table_set$Sample==i,][,10]
  temp <- unlist(strsplit(temp, split="|",fixed = TRUE))
  temp <- unlist(strsplit(temp,","))
  temp <- temp[!numbers_only(temp)]
  temp <- cbind("Num"=temp, "Sample" = i)
  temp <- data.frame(temp)
  delete_motifs2 <- rbind(temp, delete_motifs2)
}
ggplot(delete_motifs2, aes(x = Num)) +
  stat_count(geom='bar', width = 0.7, fill = "white", color = "black", alpha = 0.7) +
  facet_wrap(~Sample, nrow=7, scales = "free_y") +
  ggtitle("Deletions motifs per sample") +
  theme(plot.title = element_text(hjust = 0.5),axis.text.x = element_text(angle = 90, hjust = 1))
ggsave("../report/images_assembly_zvl/dmps.png", width = 8250, height = 3500, units = "px")
rm (delete_motifs2, temp)