Insertos

A veces es conveniente saber cual era el tamaño de los insertos entre cada una de las secuencias pareadas (ver figura), pues eso ayuda a determinar si un ensamble pobre pudo venir desde la preparación de la librería de secuenciación. En general mientras mayor sea el tamaño promedio de los insertos, mejor podrá ser el ensamble. Insertos pueden variar entre 200 y 800 bases.

Usaremos bowtie2 para hacer este ejercicio que es una alineamiento de secuencia pareadas a contigs de un genoma ensamblado previamente.

Primero revisemos que tengamos bowtie2 instalado:

 $ bowtie2 --version

Debemos tener un resultado similar a este:

/usr/bin/bowtie2-align-s version 2.3.5.1

64-bit

Indice

El primer paso consiste en crear un índice con las secuencias fasta del ensamble previamente realizado (contigs.fasta en este ejemplo). Si tenemos un ensamble llamado contigs.fasta, podemos crear el índice, que se llamara "genome" así:

$ bowtie2-build contigs.fasta genome

Mapeo de secuencias

Ahora podemos mapear las secuencias pareadas (*.R1.fastq y *.R2.fastq) con ayuda de este índice. Es importante tener las secuencias ya limpias.

$ bowtie2 --no-mixed --no-discordant --no-contain --no-overlap --maxins 1200 -x genome -p 8 -S genome.sam --no-unal -1 file.R1.fastq -2 file.R2.fastq

genome.sam será el archivo de salida con los resultados del alineamiento de las secuencias pareadas al genoma.

Insertos

Obtengamos los tamaños de los insertos en una lista:

$awk '{ if ($9 != "0" && $9 !~ /-/ && $9 !~ /genome/ && $9 != "") print $9 }' genome.sam > inserts.list

Aquí estamos extrayendo con awk del archivo generado (genome.sam) la novena column ($9, tamaño de los insertos para cada secuencia pareada) pero solo aquellos valores que NO tengan cero (!= "0"), sean negativos (!~ /-/),  las lineas que estén vacías ($9 != "") o tengan la palabra "genome" (!~ /genome/) y los estamos salvando a un archivo llamado inserts.list

Resultados

De esta lista de valores de insertos (inserts.list) podemos calcular la media, mediana y moda e incluso generar un gráfico de distribución de frecuencias.

Media:

$ awk '{ sum += $1 } END { if (NR > 0) print sum / NR }' inserts.list

Desviación estándar:
$ awk '{sum += $1; sum2 += $1^2} END { if (NR > 0) print sqrt(sum2/NR - (sum/NR)^2) }' inserts.list

Moda

$ uniq -c inserts.list | sort -nr | head -n 1 | awk '{print $2}'

Gráfica

Iniciemos R: 

$ R

Grafiquemos con R esta lista de insertos (inserts.list). una vez entrado a R o RStudio, cargar ggplot2, este comando checa si existe el paquete ggplot2 instalado y si no, lo instala (puede tardar algunos minutos en hacerlo):

if (!require(ggplot2)) {

  install.packages("ggplot2")

}

Crear un objeto leyendo el archivo de insertos:

data <- read.table("inserts.list", header = FALSE)

Crear el gráfico con ggplot:

ggplot(data, aes(x = V1)) + geom_histogram(binwidth = 5, fill = "skyblue", color = "grey") + labs(title = "Insert size distribution", x = "Insert size (bp)", y = "Frequency") + theme_minimal()

Tendremos una gráfica como la siguiente: