R/plot.R

Defines functions plotCountsPerTargetRelativeVcf plotCountsPerTargetVcf plotCountsPerSampleRelativeVcf plotCountsPerSampleVcf plotCountsPerTargetRelative plotCountsPerTarget plotCountsPerSampleRelative plotCountsPerSample plotCountsPerTargetRelativeVcf_plot plotCountsPerTargetVcf_plot plotCountsPerSampleRelativeVcf_plot plotCountsPerSampleVcf_plot plotCountsPerTargetRelative_plot plotCountsPerTargetRelative_col_insert plotCountsPerTargetRelative_insert plotCountsPerTargetRelative_col plotCountsPerTarget_plot plotCountsPerTarget_ref plotCountsPerTarget_col_insert plotCountsPerTarget_insert plotCountsPerTarget_col plotCountsPerTarget_abs plotCountsPerSampleRelative_plot plotCountsPerSampleRelative_col_insert plotCountsPerSampleRelative_insert plotCountsPerSampleRelative_col plotCountsPerSample_plot plotCountsPerSample_ref plotCountsPerSample_col_insert plotCountsPerSample_insert plotCountsPerSample_col plotCountsPerSample_abs

#Plot counts
plotCountsPerSample_abs<-function(calling,frequency_weight,long_insert,i,j,k,
                                  plot_abs){
    if(k==1){
        plot_abs[i,k]<-frequency_weight[[j]][i-long_insert,3+k]
        if(is.na(plot_abs[i,k])){
            plot_abs[i,k]<-0
        }
    }
    if(k==3){
        plot_abs[i,2]<-frequency_weight[[j]][i-long_insert,3+k]
        if(is.na(plot_abs[i,2])){
            plot_abs[i,2]<-0
        }
    }
    if(k==5){
        plot_abs[i,3]<-frequency_weight[[j]][i-long_insert,3+k]
        if(is.na(plot_abs[i,3])){
            plot_abs[i,3]<-0
        }
    }
    if(k==7){
        plot_abs[i,4]<-frequency_weight[[j]][i-long_insert,3+k]
        if(is.na(plot_abs[i,4])){
            plot_abs[i,4]<-0
        }
    }
    if(k==9){
        plot_abs[i,5]<-frequency_weight[[j]][i-long_insert,3+k]
        if(is.na(plot_abs[i,5])){
            plot_abs[i,5]<-0
        }
    }
    if(k==10){
        plot_abs[i,6]<-frequency_weight[[j]][i-long_insert,3+k]
        if(is.na(plot_abs[i,6])){
            plot_abs[i,6]<-0
        }
    }
    return(plot_abs)
}

plotCountsPerSample_col<-function(calling,frequency_weight,long_insert,
                                  qual_lower_bound,qual_upper_bound,
                                  A_col,C_col,G_col,T_col,i,j,k,plot_col){
    if(k==1&&!is.na(frequency_weight[[j]][i-long_insert,4+k])){
        if(round(frequency_weight[[j]][i-long_insert,4+k],digits=1)<qual_lower_bound){
            plot_col[i,1]<-"chartreuse"
        }
        if(round(frequency_weight[[j]][i-long_insert,4+k],digits=1)>qual_upper_bound){
            plot_col[i,1]<-"forestgreen"
        }
        if(qual_lower_bound<=round(frequency_weight[[j]][i-long_insert,4+k],digits=1)
           &&round(frequency_weight[[j]][i-long_insert,4+k],digits=1)<=qual_upper_bound){
            plot_col[i,1]<-A_col[round(frequency_weight[[j]][i-long_insert,4+k],digits=1)-qual_lower_bound+1]
        }
    }
    if(k==3&&!is.na(frequency_weight[[j]][i-long_insert,4+k])){
        if(round(frequency_weight[[j]][i-long_insert,4+k],digits=1)<qual_lower_bound){
            plot_col[i,2]<-"cyan"
        }
        if(round(frequency_weight[[j]][i-long_insert,4+k],digits=1)>qual_upper_bound){
            plot_col[i,2]<-"blue4"
        }
        if(qual_lower_bound<=round(frequency_weight[[j]][i-long_insert,4+k],digits=1)
           &&round(frequency_weight[[j]][i-long_insert,4+k],digits=1)<=qual_upper_bound){
            plot_col[i,2]<-C_col[round(frequency_weight[[j]][i-long_insert,4+k],digits=1)-qual_lower_bound+1]
        }
    }
    if(k==5&&!is.na(frequency_weight[[j]][i-long_insert,4+k])){
        if(round(frequency_weight[[j]][i-long_insert,4+k],digits=1)<qual_lower_bound){
            plot_col[i,3]<-"yellow"
        }
        if(round(frequency_weight[[j]][i-long_insert,4+k],digits=1)>qual_upper_bound){
            plot_col[i,3]<-"darkgoldenrod3"
        }
        if(qual_lower_bound<=round(frequency_weight[[j]][i-long_insert,4+k],digits=1)
           &&round(frequency_weight[[j]][i-long_insert,4+k],digits=1)<=qual_upper_bound){
            plot_col[i,3]<-G_col[round(frequency_weight[[j]][i-long_insert,4+k],digits=1)-qual_lower_bound+1]
        }
    }
    if(k==7&&!is.na(frequency_weight[[j]][i-long_insert,4+k])){
        if(round(frequency_weight[[j]][i-long_insert,4+k],digits=1)<qual_lower_bound){
            plot_col[i,4]<-"firebrick1"
        }
        if(round(frequency_weight[[j]][i-long_insert,4+k],digits=1)>qual_upper_bound){
            plot_col[i,4]<-"darkred"
        }
        if(qual_lower_bound<=round(frequency_weight[[j]][i-long_insert,4+k],digits=1)
           &&round(frequency_weight[[j]][i-long_insert,4+k],digits=1)<=qual_upper_bound){
            plot_col[i,4]<-T_col[round(frequency_weight[[j]][i-long_insert,4+k],digits=1)-qual_lower_bound+1]
        }
    }
    if(k==9&&!is.na(frequency_weight[[j]][i-long_insert,3+k])){
        plot_col[i,5]<-"black"
    }
    if(k==10&&!is.na(frequency_weight[[j]][i-long_insert,3+k])){
        plot_col[i,6]<-"purple"
    }
    return(plot_col)
}

plotCountsPerSample_insert<-function(calling,frequency_weight,long_insert,i,j,k,
                                     plot_insert){
    if(k==1){
        plot_insert[i,k]<-frequency_weight[[j]][i-long_insert,13+k]
        if(is.na(plot_insert[i,k])){
            plot_insert[i,k]<-0
        }
    }
    if(k==3){
        plot_insert[i,2]<-frequency_weight[[j]][i-long_insert,13+k]
        if(is.na(plot_insert[i,2])){
            plot_insert[i,2]<-0
        }
    }
    if(k==5){
        plot_insert[i,3]<-frequency_weight[[j]][i-long_insert,13+k]
        if(is.na(plot_insert[i,3])){
            plot_insert[i,3]<-0
        }
    }
    if(k==7){
        plot_insert[i,4]<-frequency_weight[[j]][i-long_insert,13+k]
        if(is.na(plot_insert[i,4])){
            plot_insert[i,4]<-0
        }
    } 
    return(plot_insert)
}

plotCountsPerSample_col_insert<-function(calling,frequency_weight,long_insert,
                               qual_lower_bound,qual_upper_bound,
                               A_col,C_col,G_col,T_col,i,j,k,plot_col_insert){
    if(k==1&&!is.na(frequency_weight[[j]][i-long_insert,14+k])){
        if(round(frequency_weight[[j]][i-long_insert,14+k])<qual_lower_bound){
            plot_col_insert[i,1]<-"chartreuse"
        }
        if(round(frequency_weight[[j]][i-long_insert,14+k])>qual_upper_bound){
            plot_col_insert[i,1]<-"forestgreen"
        }
        if(qual_lower_bound<=round(frequency_weight[[j]][i-long_insert,14+k])
           &&round(frequency_weight[[j]][i-long_insert,14+k])<=qual_upper_bound){
            plot_col_insert[i,1]<-A_col[round(frequency_weight[[j]][i-long_insert,14+k])-qual_lower_bound+1]
        }
    }
    if(k==3&&!is.na(frequency_weight[[j]][i-long_insert,14+k])){
        if(round(frequency_weight[[j]][i-long_insert,14+k])<qual_lower_bound){
            plot_col_insert[i,1:2]<-"cyan"
        }
        if(round(frequency_weight[[j]][i-long_insert,14+k])>qual_upper_bound){
            plot_col_insert[i,1:2]<-"blue4"
        }
        if(qual_lower_bound<=round(frequency_weight[[j]][i-long_insert,14+k])
           &&round(frequency_weight[[j]][i-long_insert,14+k])<=qual_upper_bound){
            plot_col_insert[i,1:2]<-C_col[round(frequency_weight[[j]][i-long_insert,14+k])-qual_lower_bound+1]
        }
    }
    if(k==5&&!is.na(frequency_weight[[j]][i-long_insert,14+k])){
        if(round(frequency_weight[[j]][i-long_insert,14+k])<qual_lower_bound){
            plot_col_insert[i,1:3]<-"yellow"
        }
        if(round(frequency_weight[[j]][i-long_insert,14+k])>qual_upper_bound){
            plot_col_insert[i,1:3]<-"darkgoldenrod3"
        }
        if(qual_lower_bound<=round(frequency_weight[[j]][i-long_insert,14+k])
           &&round(frequency_weight[[j]][i-long_insert,14+k])<=qual_upper_bound){
            plot_col_insert[i,1:3]<-G_col[round(frequency_weight[[j]][i-long_insert,14+k])-qual_lower_bound+1]
        }
    }
    if(k==7&&!is.na(frequency_weight[[j]][i-long_insert,14+k])){
        if(round(frequency_weight[[j]][i-long_insert,14+k])<qual_lower_bound){
            plot_col_insert[i,1:4]<-"firebrick1"
        }
        if(round(frequency_weight[[j]][i-long_insert,14+k])>qual_upper_bound){
            plot_col_insert[i,1:4]<-"darkred"
        }
        if(qual_lower_bound<=round(frequency_weight[[j]][i-long_insert,14+k])
           &&round(frequency_weight[[j]][i-long_insert,14+k])<=qual_upper_bound){
            plot_col_insert[i,1:4]<-T_col[round(frequency_weight[[j]][i-long_insert,14+k])-qual_lower_bound+1]
        }
    }  
    return(plot_col_insert)
}

plotCountsPerSample_ref<-function(calling,dbsnp_file,vcf.comp,j,i,plot_ref){
    if(dbsnp_file!=""&&length(rowRanges(vcf.comp))==1
       &&length(rowRanges(vcf.comp)$REF)>0
       &&length(rowRanges(vcf.comp)$ALT)>0
       &&nchar(as.character(rowRanges(vcf.comp)$REF))==1
       &&nchar(as.character(rowRanges(vcf.comp)$ALT[[1]]))==1){
        ref<-as.character(rowRanges(vcf.comp)$REF)
        alt<-as.character(rowRanges(vcf.comp)$ALT[[1]][1])
        alt2<-""
        if(length(vcf.comp$ALT)==2){
            alt2<-(as.character(rowRanges(vcf.comp)$ALT[[2]][1]))
        }
        if(nchar(alt2)==0){
            if(ref=="A"){
                plot_ref[i,1]<--1
            }
            if(ref=="C"){
                plot_ref[i,2]<--1
            }
            if(ref=="G"){
                plot_ref[i,3]<--1
            }
            if(ref=="T"){
                plot_ref[i,4]<--1
            }
            if(alt=="A"){
                plot_ref[i,1]<--1
            }
            if(alt=="C"){
                plot_ref[i,2]<--1
            }
            if(alt=="G"){
                plot_ref[i,3]<--1
            }
            if(alt=="T"){
                plot_ref[i,4]<--1
            }
        }
        if(nchar(alt2)>0){
            if(alt2=="A"){
                plot_ref[i,1]<--1
            }
            if(alt2=="C"){
                plot_ref[i,2]<--1
            }
            if(alt2=="G"){
                plot_ref[i,3]<--1
            }
            if(alt2=="T"){
                plot_ref[i,4]<--1
            }
        }
    }
    if(dbsnp_file==""||length(rowRanges(vcf.comp))==0
       ||nchar(as.character(rowRanges(vcf.comp)$REF))!=1
       ||nchar(as.character(rowRanges(vcf.comp)$ALT[[1]]))!=1){
        if(calling[[j]][i,3]=="A"){
            plot_ref[i,1]<--1
        }
        if(calling[[j]][i,3]=="C"){
            plot_ref[i,2]<--1
        }
        if(calling[[j]][i,3]=="G"){
            plot_ref[i,3]<--1
        }
        if(calling[[j]][i,3]=="T"){
            plot_ref[i,4]<--1
        }
    } 
    return(plot_ref)
}

plotCountsPerSample_plot<-function(directory2,samples,plot_abs,plot_col,abs_max,
                                   calling,plot_insert,plot_col_insert,A_col,
                                   C_col,G_col,T_col,plot_ref,marks,
                                   qual_lower_bound,qual_upper_bound,j){
    if(directory2!=""){
        png(filename=paste(directory2,"/",samples[j,1],".png",sep=""),width=1600,height=800)
    }
    par(mar=c(16,6,2,1)) 
    barplot(t(plot_abs[[j]]),beside=TRUE,col=t(plot_col[[j]]),
            xlim=c(-1.2*length(plot_abs[[j]][,1]),7*length(plot_abs[[j]][,1])),
            ylim=c((-1*0.1*abs_max),abs_max),ylab="Number of reads",
            mgp=c(4.5,1,-2),main=samples[j,1],
            names.arg=paste(calling[[1]][,1],calling[[1]][,2],sep=";"),
            las=2,cex.axis=2,cex.names=2,cex.lab=2,cex.main=2)
    par(new=TRUE)
    barplot((t(plot_ref[[j]]))*0.1*abs_max,beside=TRUE,
            col=c("forestgreen","blue4","darkgoldenrod3","darkred","black","purple"),
            xlim=c(-1.2*length(plot_abs[[j]][,1]),7*length(plot_abs[[j]][,1])),
            ylim=c((-1*0.1*abs_max),abs_max),yaxt='n')
    for(i in 1:length(plot_abs[[j]][,1])){
        for(k in 4:1){
            rect(xleft=(7*i-1),xright=(7*i),ybottom=0,ytop=plot_insert[[j]][i,k],
                 col=plot_col_insert[[j]][i,k])
        }
        if(plot_insert[[j]][i,4]!=0)
            rect(xleft=(7*i-1),xright=(7*i),ybottom=0,ytop=plot_insert[[j]][i,4],
                 col=NA,border="purple",lwd=4)
    }
    counter<-0
    for(i in 1:length(plot_abs[[j]][,1])){
        if(length(marks)>0){
            for(k in 1:length(marks)){
                lines(x=c((1+counter*7),(7+counter*7)),
                      y=c((sum(plot_abs[[j]][i,1:5],na.rm=TRUE)*marks[k]),
                          (sum(plot_abs[[j]][i,1:5],na.rm=TRUE)*marks[k])))
            }
            counter<-counter+1
        }
    }
    text(x=-0.5*length(plot_abs[[j]][,1]),y=c(abs_max*0.8*(-0.1)),"Reference",cex=2)
    steps<-seq(1,length(plot_abs[[j]][,1])*0.65,
               by=(length(plot_abs[[j]][,1])*0.65-1)/length(A_col))-1.2*length(plot_abs[[j]][,1])
    if(length(plot_abs[[j]][,1])==1){
        steps<-seq(1*0.65,1,by=(1-1*0.65)/length(A_col))-1.75*1
    }
    if(length(plot_abs[[j]][,1])==2){
        steps<-seq(0.7,2*0.65,by=(2*0.65-0.7)/length(A_col))-1.2*2
    }
    rect(xleft=steps[1:length(steps)-1],xright=steps[2:length(steps)],
         ybottom=abs_max*0.88,ytop=abs_max*0.93,lwd=0,col=A_col,border=NA)
    rect(xleft=steps[1:length(steps)-1],xright=steps[2:length(steps)],
         ybottom=abs_max*0.81,ytop=abs_max*0.86,lwd=0,col=C_col,border=NA)
    rect(xleft=steps[1:length(steps)-1],xright=steps[2:length(steps)],
         ybottom=abs_max*0.74,ytop=abs_max*0.79,lwd=0,col=G_col,border=NA)
    rect(xleft=steps[1:length(steps)-1],xright=steps[2:length(steps)],
         ybottom=abs_max*0.67,ytop=abs_max*0.72,lwd=0,col=T_col,border=NA)
    text(x=length(plot_abs[[j]][,1])*0.75-1.2*length(plot_abs[[j]][,1]),
         y=c(abs_max*0.90),"A",cex=2)
    text(x=length(plot_abs[[j]][,1])*0.75-1.2*length(plot_abs[[j]][,1]),
         y=c(abs_max*0.83),"C",cex=2)
    text(x=length(plot_abs[[j]][,1])*0.75-1.2*length(plot_abs[[j]][,1]),
         y=c(abs_max*0.76),"G",cex=2)
    text(x=length(plot_abs[[j]][,1])*0.75-1.2*length(plot_abs[[j]][,1]),
         y=c(abs_max*0.69),"T",cex=2)
    text(x=length(plot_abs[[j]][,1])*0.78-1.2*length(plot_abs[[j]][,1]),
         y=c(abs_max*0.62),"Del",cex=2)
    text(x=length(plot_abs[[j]][,1])*0.78-1.2*length(plot_abs[[j]][,1]),
         y=c(abs_max*0.55),"Ins",cex=2)
    text(x=length(plot_abs[[j]][,1])*1.15-1.2*length(plot_abs[[j]][,1]),
         y=c(abs_max*0.97),"Mean quality",cex=2)
    if(length(plot_abs[[j]][,1])>2){
        rect(xleft=1-1.2*length(plot_abs[[j]][,1]),
             xright=length(plot_abs[[j]][,1])*0.65-1.2*length(plot_abs[[j]][,1]),
             ybottom=abs_max*0.60,ytop=abs_max*0.65,lwd=0,col="black")
        rect(xleft=1-1.2*length(plot_abs[[j]][,1]),
             xright=length(plot_abs[[j]][,1])*0.65-1.2*length(plot_abs[[j]][,1]),
             ybottom=abs_max*0.53,ytop=abs_max*0.58,lwd=2,
             border="purple",col="white")
        text(x=1-1.2*length(plot_abs[[j]][,1]),y=c(abs_max*0.97),
             qual_lower_bound,cex=2)
        text(x=length(plot_abs[[j]][,1])*0.65-1.2*length(plot_abs[[j]][,1]),
             y=c(abs_max*0.97),qual_upper_bound,cex=2)
    }
    if(length(plot_abs[[j]][,1])==1){
        rect(xleft=-1.1,xright=-0.75,ybottom=abs_max*0.60,
             ytop=abs_max*0.65,lwd=0,col="black")
        rect(xleft=-1.1,xright=-0.75,ybottom=abs_max*0.53,
             ytop=abs_max*0.58,lwd=2,border="purple",col="white")
        text(x=-1.1,y=c(abs_max*0.97),qual_lower_bound,cex=2)
        text(x=-0.75,y=c(abs_max*0.97),qual_upper_bound,cex=2)
    }
    if(length(plot_abs[[j]][,1])==2){
        rect(xleft=-1.7,xright=-1.1,ybottom=abs_max*0.60,ytop=abs_max*0.65,
             lwd=0,col="black")
        rect(xleft=-1.7,xright=-1.1,ybottom=abs_max*0.53,ytop=abs_max*0.58,
             lwd=2,border="purple",col="white")
        text(x=-1.7,y=c(abs_max*0.97),qual_lower_bound,cex=2)
        text(x=-1.1,y=c(abs_max*0.97),qual_upper_bound,cex=2)
    }
    if(directory2!=""){
        dev.off()
    }
    return()
}


plotCountsPerSampleRelative_col<-function(calling,frequency_weight,long_insert,
                                          qual_lower_bound,qual_upper_bound,
                                          A_col,C_col,G_col,T_col,i,j,k,plot_col){
    if(k==4&&!is.na(frequency_weight[[j]][i-long_insert,1+k])){
        if(round(frequency_weight[[j]][i-long_insert,1+k],digits=1)<qual_lower_bound){
            plot_col[i,1]<-"chartreuse"
        }
        if(round(frequency_weight[[j]][i-long_insert,1+k],digits=1)>qual_upper_bound){
            plot_col[i,1]<-"forestgreen"
        }
        if(qual_lower_bound<=round(frequency_weight[[j]][i-long_insert,1+k],digits=1)&&round(frequency_weight[[j]][i-long_insert,1+k],digits=1)<=qual_upper_bound){
            plot_col[i,1]<-A_col[round(frequency_weight[[j]][i-long_insert,1+k],digits=1)-qual_lower_bound+1]
        }
    }
    if(k==5&&!is.na(frequency_weight[[j]][i-long_insert,2+k])){
        if(round(frequency_weight[[j]][i-long_insert,2+k],digits=1)<qual_lower_bound){
            plot_col[i,2]<-"cyan"
        }
        if(round(frequency_weight[[j]][i-long_insert,2+k],digits=1)>qual_upper_bound){
            plot_col[i,2]<-"blue4"
        }
        if(qual_lower_bound<=round(frequency_weight[[j]][i-long_insert,2+k],digits=1)&&round(frequency_weight[[j]][i-long_insert,2+k],digits=1)<=qual_upper_bound){
            plot_col[i,2]<-C_col[round(frequency_weight[[j]][i-long_insert,2+k],digits=1)-qual_lower_bound+1]
        }
    }
    if(k==6&&!is.na(frequency_weight[[j]][i-long_insert,3+k])){
        if(round(frequency_weight[[j]][i-long_insert,3+k],digits=1)<qual_lower_bound){
            plot_col[i,3]<-"yellow"
        }
        if(round(frequency_weight[[j]][i-long_insert,3+k],digits=1)>qual_upper_bound){
            plot_col[i,3]<-"darkgoldenrod3"
        }
        if(qual_lower_bound<=round(frequency_weight[[j]][i-long_insert,3+k],digits=1)&&round(frequency_weight[[j]][i-long_insert,3+k],digits=1)<=qual_upper_bound){
            plot_col[i,3]<-G_col[round(frequency_weight[[j]][i-long_insert,3+k],digits=1)-qual_lower_bound+1]
        }
    }
    if(k==7&&!is.na(frequency_weight[[j]][i-long_insert,4+k])){
        if(round(frequency_weight[[j]][i-long_insert,4+k],digits=1)<qual_lower_bound){
            plot_col[i,4]<-"firebrick1"
        }
        if(round(frequency_weight[[j]][i-long_insert,4+k],digits=1)>qual_upper_bound){
            plot_col[i,4]<-"darkred"
        }
        if(qual_lower_bound<=round(frequency_weight[[j]][i-long_insert,4+k],digits=1)&&round(frequency_weight[[j]][i-long_insert,4+k],digits=1)<=qual_upper_bound){
            plot_col[i,4]<-T_col[round(frequency_weight[[j]][i-long_insert,4+k],digits=1)-qual_lower_bound+1]
        }
    }
    if(k==8&&!is.na(frequency_weight[[j]][i-long_insert,4+k])){
        plot_col[i,5]<-"black"
    }
    if(k==9&&!is.na(frequency_weight[[j]][i-long_insert,4+k])){
        plot_col[i,6]<-"purple"
    }
    return(plot_col)
}

plotCountsPerSampleRelative_insert<-function(calling,frequency_weight,
                                             long_insert,i,j,k,plot_insert){
    if(k==1&&!(is.na(frequency_weight[[j]][i-long_insert,4])
               &&is.na(frequency_weight[[j]][i-long_insert,6])
               &&is.na(frequency_weight[[j]][i-long_insert,8])
               &&is.na(frequency_weight[[j]][i-long_insert,10])
               &&is.na(frequency_weight[[j]][i-long_insert,12]))){
        plot_insert[i,k]<-frequency_weight[[j]][i-long_insert,13+k]/sum(frequency_weight[[j]][i-long_insert,c(4,6,8,10,12)],na.rm=TRUE)
        if(is.na(plot_insert[i,k])){
            plot_insert[i,k]<-0
        }
    }
    if(k==3&&!(is.na(frequency_weight[[j]][i-long_insert,4])
               &&is.na(frequency_weight[[j]][i-long_insert,6])
               &&is.na(frequency_weight[[j]][i-long_insert,8])
               &&is.na(frequency_weight[[j]][i-long_insert,10])
               &&is.na(frequency_weight[[j]][i-long_insert,12]))){
        plot_insert[i,2]<-frequency_weight[[j]][i-long_insert,13+k]/sum(frequency_weight[[j]][i-long_insert,c(4,6,8,10,12)],na.rm=TRUE)
        if(is.na(plot_insert[i,2])){
            plot_insert[i,2]<-0
        }
    }
    if(k==5&&!(is.na(frequency_weight[[j]][i-long_insert,4])
               &&is.na(frequency_weight[[j]][i-long_insert,6])
               &&is.na(frequency_weight[[j]][i-long_insert,8])
               &&is.na(frequency_weight[[j]][i-long_insert,10])
               &&is.na(frequency_weight[[j]][i-long_insert,12]))){
        plot_insert[i,3]<-frequency_weight[[j]][i-long_insert,13+k]/sum(frequency_weight[[j]][i-long_insert,c(4,6,8,10,12)],na.rm=TRUE)
        if(is.na(plot_insert[i,3])){
            plot_insert[i,3]<-0
        }
    }
    if(k==7&&!(is.na(frequency_weight[[j]][i-long_insert,4])
               &&is.na(frequency_weight[[j]][i-long_insert,6])
               &&is.na(frequency_weight[[j]][i-long_insert,8])
               &&is.na(frequency_weight[[j]][i-long_insert,10])
               &&is.na(frequency_weight[[j]][i-long_insert,12]))){
        plot_insert[i,4]<-frequency_weight[[j]][i-long_insert,13+k]/sum(frequency_weight[[j]][i-long_insert,c(4,6,8,10,12)],na.rm=TRUE)
        if(is.na(plot_insert[i,4])){
            plot_insert[i,4]<-0
        }
    }
    return(plot_insert)
}

plotCountsPerSampleRelative_col_insert<-function(calling,frequency_weight,
                                                 long_insert,qual_lower_bound,
                                                 qual_upper_bound,
                                                 A_col,C_col,G_col,T_col,i,j,k,
                                                 plot_col_insert){
    if(k==1&&!is.na(frequency_weight[[j]][i-long_insert,14+k])){
        if(round(frequency_weight[[j]][i-long_insert,14+k])<qual_lower_bound){
            plot_col_insert[[j]][i,1]<-"chartreuse"
        }
        if(round(frequency_weight[[j]][i-long_insert,14+k])>qual_upper_bound){
            plot_col_insert[[j]][i,1]<-"forestgreen"
        }
        if(qual_lower_bound<=round(frequency_weight[[j]][i-long_insert,14+k])
           &&round(frequency_weight[[j]][i-long_insert,14+k])<=qual_upper_bound){
            plot_col_insert[[j]][i,1]<-A_col[round(frequency_weight[[j]][i-long_insert,14+k])-qual_lower_bound+1]
        }
    }
    if(k==3&&!is.na(frequency_weight[[j]][i-long_insert,14+k])){
        if(round(frequency_weight[[j]][i-long_insert,14+k])<qual_lower_bound){
            plot_col_insert[[j]][i,1:2]<-"cyan"
        }
        if(round(frequency_weight[[j]][i-long_insert,14+k])>qual_upper_bound){
            plot_col_insert[[j]][i,1:2]<-"blue4"
        }
        if(qual_lower_bound<=round(frequency_weight[[j]][i-long_insert,14+k])
           &&round(frequency_weight[[j]][i-long_insert,14+k])<=qual_upper_bound){
            plot_col_insert[[j]][i,1:2]<-C_col[round(frequency_weight[[j]][i-long_insert,14+k])-qual_lower_bound+1]
        }
    }
    if(k==5&&!is.na(frequency_weight[[j]][i-long_insert,14+k])){
        if(round(frequency_weight[[j]][i-long_insert,14+k])<qual_lower_bound){
            plot_col_insert[[j]][i,1:3]<-"yellow"
        }
        if(round(frequency_weight[[j]][i-long_insert,14+k])>qual_upper_bound){
            plot_col_insert[[j]][i,1:3]<-"darkgoldenrod3"
        }
        if(qual_lower_bound<=round(frequency_weight[[j]][i-long_insert,14+k])
           &&round(frequency_weight[[j]][i-long_insert,14+k])<=qual_upper_bound){
            plot_col_insert[[j]][i,1:3]<-G_col[round(frequency_weight[[j]][i-long_insert,14+k])-qual_lower_bound+1]
        }
    }
    if(k==7&&!is.na(frequency_weight[[j]][i-long_insert,14+k])){
        if(round(frequency_weight[[j]][i-long_insert,14+k])<qual_lower_bound){
            plot_col_insert[[j]][i,1:4]<-"firebrick1"
        }
        if(round(frequency_weight[[j]][i-long_insert,14+k])>qual_upper_bound){
            plot_col_insert[[j]][i,1:4]<-"darkred"
        }
        if(qual_lower_bound<=round(frequency_weight[[j]][i-long_insert,14+k])
           &&round(frequency_weight[[j]][i-long_insert,14+k])<=qual_upper_bound){
            plot_col_insert[[j]][i,1:4]<-T_col[round(frequency_weight[[j]][i-long_insert,14+k])-qual_lower_bound+1]
        }
    }
    return(plot_col_insert)
}

