Hi all! Given a BAM, I want to plot the insert size of a read pair vs average base quality of bases in a read I've actually got the insert size data point by filtering my BAM using:
samtools view -f66 mybam.bam | cut -f9 | awk '{print sqrt($0^2)}' > isizes.txt # To prevent negatives
But now I have doubts in how to fetch the mean of quality of all bases in both reads of a read pair. Could yo please give me a clue or something? Thank you very much!
Other tools/scripts are also welcomed!
•
written
by
pablosolar.r •
not fully tested, but it could be:
samtools collate -O -f in.bam TMPDIR |\
samtools view |\
awk -F '\t' 'BEGIN{N="";Q="";for(i=33;i<256;i++) ord[sprintf("%c",i)]=i-33;} {if($1==N) {Q=sprintf("%s%s",Q,$11);V=0.0; for(i=1;i<=length(Q);i++) {V+=ord[substr(Q,i,1)];} printf("%s %f\n",$1,V/length(Q)); N="";Q="";} else {N=$1;Q=$11;}}'
Hi Pierre, could you please explain that command? Thank you very much!! Is it collate really necessary?
Yes but in that case the expected results have not the same length. I mean, getting the insert sizes:
samtools view -f66 mybam.bam | cut -f9 | awk '{print sqrt($0^2)}' > isizes.txt # To prevent negatives
cat isizes.txt | wc -l
> 108853
And running the Pierre command I am getting:
samtools collate -O -f mybam.bam TMPDIR | samtools view |awk -F '\t' 'BEGIN{N="";Q="";for(i=33;i<256;i++) ord[sprintf("%c",i)]=i-33;} {if($1==N) {Q=sprintf("%s%s",Q,$11);V=0.0; for(i=1;i<=length(Q);i++) {V+=ord[substr(Q,i,1)];} printf("%s %f\n",$1,V/length(Q)); N="";Q="";} else {N=$1;Q=$11;}}' | awk '{split($0,a," "); print a[2]}' > average_bq.txt
cat average_bq.txt | wc -l
> 57604
So plot can't be performed due to the different length of data points. I am really missed with this... Thank you all!
Log in to answer this question.
See FASTQ Phred33 average base quality score for some clarification on
average base quality.Average base quality does not make a lot of sense in real world.