one-liner for insert size histogram from BAM

Hola,
ayer necesitábamos obtener rápidamente un histograma con el tamaño de los insertos de una librería de lecturas/reads paired-end. Lo logramos con este one-liner que requiere haber instalado R:

$ samtools view -q 30 -F 3916 mapped_reads.bam | cut -f 9 | \
   Rscript -e 'data=abs(scan(file="stdin")); pdf("hist.pdf"); hist(data,xlab="insert size (bp)")'

$ evince hist.pdf

Lo explico por pasos:

1) Previamente habíamos alineado/mapeado las lecturas contra una referencia y convertido el alineamiento a formato BAM. Podemos hacernos una idea qué mapeos contiene el fichero mapped_reads.bam con ayuda de samtools:

$ samtools flagstat mapped_reads.bam 

497656 + 0 in total (QC-passed reads + QC-failed reads)
0 + 0 duplicates
497409 + 0 mapped (99.95%:-nan%)
497656 + 0 paired in sequencing
248828 + 0 read1
248828 + 0 read2
492706 + 0 properly paired (99.01%:-nan%)
...

2) con -q 30 le pedimos a samtools que nos devuelve solamente lecturas con calidad de mapeo (MAPQ) >= 30

3) con -F 3916 le pedimos que ignore, según podemos averiguar aquí, los siguientes tipos de lecturas: "read unmapped, mate unmapped, first in pair, not primary alignment, read fails platform/vendor quality checks, read is PCR or optical duplicate, supplementary alignment"

NOTA1: Si el fichero BAM es muy grande se puede hacer lo mismo con una muestra al azar de las lecturas, por ejemplo el 10%, con samtools view -s 0.10.

NOTA2: Si la distribución de mapeos contiene algunos insertos anormalmente grandes el histograma por defecto puede quedar demasiado ancho. En ese caso puede ser buena idea probar algo como:

hist(data[data < quantile(data,0.99)])

Hasta luego,
Bruno