plotCountsPerSampleRelative_plot<-function(directory2,samples,plot_rel,plot_col,
                                           rel_max,calling,plot_insert,
                                           plot_col_insert,A_col,C_col,G_col,
                                           T_col,plot_ref,marks,
                                           qual_lower_bound,qual_upper_bound,j){
    if(directory2!=""){
        png(filename=paste(directory2,"/",samples[j,1],".png",sep=""),width=1600,height=800)
    }
    par(mar=c(16,6,2,1)) 
    barplot(t(plot_rel[[j]]),beside=TRUE,col=t(plot_col[[j]]),
            xlim=c(-1.2*length(plot_rel[[j]][,1]),7*length(plot_rel[[j]][,1])),
            ylim=c((-1*0.1*rel_max),rel_max),ylab="Relative frequency",
            mgp=c(4.5,1,-2),main=samples[j,1],names.arg=
                paste(calling[[1]][,1],calling[[1]][,2],sep=";"),las=2,
            cex.axis=2,cex.names=2,cex.lab=2,cex.main=2)
    par(new=TRUE)
    barplot((t(plot_ref[[j]]))*0.1*rel_max,beside=TRUE,
            col=c("forestgreen","blue4","darkgoldenrod3","darkred","black","purple"),
            xlim=c(-1.2*length(plot_rel[[j]][,1]),7*length(plot_rel[[j]][,1])),
            ylim=c((-1*0.1*rel_max),rel_max),yaxt='n')
    for(i in 1:length(plot_rel[[j]][,1])){
        for(k in 4:1){
            rect(xleft=(7*i-1),xright=(7*i),ybottom=0,
                 ytop=plot_insert[[j]][i,k],col=plot_col_insert[[j]][i,k])
        }
        if(plot_insert[[j]][i,4]!=0)
            rect(xleft=(7*i-1),xright=(7*i),ybottom=0,
                 ytop=plot_insert[[j]][i,4],col=NA,border="purple",lwd=4)
    }
    counter<-0
    for(i in 1:length(plot_rel[[j]][,1])){
        if(length(marks)>0){
            for(k in 1:length(marks)){
                lines(x=c((1+counter*7),(7+counter*7)),y=c(marks[k],marks[k]))
            }
            counter<-counter+1 
        }
    }
    text(x=-0.5*length(plot_rel[[j]][,1]),y=c(rel_max*0.8*(-0.1)),"Reference",cex=2)
    steps<-seq(1,length(plot_rel[[j]][,1])*0.65,
               by=(length(plot_rel[[j]][,1])*0.65-1)/length(A_col))-1.2*length(plot_rel[[j]][,1])
    if(length(plot_rel[[j]][,1])==1){
        steps<-seq(1*0.65,1,by=(1-1*0.65)/length(A_col))-1.75*1
    }
    if(length(plot_rel[[j]][,1])==2){
        steps<-seq(0.7,2*0.65,by=(2*0.65-0.7)/length(A_col))-1.2*2
    }
    rect(xleft=steps[1:length(steps)-1],xright=steps[2:length(steps)],
         ybottom=rel_max*0.88,ytop=rel_max*0.93,lwd=0,col=A_col,border=NA)
    rect(xleft=steps[1:length(steps)-1],xright=steps[2:length(steps)],
         ybottom=rel_max*0.81,ytop=rel_max*0.86,lwd=0,col=C_col,border=NA)
    rect(xleft=steps[1:length(steps)-1],xright=steps[2:length(steps)],
         ybottom=rel_max*0.74,ytop=rel_max*0.79,lwd=0,col=G_col,border=NA)
    rect(xleft=steps[1:length(steps)-1],xright=steps[2:length(steps)],
         ybottom=rel_max*0.67,ytop=rel_max*0.72,lwd=0,col=T_col,border=NA)
    text(x=length(plot_rel[[j]][,1])*0.75-1.2*length(plot_rel[[j]][,1]),
         y=c(rel_max*0.90),"A",cex=2)
    text(x=length(plot_rel[[j]][,1])*0.75-1.2*length(plot_rel[[j]][,1]),
         y=c(rel_max*0.83),"C",cex=2)
    text(x=length(plot_rel[[j]][,1])*0.75-1.2*length(plot_rel[[j]][,1]),
         y=c(rel_max*0.76),"G",cex=2)
    text(x=length(plot_rel[[j]][,1])*0.75-1.2*length(plot_rel[[j]][,1]),
         y=c(rel_max*0.69),"T",cex=2)
    text(x=length(plot_rel[[j]][,1])*0.78-1.2*length(plot_rel[[j]][,1]),
         y=c(rel_max*0.62),"Del",cex=2)
    text(x=length(plot_rel[[j]][,1])*0.78-1.2*length(plot_rel[[j]][,1]),
         y=c(rel_max*0.55),"Ins",cex=2)
    text(x=length(plot_rel[[j]][,1])*1.15-1.2*length(plot_rel[[j]][,1]),
         y=c(rel_max*0.97),"Mean quality",cex=2)
    if(length(plot_rel[[j]][,1])>2){
        rect(xleft=1-1.2*length(plot_rel[[j]][,1]),
             xright=length(plot_rel[[j]][,1])*0.65-1.2*length(plot_rel[[j]][,1]),
             ybottom=rel_max*0.60,ytop=rel_max*0.65,lwd=0,col="black")
        rect(xleft=1-1.2*length(plot_rel[[j]][,1]),
             xright=length(plot_rel[[j]][,1])*0.65-1.2*length(plot_rel[[j]][,1]),
             ybottom=rel_max*0.53,ytop=rel_max*0.58,lwd=2,border="purple",
             col="white")
        text(x=1-1.2*length(plot_rel[[j]][,1]),y=c(rel_max*0.97),
             qual_lower_bound,cex=2)
        text(x=length(plot_rel[[j]][,1])*0.65-1.2*length(plot_rel[[j]][,1]),
             y=c(rel_max*0.97),qual_upper_bound,cex=2)
    }
    if(length(plot_rel[[j]][,1])==1){
        rect(xleft=-1.1,xright=-0.75,ybottom=rel_max*0.60,ytop=rel_max*0.65,
             lwd=0,col="black")
        rect(xleft=-1.1,xright=-0.75,ybottom=rel_max*0.53,ytop=rel_max*0.58,
             lwd=2,border="purple",col="white")
        text(x=-1.1,y=c(rel_max*0.97),qual_lower_bound,cex=2)
        text(x=-0.75,y=c(rel_max*0.97),qual_upper_bound,cex=2)
    }
    if(length(plot_rel[[j]][,1])==2){
        rect(xleft=-1.7,xright=-1.1,ybottom=rel_max*0.60,ytop=rel_max*0.65,
             lwd=0,col="black")
        rect(xleft=-1.7,xright=-1.1,ybottom=rel_max*0.53,ytop=rel_max*0.58,
             lwd=2,border="purple",col="white")
        text(x=-1.7,y=c(rel_max*0.97),qual_lower_bound,cex=2)
        text(x=-1.1,y=c(rel_max*0.97),qual_upper_bound,cex=2)
    }
    if(directory2!=""){
        dev.off()
    }
}


plotCountsPerTarget_abs<-function(samples,frequency_weight,long_insert,i,j,k,
                                  plot_abs){
    if(k==1){
        plot_abs[i,k]<-frequency_weight[[i]][j-long_insert[i],3+k]
        if(is.na(plot_abs[i,k])){
            plot_abs[i,k]<-0
        }
    }
    if(k==3){
        plot_abs[i,2]<-frequency_weight[[i]][j-long_insert[i],3+k]
        if(is.na(plot_abs[i,2])){
            plot_abs[i,2]<-0
        }
    }
    if(k==5){
        plot_abs[i,3]<-frequency_weight[[i]][j-long_insert[i],3+k]
        if(is.na(plot_abs[i,3])){
            plot_abs[i,3]<-0
        }
    }
    if(k==7){
        plot_abs[i,4]<-frequency_weight[[i]][j-long_insert[i],3+k]
        if(is.na(plot_abs[i,4])){
            plot_abs[i,4]<-0
        }
    }
    if(k==9){
        plot_abs[i,5]<-frequency_weight[[i]][j-long_insert[i],3+k]
        if(is.na(plot_abs[i,5])){
            plot_abs[i,5]<-0
        }
    }
    if(k==10){
        plot_abs[i,6]<-frequency_weight[[i]][j-long_insert[i],3+k]
        if(is.na(plot_abs[i,6])){
            plot_abs[i,6]<-0
        }
    }
    return(plot_abs)
}

plotCountsPerTarget_col<-function(samples,frequency_weight,long_insert,
                                  qual_lower_bound,qual_upper_bound,
                                  A_col,C_col,G_col,T_col,i,j,k,plot_col){
    if(k==1&&!is.na(frequency_weight[[i]][j-long_insert[i],4+k])){
        if(round(frequency_weight[[i]][j-long_insert[i],4+k])<qual_lower_bound){
            plot_col[i,1]<-"chartreuse"
        }
        if(round(frequency_weight[[i]][j-long_insert[i],4+k])>qual_upper_bound){
            plot_col[i,1]<-"forestgreen"
        }
        if(qual_lower_bound<=round(frequency_weight[[i]][j-long_insert[i],4+k])
           &&round(frequency_weight[[i]][j-long_insert[i],4+k])<=qual_upper_bound){
            plot_col[i,1]<-A_col[round(frequency_weight[[i]][j-long_insert[i],4+k])-qual_lower_bound+1]
        }
    }
    if(k==3&&!is.na(frequency_weight[[i]][j-long_insert[i],4+k])){
        if(round(frequency_weight[[i]][j-long_insert[i],4+k])<qual_lower_bound){
            plot_col[i,2]<-"cyan"
        }
        if(round(frequency_weight[[i]][j-long_insert[i],4+k])>qual_upper_bound){
            plot_col[i,2]<-"blue4"
        }
        if(qual_lower_bound<=round(frequency_weight[[i]][j-long_insert[i],4+k])
           &&round(frequency_weight[[i]][j-long_insert[i],4+k])<=qual_upper_bound){
            plot_col[i,2]<-C_col[round(frequency_weight[[i]][j-long_insert[i],4+k])-qual_lower_bound+1]
        }
    }
    if(k==5&&!is.na(frequency_weight[[i]][j-long_insert[i],4+k])){
        if(round(frequency_weight[[i]][j-long_insert[i],4+k])<qual_lower_bound){
            plot_col[i,3]<-"yellow"
        }
        if(round(frequency_weight[[i]][j-long_insert[i],4+k])>qual_upper_bound){
            plot_col[i,3]<-"darkgoldenrod3"
        }
        if(qual_lower_bound<=round(frequency_weight[[i]][j-long_insert[i],4+k])
           &&round(frequency_weight[[i]][j-long_insert[i],4+k])<=qual_upper_bound){
            plot_col[i,3]<-G_col[round(frequency_weight[[i]][j-long_insert[i],4+k])-qual_lower_bound+1]
        }
    }
    if(k==7&&!is.na(frequency_weight[[i]][j-long_insert[i],4+k])){
        if(round(frequency_weight[[i]][j-long_insert[i],4+k])<qual_lower_bound){
            plot_col[i,4]<-"firebrick1"
        }
        if(round(frequency_weight[[i]][j-long_insert[i],4+k])>qual_upper_bound){
            plot_col[i,4]<-"darkred"
        }
        if(qual_lower_bound<=round(frequency_weight[[i]][j-long_insert[i],4+k])
           &&round(frequency_weight[[i]][j-long_insert[i],4+k])<=qual_upper_bound){
            plot_col[i,4]<-T_col[round(frequency_weight[[i]][j-long_insert[i],4+k])-qual_lower_bound+1]
        }
    }
    if(k==9&&!is.na(frequency_weight[[i]][j-long_insert[i],3+k])){
        plot_col[i,5]<-"black"
    }
    if(k==10&&!is.na(frequency_weight[[i]][j-long_insert[i],3+k])){
        plot_col[i,6]<-"purple"
    } 
    return(plot_col)
}

plotCountsPerTarget_insert<-function(samples,frequency_weight,long_insert,i,j,k,
                                     plot_insert){
    if(k==1){
        plot_insert[i,k]<-frequency_weight[[i]][j-long_insert[i],13+k]
        if(is.na(plot_insert[i,k])){
            plot_insert[i,k]<-0
        }
    }
    if(k==3){
        plot_insert[i,2]<-frequency_weight[[i]][j-long_insert[i],13+k]
        if(is.na(plot_insert[i,2])){
            plot_insert[i,2]<-0
        }
    }
    if(k==5){
        plot_insert[i,3]<-frequency_weight[[i]][j-long_insert[i],13+k]
        if(is.na(plot_insert[i,3])){
            plot_insert[i,3]<-0
        }
    }
    if(k==7){
        plot_insert[i,4]<-frequency_weight[[i]][j-long_insert[i],13+k]
        if(is.na(plot_insert[i,4])){
            plot_insert[i,4]<-0
        }
    }
    return(plot_insert)
}

plotCountsPerTarget_col_insert<-function(samples,frequency_weight,long_insert,
                                         qual_lower_bound,qual_upper_bound,
                                         A_col,C_col,G_col,T_col,i,j,k,
                                         plot_col_insert){
    if(k==1&&!is.na(frequency_weight[[i]][j-long_insert[i],14+k])){
        if(round(frequency_weight[[i]][j-long_insert[i],14+k])<qual_lower_bound){
            plot_col_insert[i,1]<-"chartreuse"
        }
        if(round(frequency_weight[[i]][j-long_insert[i],14+k])>qual_upper_bound){
            plot_col_insert[i,1]<-"forestgreen"
        }
        if(qual_lower_bound<=round(frequency_weight[[i]][j-long_insert[i],14+k])
           &&round(frequency_weight[[i]][j-long_insert[i],14+k])<=qual_upper_bound){
            plot_col_insert[i,1]<-A_col[round(frequency_weight[[i]][j-long_insert[i],14+k])-qual_lower_bound+1]
        }
    }
    if(k==3&&!is.na(frequency_weight[[i]][j-long_insert[i],14+k])){
        if(round(frequency_weight[[i]][j-long_insert[i],14+k])<qual_lower_bound){
            plot_col_insert[i,1:2]<-"cyan"
        }
        if(round(frequency_weight[[i]][j-long_insert[i],14+k])>qual_upper_bound){
            plot_col_insert[i,1:2]<-"blue4"
        }
        if(qual_lower_bound<=round(frequency_weight[[i]][j-long_insert[i],14+k])
           &&round(frequency_weight[[i]][j-long_insert[i],14+k])<=qual_upper_bound){
            plot_col_insert[i,1:2]<-C_col[round(frequency_weight[[i]][j-long_insert[i],14+k])-qual_lower_bound+1]
        }
    }
    if(k==5&&!is.na(frequency_weight[[i]][j-long_insert[i],14+k])){
        if(round(frequency_weight[[i]][j-long_insert[i],14+k])<qual_lower_bound){
            plot_col_insert[i,1:3]<-"yellow"
        }
        if(round(frequency_weight[[i]][j-long_insert[i],14+k])>qual_upper_bound){
            plot_col_insert[i,1:3]<-"darkgoldenrod3"
        }
        if(qual_lower_bound<=round(frequency_weight[[i]][j-long_insert[i],14+k])
           &&round(frequency_weight[[i]][j-long_insert[i],14+k])<=qual_upper_bound){
            plot_col_insert[i,1:3]<-G_col[round(frequency_weight[[i]][j-long_insert[i],14+k])-qual_lower_bound+1]
        }
    }
    if(k==7&&!is.na(frequency_weight[[i]][j-long_insert[i],14+k])){
        if(round(frequency_weight[[i]][j-long_insert[i],14+k])<qual_lower_bound){
            plot_col_insert[i,1:4]<-"firebrick1"
        }
        if(round(frequency_weight[[i]][j-long_insert[i],14+k])>qual_upper_bound){
            plot_col_insert[i,1:4]<-"darkred"
        }
        if(qual_lower_bound<=round(frequency_weight[[i]][j-long_insert[i],14+k])
           &&round(frequency_weight[[i]][j-long_insert[i],14+k])<=qual_upper_bound){
            plot_col_insert[i,1:4]<-T_col[round(frequency_weight[[i]][j-long_insert[i],14+k])-qual_lower_bound+1]
        }
    } 
    return(plot_col_insert)
}

plotCountsPerTarget_ref<-function(samples,calling,dbsnp_file,vcf.comp,j,i,
                                  plot_ref){
    if(dbsnp_file!=""&&length(rowRanges(vcf.comp))==1
       &&length(rowRanges(vcf.comp)$REF)>0
       &&length(rowRanges(vcf.comp)$ALT)>0
       &&nchar(as.character(rowRanges(vcf.comp)$REF))==1
       &&nchar(as.character(rowRanges(vcf.comp)$ALT[[1]]))==1){
        ref<-as.character(rowRanges(vcf.comp)$REF)
        alt<-as.character(rowRanges(vcf.comp)$ALT[[1]][1])
        alt2<-""
        if(length(vcf.comp$ALT)==2){
            alt2<-(as.character(rowRanges(vcf.comp)$ALT[[2]][1]))
        }
        if(nchar(alt2)==0){
            if(ref=="A"){
                plot_ref[i,1]<--1
            }
            if(ref=="C"){
                plot_ref[i,2]<--1
            }
            if(ref=="G"){
                plot_ref[i,3]<--1
            }
            if(ref=="T"){
                plot_ref[i,4]<--1
            }
            if(alt=="A"){
                plot_ref[i,1]<--1
            }
            if(alt=="C"){
                plot_ref[i,2]<--1
            }
            if(alt=="G"){
                plot_ref[i,3]<--1
            }
            if(alt=="T"){
                plot_ref[i,4]<--1
            }
        }
        if(nchar(alt2)>0){
            if(alt2=="A"){
                plot_ref[i,1]<--1
            }
            if(alt2=="C"){
                plot_ref[i,2]<--1
            }
            if(alt2=="G"){
                plot_ref[i,3]<--1
            }
            if(alt2=="T"){
                plot_ref[i,4]<--1
            }
        }
    }
    if(dbsnp_file==""||length(rowRanges(vcf.comp))==0
       ||nchar(as.character(rowRanges(vcf.comp)$REF))!=1
       ||nchar(as.character(rowRanges(vcf.comp)$ALT[[1]]))!=1){
        if(calling[[i]][j,3]=="A"){
            plot_ref[i,1]<--1
        }
        if(calling[[i]][j,3]=="C"){
            plot_ref[i,2]<--1
        }
        if(calling[[i]][j,3]=="G"){
            plot_ref[i,3]<--1
        }
        if(calling[[i]][j,3]=="T"){
            plot_ref[i,4]<--1
        }
    }  
    return(plot_ref)
}

plotCountsPerTarget_plot<-function(directory2,samples,plot_abs,plot_col,abs_max,
                                   calling,plot_insert,plot_col_insert,A_col,
                                   C_col,G_col,T_col,plot_ref,marks,
                                   qual_lower_bound,qual_upper_bound,counter,
                                   j){
    if(directory2!=""){
        if(j>1&&calling[[1]][j-1,1]==calling[[1]][j,1]
           &&calling[[1]][j-1,2]==calling[[1]][j,2]){
            counter<-counter+1
            png(filename=paste(directory2,"/",calling[[1]][j,1],";",
                           calling[[1]][j,2],"_",counter,".png",sep=""),
                width=1600,height=600)
        }
        if(j==1||calling[[1]][j-1,1]!=calling[[1]][j,1]
           ||calling[[1]][j-1,2]!=calling[[1]][j,2]){
            png(filename=paste(directory2,"/",calling[[1]][j,1],";",
                           calling[[1]][j,2],".png",sep=""),
                width=1600,height=600)
            counter<-0
        }
    }
    par(mar=c(16,6,2,1)) 
    if(counter==0){
        barplot(t(plot_abs[[j]]),beside=TRUE,col=t(plot_col[[j]]),
                xlim=c(-1.2*length(plot_abs[[j]][,1]),7*length(plot_abs[[j]][,1])),
                ylim=c((-1*0.1*abs_max),abs_max),ylab="Number of reads",
                main=paste(calling[[1]][j,1],";",calling[[1]][j,2],sep=""),
                names.arg=t(samples[1]),las=2,cex.axis=2,cex.names=2,cex.lab=2,
                cex.main=2)
    }
    if(counter>0){
        barplot(t(plot_abs[[j]]),beside=TRUE,col=t(plot_col[[j]]),
                xlim=c(-1.2*length(plot_abs[[j]][,1]),7*length(plot_abs[[j]][,1])),
                ylim=c((-1*0.1*abs_max),abs_max),ylab="Number of reads",
                main=paste(calling[[1]][j,1],";",calling[[1]][j,2],"_",counter,sep=""),
                names.arg=t(samples[1]),las=2,cex.axis=2,cex.names=2,
                cex.lab=2,cex.main=2)
    }
    par(new=TRUE)
    barplot((t(plot_ref[[j]]))*0.1*abs_max,beside=TRUE,
            col=c("forestgreen","blue4","darkgoldenrod3","darkred","black","purple"),
            xlim=c(-1.2*length(plot_abs[[j]][,1]),7*length(plot_abs[[j]][,1])),
            ylim=c((-1*0.1*abs_max),abs_max),yaxt='n')
    for(i in 1:length(plot_abs[[j]][,1])){
        for(k in 4:1){
            rect(xleft=(7*i-1),xright=(7*i),ybottom=0,
                 ytop=plot_insert[[j]][i,k],col=plot_col_insert[[j]][i,k])
        }
        if(plot_insert[[j]][i,4]!=0)
            rect(xleft=(7*i-1),xright=(7*i),ybottom=0,
                 ytop=plot_insert[[j]][i,4],col=NA,border="purple",lwd=4)
    }
    counter2<-0
    for(i in 1:length(plot_abs[[j]][,1])){
        if(length(marks)>0){
            for(k in 1:length(marks)){
                lines(x=c((1+counter2*7),(7+counter2*7)),
                      y=c((sum(plot_abs[[j]][i,1:5],na.rm=TRUE)*marks[k]),
                          (sum(plot_abs[[j]][i,1:5],na.rm=TRUE)*marks[k])))
            }
            counter2<-counter2+1
        }
    }
    text(x=-0.5*length(plot_abs[[j]][,1]),y=c(abs_max*0.8*(-0.1)),"Reference",cex=2)
    steps<-seq(1,length(plot_abs[[j]][,1])*0.65,
               by=(length(plot_abs[[j]][,1])*0.65-1)/length(A_col))-1.2*length(plot_abs[[j]][,1])
    if(length(plot_abs[[j]][,1])==1){
        steps<-seq(1*0.65,1,by=(1-1*0.65)/length(A_col))-1.75*1
    }
    if(length(plot_abs[[j]][,1])==2){
        steps<-seq(0.7,2*0.65,by=(2*0.65-0.7)/length(A_col))-1.2*2
    }   
    rect(xleft=steps[1:length(steps)-1],xright=steps[2:length(steps)],
         ybottom=abs_max*0.88,ytop=abs_max*0.93,lwd=0,col=A_col,border=NA)
    rect(xleft=steps[1:length(steps)-1],xright=steps[2:length(steps)],
         ybottom=abs_max*0.81,ytop=abs_max*0.86,lwd=0,col=C_col,border=NA)
    rect(xleft=steps[1:length(steps)-1],xright=steps[2:length(steps)],
         ybottom=abs_max*0.74,ytop=abs_max*0.79,lwd=0,col=G_col,border=NA)
    rect(xleft=steps[1:length(steps)-1],xright=steps[2:length(steps)],
         ybottom=abs_max*0.67,ytop=abs_max*0.72,lwd=0,col=T_col,border=NA)
    text(x=length(plot_abs[[j]][,1])*0.75-1.2*length(plot_abs[[j]][,1]),
         y=c(abs_max*0.90),"A",cex=2)
    text(x=length(plot_abs[[j]][,1])*0.75-1.2*length(plot_abs[[j]][,1]),
         y=c(abs_max*0.83),"C",cex=2)
    text(x=length(plot_abs[[j]][,1])*0.75-1.2*length(plot_abs[[j]][,1]),
         y=c(abs_max*0.76),"G",cex=2)
    text(x=length(plot_abs[[j]][,1])*0.75-1.2*length(plot_abs[[j]][,1]),
         y=c(abs_max*0.69),"T",cex=2)
    text(x=length(plot_abs[[j]][,1])*0.78-1.2*length(plot_abs[[j]][,1]),
         y=c(abs_max*0.62),"Del",cex=2)
    text(x=length(plot_abs[[j]][,1])*0.78-1.2*length(plot_abs[[j]][,1]),
         y=c(abs_max*0.55),"Ins",cex=2)
    text(x=length(plot_abs[[j]][,1])*1.15-1.2*length(plot_abs[[j]][,1]),
         y=c(abs_max*0.97),"Mean quality",cex=2)
    if(length(plot_abs[[j]][,1])>2){
        rect(xleft=1-1.2*length(plot_abs[[j]][,1]),
             xright=length(plot_abs[[j]][,1])*0.65-1.2*length(plot_abs[[j]][,1]),
             ybottom=abs_max*0.60,ytop=abs_max*0.65,lwd=0,col="black")
        rect(xleft=1-1.2*length(plot_abs[[j]][,1]),
             xright=length(plot_abs[[j]][,1])*0.65-1.2*length(plot_abs[[j]][,1]),
             ybottom=abs_max*0.53,ytop=abs_max*0.58,lwd=2,
             border="purple",col="white")
        text(x=1-1.2*length(plot_abs[[j]][,1]),y=c(abs_max*0.97),
             qual_lower_bound,cex=2)
        text(x=length(plot_abs[[j]][,1])*0.65-1.2*length(plot_abs[[j]][,1]),
             y=c(abs_max*0.97),qual_upper_bound,cex=2)
    }
    if(length(plot_abs[[j]][,1])==1){
        rect(xleft=-1.1,xright=-0.75,ybottom=abs_max*0.60,ytop=abs_max*0.65,
             lwd=0,col="black")
        rect(xleft=-1.1,xright=-0.75,ybottom=abs_max*0.53,ytop=abs_max*0.58,
             lwd=2,border="purple",col="white")
        text(x=-1.1,y=c(abs_max*0.97),qual_lower_bound,cex=2)
        text(x=-0.75,y=c(abs_max*0.97),qual_upper_bound,cex=2)
    }
    if(length(plot_abs[[j]][,1])==2){
        rect(xleft=-1.7,xright=-1.1,ybottom=abs_max*0.60,ytop=abs_max*0.65,
             lwd=0,col="black")
        rect(xleft=-1.7,xright=-1.1,ybottom=abs_max*0.53,ytop=abs_max*0.58,
             lwd=2,border="purple",col="white")
        text(x=-1.7,y=c(abs_max*0.97),qual_lower_bound,cex=2)
        text(x=-1.1,y=c(abs_max*0.97),qual_upper_bound,cex=2)
    }
    if(directory2!=""){
        dev.off()
    }
    return(counter)
}


plotCountsPerTargetRelative_col<-function(samples,frequency_weight,long_insert,
                                          qual_lower_bound,qual_upper_bound,
                                          A_col,C_col,G_col,T_col,i,j,k,plot_col){
    if(k==4&&!is.na(frequency_weight[[i]][j-long_insert[i],1+k])){
        if(round(frequency_weight[[i]][j-long_insert[i],1+k])<qual_lower_bound){
            plot_col[i,1]<-"chartreuse"
        }
        if(round(frequency_weight[[i]][j-long_insert[i],1+k])>qual_upper_bound){
            plot_col[i,1]<-"forestgreen"
        }
        if(qual_lower_bound<=round(frequency_weight[[i]][j-long_insert[i],1+k])&&round(frequency_weight[[i]][j-long_insert[i],1+k])<=qual_upper_bound){
            plot_col[i,1]<-A_col[round(frequency_weight[[i]][j-long_insert[i],1+k])-qual_lower_bound+1]
        }
    }
    if(k==5&&!is.na(frequency_weight[[i]][j-long_insert[i],2+k])){
        if(round(frequency_weight[[i]][j-long_insert[i],2+k])<qual_lower_bound){
            plot_col[i,2]<-"cyan"
        }
        if(round(frequency_weight[[i]][j-long_insert[i],2+k])>qual_upper_bound){
            plot_col[i,2]<-"blue4"
        }
        if(qual_lower_bound<=round(frequency_weight[[i]][j-long_insert[i],2+k])&&round(frequency_weight[[i]][j-long_insert[i],2+k])<=qual_upper_bound){
            plot_col[i,2]<-C_col[round(frequency_weight[[i]][j-long_insert[i],2+k])-qual_lower_bound+1]
        }
    }
    if(k==6&&!is.na(frequency_weight[[i]][j-long_insert[i],3+k])){
        if(round(frequency_weight[[i]][j-long_insert[i],3+k])<qual_lower_bound){
            plot_col[i,3]<-"yellow"
        }
        if(round(frequency_weight[[i]][j-long_insert[i],3+k])>qual_upper_bound){
            plot_col[i,3]<-"darkgoldenrod3"
        }
        if(qual_lower_bound<=round(frequency_weight[[i]][j-long_insert[i],3+k])&&round(frequency_weight[[i]][j-long_insert[i],3+k])<=qual_upper_bound){
            plot_col[i,3]<-G_col[round(frequency_weight[[i]][j-long_insert[i],3+k])-qual_lower_bound+1]
        }
    }
    if(k==7&&!is.na(frequency_weight[[i]][j-long_insert[i],4+k])){
        if(round(frequency_weight[[i]][j-long_insert[i],4+k])<qual_lower_bound){
            plot_col[i,4]<-"firebrick1"
        }
        if(round(frequency_weight[[i]][j-long_insert[i],4+k])>qual_upper_bound){
            plot_col[i,4]<-"darkred"
        }
        if(qual_lower_bound<=round(frequency_weight[[i]][j-long_insert[i],4+k])&&round(frequency_weight[[i]][j-long_insert[i],4+k])<=qual_upper_bound){
            plot_col[i,4]<-T_col[round(frequency_weight[[i]][j-long_insert[i],4+k])-qual_lower_bound+1]
        }
    }
    if(k==8&&!is.na(frequency_weight[[i]][j-long_insert[i],4+k])){
        plot_col[i,5]<-"black"
    }
    if(k==9&&!is.na(frequency_weight[[i]][j-long_insert[i],4+k])){
        plot_col[i,6]<-"purple"
    }
    return(plot_col)
}

plotCountsPerTargetRelative_insert<-function(samples,frequency_weight,
                                             long_insert,i,j,k,plot_insert){
    if(k==1&&!(is.na(frequency_weight[[i]][j-long_insert[i],4])
               &&is.na(frequency_weight[[i]][j-long_insert[i],6])
               &&is.na(frequency_weight[[i]][j-long_insert[i],8])
               &&is.na(frequency_weight[[i]][j-long_insert[i],10])
               &&is.na(frequency_weight[[i]][j-long_insert[i],12]))){
        plot_insert[i,k]<-frequency_weight[[i]][j-long_insert[i],13+k]/sum(frequency_weight[[i]][j-long_insert[i],c(4,6,8,10,12)],na.rm=TRUE)
        if(is.na(plot_insert[i,k])){
            plot_insert[i,k]<-0
        }
    }
    if(k==3&&!(is.na(frequency_weight[[i]][j-long_insert[i],4])
               &&is.na(frequency_weight[[i]][j-long_insert[i],6])
               &&is.na(frequency_weight[[i]][j-long_insert[i],8])
               &&is.na(frequency_weight[[i]][j-long_insert[i],10])
               &&is.na(frequency_weight[[i]][j-long_insert[i],12]))){
        plot_insert[i,2]<-frequency_weight[[i]][j-long_insert[i],13+k]/sum(frequency_weight[[i]][j-long_insert[i],c(4,6,8,10,12)],na.rm=TRUE)
        if(is.na(plot_insert[i,2])){
            plot_insert[i,2]<-0
        }
    }
    if(k==5&&!(is.na(frequency_weight[[i]][j-long_insert[i],4])
               &&is.na(frequency_weight[[i]][j-long_insert[i],6])
               &&is.na(frequency_weight[[i]][j-long_insert[i],8])
               &&is.na(frequency_weight[[i]][j-long_insert[i],10])
               &&is.na(frequency_weight[[i]][j-long_insert[i],12]))){
        plot_insert[i,3]<-frequency_weight[[i]][j-long_insert[i],13+k]/sum(frequency_weight[[i]][j-long_insert[i],c(4,6,8,10,12)],na.rm=TRUE)
        if(is.na(plot_insert[i,3])){
            plot_insert[i,3]<-0
        }
    }
    if(k==7&&!(is.na(frequency_weight[[i]][j-long_insert[i],4])
               &&is.na(frequency_weight[[i]][j-long_insert[i],6])
               &&is.na(frequency_weight[[i]][j-long_insert[i],8])
               &&is.na(frequency_weight[[i]][j-long_insert[i],10])
               &&is.na(frequency_weight[[i]][j-long_insert[i],12]))){
        plot_insert[i,4]<-frequency_weight[[i]][j-long_insert[i],13+k]/sum(frequency_weight[[i]][j-long_insert[i],c(4,6,8,10,12)],na.rm=TRUE)
        if(is.na(plot_insert[i,4])){
            plot_insert[i,4]<-0
        }
    }
    return(plot_insert)
}

plotCountsPerTargetRelative_col_insert<-function(samples,frequency_weight,
                                                 long_insert,qual_lower_bound,
                                                 qual_upper_bound,
                                                 A_col,C_col,G_col,T_col,i,j,k,
                                                 plot_col_insert){
    if(k==1&&!is.na(frequency_weight[[i]][j-long_insert[i],14+k])){
        if(round(frequency_weight[[i]][j-long_insert[i],14+k])<qual_lower_bound){
            plot_col_insert[i,1]<-"chartreuse"
        }
        if(round(frequency_weight[[i]][j-long_insert[i],14+k])>qual_upper_bound){
            plot_col_insert[i,1]<-"forestgreen"
        }
        if(qual_lower_bound<=round(frequency_weight[[i]][j-long_insert[i],14+k])
           &&round(frequency_weight[[i]][j-long_insert[i],14+k])<=qual_upper_bound){
            plot_col_insert[i,1]<-A_col[round(frequency_weight[[i]][j-long_insert[i],14+k])-qual_lower_bound+1]
        }
    }
    if(k==3&&!is.na(frequency_weight[[i]][j-long_insert[i],14+k])){
        if(round(frequency_weight[[i]][j-long_insert[i],14+k])<qual_lower_bound){
            plot_col_insert[i,1:2]<-"cyan"
        }
        if(round(frequency_weight[[i]][j-long_insert[i],14+k])>qual_upper_bound){
            plot_col_insert[i,1:2]<-"blue4"
        }
        if(qual_lower_bound<=round(frequency_weight[[i]][j-long_insert[i],14+k])
           &&round(frequency_weight[[i]][j-long_insert[i],14+k])<=qual_upper_bound){
            plot_col_insert[i,1:2]<-C_col[round(frequency_weight[[i]][j-long_insert[i],14+k])-qual_lower_bound+1]
        }
    }
    if(k==5&&!is.na(frequency_weight[[i]][j-long_insert[i],14+k])){
        if(round(frequency_weight[[i]][j-long_insert[i],14+k])<qual_lower_bound){
            plot_col_insert[i,1:3]<-"yellow"
        }
        if(round(frequency_weight[[i]][j-long_insert[i],14+k])>qual_upper_bound){
            plot_col_insert[i,1:3]<-"darkgoldenrod3"
        }
        if(qual_lower_bound<=round(frequency_weight[[i]][j-long_insert[i],14+k])
           &&round(frequency_weight[[i]][j-long_insert[i],14+k])<=qual_upper_bound){
            plot_col_insert[i,1:3]<-G_col[round(frequency_weight[[i]][j-long_insert[i],14+k])-qual_lower_bound+1]
        }
    }
    if(k==7&&!is.na(frequency_weight[[i]][j-long_insert[i],14+k])){
        if(round(frequency_weight[[i]][j-long_insert[i],14+k])<qual_lower_bound){
            plot_col_insert[i,1:4]<-"firebrick1"
        }
        if(round(frequency_weight[[i]][j-long_insert[i],14+k])>qual_upper_bound){
            plot_col_insert[i,1:4]<-"darkred"
        }
        if(qual_lower_bound<=round(frequency_weight[[i]][j-long_insert[i],14+k])
           &&round(frequency_weight[[i]][j-long_insert[i],14+k])<=qual_upper_bound){
            plot_col_insert[i,1:4]<-T_col[round(frequency_weight[[i]][j-long_insert[i],14+k])-qual_lower_bound+1]
        }
    }  
    return(plot_col_insert)
}

plotCountsPerTargetRelative_plot<-function(directory2,samples,plot_rel,plot_col,
                                           rel_max,calling,plot_insert,
                                           plot_col_insert,A_col,C_col,G_col,
                                           T_col,plot_ref,marks,
                                           qual_lower_bound,qual_upper_bound,
                                           counter,j){
    if(directory2!=""){
        if(j>1&&calling[[1]][j-1,1]==calling[[1]][j,1]
           &&calling[[1]][j-1,2]==calling[[1]][j,2]){
            counter<-counter+1
            png(filename=paste(directory2,"/",calling[[1]][j,1],";",calling[[1]][j,2],"_",counter,".png",sep=""),
                width=1600,height=600)
        }
        if(j==1||calling[[1]][j-1,1]!=calling[[1]][j,1]||calling[[1]][j-1,2]!=calling[[1]][j,2]){
            png(filename=paste(directory2,"/",calling[[1]][j,1],";",calling[[1]][j,2],".png",sep=""),
                width=1600,height=600)
            counter<-0
        }
    }
    par(mar=c(16,6,2,1)) 
    if(counter==0){
        barplot(t(plot_rel[[j]]),beside=TRUE,col=t(plot_col[[j]]),
                xlim=c(-1.2*length(plot_rel[[j]][,1]),7*length(plot_rel[[j]][,1])),
                ylim=c((-1*0.1*rel_max),rel_max),ylab="Relative frequency",
                main=paste(calling[[1]][j,1],";",calling[[1]][j,2],sep=""),
                names.arg=t(samples[1]),las=2,cex.axis=2,cex.names=2,cex.lab=2,
                cex.main=2)
    }
    if(counter>0){
        barplot(t(plot_rel[[j]]),beside=TRUE,col=t(plot_col[[j]]),
                xlim=c(-1.2*length(plot_rel[[j]][,1]),7*length(plot_rel[[j]][,1])),
                ylim=c((-1*0.1*rel_max),rel_max),ylab="Relative frequency",
                main=paste(calling[[1]][j,1],";",calling[[1]][j,2],"_",counter,sep=""),
                names.arg=t(samples[1]),las=2,cex.axis=2,cex.names=2,cex.lab=2,
                cex.main=2)
    }
    par(new=TRUE)
    barplot((t(plot_ref[[j]]))*0.1*rel_max,beside=TRUE,
            col=c("forestgreen","blue4","darkgoldenrod3","darkred","black","purple"),
            xlim=c(-1.2*length(plot_rel[[j]][,1]),7*length(plot_rel[[j]][,1])),
            ylim=c((-1*0.1*rel_max),rel_max),yaxt='n')
    for(i in 1:length(plot_rel[[j]][,1])){
        for(k in 4:1){
            rect(xleft=(7*i-1),xright=(7*i),ybottom=0,
                 ytop=plot_insert[[j]][i,k],col=plot_col_insert[[j]][i,k])
        }
        if(plot_insert[[j]][i,4]!=0)
            rect(xleft=(7*i-1),xright=(7*i),ybottom=0,
                 ytop=plot_insert[[j]][i,4],col=NA,border="purple",lwd=4)
    }
    counter2<-0
    for(i in 1:length(plot_rel[[j]][,1])){
        if(length(marks)>0){
            for(k in 1:length(marks)){
                lines(x=c((1+counter2*7),(7+counter2*7)),y=c(marks[k],marks[k]))
            }
            counter2<-counter2+1
        }
    }
    text(x=-0.5*length(plot_rel[[j]][,1]),y=c(rel_max*0.8*(-0.1)),"Reference",cex=2)
    steps<-seq(1,length(plot_rel[[j]][,1])*0.65,
               by=(length(plot_rel[[j]][,1])*0.65-1)/length(A_col))-1.2*length(plot_rel[[j]][,1])
    if(length(plot_rel[[j]][,1])==1){
        steps<-seq(1*0.65,1,by=(1-1*0.65)/length(A_col))-1.75*1
    }
    if(length(plot_rel[[j]][,1])==2){
        steps<-seq(0.7,2*0.65,by=(2*0.65-0.7)/length(A_col))-1.2*2
    }    
    rect(xleft=steps[1:length(steps)-1],xright=steps[2:length(steps)],
         ybottom=rel_max*0.88,ytop=rel_max*0.93,lwd=0,col=A_col,border=NA)
    rect(xleft=steps[1:length(steps)-1],xright=steps[2:length(steps)],
         ybottom=rel_max*0.81,ytop=rel_max*0.86,lwd=0,col=C_col,border=NA)
    rect(xleft=steps[1:length(steps)-1],xright=steps[2:length(steps)],
         ybottom=rel_max*0.74,ytop=rel_max*0.79,lwd=0,col=G_col,border=NA)
    rect(xleft=steps[1:length(steps)-1],xright=steps[2:length(steps)],
         ybottom=rel_max*0.67,ytop=rel_max*0.72,lwd=0,col=T_col,border=NA)
    text(x=length(plot_rel[[j]][,1])*0.75-1.2*length(plot_rel[[j]][,1]),
         y=c(rel_max*0.90),"A",cex=2)
    text(x=length(plot_rel[[j]][,1])*0.75-1.2*length(plot_rel[[j]][,1]),
         y=c(rel_max*0.83),"C",cex=2)
    text(x=length(plot_rel[[j]][,1])*0.75-1.2*length(plot_rel[[j]][,1]),
         y=c(rel_max*0.76),"G",cex=2)
    text(x=length(plot_rel[[j]][,1])*0.75-1.2*length(plot_rel[[j]][,1]),
         y=c(rel_max*0.69),"T",cex=2)
    text(x=length(plot_rel[[j]][,1])*0.78-1.2*length(plot_rel[[j]][,1]),
         y=c(rel_max*0.62),"Del",cex=2)
    text(x=length(plot_rel[[j]][,1])*0.78-1.2*length(plot_rel[[j]][,1]),
         y=c(rel_max*0.55),"Ins",cex=2)
    text(x=length(plot_rel[[j]][,1])*1.15-1.2*length(plot_rel[[j]][,1]),
         y=c(rel_max*0.97),"Mean quality",cex=2)
    if(length(plot_rel[[j]][,1])>2){
        rect(xleft=1-1.2*length(plot_rel[[j]][,1]),
             xright=length(plot_rel[[j]][,1])*0.65-1.2*length(plot_rel[[j]][,1]),
             ybottom=rel_max*0.60,ytop=rel_max*0.65,lwd=0,col="black")
        rect(xleft=1-1.2*length(plot_rel[[j]][,1]),
             xright=length(plot_rel[[j]][,1])*0.65-1.2*length(plot_rel[[j]][,1]),
             ybottom=rel_max*0.53,ytop=rel_max*0.58,lwd=2,border="purple",
             col="white")
        text(x=1-1.2*length(plot_rel[[j]][,1]),y=c(rel_max*0.97),
             qual_lower_bound,cex=2)
        text(x=length(plot_rel[[j]][,1])*0.65-1.2*length(plot_rel[[j]][,1]),
             y=c(rel_max*0.97),qual_upper_bound,cex=2)
    }
    if(length(plot_rel[[j]][,1])==1){
        rect(xleft=-1.1,xright=-0.75,ybottom=rel_max*0.60,ytop=rel_max*0.65,
             lwd=0,col="black")
        rect(xleft=-1.1,xright=-0.75,ybottom=rel_max*0.53,ytop=rel_max*0.58,
             lwd=2,border="purple",col="white")
        text(x=-1.1,y=c(rel_max*0.97),qual_lower_bound,cex=2)
        text(x=-0.75,y=c(rel_max*0.97),qual_upper_bound,cex=2)
    }
    if(length(plot_rel[[j]][,1])==2){
        rect(xleft=-1.7,xright=-1.1,ybottom=rel_max*0.60,ytop=rel_max*0.65,
             lwd=0,col="black")
        rect(xleft=-1.7,xright=-1.1,ybottom=rel_max*0.53,ytop=rel_max*0.58,
             lwd=2,border="purple",col="white")
        text(x=-1.7,y=c(rel_max*0.97),qual_lower_bound,cex=2)
        text(x=-1.1,y=c(rel_max*0.97),qual_upper_bound,cex=2)
    }
    if(directory2!=""){
        dev.off()
    }
    return(counter)
}


plotCountsPerSampleVcf_plot<-function(directory2,samples,plot_abs,plot_col,
                                      abs_max,calling,plot_insert,
                                      plot_col_insert,A_col,C_col,G_col,T_col,
                                      plot_ref,marks,qual_lower_bound,
                                      qual_upper_bound,j){
    if(directory2!=""){
        png(filename=paste(directory2,"/",samples[j,1],".png",sep=""),width=1600,
            height=800)   
    }
    par(mar=c(16,6,2,1)) 
    barplot(t(plot_abs[[j]]),beside=TRUE,col=t(plot_col[[j]]),
            xlim=c(-1.2*length(plot_abs[[j]][,1]),7*length(plot_abs[[j]][,1])),
            ylim=c((-1*0.1*abs_max),abs_max),ylab="Number of reads",
            mgp=c(4.5,1,-2),main=samples[j,1],
            names.arg=paste(calling[[1]][,1],calling[[1]][,2],sep=";"),
            las=2,cex.axis=2,cex.names=2,cex.lab=2,cex.main=2)
    par(new=TRUE)
    barplot((t(plot_ref[[j]]))*0.1*abs_max,beside=TRUE,
            col=c("forestgreen","blue4","darkgoldenrod3","darkred","black","purple"),
            xlim=c(-1.2*length(plot_abs[[j]][,1]),7*length(plot_abs[[j]][,1])),
            ylim=c((-1*0.1*abs_max),abs_max),yaxt='n')
    for(i in 1:length(plot_abs[[j]][,1])){
        for(k in 4:1){
            rect(xleft=(7*i-1),xright=(7*i),ybottom=0,
                 ytop=plot_insert[[j]][i,k],col=plot_col_insert[[j]][i,k])
        }
        if(plot_insert[[j]][i,4]!=0)
            rect(xleft=(7*i-1),xright=(7*i),ybottom=0,
                 ytop=plot_insert[[j]][i,4],col=NA,border="purple",lwd=4)
    }
    counter<-0
    for(i in 1:length(plot_abs[[j]][,1])){
        if(length(marks)){
            for(k in 1:length(marks)){
                lines(x=c((1+counter*7),(7+counter*7)),
                      y=c((sum(plot_abs[[j]][i,1:5],na.rm=TRUE)*marks[k]),
                          (sum(plot_abs[[j]][i,1:5],na.rm=TRUE)*marks[k])))
            }
            counter<-counter+1 
        }
    }  
    for(i in 1:length(calling[[j]][,1])){
        if(is.na(calling[[j]][i,16])&&sum(plot_abs[[j]][i,1:5],na.rm=TRUE)!=0){
            if(calling[[j]][i,3]=="A"){
                rect(xleft=1+(i-1)*7,xright=2+(i-1)*7,ybottom=0,
                     ytop=sum(plot_abs[[j]][i,1:5],na.rm=TRUE),border="grey50",
                     lty=3,lwd=4)
            }
            if(calling[[j]][i,3]=="C"){
                rect(xleft=2+(i-1)*7,xright=3+(i-1)*7,ybottom=0,
                     ytop=sum(plot_abs[[j]][i,1:5],na.rm=TRUE),border="grey50",
                     lty=3,lwd=4)
            }
            if(calling[[j]][i,3]=="G"){
                rect(xleft=3+(i-1)*7,xright=4+(i-1)*7,ybottom=0,
                     ytop=sum(plot_abs[[j]][i,1:5],na.rm=TRUE),border="grey50",
                     lty=3,lwd=4)
            }
            if(calling[[j]][i,3]=="T"){
                rect(xleft=4+(i-1)*7,xright=5+(i-1)*7,ybottom=0,
                     ytop=sum(plot_abs[[j]][i,1:5],na.rm=TRUE),border="grey50",
                     lty=3,lwd=4)
            }
        }
        if(!is.na(calling[[j]][i,16])){
            if(calling[[j]][i,16]==calling[[j]][i,17]&&is.na(calling[[j]][i,18])){
                if(calling[[j]][i,16]=="A"){
                    rect(xleft=1+(i-1)*7,xright=2+(i-1)*7,ybottom=0,
                         ytop=sum(plot_abs[[j]][i,1:5],na.rm=TRUE),
                         border="grey50",lty=3,lwd=4)
                }
                if(calling[[j]][i,16]=="C"){
                    rect(xleft=2+(i-1)*7,xright=3+(i-1)*7,ybottom=0,
                         ytop=sum(plot_abs[[j]][i,1:5],na.rm=TRUE),
                         border="grey50",lty=3,lwd=4)
                }
                if(calling[[j]][i,16]=="G"){
                    rect(xleft=3+(i-1)*7,xright=4+(i-1)*7,ybottom=0,
                         ytop=sum(plot_abs[[j]][i,1:5],na.rm=TRUE),
                         border="grey50",lty=3,lwd=4)
                }
                if(calling[[j]][i,16]=="T"){
                    rect(xleft=4+(i-1)*7,xright=5+(i-1)*7,ybottom=0,
                         ytop=sum(plot_abs[[j]][i,1:5],na.rm=TRUE),
                         border="grey50",lty=3,lwd=4)
                } 
                if(calling[[j]][i,16]==""){
                    rect(xleft=5+(i-1)*7,xright=6+(i-1)*7,ybottom=0,
                         ytop=sum(plot_abs[[j]][i,1:5],na.rm=TRUE),
                         border="grey50",lty=3,lwd=4)
                } 
            }
            if(calling[[j]][i,16]!=calling[[j]][i,17]
               &&is.na(calling[[j]][i,18])){
                if(calling[[j]][i,16]=="A"){
                    rect(xleft=1+(i-1)*7,xright=2+(i-1)*7,ybottom=0,
                         ytop=(sum(plot_abs[[j]][i,1:5],na.rm=TRUE)/2),
                         border="grey50",lty=3,lwd=4)
                }
                if(calling[[j]][i,16]=="C"){
                    rect(xleft=2+(i-1)*7,xright=3+(i-1)*7,ybottom=0,
                         ytop=(sum(plot_abs[[j]][i,1:5],na.rm=TRUE)/2),
                         border="grey50",lty=3,lwd=4)
                }
                if(calling[[j]][i,16]=="G"){
                    rect(xleft=3+(i-1)*7,xright=4+(i-1)*7,ybottom=0,
                         ytop=(sum(plot_abs[[j]][i,1:5],na.rm=TRUE)/2),
                         border="grey50",lty=3,lwd=4)
                }
                if(calling[[j]][i,16]=="T"){
                    rect(xleft=4+(i-1)*7,xright=5+(i-1)*7,ybottom=0,
                         ytop=(sum(plot_abs[[j]][i,1:5],na.rm=TRUE)/2),
                         border="grey50",lty=3,lwd=4)
                }
                if(calling[[j]][i,16]==""){
                    rect(xleft=5+(i-1)*7,xright=6+(i-1)*7,ybottom=0,
                         ytop=(sum(plot_abs[[j]][i,1:5],na.rm=TRUE)/2),
                         border="grey50",lty=3,lwd=4)
                }
                if(calling[[j]][i,17]=="A"){
                    rect(xleft=1+(i-1)*7,xright=2+(i-1)*7,ybottom=0,
                         ytop=(sum(plot_abs[[j]][i,1:5],na.rm=TRUE)/2),
                         border="grey50",lty=3,lwd=4)
                }
                if(calling[[j]][i,17]=="C"){
                    rect(xleft=2+(i-1)*7,xright=3+(i-1)*7,ybottom=0,
                         ytop=(sum(plot_abs[[j]][i,1:5],na.rm=TRUE)/2),
                         border="grey50",lty=3,lwd=4)
                }
                if(calling[[j]][i,17]=="G"){
                    rect(xleft=3+(i-1)*7,xright=4+(i-1)*7,ybottom=0,
                         ytop=(sum(plot_abs[[j]][i,1:5],na.rm=TRUE)/2),
                         border="grey50",lty=3,lwd=4)
                }
                if(calling[[j]][i,17]=="T"){
                    rect(xleft=4+(i-1)*7,xright=5+(i-1)*7,ybottom=0,
                         ytop=(sum(plot_abs[[j]][i,1:5],na.rm=TRUE)/2),
                         border="grey50",lty=3,lwd=4)
                }
                if(calling[[j]][i,17]==""){
                    rect(xleft=5+(i-1)*7,xright=6+(i-1)*7,ybottom=0,
                         ytop=(sum(plot_abs[[j]][i,1:5],na.rm=TRUE)/2),
                         border="grey50",lty=3,lwd=4)
                }
            }
            if(is.na(calling[[j]][i,19])&&!is.na(calling[[j]][i,18])){
                rect(xleft=6+(i-1)*7,xright=7+(i-1)*7,ybottom=0,
                     ytop=sum(plot_abs[[j]][i,1:5],na.rm=TRUE),
                     border="grey50",lty=3,lwd=4)
                if(calling[[j]][i,16]=="A"){
                    rect(xleft=1+(i-1)*7,xright=2+(i-1)*7,ybottom=0,
                         ytop=sum(plot_abs[[j]][i,1:5],na.rm=TRUE),
                         border="grey50",lty=3,lwd=4)
                }
                if(calling[[j]][i,16]=="C"){
                    rect(xleft=2+(i-1)*7,xright=3+(i-1)*7,ybottom=0,
                         ytop=sum(plot_abs[[j]][i,1:5],na.rm=TRUE),
                         border="grey50",lty=3,lwd=4)
                }
                if(calling[[j]][i,16]=="G"){
                    rect(xleft=3+(i-1)*7,xright=4+(i-1)*7,ybottom=0,
                         ytop=sum(plot_abs[[j]][i,1:5],na.rm=TRUE),
                         border="grey50",lty=3,lwd=4)
                }
                if(calling[[j]][i,16]=="T"){
                    rect(xleft=4+(i-1)*7,xright=5+(i-1)*7,ybottom=0,
                         ytop=sum(plot_abs[[j]][i,1:5],na.rm=TRUE),
                         border="grey50",lty=3,lwd=4)
                } 
            }
            if(!is.na(calling[[j]][i,19])&&!is.na(calling[[j]][i,18])){
                rect(xleft=6+(i-1)*7,xright=7+(i-1)*7,ybottom=0,
                     ytop=(sum(plot_abs[[j]][i,1:5],na.rm=TRUE)/2),
                     border="grey50",lty=3,lwd=4)
                if(calling[[j]][i,16]==calling[[j]][i,17]){
                    if(calling[[j]][i,16]=="A"){
                        rect(xleft=1+(i-1)*7,xright=2+(i-1)*7,ybottom=0,
                             ytop=sum(plot_abs[[j]][i,1:5],na.rm=TRUE),
                             border="grey50",lty=3,lwd=4)
                    }
                    if(calling[[j]][i,16]=="C"){
                        rect(xleft=2+(i-1)*7,xright=3+(i-1)*7,ybottom=0,
                             ytop=sum(plot_abs[[j]][i,1:5],na.rm=TRUE),
                             border="grey50",lty=3,lwd=4)
                    }
                    if(calling[[j]][i,16]=="G"){
                        rect(xleft=3+(i-1)*7,xright=4+(i-1)*7,ybottom=0,
                             ytop=sum(plot_abs[[j]][i,1:5],na.rm=TRUE),
                             border="grey50",lty=3,lwd=4)
                    }
                    if(calling[[j]][i,16]=="T"){
                        rect(xleft=4+(i-1)*7,xright=5+(i-1)*7,ybottom=0,
                             ytop=sum(plot_abs[[j]][i,1:5],na.rm=TRUE),
                             border="grey50",lty=3,lwd=4)
                    } 
                }
                if(calling[[j]][i,16]!=calling[[j]][i,17]){
                    if(calling[[j]][i,16]=="A"){
                        rect(xleft=1+(i-1)*7,xright=2+(i-1)*7,ybottom=0,
                             ytop=(sum(plot_abs[[j]][i,1:5],na.rm=TRUE)/2),
                             border="grey50",lty=3,lwd=4)
                    }
                    if(calling[[j]][i,16]=="C"){
                        rect(xleft=2+(i-1)*7,xright=3+(i-1)*7,ybottom=0,
                             ytop=(sum(plot_abs[[j]][i,1:5],na.rm=TRUE)/2),
                             border="grey50",lty=3,lwd=4)
                    }
                    if(calling[[j]][i,16]=="G"){
                        rect(xleft=3+(i-1)*7,xright=4+(i-1)*7,ybottom=0,
                             ytop=(sum(plot_abs[[j]][i,1:5],na.rm=TRUE)/2),
                             border="grey50",lty=3,lwd=4)
                    }
                    if(calling[[j]][i,16]=="T"){
                        rect(xleft=4+(i-1)*7,xright=5+(i-1)*7,ybottom=0,
                             ytop=(sum(plot_abs[[j]][i,1:5],na.rm=TRUE)/2),
                             border="grey50",lty=3,lwd=4)
                    }
                    if(calling[[j]][i,16]==""){
                        rect(xleft=5+(i-1)*7,xright=6+(i-1)*7,ybottom=0,
                             ytop=(sum(plot_abs[[j]][i,1:5],na.rm=TRUE)/2),
                             border="grey50",lty=3,lwd=4)
                    }
                    if(calling[[j]][i,17]=="A"){
                        rect(xleft=1+(i-1)*7,xright=2+(i-1)*7,ybottom=0,
                             ytop=(sum(plot_abs[[j]][i,1:5],na.rm=TRUE)/2),
                             border="grey50",lty=3,lwd=4)
                    }
                    if(calling[[j]][i,17]=="C"){
                        rect(xleft=2+(i-1)*7,xright=3+(i-1)*7,ybottom=0,
                             ytop=(sum(plot_abs[[j]][i,1:5],na.rm=TRUE)/2),
                             border="grey50",lty=3,lwd=4)
                    }
                    if(calling[[j]][i,17]=="G"){
                        rect(xleft=3+(i-1)*7,xright=4+(i-1)*7,ybottom=0,
                             ytop=(sum(plot_abs[[j]][i,1:5],na.rm=TRUE)/2),
                             border="grey50",lty=3,lwd=4)
                    }
                    if(calling[[j]][i,17]=="T"){
                        rect(xleft=4+(i-1)*7,xright=5+(i-1)*7,ybottom=0,
                             ytop=(sum(plot_abs[[j]][i,1:5],na.rm=TRUE)/2),
                             border="grey50",lty=3,lwd=4)
                    }
                    if(calling[[j]][i,17]==""){
                        rect(xleft=5+(i-1)*7,xright=6+(i-1)*7,ybottom=0,
                             ytop=(sum(plot_abs[[j]][i,1:5],na.rm=TRUE)/2),
                             border="grey50",lty=3,lwd=4)
                    }
                }
            }
        }
    }
    
    text(x=-0.5*length(plot_abs[[j]][,1]),y=c(abs_max*0.8*(-0.1)),"Reference",cex=2)
    steps<-seq(1,length(plot_abs[[j]][,1])*0.65,
               by=(length(plot_abs[[j]][,1])*0.65-1)/length(A_col))-1.2*length(plot_abs[[j]][,1])
    if(length(plot_abs[[j]][,1])==1){
        steps<-seq(1*0.65,1,by=(1-1*0.65)/length(A_col))-1.75*1
    }
    if(length(plot_abs[[j]][,1])==2){
        steps<-seq(0.7,2*0.65,by=(2*0.65-0.7)/length(A_col))-1.2*2
    }
    rect(xleft=steps[1:length(steps)-1],xright=steps[2:length(steps)],
         ybottom=abs_max*0.88,ytop=abs_max*0.93,lwd=0,col=A_col,border=NA)
    rect(xleft=steps[1:length(steps)-1],xright=steps[2:length(steps)],
         ybottom=abs_max*0.81,ytop=abs_max*0.86,lwd=0,col=C_col,border=NA)
    rect(xleft=steps[1:length(steps)-1],xright=steps[2:length(steps)],
         ybottom=abs_max*0.74,ytop=abs_max*0.79,lwd=0,col=G_col,border=NA)
    rect(xleft=steps[1:length(steps)-1],xright=steps[2:length(steps)],
         ybottom=abs_max*0.67,ytop=abs_max*0.72,lwd=0,col=T_col,border=NA)
    text(x=length(plot_abs[[j]][,1])*0.75-1.2*length(plot_abs[[j]][,1]),
         y=c(abs_max*0.90),"A",cex=2)
    text(x=length(plot_abs[[j]][,1])*0.75-1.2*length(plot_abs[[j]][,1]),
         y=c(abs_max*0.83),"C",cex=2)
    text(x=length(plot_abs[[j]][,1])*0.75-1.2*length(plot_abs[[j]][,1]),
         y=c(abs_max*0.76),"G",cex=2)
    text(x=length(plot_abs[[j]][,1])*0.75-1.2*length(plot_abs[[j]][,1]),
         y=c(abs_max*0.69),"T",cex=2)
    text(x=length(plot_abs[[j]][,1])*0.78-1.2*length(plot_abs[[j]][,1]),
         y=c(abs_max*0.62),"Del",cex=2)
    text(x=length(plot_abs[[j]][,1])*0.78-1.2*length(plot_abs[[j]][,1]),
         y=c(abs_max*0.55),"Ins",cex=2)
    text(x=length(plot_abs[[j]][,1])*0.98-1.2*length(plot_abs[[j]][,1]),
         y=c(abs_max*0.48),"Expected",cex=2)
    text(x=length(plot_abs[[j]][,1])*1.15-1.2*length(plot_abs[[j]][,1]),
         y=c(abs_max*0.97),"Mean quality",cex=2)
    if(length(plot_abs[[j]][,1])>2){
        rect(xleft=1-1.2*length(plot_abs[[j]][,1]),
             xright=length(plot_abs[[j]][,1])*0.65-1.2*length(plot_abs[[j]][,1]),
             ybottom=abs_max*0.60,ytop=abs_max*0.65,lwd=0,col="black")
        rect(xleft=1-1.2*length(plot_abs[[j]][,1]),
             xright=length(plot_abs[[j]][,1])*0.65-1.2*length(plot_abs[[j]][,1]),
             ybottom=abs_max*0.53,ytop=abs_max*0.58,lwd=2,border="purple",
             col="white")
        rect(xleft=1-1.2*length(plot_abs[[j]][,1]),
             xright=length(plot_abs[[j]][,1])*0.65-1.2*length(plot_abs[[j]][,1]),
             ybottom=abs_max*0.46,ytop=abs_max*0.51,lwd=2,border="grey50",
             col="white",lty=3)
        text(x=1-1.2*length(plot_abs[[j]][,1]),y=c(abs_max*0.97),
             qual_lower_bound,cex=2)
        text(x=length(plot_abs[[j]][,1])*0.65-1.2*length(plot_abs[[j]][,1]),
             y=c(abs_max*0.97),qual_upper_bound,cex=2)
    }
    if(length(plot_abs[[j]][,1])==1){
        rect(xleft=-1.1,xright=-0.75,ybottom=abs_max*0.60,ytop=abs_max*0.65,
             lwd=0,col="black")
        rect(xleft=-1.1,xright=-0.75,ybottom=abs_max*0.53,ytop=abs_max*0.58,
             lwd=2,border="purple",col="white")
        rect(xleft=-1.1,xright=-0.75,ybottom=abs_max*0.46,ytop=abs_max*0.51,
             lwd=2,border="grey50",col="white",lty=3)
        text(x=-1.1,y=c(abs_max*0.97),qual_lower_bound,cex=2)
        text(x=-0.75,y=c(abs_max*0.97),qual_upper_bound,cex=2)
    }
    if(length(plot_abs[[j]][,1])==2){
        rect(xleft=-1.7,xright=-1.1,ybottom=abs_max*0.60,ytop=abs_max*0.65,
             lwd=0,col="black")
        rect(xleft=-1.7,xright=-1.1,ybottom=abs_max*0.53,ytop=abs_max*0.58,
             lwd=2,border="purple",col="white")
        rect(xleft=-1.7,xright=-1.1,ybottom=abs_max*0.46,ytop=abs_max*0.51,
             lwd=2,border="grey50",col="white",lty=3)
        text(x=-1.7,y=c(abs_max*0.97),qual_lower_bound,cex=2)
        text(x=-1.1,y=c(abs_max*0.97),qual_upper_bound,cex=2)
    }
    if(directory2!=""){
        dev.off()   
    }
    return()
}

