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((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
        }
    }
    for(i in 1:length(calling)){
        if(is.na(calling[[i]][j,16])&&sum(plot_abs[[j]][i,1:5],na.rm=TRUE)!=0){
            if(calling[[i]][j,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[[i]][j,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[[i]][j,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[[i]][j,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[[i]][j,16])){
            if(calling[[i]][j,16]==calling[[i]][j,17]&&is.na(calling[[i]][j,18])){
                if(calling[[i]][j,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[[i]][j,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[[i]][j,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[[i]][j,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[[i]][j,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[[i]][j,16]!=calling[[i]][j,17]&&is.na(calling[[i]][j,18])){
                if(calling[[i]][j,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[[i]][j,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[[i]][j,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[[i]][j,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[[i]][j,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[[i]][j,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[[i]][j,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[[i]][j,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[[i]][j,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[[i]][j,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[[i]][j,19])&&!is.na(calling[[i]][j,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[[i]][j,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[[i]][j,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[[i]][j,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[[i]][j,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[[i]][j,19])&&!is.na(calling[[i]][j,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[[i]][j,16]==calling[[i]][j,17]){
                    if(calling[[i]][j,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[[i]][j,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[[i]][j,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[[i]][j,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[[i]][j,16]!=calling[[i]][j,17]){
                    if(calling[[i]][j,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[[i]][j,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[[i]][j,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[[i]][j,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[[i]][j,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[[i]][j,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[[i]][j,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[[i]][j,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[[i]][j,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[[i]][j,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(counter)
}

plotCountsPerTargetRelativeVcf_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
        }
    }
    for(i in 1:length(calling)){
        if(is.na(calling[[i]][j,16])&&sum(plot_rel[[j]][i,1:5],na.rm=TRUE)!=0){
            if(calling[[i]][j,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[[i]][j,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[[i]][j,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[[i]][j,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[[i]][j,16])){
            if(calling[[i]][j,16]==calling[[i]][j,17]
               &&is.na(calling[[i]][j,18])){
                if(calling[[i]][j,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[[i]][j,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[[i]][j,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[[i]][j,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[[i]][j,16]==""){
                    rect(xleft=5+(i-1)*7,xright=6+(i-1)*7,ybottom=0,ytop=1,
                         border="grey50",lty=3,lwd=4)
                } 
            }
            if(calling[[i]][j,16]!=calling[[i]][j,17]&&is.na(calling[[i]][j,18])){
                if(calling[[i]][j,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[[i]][j,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[[i]][j,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[[i]][j,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[[i]][j,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[[i]][j,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[[i]][j,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[[i]][j,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[[i]][j,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[[i]][j,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[[i]][j,19])&&!is.na(calling[[i]][j,18])){
                rect(xleft=6+(i-1)*7,xright=7+(i-1)*7,ybottom=0,ytop=1,
                     border="grey50",lty=3,lwd=4)
                if(calling[[i]][j,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[[i]][j,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[[i]][j,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[[i]][j,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[[i]][j,19])&&!is.na(calling[[i]][j,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[[i]][j,16]==calling[[i]][j,17]){
                    if(calling[[i]][j,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[[i]][j,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[[i]][j,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[[i]][j,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[[i]][j,16]!=calling[[i]][j,17]){
                    if(calling[[i]][j,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[[i]][j,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[[i]][j,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[[i]][j,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[[i]][j,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[[i]][j,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[[i]][j,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[[i]][j,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[[i]][j,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[[i]][j,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(counter)
}



#plot absolute (number of reads), plot frequency threshold, 
#quality by color (intensity) (define upper and lower bound for color-coding)
plotCountsPerSample<-function(samples,
                              directory2,
                              calling,
                              dbsnp_file,
                              frequency_weight,
                              qual_lower_bound,
                              qual_upper_bound,
                              marks){
    plot_abs<-list()
    plot_ref<-list()
    plot_col<-list()
    plot_insert<-list()
    plot_col_insert<-list()
    A_Pal<-colorRampPalette(c('chartreuse','forestgreen'))
    A_col<-A_Pal(10*(qual_upper_bound-qual_lower_bound)+1)
    C_Pal<-colorRampPalette(c('cyan','blue4'))
    C_col<-C_Pal(10*(qual_upper_bound-qual_lower_bound)+1)
    G_Pal<-colorRampPalette(c('yellow','darkgoldenrod3'))
    G_col<-G_Pal(10*(qual_upper_bound-qual_lower_bound)+1)
    T_Pal<-colorRampPalette(c('firebrick1','darkred'))
    T_col<-T_Pal(10*(qual_upper_bound-qual_lower_bound)+1)
    if((qual_upper_bound-qual_lower_bound)==0){
        A_col<-"forestgreen"
        C_col<-"blue4"
        G_col<-"darkgoldenrod3"
        T_col<-"darkred"
    }   
    for(j in 1:length(samples[,1])){
        plot_abs[[j]]<-matrix(rep(0,6*length(calling[[j]][,1])),ncol=6)
        plot_ref[[j]]<-matrix(rep(0,6*length(calling[[j]][,1])),ncol=6)
        plot_col[[j]]<-matrix(rep(0,6*length(calling[[j]][,1])),ncol=6)
        plot_insert[[j]]<-matrix(rep(0,4*length(calling[[j]][,1])),ncol=4)
        plot_col_insert[[j]]<-matrix(rep(NA,4*length(calling[[j]][,1])),ncol=4)
    }
    abs_max<-0
    for(j in 1:length(samples[,1])){
        long_insert<-0
        for(i in 1:length(calling[[j]][,1])){
            if(!is.na(frequency_weight[[j]][i-long_insert,1])
               &&paste(calling[[j]][i,1],calling[[j]][i,2])
               !=paste(frequency_weight[[j]][i-long_insert,1],frequency_weight[[j]][i-long_insert,2])){
                long_insert<-long_insert+1
            }
            if(is.na(frequency_weight[[j]][i-long_insert,1])){
                long_insert<-long_insert+1
            }
            for(k in c(1,3,5,7,9,10)){
                plot_abs[[j]]<-plotCountsPerSample_abs(calling,
                                                       frequency_weight,
                                                       long_insert,
                                                       i,j,k,plot_abs[[j]])
                plot_col[[j]]<-plotCountsPerSample_col(calling,
                                                       frequency_weight,
                                                       long_insert,
                                                       qual_lower_bound,
                                                       qual_upper_bound,
                                                       A_col,C_col,G_col,T_col,
                                                       i,j,k,plot_col[[j]])
   
            }
            if(sum(plot_abs[[j]][i,1:5],na.rm=TRUE)>abs_max){
                abs_max<-sum(plot_abs[[j]][i,1:5],na.rm=TRUE)
            }
            
            for(k in c(1,3,5,7)){
                plot_insert[[j]]<-plotCountsPerSample_insert(calling,
                                                             frequency_weight,
                                                             long_insert,
                                                             i,j,k,plot_insert[[j]])
                plot_col_insert[[j]]<-plotCountsPerSample_col(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[[j]])
            }
            
            plot_insert[[j]][i,2]<-sum(plot_insert[[j]][i,1],plot_insert[[j]][i,2])
            plot_insert[[j]][i,3]<-sum(plot_insert[[j]][i,2],plot_insert[[j]][i,3])
            plot_insert[[j]][i,4]<-sum(plot_insert[[j]][i,3],plot_insert[[j]][i,4])
            
            chr<-substr(calling[[j]][i,1],4,nchar(calling[[j]][i,1]))
            pos<-calling[[j]][i,2]
            if(dbsnp_file!=""){
                vcfFile<-system.file(dbsnp_file, package="VariantAnnotation")
                param = ScanVcfParam(fixed=character(),info=character(), geno=character(), 
                                     samples=character(),which=GRanges(seqnames=chr,ranges=IRanges(pos,pos)))
                vcf.comp<-readVcf(dbsnp_file,"hg19",param=param)
            }
            plot_ref[[j]]<-plotCountsPerSample_ref(calling,
                                                   dbsnp_file,
                                                   vcf.comp,
                                                   j,i,plot_ref[[j]])  
        }
    }
    if(abs_max==0){
        abs_max<-1
    }
    for(j in 1:length(plot_abs)){
        plotCountsPerSample_plot(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)
    }
    return()
}


#plot relative (number of reads), plot frequency threshold, 
#quality by color (intensity) (define upper and lower bound for color-coding)
plotCountsPerSampleRelative<-function(samples,
                                      directory2,
                                      calling,
                                      dbsnp_file,
                                      frequency_weight,
                                      qual_lower_bound,
                                      qual_upper_bound,
                                      marks){
    plot_rel<-list()
    plot_ref<-list()
    plot_col<-list()
    plot_insert<-list()
    plot_col_insert<-list()
    A_Pal<-colorRampPalette(c('chartreuse','forestgreen'))
    A_col<-A_Pal(10*(qual_upper_bound-qual_lower_bound)+1)
    C_Pal<-colorRampPalette(c('cyan','blue4'))
    C_col<-C_Pal(10*(qual_upper_bound-qual_lower_bound)+1)
    G_Pal<-colorRampPalette(c('yellow','darkgoldenrod3'))
    G_col<-G_Pal(10*(qual_upper_bound-qual_lower_bound)+1)
    T_Pal<-colorRampPalette(c('firebrick1','darkred'))
    T_col<-T_Pal(10*(qual_upper_bound-qual_lower_bound)+1)
    if((qual_upper_bound-qual_lower_bound)==0){
        A_col<-"forestgreen"
        C_col<-"blue4"
        G_col<-"darkgoldenrod3"
        T_col<-"darkred"
    }
    
    for(j in 1:length(samples[,1])){
        plot_rel[[j]]<-matrix(rep(0,6*length(calling[[j]][,1])),ncol=6)
        plot_ref[[j]]<-matrix(rep(0,6*length(calling[[j]][,1])),ncol=6)
        plot_col[[j]]<-matrix(rep(0,6*length(calling[[j]][,1])),ncol=6)
        plot_insert[[j]]<-matrix(rep(0,4*length(calling[[j]][,1])),ncol=4)
        plot_col_insert[[j]]<-matrix(rep(NA,4*length(calling[[j]][,1])),ncol=4)
    }
    rel_max<-1
    for(j in 1:length(samples[,1])){
        long_insert<-0
        for(i in 1:length(calling[[j]][,1])){
            if(!is.na(frequency_weight[[j]][i-long_insert,1])
               &&paste(calling[[j]][i,1],calling[[j]][i,2])
               !=paste(frequency_weight[[j]][i-long_insert,1],frequency_weight[[j]][i-long_insert,2])){
                long_insert<-long_insert+1
            }
            if(is.na(frequency_weight[[j]][i-long_insert,1])){
                long_insert<-long_insert+1
            }
            for(k in 4:9){
                plot_rel[[j]][i,k-3]<-calling[[j]][i,k]
                if(is.na(plot_rel[[j]][i,k-3])){
                    plot_rel[[j]][i,k-3]<-0
                }
                plot_col[[j]]<-plotCountsPerSampleRelative_col(calling,
                                                               frequency_weight,
                                                               long_insert,
                                                               qual_lower_bound,
                                                               qual_upper_bound,
                                                               A_col,C_col,
                                                               G_col,T_col,
                                                               i,j,k,plot_col[[j]])   
            }
            
            for(k in c(1,3,5,7)){
                plot_insert[[j]]<-plotCountsPerSampleRelative_insert(calling,
                                                                     frequency_weight,
                                                                     long_insert,
                                                                     i,j,k,
                                                                     plot_insert[[j]])
                plot_col_insert[[j]]<-plotCountsPerSampleRelative_col(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[[j]])
            }
            plot_insert[[j]][i,2]<-sum(plot_insert[[j]][i,1],plot_insert[[j]][i,2])
            plot_insert[[j]][i,3]<-sum(plot_insert[[j]][i,2],plot_insert[[j]][i,3])
            plot_insert[[j]][i,4]<-sum(plot_insert[[j]][i,3],plot_insert[[j]][i,4])
            
            chr<-substr(calling[[j]][i,1],4,nchar(calling[[j]][i,1]))
            pos<-calling[[j]][i,2]
            if(dbsnp_file!=""){
                vcfFile<-system.file(dbsnp_file, package="VariantAnnotation")
                param = ScanVcfParam(fixed=character(),info=character(), geno=character(), 
                                     samples=character(),which=GRanges(seqnames=chr,ranges=IRanges(pos,pos)))
                vcf.comp<-readVcf(dbsnp_file,"hg19",param=param)
            }
            plot_ref[[j]]<-plotCountsPerSample_ref(calling,
                                                   dbsnp_file,
                                                   vcf.comp,
                                                   j,i,plot_ref[[j]])  
        }
    }
    
    for(j in 1:length(plot_rel)){
        plotCountsPerSampleRelative_plot(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)
    }
    return()
}

#Plot counts per Target
plotCountsPerTarget<-function(samples,
                              directory2,
                              calling,
                              dbsnp_file,
                              frequency_weight,
                              qual_lower_bound,
                              qual_upper_bound,
                              marks){
    plot_abs<-list()
    plot_ref<-list()
    plot_col<-list()
    plot_insert<-list()
    plot_col_insert<-list()
    A_Pal<-colorRampPalette(c('chartreuse','forestgreen'))
    A_col<-A_Pal(10*(qual_upper_bound-qual_lower_bound)+1)
    C_Pal<-colorRampPalette(c('cyan','blue4'))
    C_col<-C_Pal(10*(qual_upper_bound-qual_lower_bound)+1)
    G_Pal<-colorRampPalette(c('yellow','darkgoldenrod3'))
    G_col<-G_Pal(10*(qual_upper_bound-qual_lower_bound)+1)
    T_Pal<-colorRampPalette(c('firebrick1','darkred'))
    T_col<-T_Pal(10*(qual_upper_bound-qual_lower_bound)+1)  
    if((qual_upper_bound-qual_lower_bound)==0){
        A_col<-"forestgreen"
        C_col<-"blue4"
        G_col<-"darkgoldenrod3"
        T_col<-"darkred"
    }
    
    for(j in 1:length(calling[[1]][,1])){
        plot_abs[[j]]<-matrix(rep(0,6*length(samples[,1])),ncol=6)
        plot_ref[[j]]<-matrix(rep(0,6*length(samples[,1])),ncol=6)
        plot_col[[j]]<-matrix(rep(0,6*length(samples[,1])),ncol=6)
        plot_insert[[j]]<-matrix(rep(0,4*length(samples[,1])),ncol=4)
        plot_col_insert[[j]]<-matrix(rep(NA,4*length(samples[,1])),ncol=4)
    }
    abs_max<-0
    long_insert<-rep(0,length(samples[,1]))
    
    for(j in 1:length(calling[[1]][,1])){
        for(i in 1:length(samples[,1])){
            if(!is.na(frequency_weight[[i]][j-long_insert[i],1])
               &&paste(calling[[i]][j,1],calling[[i]][j,2])
               !=paste(frequency_weight[[i]][j-long_insert[i],1],frequency_weight[[i]][j-long_insert[i],2])){
                long_insert[i]<-long_insert[i]+1
            }
            if(is.na(frequency_weight[[i]][j-long_insert[i],1])){
                long_insert[i]<-long_insert[i]+1
            }
            for(k in c(1,3,5,7,9,10)){
                plot_abs[[j]]<-plotCountsPerTarget_abs(samples,
                                                       frequency_weight,
                                                       long_insert,
                                                       i,j,k,plot_abs[[j]])
                plot_col[[j]]<-plotCountsPerTarget_col(samples,
                                                       frequency_weight,
                                                       long_insert,
                                                       qual_lower_bound,
                                                       qual_upper_bound,
                                                       A_col,C_col,G_col,T_col,
                                                       i,j,k,plot_col[[j]])
            }
            if(sum(plot_abs[[j]][i,1:5],na.rm=TRUE)>abs_max){
                abs_max<-sum(plot_abs[[j]][i,1:5],na.rm=TRUE)
            }
            for(k in c(1,3,5,7)){
                plot_insert[[j]]<-plotCountsPerTarget_insert(samples,
                                                             frequency_weight,
                                                             long_insert,
                                                             i,j,k,plot_insert[[j]])
                plot_col_insert[[j]]<-plotCountsPerTarget_col_insert(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[[j]])
            }
            plot_insert[[j]][i,2]<-sum(plot_insert[[j]][i,1],plot_insert[[j]][i,2])
            plot_insert[[j]][i,3]<-sum(plot_insert[[j]][i,2],plot_insert[[j]][i,3])
            plot_insert[[j]][i,4]<-sum(plot_insert[[j]][i,3],plot_insert[[j]][i,4])
            
            chr<-substr(calling[[i]][j,1],4,nchar(calling[[i]][j,1]))
            pos<-calling[[i]][j,2]
            if(dbsnp_file!=""){
                vcfFile<-system.file(dbsnp_file, package="VariantAnnotation")
                param = ScanVcfParam(fixed=character(),info=character(), geno=character(), 
                                     samples=character(),which=GRanges(seqnames=chr,ranges=IRanges(pos,pos)))
                vcf.comp<-readVcf(dbsnp_file,"hg19",param=param)
            }
            plot_ref[[j]]<-plotCountsPerTarget_ref(samples,
                                                   calling,
                                                   dbsnp_file,
                                                   vcf.comp,
                                                   j,i,plot_ref[[j]])
         }
    }
    if(abs_max==0){
        abs_max<-1
    }
    counter2<-0
    for(j in 1:length(calling[[1]][,1])){
        counter<-counter2
        counter2<-plotCountsPerTarget_plot(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)
    }
    return()
}

#Plot counts per Target
plotCountsPerTargetRelative<-function(samples,
                                      directory2,
                                      calling,
                                      dbsnp_file,
                                      frequency_weight,
                                      qual_lower_bound,
                                      qual_upper_bound,
                                      marks){
    plot_rel<-list()
    plot_ref<-list()
    plot_col<-list()
    plot_insert<-list()
    plot_col_insert<-list()
    A_Pal<-colorRampPalette(c('chartreuse','forestgreen'))
    A_col<-A_Pal(10*(qual_upper_bound-qual_lower_bound)+1)
    C_Pal<-colorRampPalette(c('cyan','blue4'))
    C_col<-C_Pal(10*(qual_upper_bound-qual_lower_bound)+1)
    G_Pal<-colorRampPalette(c('yellow','darkgoldenrod3'))
    G_col<-G_Pal(10*(qual_upper_bound-qual_lower_bound)+1)
    T_Pal<-colorRampPalette(c('firebrick1','darkred'))
    T_col<-T_Pal(10*(qual_upper_bound-qual_lower_bound)+1)
    if((qual_upper_bound-qual_lower_bound)==0){
        A_col<-"forestgreen"
        C_col<-"blue4"
        G_col<-"darkgoldenrod3"
        T_col<-"darkred"
    }
    
    for(j in 1:length(calling[[1]][,1])){
        plot_rel[[j]]<-matrix(rep(0,6*length(samples[,1])),ncol=6)
        plot_ref[[j]]<-matrix(rep(0,6*length(samples[,1])),ncol=6)
        plot_col[[j]]<-matrix(rep(0,6*length(samples[,1])),ncol=6)
        plot_insert[[j]]<-matrix(rep(0,4*length(samples[,1])),ncol=4)
        plot_col_insert[[j]]<-matrix(rep(NA,4*length(samples[,1])),ncol=4)
    }
    rel_max<-1
    long_insert<-rep(0,length(samples[,1]))
    
    for(j in 1:length(calling[[1]][,1])){
        for(i in 1:length(samples[,1])){
            if(!is.na(frequency_weight[[i]][j-long_insert[i],1])
               &&paste(calling[[i]][j,1],calling[[i]][j,2])
               !=paste(frequency_weight[[i]][j-long_insert[i],1],
                       frequency_weight[[i]][j-long_insert[i],2])){
                long_insert[i]<-long_insert[i]+1
            }
            if(is.na(frequency_weight[[i]][j-long_insert[i],1])){
                long_insert[i]<-long_insert[i]+1
            }
            for(k in 4:9){
                plot_rel[[j]][i,k-3]<-calling[[i]][j,k]
                if(is.na(plot_rel[[j]][i,k-3])){
                    plot_rel[[j]][i,k-3]<-0
                }
                plot_col[[j]]<-plotCountsPerTargetRelative_col(samples,
                                                               frequency_weight,
                                                               long_insert,
                                                               qual_lower_bound,
                                                               qual_upper_bound,
                                                               A_col,C_col,
                                                               G_col,T_col,
                                                               i,j,k,plot_col[[j]])    
            }            
            for(k in c(1,3,5,7)){
                plot_insert[[j]]<-plotCountsPerTargetRelative_insert(samples,
                                                                     frequency_weight,
                                                                     long_insert,
                                                                     i,j,k,
                                                                     plot_insert[[j]])
                plot_col_insert[[j]]<-plotCountsPerTargetRelative_col_insert(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[[j]]) 
            }
            plot_insert[[j]][i,2]<-sum(plot_insert[[j]][i,1],plot_insert[[j]][i,2])
            plot_insert[[j]][i,3]<-sum(plot_insert[[j]][i,2],plot_insert[[j]][i,3])
            plot_insert[[j]][i,4]<-sum(plot_insert[[j]][i,3],plot_insert[[j]][i,4])
            
            chr<-substr(calling[[i]][j,1],4,nchar(calling[[i]][j,1]))
            pos<-calling[[i]][j,2]
            if(dbsnp_file!=""){
                vcfFile<-system.file(dbsnp_file, package="VariantAnnotation")
                param = ScanVcfParam(fixed=character(),info=character(), geno=character(), 
                                     samples=character(),which=GRanges(seqnames=chr,ranges=IRanges(pos,pos)))
                vcf.comp<-readVcf(dbsnp_file,"hg19",param=param)
            }
            plot_ref[[j]]<-plotCountsPerTarget_ref(samples,
                                                   calling,
                                                   dbsnp_file,
                                                   vcf.comp,
                                                   j,i,plot_ref[[j]])  
        }
    }
    counter2<-0
    for(j in 1:length(calling[[1]][,1])){
        counter<-counter2
        counter2<-plotCountsPerTargetRelative_plot(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)  
    }
    return()
}





#plot absolute (number of reads), plot frequency threshold, quality by color 
#(intensity) (define upper and lower bound for color-coding)
plotCountsPerSampleVcf<-function(samples,
                                 directory2,
                                 calling,
                                 dbsnp_file,
                                 frequency_weight,
                                 qual_lower_bound,
                                 qual_upper_bound,
                                 marks){
    plot_abs<-list()
    plot_ref<-list()
    plot_col<-list()
    plot_insert<-list()
    plot_col_insert<-list()
    A_Pal<-colorRampPalette(c('chartreuse','forestgreen'))
    A_col<-A_Pal(10*(qual_upper_bound-qual_lower_bound)+1)
    C_Pal<-colorRampPalette(c('cyan','blue4'))
    C_col<-C_Pal(10*(qual_upper_bound-qual_lower_bound)+1)
    G_Pal<-colorRampPalette(c('yellow','darkgoldenrod3'))
    G_col<-G_Pal(10*(qual_upper_bound-qual_lower_bound)+1)
    T_Pal<-colorRampPalette(c('firebrick1','darkred'))
    T_col<-T_Pal(10*(qual_upper_bound-qual_lower_bound)+1)
    if((qual_upper_bound-qual_lower_bound)==0){
        A_col<-"forestgreen"
        C_col<-"blue4"
        G_col<-"darkgoldenrod3"
        T_col<-"darkred"
    }
    
    for(j in 1:length(samples[,1])){
        plot_abs[[j]]<-matrix(rep(0,6*length(calling[[j]][,1])),ncol=6)
        plot_ref[[j]]<-matrix(rep(0,6*length(calling[[j]][,1])),ncol=6)
        plot_col[[j]]<-matrix(rep(0,6*length(calling[[j]][,1])),ncol=6)
        plot_insert[[j]]<-matrix(rep(0,4*length(calling[[j]][,1])),ncol=4)
        plot_col_insert[[j]]<-matrix(rep(NA,4*length(calling[[j]][,1])),ncol=4)
    }
    abs_max<-0
    for(j in 1:length(samples[,1])){
        long_insert<-0
        for(i in 1:length(calling[[j]][,1])){
            if(!is.na(frequency_weight[[j]][i-long_insert,1])
               &&paste(calling[[j]][i,1],calling[[j]][i,2])
               !=paste(frequency_weight[[j]][i-long_insert,1],frequency_weight[[j]][i-long_insert,2])){
                long_insert<-long_insert+1
            }
            if(is.na(frequency_weight[[j]][i-long_insert,1])){
                long_insert<-long_insert+1
            }
            for(k in c(1,3,5,7,9,10)){
                plot_abs[[j]]<-plotCountsPerSample_abs(calling,
                                                       frequency_weight,
                                                       long_insert,
                                                       i,j,k,plot_abs[[j]])
                plot_col[[j]]<-plotCountsPerSample_col(calling,
                                                       frequency_weight,
                                                       long_insert,
                                                       qual_lower_bound,
                                                       qual_upper_bound,
                                                       A_col,C_col,G_col,T_col,
                                                       i,j,k,plot_col[[j]])   
            }
            if(sum(plot_abs[[j]][i,1:5],na.rm=TRUE)>abs_max){
                abs_max<-sum(plot_abs[[j]][i,1:5],na.rm=TRUE)
            }
            
            for(k in c(1,3,5,7)){
                plot_insert[[j]]<-plotCountsPerSample_insert(calling,
                                                             frequency_weight,
                                                             long_insert,
                                                             i,j,k,plot_insert[[j]])
                plot_col_insert[[j]]<-plotCountsPerSample_col_insert(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[[j]]) 
            }
            plot_insert[[j]][i,2]<-sum(plot_insert[[j]][i,1],plot_insert[[j]][i,2])
            plot_insert[[j]][i,3]<-sum(plot_insert[[j]][i,2],plot_insert[[j]][i,3])
            plot_insert[[j]][i,4]<-sum(plot_insert[[j]][i,3],plot_insert[[j]][i,4])
            
            chr<-substr(calling[[j]][i,1],4,nchar(calling[[j]][i,1]))
            pos<-calling[[j]][i,2]
            if(dbsnp_file!=""){
                vcfFile<-system.file(dbsnp_file, package="VariantAnnotation")
                param = ScanVcfParam(fixed=character(),info=character(), geno=character(), 
                                     samples=character(),which=GRanges(seqnames=chr,ranges=IRanges(pos,pos)))
                vcf.comp<-readVcf(dbsnp_file,"hg19",param=param)
            }
            plot_ref<-plotCountsPerSample_ref(calling,
                                              dbsnp_file,
                                              vcf.comp,
                                              j,i,plot_ref[[j]])
            
        }
    }
    if(abs_max==0){
        abs_max<-1
    }
    for(j in 1:length(plot_abs)){
        plotCountsPerSampleVcf_plot(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)
    }
    return()
}

#plot relative (number of reads), plot frequency threshold, quality by color 
#(intensity) (define upper and lower bound for color-coding)
plotCountsPerSampleRelativeVcf<-function(samples,
                                         directory2,
                                         calling,
                                         dbsnp_file,
                                         frequency_weight,
                                         qual_lower_bound,
                                         qual_upper_bound,
                                         marks){
    plot_rel<-list()
    plot_ref<-list()
    plot_col<-list()
    plot_insert<-list()
    plot_col_insert<-list()
    A_Pal<-colorRampPalette(c('chartreuse','forestgreen'))
    A_col<-A_Pal(10*(qual_upper_bound-qual_lower_bound)+1)
    C_Pal<-colorRampPalette(c('cyan','blue4'))
    C_col<-C_Pal(10*(qual_upper_bound-qual_lower_bound)+1)
    G_Pal<-colorRampPalette(c('yellow','darkgoldenrod3'))
    G_col<-G_Pal(10*(qual_upper_bound-qual_lower_bound)+1)
    T_Pal<-colorRampPalette(c('firebrick1','darkred'))
    T_col<-T_Pal(10*(qual_upper_bound-qual_lower_bound)+1)
    if((qual_upper_bound-qual_lower_bound)==0){
        A_col<-"forestgreen"
        C_col<-"blue4"
        G_col<-"darkgoldenrod3"
        T_col<-"darkred"
    }
    
    for(j in 1:length(samples[,1])){
        plot_rel[[j]]<-matrix(rep(0,6*length(calling[[j]][,1])),ncol=6)
        plot_ref[[j]]<-matrix(rep(0,6*length(calling[[j]][,1])),ncol=6)
        plot_col[[j]]<-matrix(rep(0,6*length(calling[[j]][,1])),ncol=6)
        plot_insert[[j]]<-matrix(rep(0,4*length(calling[[j]][,1])),ncol=4)
        plot_col_insert[[j]]<-matrix(rep(NA,4*length(calling[[j]][,1])),ncol=4)
    }
    rel_max<-1
    for(j in 1:length(samples[,1])){
        long_insert<-0
        for(i in 1:length(calling[[j]][,1])){
            if(!is.na(frequency_weight[[j]][i-long_insert,1])
               &&paste(calling[[j]][i,1],calling[[j]][i,2])
               !=paste(frequency_weight[[j]][i-long_insert,1],frequency_weight[[j]][i-long_insert,2])){
                long_insert<-long_insert+1
            }
            if(is.na(frequency_weight[[j]][i-long_insert,1])){
                long_insert<-long_insert+1
            }
            for(k in 4:9){
                plot_rel[[j]][i,k-3]<-calling[[j]][i,k]
                if(is.na(plot_rel[[j]][i,k-3])){
                    plot_rel[[j]][i,k-3]<-0
                }
                
                plot_col[[j]]<-plotCountsPerSampleRelative_col(calling,
                                                          frequency_weight,
                                                          long_insert,
                                                          qual_lower_bound,
                                                          qual_upper_bound,
                                                          A_col,C_col,
                                                          G_col,T_col,
                                                          i,j,k,plot_col[[j]])
            }
            for(k in c(1,3,5,7)){
                plot_insert[[j]]<-plotCountsPerSampleRelative_insert(calling,
                                                                     frequency_weight,
                                                                     long_insert,
                                                                     i,j,k,
                                                                     plot_insert[[j]])
                plot_col_insert[[j]]<-plotCountsPerSampleRelative_col_insert(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[[j]])
            }
            plot_insert[[j]][i,2]<-sum(plot_insert[[j]][i,1],plot_insert[[j]][i,2])
            plot_insert[[j]][i,3]<-sum(plot_insert[[j]][i,2],plot_insert[[j]][i,3])
            plot_insert[[j]][i,4]<-sum(plot_insert[[j]][i,3],plot_insert[[j]][i,4])
            
            chr<-substr(calling[[j]][i,1],4,nchar(calling[[j]][i,1]))
            pos<-calling[[j]][i,2]
            if(dbsnp_file!=""){
                vcfFile<-system.file(dbsnp_file, package="VariantAnnotation")
                param = ScanVcfParam(fixed=character(),info=character(), geno=character(), 
                                     samples=character(),which=GRanges(seqnames=chr,ranges=IRanges(pos,pos)))
                vcf.comp<-readVcf(dbsnp_file,"hg19",param=param)
            }
            plot_ref[[j]]<-plotCountsPerSample_ref(calling,dbsnp_file,vcf.comp,j,i,
                                                   plot_ref[[j]])
        }
    }
    for(j in 1:length(plot_rel)){
        plotCountsPerSampleRelativeVcf_plot(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)
    }
    return()
}

#Plot counts per Target
plotCountsPerTargetVcf<-function(samples,
                                 directory2,
                                 calling,
                                 dbsnp_file,
                                 frequency_weight,
                                 qual_lower_bound,
                                 qual_upper_bound,
                                 marks){
    plot_abs<-list()
    plot_ref<-list()
    plot_col<-list()
    plot_insert<-list()
    plot_col_insert<-list()
    A_Pal<-colorRampPalette(c('chartreuse','forestgreen'))
    A_col<-A_Pal(10*(qual_upper_bound-qual_lower_bound)+1)
    C_Pal<-colorRampPalette(c('cyan','blue4'))
    C_col<-C_Pal(10*(qual_upper_bound-qual_lower_bound)+1)
    G_Pal<-colorRampPalette(c('yellow','darkgoldenrod3'))
    G_col<-G_Pal(10*(qual_upper_bound-qual_lower_bound)+1)
    T_Pal<-colorRampPalette(c('firebrick1','darkred'))
    T_col<-T_Pal(10*(qual_upper_bound-qual_lower_bound)+1)
    if((qual_upper_bound-qual_lower_bound)==0){
        A_col<-"forestgreen"
        C_col<-"blue4"
        G_col<-"darkgoldenrod3"
        T_col<-"darkred"
    }
    
    for(j in 1:length(calling[[1]][,1])){
        plot_abs[[j]]<-matrix(rep(0,6*length(samples[,1])),ncol=6)
        plot_ref[[j]]<-matrix(rep(0,6*length(samples[,1])),ncol=6)
        plot_col[[j]]<-matrix(rep(0,6*length(samples[,1])),ncol=6)
        plot_insert[[j]]<-matrix(rep(0,4*length(samples[,1])),ncol=4)
        plot_col_insert[[j]]<-matrix(rep(NA,4*length(samples[,1])),ncol=4)
    }
    abs_max<-0
    long_insert<-rep(0,length(samples[,1]))
    
    for(j in 1:length(calling[[1]][,1])){
        for(i in 1:length(samples[,1])){
            if(!is.na(frequency_weight[[i]][j-long_insert[i],1])
               &&paste(calling[[i]][j,1],calling[[i]][j,2])
               !=paste(frequency_weight[[i]][j-long_insert[i],1],frequency_weight[[i]][j-long_insert[i],2])){
                long_insert[i]<-long_insert[i]+1
            }
            if(is.na(frequency_weight[[i]][j-long_insert[i],1])){
                long_insert[i]<-long_insert[i]+1
            }
            for(k in c(1,3,5,7,9,10)){
                plot_abs[[j]]<-plotCountsPerTarget_abs(samples,
                                                       frequency_weight,
                                                       long_insert,
                                                       i,j,k,plot_abs[[j]])
                plot_col[[j]]<-plotCountsPerTarget_col(samples,
                                                       frequency_weight,
                                                       long_insert,
                                                       qual_lower_bound,
                                                       qual_upper_bound,
                                                       A_col,C_col,G_col,T_col,
                                                       i,j,k,plot_col[[j]])  
            }
            if(sum(plot_abs[[j]][i,1:5],na.rm=TRUE)>abs_max){
                abs_max<-sum(plot_abs[[j]][i,1:5],na.rm=TRUE)
            }
            for(k in c(1,3,5,7)){
                plot_insert[[j]]<-plotCountsPerTarget_insert(samples,
                                                             frequency_weight,
                                                             long_insert,
                                                             i,j,k,plot_insert[[j]])
                plot_col_insert[[j]]<-plotCountsPerTarget_col_insert(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[[j]])
            }
            plot_insert[[j]][i,2]<-sum(plot_insert[[j]][i,1],plot_insert[[j]][i,2])
            plot_insert[[j]][i,3]<-sum(plot_insert[[j]][i,2],plot_insert[[j]][i,3])
            plot_insert[[j]][i,4]<-sum(plot_insert[[j]][i,3],plot_insert[[j]][i,4])
            
            chr<-substr(calling[[i]][j,1],4,nchar(calling[[i]][j,1]))
            pos<-calling[[i]][j,2]
            if(dbsnp_file!=""){
                vcfFile<-system.file(dbsnp_file, package="VariantAnnotation")
                param = ScanVcfParam(fixed=character(),info=character(), geno=character(), 
                                     samples=character(),which=GRanges(seqnames=chr,ranges=IRanges(pos,pos)))
                vcf.comp<-readVcf(dbsnp_file,"hg19",param=param)
            }
            plot_ref[[j]]<-plotCountsPerTarget_ref(samples,calling,
                                                   dbsnp_file,vcf.comp,j,i,
                                                   plot_ref[[j]])
        }
    }
    if(abs_max==0){
        abs_max<-1
    }
    counter2<-0
    for(j in 1:length(calling[[1]][,1])){
        counter<-counter2
        counter2<-plotCountsPerTargetVcf_plot(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)
    }
    return()
}

#Plot counts per Target
plotCountsPerTargetRelativeVcf<-function(samples,
                                         directory2,
                                         calling,
                                         dbsnp_file,
                                         frequency_weight,
                                         qual_lower_bound,
                                         qual_upper_bound,
                                         marks){
    plot_rel<-list()
    plot_ref<-list()
    plot_col<-list()
    plot_insert<-list()
    plot_col_insert<-list()
    A_Pal<-colorRampPalette(c('chartreuse','forestgreen'))
    A_col<-A_Pal(10*(qual_upper_bound-qual_lower_bound)+1)
    C_Pal<-colorRampPalette(c('cyan','blue4'))
    C_col<-C_Pal(10*(qual_upper_bound-qual_lower_bound)+1)
    G_Pal<-colorRampPalette(c('yellow','darkgoldenrod3'))
    G_col<-G_Pal(10*(qual_upper_bound-qual_lower_bound)+1)
    T_Pal<-colorRampPalette(c('firebrick1','darkred'))
    T_col<-T_Pal(10*(qual_upper_bound-qual_lower_bound)+1)
    if((qual_upper_bound-qual_lower_bound)==0){
        A_col<-"forestgreen"
        C_col<-"blue4"
        G_col<-"darkgoldenrod3"
        T_col<-"darkred"
    }
    
    for(j in 1:length(calling[[1]][,1])){
        plot_rel[[j]]<-matrix(rep(0,6*length(samples[,1])),ncol=6)
        plot_ref[[j]]<-matrix(rep(0,6*length(samples[,1])),ncol=6)
        plot_col[[j]]<-matrix(rep(0,6*length(samples[,1])),ncol=6)
        plot_insert[[j]]<-matrix(rep(0,4*length(samples[,1])),ncol=4)
        plot_col_insert[[j]]<-matrix(rep(NA,4*length(samples[,1])),ncol=4)
    }
    rel_max<-1
    long_insert<-rep(0,length(samples[,1]))
    
    for(j in 1:length(calling[[1]][,1])){
        for(i in 1:length(samples[,1])){
            if(!is.na(frequency_weight[[i]][j-long_insert[i],1])
               &&paste(calling[[i]][j,1],calling[[i]][j,2])
               !=paste(frequency_weight[[i]][j-long_insert[i],1],frequency_weight[[i]][j-long_insert[i],2])){
                long_insert[i]<-long_insert[i]+1
            }
            if(is.na(frequency_weight[[i]][j-long_insert[i],1])){
                long_insert[i]<-long_insert[i]+1
            }
            for(k in 4:9){
                plot_rel[[j]][i,k-3]<-calling[[i]][j,k]
                if(is.na(plot_rel[[j]][i,k-3])){
                    plot_rel[[j]][i,k-3]<-0
                }
                plot_col[[j]]<-plotCountsPerTargetRelative_col(samples,
                                                               frequency_weight,
                                                               long_insert,
                                                               qual_lower_bound,
                                                               qual_upper_bound,
                                                               A_col,C_col,
                                                               G_col,T_col,
                                                               i,j,k,plot_col[[j]]) 
            }
            for(k in c(1,3,5,7)){
                plot_insert[[j]]<-plotCountsPerTargetRelative_insert(samples,
                                                                     frequency_weight,
                                                                     long_insert,
                                                                     i,j,k,
                                                                     plot_col[[j]])
                plot_col_insert[[j]]<-plotCountsPerTargetRelative_col_insert(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[[j]]) 
            }
            plot_insert[[j]][i,2]<-sum(plot_insert[[j]][i,1],plot_insert[[j]][i,2])
            plot_insert[[j]][i,3]<-sum(plot_insert[[j]][i,2],plot_insert[[j]][i,3])
            plot_insert[[j]][i,4]<-sum(plot_insert[[j]][i,3],plot_insert[[j]][i,4])
            
            chr<-substr(calling[[i]][j,1],4,nchar(calling[[i]][j,1]))
            pos<-calling[[i]][j,2]
            if(dbsnp_file!=""){
                vcfFile<-system.file(dbsnp_file, package="VariantAnnotation")
                param = ScanVcfParam(fixed=character(),info=character(), geno=character(), 
                                     samples=character(),which=GRanges(seqnames=chr,ranges=IRanges(pos,pos)))
                vcf.comp<-readVcf(dbsnp_file,"hg19",param=param)
            }
            plot_ref[[j]]<-plotCountsPerTarget_ref(samples,calling,
                                                   dbsnp_file,vcf.comp,j,i,
                                                   plot_ref[[j]]) 
        }
    }
    counter2<-0
    for(j in 1:length(calling[[1]][,1])){
        counter<-counter2
        counter2<-plotCountsPerTargetRelativeVcf_plot(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)
    }
    return()
}

Try the BBCAnalyzer package in your browser

Any scripts or data that you put into this service are public.

BBCAnalyzer documentation built on Nov. 8, 2020, 5:05 p.m.