require(reshape);
require(ggplot2)

plotQualityHeatmap=function(assignment, cells, minlim=3, maxlim=150) {
  if(sum(assignment$QUAL>150)>0) {
    assignment[assignment$QUAL>150,]$QUAL = 150
  }
  qual=melt(assignment[assignment$Barcode %in% cells,][,c(1,3)],id.vars="Barcode")
  qual$value=as.numeric(qual$value)
  qual=qual[order(qual$value,decreasing=T),]
  qual_reorder=qual
  qual_reorder$Barcode = c(1:nrow(qual_reorder))
  p = ggplot(qual_reorder, aes(x = Barcode,y=variable)) +
    geom_raster(aes(fill = value), interpolate=TRUE) +
    #scale_fill_gradient2(low="grey26", mid="grey85", high="red", 
    #                     midpoint=10, limits=c(minlim, maxlim)) +
    scale_fill_gradientn(colours = c("#CEEBFA","#EEDCC5","#F4A473","#E67A56","#CD473E","#AE1D3E"), limits=c(minlim, maxlim))+
    theme_classic()
  print(p)
}
convertCellVariantMatrixSortCellByQual = function(data,cells,variant_range=NA,assignment){
  qual=assignment[assignment$Barcode %in% cells,][,c(1,3)]
  order = order(qual$QUAL,decreasing=T)
  
  if (!is.na(variant_range)) {
    matrix = data[variant_range,][,colnames(data) %in% cells]
    matrix = cbind(data[variant_range,]$variant,matrix[,order])
  } else {
    matrix = data[,colnames(data) %in% cells]
    matrix = cbind(data$variant,matrix[,order])
  }
  #print(head(matrix[,c(1:5)]))
  colnames(matrix) = c("variant", 1:(ncol(matrix)-1))
  matrix$variant = c(1:nrow(matrix))
  #print(head(matrix[,c(1:5)]))
  df_heatmap=melt(matrix,id.vars="variant")
  #print(head(df_heatmap))
  df_heatmap$variant = as.numeric(df_heatmap$variant)
  df_heatmap$variable = as.numeric(df_heatmap$variable)
  df_heatmap$value=as.numeric(df_heatmap$value)
  #print(head(df_heatmap))
  return(df_heatmap)
}
convertCellVariantMatrix = function(data,cells,variant_range=NA){
  if (!is.na(variant_range)) {
    matrix = data[variant_range,][,colnames(data) %in% cells]
    matrix = cbind(data[variant_range,]$variant,matrix)
  } else {
    matrix = data[,colnames(data) %in% cells]
    matrix = cbind(data$variant,matrix)
  }
  #print(head(matrix[,c(1:5)]))
  colnames(matrix) = c("variant", 1:(ncol(matrix)-1))
  matrix$variant = c(1:nrow(matrix))
  #print(head(matrix[,c(1:5)]))
  df_heatmap=melt(matrix,id.vars="variant")
  #print(head(df_heatmap))
  df_heatmap$variant = as.numeric(df_heatmap$variant)
  df_heatmap$variable = as.numeric(df_heatmap$variable)
  df_heatmap$value=as.numeric(df_heatmap$value)
  #print(head(df_heatmap))
  return(df_heatmap)
}

sample = "pre"

data = read.table(paste("../",sample,".scaf.tsv",sep=""),header=T)

assignment = read.table(paste("../expected_results/",sample,".assign_results.tsv",sep=""),header = T)
assignment$Barcode = paste("X",assignment$Barcode,sep="")
rownames(assignment) = assignment$Barcode

assignment$QUAL[assignment$ASIG=="UNASSIGN"] = NA
assignment$QUAL=as.numeric(as.character(assignment$QUAL))
assignment$QUAL[is.na(assignment$QUAL)] = 0

if (sum(assignment$QUAL=="Inf")>0){
  assignment[assignment$QUAL == Inf,]$QUAL = max(assignment[assignment$QUAL!=Inf,]$QUAL)
}
max(assignment$QUAL)
head(assignment)