plotCountsPerSampleRelativeVcf_plot<-function(directory2,samples,plot_rel,
                                              plot_col,rel_max,calling,
                                              plot_insert,plot_col_insert,
                                              A_col,C_col,G_col,T_col,plot_ref,
                                              marks,qual_lower_bound,
                                              qual_upper_bound,j){
    if(directory2!=""){
        png(filename=paste(directory2,"/",samples[j,1],".png",sep=""),width=1600,
            height=800)
    }
    par(mar=c(16,6,2,1)) 
    barplot(t(plot_rel[[j]]),beside=TRUE,col=t(plot_col[[j]]),
            xlim=c(-1.2*length(plot_rel[[j]][,1]),7*length(plot_rel[[j]][,1])),
            ylim=c((-1*0.1*rel_max),rel_max),ylab="Relative frequency",
            mgp=c(4.5,1,-2),main=samples[j,1],
            names.arg=paste(calling[[1]][,1],calling[[1]][,2],sep=";"),
            las=2,cex.axis=2,cex.names=2,cex.lab=2,cex.main=2)
    par(new=TRUE)
    barplot((t(plot_ref[[j]]))*0.1*rel_max,beside=TRUE,
            col=c("forestgreen","blue4","darkgoldenrod3","darkred","black","purple"),
            xlim=c(-1.2*length(plot_rel[[j]][,1]),7*length(plot_rel[[j]][,1])),
            ylim=c((-1*0.1*rel_max),rel_max),yaxt='n')
    for(i in 1:length(plot_rel[[j]][,1])){
        for(k in 4:1){
            rect(xleft=(7*i-1),xright=(7*i),ybottom=0,
                 ytop=plot_insert[[j]][i,k],col=plot_col_insert[[j]][i,k])
        }
        if(plot_insert[[j]][i,4]!=0)
            rect(xleft=(7*i-1),xright=(7*i),ybottom=0,
                 ytop=plot_insert[[j]][i,4],col=NA,border="purple",lwd=4)
    }
    counter<-0
    for(i in 1:length(plot_rel[[j]][,1])){
        if(length(marks)>0){
            for(k in 1:length(marks)){
                lines(x=c((1+counter*7),(7+counter*7)),y=c(marks[k],marks[k]))
            }
            counter<-counter+1 
        }
    }  
    for(i in 1:length(calling[[j]][,1])){
        if(is.na(calling[[j]][i,16])&&sum(plot_rel[[j]][i,1:5],na.rm=TRUE)!=0){
            if(calling[[j]][i,3]=="A"){
                rect(xleft=1+(i-1)*7,xright=2+(i-1)*7,ybottom=0,ytop=1,
                     border="grey50",lty=3,lwd=4)
            }
            if(calling[[j]][i,3]=="C"){
                rect(xleft=2+(i-1)*7,xright=3+(i-1)*7,ybottom=0,ytop=1,
                     border="grey50",lty=3,lwd=4)
            }
            if(calling[[j]][i,3]=="G"){
                rect(xleft=3+(i-1)*7,xright=4+(i-1)*7,ybottom=0,ytop=1,
                     border="grey50",lty=3,lwd=4)
            }
            if(calling[[j]][i,3]=="T"){
                rect(xleft=4+(i-1)*7,xright=5+(i-1)*7,ybottom=0,ytop=1,
                     border="grey50",lty=3,lwd=4)
            }
        }
        if(!is.na(calling[[j]][i,16])){
            if(calling[[j]][i,16]==calling[[j]][i,17]&&is.na(calling[[j]][i,18])){
                if(calling[[j]][i,16]=="A"){
                    rect(xleft=1+(i-1)*7,xright=2+(i-1)*7,ybottom=0,ytop=1,
                         border="grey50",lty=3,lwd=4)
                }
                if(calling[[j]][i,16]=="C"){
                    rect(xleft=2+(i-1)*7,xright=3+(i-1)*7,ybottom=0,ytop=1,
                         border="grey50",lty=3,lwd=4)
                }
                if(calling[[j]][i,16]=="G"){
                    rect(xleft=3+(i-1)*7,xright=4+(i-1)*7,ybottom=0,ytop=1,
                         border="grey50",lty=3,lwd=4)
                }
                if(calling[[j]][i,16]=="T"){
                    rect(xleft=4+(i-1)*7,xright=5+(i-1)*7,ybottom=0,ytop=1,
                         border="grey50",lty=3,lwd=4)
                } 
                if(calling[[j]][i,16]==""){
                    rect(xleft=5+(i-1)*7,xright=6+(i-1)*7,ybottom=0,ytop=1,
                         border="grey50",lty=3,lwd=4)
                } 
            }
            if(calling[[j]][i,16]!=calling[[j]][i,17]&&is.na(calling[[j]][i,18])){
                if(calling[[j]][i,16]=="A"){
                    rect(xleft=1+(i-1)*7,xright=2+(i-1)*7,ybottom=0,ytop=(1/2),
                         border="grey50",lty=3,lwd=4)
                }
                if(calling[[j]][i,16]=="C"){
                    rect(xleft=2+(i-1)*7,xright=3+(i-1)*7,ybottom=0,ytop=(1/2),
                         border="grey50",lty=3,lwd=4)
                }
                if(calling[[j]][i,16]=="G"){
                    rect(xleft=3+(i-1)*7,xright=4+(i-1)*7,ybottom=0,ytop=(1/2),
                         border="grey50",lty=3,lwd=4)
                }
                if(calling[[j]][i,16]=="T"){
                    rect(xleft=4+(i-1)*7,xright=5+(i-1)*7,ybottom=0,ytop=(1/2),
                         border="grey50",lty=3,lwd=4)
                }
                if(calling[[j]][i,16]==""){
                    rect(xleft=5+(i-1)*7,xright=6+(i-1)*7,ybottom=0,ytop=(1/2),
                         border="grey50",lty=3,lwd=4)
                }
                if(calling[[j]][i,17]=="A"){
                    rect(xleft=1+(i-1)*7,xright=2+(i-1)*7,ybottom=0,ytop=(1/2),
                         border="grey50",lty=3,lwd=4)
                }
                if(calling[[j]][i,17]=="C"){
                    rect(xleft=2+(i-1)*7,xright=3+(i-1)*7,ybottom=0,ytop=(1/2),
                         border="grey50",lty=3,lwd=4)
                }
                if(calling[[j]][i,17]=="G"){
                    rect(xleft=3+(i-1)*7,xright=4+(i-1)*7,ybottom=0,ytop=(1/2),
                         border="grey50",lty=3,lwd=4)
                }
                if(calling[[j]][i,17]=="T"){
                    rect(xleft=4+(i-1)*7,xright=5+(i-1)*7,ybottom=0,ytop=(1/2),
                         border="grey50",lty=3,lwd=4)
                }
                if(calling[[j]][i,17]==""){
                    rect(xleft=5+(i-1)*7,xright=6+(i-1)*7,ybottom=0,ytop=(1/2),
                         border="grey50",lty=3,lwd=4)
                }
            }
            if(is.na(calling[[j]][i,19])&&!is.na(calling[[j]][i,18])){
                rect(xleft=6+(i-1)*7,xright=7+(i-1)*7,ybottom=0,ytop=1,
                     border="grey50",lty=3,lwd=4)
                if(calling[[j]][i,16]=="A"){
                    rect(xleft=1+(i-1)*7,xright=2+(i-1)*7,ybottom=0,ytop=1,
                         border="grey50",lty=3,lwd=4)
                }
                if(calling[[j]][i,16]=="C"){
                    rect(xleft=2+(i-1)*7,xright=3+(i-1)*7,ybottom=0,ytop=1,
                         border="grey50",lty=3,lwd=4)
                }
                if(calling[[j]][i,16]=="G"){
                    rect(xleft=3+(i-1)*7,xright=4+(i-1)*7,ybottom=0,ytop=1,
                         border="grey50",lty=3,lwd=4)
                }
                if(calling[[j]][i,16]=="T"){
                    rect(xleft=4+(i-1)*7,xright=5+(i-1)*7,ybottom=0,ytop=1,
                         border="grey50",lty=3,lwd=4)
                } 
            }
            if(!is.na(calling[[j]][i,19])&&!is.na(calling[[j]][i,18])){
                rect(xleft=6+(i-1)*7,xright=7+(i-1)*7,ybottom=0,ytop=(1/2),
                     border="grey50",lty=3,lwd=4)
                if(calling[[j]][i,16]==calling[[j]][i,17]){
                    if(calling[[j]][i,16]=="A"){
                        rect(xleft=1+(i-1)*7,xright=2+(i-1)*7,ybottom=0,ytop=1,
                             border="grey50",lty=3,lwd=4)
                    }
                    if(calling[[j]][i,16]=="C"){
                        rect(xleft=2+(i-1)*7,xright=3+(i-1)*7,ybottom=0,ytop=1,
                             border="grey50",lty=3,lwd=4)
                    }
                    if(calling[[j]][i,16]=="G"){
                        rect(xleft=3+(i-1)*7,xright=4+(i-1)*7,ybottom=0,ytop=1,
                             border="grey50",lty=3,lwd=4)
                    }
                    if(calling[[j]][i,16]=="T"){
                        rect(xleft=4+(i-1)*7,xright=5+(i-1)*7,ybottom=0,ytop=1,
                             border="grey50",lty=3,lwd=4)
                    } 
                }
                if(calling[[j]][i,16]!=calling[[j]][i,17]){
                    if(calling[[j]][i,16]=="A"){
                        rect(xleft=1+(i-1)*7,xright=2+(i-1)*7,ybottom=0,
                             ytop=(1/2),border="grey50",lty=3,lwd=4)
                    }
                    if(calling[[j]][i,16]=="C"){
                        rect(xleft=2+(i-1)*7,xright=3+(i-1)*7,ybottom=0,
                             ytop=(1/2),border="grey50",lty=3,lwd=4)
                    }
                    if(calling[[j]][i,16]=="G"){
                        rect(xleft=3+(i-1)*7,xright=4+(i-1)*7,ybottom=0,
                             ytop=(1/2),border="grey50",lty=3,lwd=4)
                    }
                    if(calling[[j]][i,16]=="T"){
                        rect(xleft=4+(i-1)*7,xright=5+(i-1)*7,ybottom=0,
                             ytop=(1/2),border="grey50",lty=3,lwd=4)
                    }
                    if(calling[[j]][i,16]==""){
                        rect(xleft=5+(i-1)*7,xright=6+(i-1)*7,ybottom=0,
                             ytop=(1/2),border="grey50",lty=3,lwd=4)
                    }
                    if(calling[[j]][i,17]=="A"){
                        rect(xleft=1+(i-1)*7,xright=2+(i-1)*7,ybottom=0,
                             ytop=(1/2),border="grey50",lty=3,lwd=4)
                    }
                    if(calling[[j]][i,17]=="C"){
                        rect(xleft=2+(i-1)*7,xright=3+(i-1)*7,ybottom=0,
                             ytop=(1/2),border="grey50",lty=3,lwd=4)
                    }
                    if(calling[[j]][i,17]=="G"){
                        rect(xleft=3+(i-1)*7,xright=4+(i-1)*7,ybottom=0,
                             ytop=(1/2),border="grey50",lty=3,lwd=4)
                    }
                    if(calling[[j]][i,17]=="T"){
                        rect(xleft=4+(i-1)*7,xright=5+(i-1)*7,ybottom=0,
                             ytop=(1/2),border="grey50",lty=3,lwd=4)
                    }
                    if(calling[[j]][i,17]==""){
                        rect(xleft=5+(i-1)*7,xright=6+(i-1)*7,ybottom=0,
                             ytop=(1/2),border="grey50",lty=3,lwd=4)
                    }
                }
            }
        }
    }
    
    text(x=-0.5*length(plot_rel[[j]][,1]),y=c(rel_max*0.8*(-0.1)),"Reference",cex=2)
    steps<-seq(1,length(plot_rel[[j]][,1])*0.65,
               by=(length(plot_rel[[j]][,1])*0.65-1)/length(A_col))-1.2*length(plot_rel[[j]][,1])
    if(length(plot_rel[[j]][,1])==1){
        steps<-seq(1*0.65,1,by=(1-1*0.65)/length(A_col))-1.75*1
    }
    if(length(plot_rel[[j]][,1])==2){
        steps<-seq(0.7,2*0.65,by=(2*0.65-0.7)/length(A_col))-1.2*2
    }
    rect(xleft=steps[1:length(steps)-1],xright=steps[2:length(steps)],
         ybottom=rel_max*0.88,ytop=rel_max*0.93,lwd=0,col=A_col,border=NA)
    rect(xleft=steps[1:length(steps)-1],xright=steps[2:length(steps)],
         ybottom=rel_max*0.81,ytop=rel_max*0.86,lwd=0,col=C_col,border=NA)
    rect(xleft=steps[1:length(steps)-1],xright=steps[2:length(steps)],
         ybottom=rel_max*0.74,ytop=rel_max*0.79,lwd=0,col=G_col,border=NA)
    rect(xleft=steps[1:length(steps)-1],xright=steps[2:length(steps)],
         ybottom=rel_max*0.67,ytop=rel_max*0.72,lwd=0,col=T_col,border=NA)
    text(x=length(plot_rel[[j]][,1])*0.75-1.2*length(plot_rel[[j]][,1]),
         y=c(rel_max*0.90),"A",cex=2)
    text(x=length(plot_rel[[j]][,1])*0.75-1.2*length(plot_rel[[j]][,1]),
         y=c(rel_max*0.83),"C",cex=2)
    text(x=length(plot_rel[[j]][,1])*0.75-1.2*length(plot_rel[[j]][,1]),
         y=c(rel_max*0.76),"G",cex=2)
    text(x=length(plot_rel[[j]][,1])*0.75-1.2*length(plot_rel[[j]][,1]),
         y=c(rel_max*0.69),"T",cex=2)
    text(x=length(plot_rel[[j]][,1])*0.78-1.2*length(plot_rel[[j]][,1]),
         y=c(rel_max*0.62),"Del",cex=2)
    text(x=length(plot_rel[[j]][,1])*0.78-1.2*length(plot_rel[[j]][,1]),
         y=c(rel_max*0.55),"Ins",cex=2)
    text(x=length(plot_rel[[j]][,1])*0.98-1.2*length(plot_rel[[j]][,1]),
         y=c(rel_max*0.48),"Expected",cex=2)
    text(x=length(plot_rel[[j]][,1])*1.15-1.2*length(plot_rel[[j]][,1]),
         y=c(rel_max*0.97),"Mean quality",cex=2)
    if(length(plot_rel[[j]][,1])>2){
        rect(xleft=1-1.2*length(plot_rel[[j]][,1]),
             xright=length(plot_rel[[j]][,1])*0.65-1.2*length(plot_rel[[j]][,1]),
             ybottom=rel_max*0.60,ytop=rel_max*0.65,lwd=0,col="black")
        rect(xleft=1-1.2*length(plot_rel[[j]][,1]),
             xright=length(plot_rel[[j]][,1])*0.65-1.2*length(plot_rel[[j]][,1]),
             ybottom=rel_max*0.53,ytop=rel_max*0.58,lwd=2,border="purple",
             col="white")
        rect(xleft=1-1.2*length(plot_rel[[j]][,1]),
             xright=length(plot_rel[[j]][,1])*0.65-1.2*length(plot_rel[[j]][,1]),
             ybottom=rel_max*0.46,ytop=rel_max*0.51,lwd=2,border="grey50",
             col="white",lty=3)
        text(x=1-1.2*length(plot_rel[[j]][,1]),y=c(rel_max*0.97),
             qual_lower_bound,cex=2)
        text(x=length(plot_rel[[j]][,1])*0.65-1.2*length(plot_rel[[j]][,1]),
             y=c(rel_max*0.97),qual_upper_bound,cex=2)
    }
    if(length(plot_rel[[j]][,1])==1){
        rect(xleft=-1.1,xright=-0.75,ybottom=rel_max*0.60,ytop=rel_max*0.65,
             lwd=0,col="black")
        rect(xleft=-1.1,xright=-0.75,ybottom=rel_max*0.53,ytop=rel_max*0.58,
             lwd=2,border="purple",col="white")
        rect(xleft=-1.1,xright=-0.75,ybottom=rel_max*0.46,ytop=rel_max*0.51,
             lwd=2,border="grey50",col="white",lty=3)
        text(x=-1.1,y=c(rel_max*0.97),qual_lower_bound,cex=2)
        text(x=-0.75,y=c(rel_max*0.97),qual_upper_bound,cex=2)
    }
    if(length(plot_rel[[j]][,1])==2){
        rect(xleft=-1.7,xright=-1.1,ybottom=rel_max*0.60,ytop=rel_max*0.65,
             lwd=0,col="black")
        rect(xleft=-1.7,xright=-1.1,ybottom=rel_max*0.53,ytop=rel_max*0.58,
             lwd=2,border="purple",col="white")
        rect(xleft=-1.7,xright=-1.1,ybottom=rel_max*0.46,ytop=rel_max*0.51,
             lwd=2,border="grey50",col="white",lty=3)
        text(x=-1.7,y=c(rel_max*0.97),qual_lower_bound,cex=2)
        text(x=-1.1,y=c(rel_max*0.97),qual_upper_bound,cex=2)
    }
    if(directory2!=""){
        dev.off()
    }
    return()
}

