Nothing
#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()
}
Any scripts or data that you put into this service are public.
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.