## genotyping assigned cells ##
cells = rownames(assignment[assignment$ASIG == "SC1",])
length(cells) #0

cells = rownames(assignment[assignment$ASIG == "SC2",])
length(cells)
df2 = convertCellVariantMatrixSortCellByQual(data, cells, assignment=assignment)
df2$variable = df2$variable
max(df2$variable)
pdf("pretreatment_SC2_cell_qual_sort.pdf",width=4,height=2)
plotQualityHeatmap(assignment, cells)
dev.off()

cells = rownames(assignment[assignment$ASIG == "SC4",])
length(cells)
df3 = convertCellVariantMatrixSortCellByQual(data, cells, assignment=assignment)
df3$variable = df3$variable + max(df2$variable)
max(df3$variable)
pdf("pretreatment_SC4_cell_qual_sort.pdf",width=4,height=2)
plotQualityHeatmap(assignment, cells)
dev.off()


cells = rownames(assignment[assignment$ASIG == "SC3",])
length(cells)
df4 = convertCellVariantMatrixSortCellByQual(data, cells, assignment=assignment)
df4$variable = df4$variable + max(df3$variable)
max(df4$variable)
pdf("pretreatment_SC3_cell_qual_sort.pdf",width=4,height=2)
plotQualityHeatmap(assignment, cells)
dev.off()


cells = rownames(assignment[assignment$ASIG == "SC5",])
length(cells)
df5 = convertCellVariantMatrixSortCellByQual(data, cells, assignment=assignment)
df5$variable = df5$variable + max(df4$variable)
max(df5$variable)
pdf("pretreatment_SC5_cell_qual_sort.pdf",width=4,height=2)
plotQualityHeatmap(assignment, cells)
dev.off()

cells = rownames(assignment[assignment$ASIG == "UNASSIGN",])
length(cells)
df6 = convertCellVariantMatrix(data,cells)
df6$variable = df6$variable + max(df5$variable)
max(df6$variable)

##### combine plot ######
pdf("pretreatment cell genotype_sort_qual.pdf",width=6,height=8)
par(mfrow=c(1,1))
df_combined_1 = rbind(df2,df3,df4,df5,df6)
b = df_combined_1[df_combined_1$value==0,]
c = df_combined_1[df_combined_1$value>0,]
c1_num = 4043
c2_num = 663
c3_num = 628
c4_num = 958
plot(b$variable, b$variant,pch=16, cex=0.7,col="#00FF0080",ylim=c(1,max(df_combined_1$variant)),xlim=c(1,max(df_combined_1$variable)),
     xlab="Cell index",ylab="Variant index", main="Genotype of pretreatment cells")
points(c$variable, c$variant,pch=16, cex=0.7,col="#0000FF80",ylim=c(1,max(df_combined_1$variant)),xlim=c(1,max(df_combined_1$variable)))
abline(v=max(df2$variable) + 0.5)
abline(v=max(df3$variable) + 0.5)
abline(v=max(df4$variable) + 0.5)
abline(v=max(df5$variable) + 0.5)
abline(h = c1_num +0.5)
abline(h = c1_num+c2_num+0.5)
abline(h = c1_num+c2_num+c4_num+0.5)
abline(h = c1_num+c2_num+c4_num+c3_num+0.5)
dev.off()



#################################################
#################################################
#################################################
sample = "post"

data = read.table(paste("../",sample,".scaf.tsv",sep=""),header=T)

assignment = read.table(paste("../expected_results/",sample,".assign_results.tsv",sep=""),header = T)
assignment$Barcode = paste("X",assignment$Barcode,sep="")
rownames(assignment) = assignment$Barcode

assignment$QUAL[assignment$ASIG=="UNASSIGN"] = NA
assignment$QUAL=as.numeric(as.character(assignment$QUAL))
assignment$QUAL[is.na(assignment$QUAL)] = 0