plotCountsPerTargetVcf_plot<-function(directory2,samples,plot_abs,plot_col,
                                      abs_max,calling,plot_insert,
                                      plot_col_insert,A_col,C_col,G_col,T_col,
                                      plot_ref,marks,qual_lower_bound,
                                      qual_upper_bound,counter,j){
    if(directory2!=""){
        if(j>1&&calling[[1]][j-1,1]==calling[[1]][j,1]
           &&calling[[1]][j-1,2]==calling[[1]][j,2]){
            counter<-counter+1
            png(filename=paste(directory2,"/",calling[[1]][j,1],";",calling[[1]][j,2],"_",counter,".png",sep=""),
                width=1600,height=600)
        }
        if(j==1||calling[[1]][j-1,1]!=calling[[1]][j,1]||calling[[1]][j-1,2]!=calling[[1]][j,2]){
            png(filename=paste(directory2,"/",calling[[1]][j,1],";",calling[[1]][j,2],".png",sep=""),
                width=1600,height=600)
            counter<-0
        }
    }
    par(mar=c(16,6,2,1)) 
    if(counter==0){
        barplot(t(plot_abs[[j]]),beside=TRUE,col=t(plot_col[[j]]),
                xlim=c(-1.2*length(plot_abs[[j]][,1]),7*length(plot_abs[[j]][,1])),
                ylim=c((-1*0.1*abs_max),abs_max),ylab="Number of reads",
                main=paste(calling[[1]][j,1],";",calling[[1]][j,2],sep=""),
                names.arg=t(samples[1]),las=2,cex.axis=2,cex.names=2,cex.lab=2,
                cex.main=2)
    }
    if(counter>0){
        barplot(t(plot_abs[[j]]),beside=TRUE,col=t(plot_col[[j]]),
                xlim=c(-1.2*length(plot_abs[[j]][,1]),7*length(plot_abs[[j]][,1])),
                ylim=c((-1*0.1*abs_max),abs_max),ylab="Number of reads",
                main=paste(calling[[1]][j,1],";",calling[[1]][j,2],"_",counter,sep=""),
                names.arg=t(samples[1]),las=2,cex.axis=2,cex.names=2,cex.lab=2,
                cex.main=2)
    }
    par(new=TRUE)
    barplot((t(plot_ref[[j]]))*0.1*abs_max,beside=TRUE,
            col=c("forestgreen","blue4","darkgoldenrod3","darkred","black","purple"),
            xlim=c(-1.2*length(plot_abs[[j]][,1]),7*length(plot_abs[[j]][,1])),
            ylim=c((-1*0.1*abs_max),abs_max),yaxt='n')
    for(i in 1:length(plot_abs[[j]][,1])){
        for(k in 4:1){
            rect(xleft=(7*i-1),xright=(7*i),ybottom=0,
                 ytop=plot_insert[[j]][i,k],col=plot_col_insert[[j]][i,k])
        }
        if(plot_insert[[j]][i,4]!=0)
            rect(xleft=(7*i-1),xright=(7*i),ybottom=0,
                 ytop=plot_insert[[j]][i,4],col=NA,border="purple",lwd=4)
    }
    counter2<-0
    for(i in 1:length(plot_abs[[j]][,1])){
        if(length(marks)>0){
            for(k in 1:length(marks)){
                lines(x=c((