if (sum(assignment$QUAL=="Inf")>0){
  assignment[assignment$QUAL == Inf,]$QUAL = max(assignment[assignment$QUAL!=Inf,]$QUAL)
}
max(assignment$QUAL)
head(assignment)

## genotyping assigned cells ##
cells = rownames(assignment[assignment$ASIG == "SC1",])
length(cells) #0


cells = rownames(assignment[assignment$ASIG == "SC2",])
length(cells) #0


cells = rownames(assignment[assignment$ASIG == "SC4",])
length(cells)
df3 = convertCellVariantMatrixSortCellByQual(data, cells, assignment=assignment)
df3$variable = df3$variable
max(df3$variable)
pdf("posttreatment_SC4_cell_qual_sort.pdf",width=4,height=2)
plotQualityHeatmap(assignment, cells)
dev.off()

cells = rownames(assignment[assignment$ASIG == "SC3",])
length(cells) 
x = data[,colnames(data) %in% cells]
head(x)
a = rep(1,nrow(data))
df4 = cbind(c(1:nrow(data)),a,x)
colnames(df4) = c("variant", "variable","value")
head(df4)
df4 = as.data.frame(df4) 
df4$variable = df4$variable + max(df3$variable)
max(df4$variable)
pdf("posttreatment_SC3_cell_qual_sort.pdf",width=4,height=2)
plotQualityHeatmap(assignment, cells)
dev.off()

cells = rownames(assignment[assignment$ASIG == "SC5",])
length(cells)
df5 = convertCellVariantMatrixSortCellByQual(data, cells, assignment=assignment)
df5$variable = df5$variable + max(df4$variable)
max(df5$variable)
pdf("posttreatment_SC5_cell_qual_sort.pdf",width=4,height=2)
plotQualityHeatmap(assignment, cells)
dev.off()

cells = rownames(assignment[assignment$ASIG == "normal",])
length(cells)
df6 = convertCellVariantMatrixSortCellByQual(data, cells, assignment=assignment)
df6$variable = df6$variable + max(df5$variable)
max(df6$variable)
pdf("posttreatment_normal_cell_qual_sort.pdf",width=4,height=2)
plotQualityHeatmap(assignment, cells)
dev.off()

cells = rownames(assignment[assignment$ASIG == "UNASSIGN",])
length(cells)
x = data[,colnames(data) %in% cells]
head(x)
a = rep(1,nrow(data))
df7 = cbind(c(1:nrow(data)),a,x)
colnames(df7) = c("variant", "variable","value")
head(df7)
df7 = as.data.frame(df7)
df7$variable= df7$variable+ max(df6$variable)
max(df7$variable)


##### combine plot ######
pdf("posttreatment cell genotype_sort_qual.pdf",width=6,height=8)
par(mfrow=c(1,1))
df_combined_1 = rbind(df3,df4,df5,df6,df7)
b = df_combined_1[df_combined_1$value==0,]
c = df_combined_1[df_combined_1$value>0,]
c1_num = 4043
c2_num = 663
c3_num = 628
c4_num = 958
plot(b$variable, b$variant,pch=16, cex=0.7,col="#00FF0080",ylim=c(1,max(df_combined_1$variant)),xlim=c(1,max(df_combined_1$variable)),
     xlab="Cell index",ylab="Variant index", main="Genotype of posttreatment cells")
points(c$variable, c$variant,pch=16, cex=0.7,col="#0000FF80",ylim=c(1,max(df_combined_1$variant)),xlim=c(1,max(df_combined_1$variable)))
abline(v=max(df3$variable) + 0.5)
abline(v=max(df4$variable) + 0.5)
abline(v=max(df5$variable) + 0.5)
abline(v=max(df6$variable) + 0.5)
abline(h = c1_num +0.5)
abline(h = c1_num+c2_num+0.5)
abline(h = c1_num+c2_num+c4_num+0.5)
abline(h = c1_num+c2_num+c4_num+c3_num+0.5)
dev.off()


