Thursday, January 14, 2010
Find nnnn in genome
perl -e'open(FILE,"chr1.fa");my $content=""; while(){chomp $_; $content.=$_;} while($content=~/(NNNNNNN.*?N[A|C|G|T|a|c|g|t])/g){print length($1),"\t",pos($content),"\n";}'
Monday, January 4, 2010
maq
esión 2 - Mapeo de secuencias cortas a un genoma con MAQ
Introducción
En esta sesión vamos a aprender a utilizar MAQ (Mapping and Assembly with Qualities) que se describe en el artículo: “Mapping short DNA sequencing reads and calling variants using mapping quality scores” (pdf)
Vamos a utilizar datos como los que obtuvimos al final de la práctica pasada, es decir, archivos en formato fastq ya depurados y recortados. Para que todos podamos trabajar con los mismos datos, les estoy incluyendo la liga a una carpeta.
Se necesita lo siguiente para realizar esta práctica
maq Instalado (src)
maqview No instalado (src)
Archivos en formato fastq carpeta
Genoma de referencia Salmonella CT18
El manual de maq viene en el src, y está disponible en pdf por si necesitan consultarlo.
Una vez más, durante esta página web, secciones de ejemplo de código, que pueden copiar directamente serán mostrados así:
echo "Hello world!"
Los resultados o salidas, se indicarán así:
Hello world
Funciones de MAQ y procesamiento de archivos de entrada
Dentro de un sólo programa maq tenemos muchas funciones, que van generando archivos con distintos propósitos. El siguiente diagrama de flujo muestra estas funciones, así como sus entradas y salidas.
Como se ve en la figura, maq utiliza como entradas un genoma de referencia en formato fasta y archivos de secuencia en formato fastq. Sin embargo, utiliza formatos de archivo binarios, para ser más eficiente. Lo primero que necesitamos hacer es convertir todos nuestros archivos de secuencia a su formato binario correspondiente.
maq fasta2bfa genome.fa genome.bfa # Para genomas de referencia, formato fasta
maq fastq2bfq reads.fq reads.bfq # Para archivos de secuencias, formato fastq
Con un ejemplo de los archivos que estamos usando:
gzip -dc Salmonella_CT18.fa.gz | maq fasta2bfa - ref.bfa
-- 3 sequences have been converted.
gzip -dc Typhi_E022759_solexa_3.a.filtered.fq.gz | maq fastq2bfq - solexa_3.a.bfq
-- finish writing file 'solexa_3.a.bfq'
-- 248010 sequences were loaded.
Mapeo de secuencias al genoma de referencia
El siguiente paso es mapear las secuencias al genoma de referencia:
maq map # Dejando las opciones que nos interesan
Usage: maq map [options]
Options: -1 INT length of the first read (<=127) [0]
-e INT maximum allowed sum of qualities of mismatches [70]
-n INT number of mismatches in the first 24bp [2]
-u FILE dump unmapped and poorly aligned reads to FILE [null]
Debido a que ya tenemos nuestras secuencias recortadas, no necesitamos la opcion -1. Las siguientes dos opciones van a controlar la cantidad de secuencias que podremos mapear. Si por alguna razón esperamos que nuestras secuencias tengan más de 2 errores, y que estas posiciones sean de alta calidad podríamos incrementarlas un poco. Sin embargo, es poco probable que haya más de un SNP dentro de una secuencia corta. Si las incrementamos, encontraremos más SNPs, pero también tendremos más falsos positivos. El incrementar -n ocasionará tambien que maq use más memoria y sea más lento. Si queremos guardar las secuencias que no pueden ser mapeadas (para posteriores análisis, por ejemplo para tratar de ensamblarlas) nos será útil la opción -u. Podemos intentar entonces:
maq map -n 2 -e 70 -u solexa_3.a.unmap solexa_3.a.map ref.bfa solexa_3.a.bfq
-- maq-0.7.1
[ma_load_reads] loading reads...
...
# tardó 2:30 minutos en palenque
Las secuencias que no pudieron ser mapeadas con los parámetros que elegimos las podemos ver directamente, pues es un archivo de texto:
more solexa_3.a.unmap
Si queremos a estas alturas ver el resultado del alineamiento, tenemos que extraerlo con el siguiente comando:
maq mapview solexa_3.a.map > solexa_3.a.map.txt
head solexa_3.a.map.txt
IL2_30_2_3_400_762 chr 121 - 0 0 85 85 85 0 0 1 0 30 cGgGCAGATACTTTAACCAATATAGGAATA 0I2IIIIIIIIIIIIIIIIIIIIIIIIIII
IL2_30_2_17_1003_180 chr 168 + 0 0 76 76 76 0 0 1 0 30 AATGACAGAGTACACAAcAtcCaTgaAcCg IIIIIIIIIIIIIIIII5I0*I0I.,I8G%
Las columnas de este archivo se describen a continuación:
1 identificador del read
2 nombre de la secuencia de referencia
3 posición en la secuencia de referencia
4 cadena de la secuencia de referencia
5-6 información para mapeo de secuencias apareadas
7-9 calidad de mapeo, idénticas por no ser secuencias apareadas
10 número de mismatches del mejor hit
11 suma de calidades de los mismatches del mejor hit
12 número de hits con 0-mismatches en las primeras 24 bases
13 número de hits con 1-mismatches en las primeras 24 bases
14 longitud del read
15 secuencia del read
16 calidades del read
El archivo .map es uno de los centrales del proceso de maq. Cuando tengamos todos los archivos .map correspondientes a los archivos fastq iniciales, podemos combinarlos usando:
maq mapmerge all.map solexa_3.a.map solexa_3.b.map solexa_3.c.map ...
Otras funciones para extraer información del archivo de mapeo:
maq mapcheck -s ref.bfa all.map
Number of reference sequences: 3
Length of reference sequences exlcuding gaps: 5133713
Length of gaps in the reference sequences: 0
Length of non-gap regions covered by reads: 4845450
Length of 24bp unique regions of the reference: 4694282
Reference nucleotide composition: A: 24.01%, C: 25.90%, G: 25.97%, T: 24.12%
Reads nucleotide composition: A: 23.58%, C: 26.32%, G: 26.37%, T: 23.73%
Average depth across all non-gap regions: 33.773
Average depth across 24bp unique regions: 33.854
A C G T : AC AG AT CA CG CT GA GC GT TA TC TG : 0? 1? 2? 3? 4? : 0? 1? 2? 3? 4?
1 25.3 25.4 25.0 24.3 : 0 0 0 3 0 1 2 1 0 2 1 1 : 4 12 20 30 933 : 248 169 0 0 0
2 24.6 24.8 25.9 24.7 : 1 1 0 2 0 0 0 0 0 1 1 2 : 5 15 25 35 919 : 200 33 7 3 1
3 24.6 25.8 25.0 24.6 : 1 1 0 4 0 1 0 0 0 1 1 2 : 6 18 29 40 906 : 191 30 8 4 1
4 23.9 25.7 26.7 23.7 : 1 1 0 3 0 0 0 0 0 1 1 3 : 6 17 28 38 909 : 204 34 9 4 1
5 25.1 24.2 26.6 24.1 : 2 2 1 5 1 1 2 1 1 2 2 4 : 10 17 28 37 907 : 405 41 8 4 1
...
25 23.6 26.6 26.4 23.5 : 15 5 2 9 5 2 1 3 3 3 7 8 : 126 206 161 120 385 : 99 8 2 1 1
26 23.6 26.0 26.9 23.6 : 13 9 3 12 13 4 1 2 3 3 6 13 : 149 209 159 117 364 : 114 10 3 2 1
27 23.2 26.7 26.8 23.4 : 23 9 3 12 10 4 1 5 5 4 10 13 : 178 238 166 114 304 : 117 10 3 2 1
28 23.4 26.0 26.9 23.7 : 21 11 6 16 14 8 2 4 8 4 8 12 : 203 248 165 109 273 : 121 12 3 2 1
29 22.7 26.7 27.1 23.5 : 46 19 9 20 20 10 3 10 10 7 18 20 : 245 265 161 100 228 : 169 18 4 2 1
30 22.9 26.1 27.0 24.0 : 39 20 16 27 25 20 4 11 16 8 16 20 : 278 274 156 93 199 : 173 25 0 0 0
<-contenido de bases-> <-frecuencia de cambios observados-> <-calidades reads -> <-cambios/calidad->
maq mapstat all.map
-- Total number of reads: 5791504
-- Sum of read length: 173745120
-- Error rate: 0.016155
-- Average read length: 30.00
-- Mismatch statistics:
-- MM 0 4051603
-- MM 1 1036434
-- MM 2 423348
-- MM 3 206417
-- MM 4 64118
-- MM 5 9584
-- Mapping quality statistics:
-- MQ 00-09 146191 146191
-- MQ 10-19 150298 150298
-- MQ 20-29 18452 18452
-- MQ 30-39 240481 240481
-- MQ 40-49 389043 389043
-- MQ 50-59 638169 638169
-- MQ 60-69 475203 475203
-- MQ 70-79 1216844 1216844
-- MQ 80-89 2312666 2312666
-- MQ 90-99 204157 204157
Para obtener un resúmen por cada base de nuestra referencia, sobre la cobertura y las secuencias que ahí maapean, podemos generar un archivo llamado pileup.
maq pileup -v ref.bfa all.map | head -n 20
chr 1 A 0 @ @ @
chr 2 G 0 @ @ @
chr 3 A 0 @ @ @
chr 4 G 0 @ @ @
chr 5 A 0 @ @ @
chr 6 T 1 @, @I @v
chr 7 T 1 @, @I @v
chr 8 A 2 @,g @I$ @vy
chr 9 C 3 @,., @III @vyh
chr 10 G 4 @,.,. @IIII @vyhz
chr 11 T 4 @,.,. @IIII @vyhz
chr 12 C 5 @,.,.. @II8II @vyhz{
chr 13 T 6 @,.,.., @I4II1I @vyhz{{
chr 14 G 9 @,.,..,,,. @IIIIIIII9 @vyhz{{uv{
chr 15 G 9 @,.,..,,,. @IIIIIIII< @vyhz{{uv{
chr 16 T 9 @,.,..,,,. @HIIIIIIII @vyhz{{uv{
chr 17 T 10 @,.,..,,,., @I?IIIIIIII @vyhz{{uv{v
chr 18 G 11 @,.,..,,,.,, @IIIIIIIIIII @vyhz{{uv{vw
chr 19 C 11 @,.,..,,,.,, @IIHIIIIIIII @vyhz{{uv{vw
chr 20 A 11 @,.,..,,,.,, @IIIIIIIIIII @vyhz{{uv{vw
El cual incluye el cromosoma, posición, base de referencia, cobertura en esa posición y 3 columnas codificadas: secuencias de los reads, calidad de los reads, calidad de mapeo. Estas últimas columnas siempre empiezan con "@". La secuencia de los reads se muestra en "." o "," cuando la base es igual a la de referencia si no, se muestra la base. Para indicar la cadena "+" se usa la "," o una letra mayúscula, para la cadena "-" se usa el "." o una letra minúscula.
Generación del consenso de mapeo y predicción de SNPs
El siguiente archivo central del proceso maq es el consenso (.cns). Lo obtenemos así:
maq assemble -N 1 -s consensus.cns ref.bfa all.map
A partir del consenso, existen varias funciones para extraer archivos con distinto tipo de información. Podemos por ejemplo extraer el consenso completo, es decir la secuencia del nuevo genoma de acuerdo al consenso de mapeo, junto con sus valores de calidad de cada base, en formato fastq:
maq cns2fq consensus.cns > cns.fq
head(cns.fq)
@chr
nnnnnttaCGTCTGGTTGCAAGAGATCATAACAGGGGAAATTGATTGAAAATAAATATAT
CGCCAGCAGCACATGAACAAGTTTCGGAATGTGATCAATTTAAAAATTTATTGACTTAGG
CGGGCAGATACTTTAACCAATATAGGAATACAAGACAGACAAATAAAAATGACAGAGTAC
ACAACATCCATGAACCGCATCAGCACCACCACCATTACCACCATCACCATTACCACAGGT
AACGGTGCGGGCTGACGCGTACAGGAAACACAGAAAAAAGCCCGCACCTGAACAGTGCGG
GCTTTTTTTTCGACCAGAGATCACGAGGTAACAACCATGCGAGTGTTGAAGTTCGGCGGT
ACATCAGTGGCAAATGCAGAACGTTTTCTGCGTGTTGCCGATATTCTGGAAAGCAATTCC
AGGCAAGGGCAGGTAGCGACCGTACTTTCCGCCCCCGCGAAAATTACCAACCATCTGGTG
GCGATGATTGAAAAAACTATCGGCGGCCAGGATGCTTTGCCGAATATCAGCGATGCCGAA
El objetivo de muchos estudios de resecuenciación es encontrar los SNPs. Para obtener todas las posiciones que varían entre el consenso y el genoma de referencia usamos cns2snp:
maq cns2snp consensus.cns > cns.snp
head cns.snp
chr 12478 C T 255 39 1.00 63 62 C 255 Y
chr 17187 C T 255 33 1.00 63 62 N 255 N
chr 22623 T C 255 36 1.00 63 62 A 255 M
los campos de esta salida se describen a continuación:
1 nombre del cromosoma de referencia
2 posición en el cromosoma
3 base de referencia
4 base consenso
5 calidad del consenso
6 profundidad de la cobertura
7 número promedio de hits de los reads que cubren esta posición
8 calidad de mapeo más alta de los reads que cubren esta posición
9 calidad mínima del consenso de las 3 bases de cada lado del sitio
10 segunda base consenso
11 log likelihood ratio de la segunda y tercera base
12 tercera base consenso
En particular el campo 5 nos define la calidad del consenso y por lo tanto la confianza que tenemos de que es un SNP verdadero. También es bueno considerar la calidad de las bases vecinas, definidas en el campo 9. Una región de baja calidad puede implicar un problema de mapeo, no variabilidad puntual. La columna 7 no da información sobre la redundancia de esta posición. Si el número es cercano a 1, esta región es única, pero si es mayor puede implicar una duplicación o una región repetitiva.
Debido que estos SNP potenciales pueden contener muchos falsos positivos, conviene filtrarlos para elegir los mejores candidatos. Para esto maq incluye un script de perl que nos ayuda:
maq.pl SNPfilter
Usage: maq.pl SNPfilter [options]
Options: -d INT minimum depth to call a SNP [3]
-D INT maximum depth (<=254), otherwise ignored [256]
-n INT minimum neighbouring quality [20]
-Q INT required max mapping quality of the reads covering the SNP [40]
-q INT minimum consensus quality [20]
-w INT size of the window in which SNPs should be filtered out [5] (cerca de un indel)
-S FILE splitread output [null]
-F FILE indelpe output [null]
-f FILE indelsoa output [null]
-s INT indelsoa score (= left_clip + right_clip - across) [3]
-m INT indelsoa: max number of reads mapped across the indel [1]
-W INT window size for filtering dense SNPs [10]
-N INT maximum number of SNPs in a window [2]
-a alternative filter for single end reads
Recuerden que no podemos usar información de indeles porque no contamos con secuencias apareadas. Podemos usar algunas de las reglas que sugieren en el artículo de maq:
maq.pl SNPfilter -d 4 -n 20 -Q 40 -q 40 -W 10 -N 3 cns.snp > filtered.snp
wc -l cns.snp filtered.snp
7489 cns.snp
256 filtered.snp
Como podrán ver, se puede eliminar un gran numero de SNPs de esta manera.
Visualización del mapeo con maqview
En un momento dado vamos a querer explorar manualmente ciertas regiones de mapeo, para visualizar mejor la cobertura, o la calidad del mapeo al rededor de un SNP potencial. Para esto nos sirven programas como maqview.
Primero debemos indexar los archivos de mapeo y consenso:
maqindex -i -c consensus.cns all.map
Maqview es un programa que utliza OpenGL para darnos una interfaz gráfica al mapeo. Nos permite navegar rápidamente por nuestro genoma de referencia, observando las secuencias que han sido mapeadas, posiciones de error, calidad de secuencia, calidad de mapeo, SNPs, etc.
maqview -c consensus.cns all.map
Un ejemplo mostrando un SNP
Una de las modalidades gráficas, los colores indican la cadena de mapeo
Ejercicios:
Escriban un script para poder trabajar fácilmente con maq y múltiples archivos fastq de entrada.
Procesen 10 de los archivos fastq, hasta el punto de predecir SNPs.
Comparen sus predicciones con las de este archivo: filtered.snp
Si asumen que esos son los verdaderos SNPs, qué pueden decir sobre falsos/verdaderos positivos/negativos?
Prueben algunas combinaciones de filtrado, para ver cómo cambia esto.
Introducción
En esta sesión vamos a aprender a utilizar MAQ (Mapping and Assembly with Qualities) que se describe en el artículo: “Mapping short DNA sequencing reads and calling variants using mapping quality scores” (pdf)
Vamos a utilizar datos como los que obtuvimos al final de la práctica pasada, es decir, archivos en formato fastq ya depurados y recortados. Para que todos podamos trabajar con los mismos datos, les estoy incluyendo la liga a una carpeta.
Se necesita lo siguiente para realizar esta práctica
maq Instalado (src)
maqview No instalado (src)
Archivos en formato fastq carpeta
Genoma de referencia Salmonella CT18
El manual de maq viene en el src, y está disponible en pdf por si necesitan consultarlo.
Una vez más, durante esta página web, secciones de ejemplo de código, que pueden copiar directamente serán mostrados así:
echo "Hello world!"
Los resultados o salidas, se indicarán así:
Hello world
Funciones de MAQ y procesamiento de archivos de entrada
Dentro de un sólo programa maq tenemos muchas funciones, que van generando archivos con distintos propósitos. El siguiente diagrama de flujo muestra estas funciones, así como sus entradas y salidas.
Como se ve en la figura, maq utiliza como entradas un genoma de referencia en formato fasta y archivos de secuencia en formato fastq. Sin embargo, utiliza formatos de archivo binarios, para ser más eficiente. Lo primero que necesitamos hacer es convertir todos nuestros archivos de secuencia a su formato binario correspondiente.
maq fasta2bfa genome.fa genome.bfa # Para genomas de referencia, formato fasta
maq fastq2bfq reads.fq reads.bfq # Para archivos de secuencias, formato fastq
Con un ejemplo de los archivos que estamos usando:
gzip -dc Salmonella_CT18.fa.gz | maq fasta2bfa - ref.bfa
-- 3 sequences have been converted.
gzip -dc Typhi_E022759_solexa_3.a.filtered.fq.gz | maq fastq2bfq - solexa_3.a.bfq
-- finish writing file 'solexa_3.a.bfq'
-- 248010 sequences were loaded.
Mapeo de secuencias al genoma de referencia
El siguiente paso es mapear las secuencias al genoma de referencia:
maq map # Dejando las opciones que nos interesan
Usage: maq map [options]
Options: -1 INT length of the first read (<=127) [0]
-e INT maximum allowed sum of qualities of mismatches [70]
-n INT number of mismatches in the first 24bp [2]
-u FILE dump unmapped and poorly aligned reads to FILE [null]
Debido a que ya tenemos nuestras secuencias recortadas, no necesitamos la opcion -1. Las siguientes dos opciones van a controlar la cantidad de secuencias que podremos mapear. Si por alguna razón esperamos que nuestras secuencias tengan más de 2 errores, y que estas posiciones sean de alta calidad podríamos incrementarlas un poco. Sin embargo, es poco probable que haya más de un SNP dentro de una secuencia corta. Si las incrementamos, encontraremos más SNPs, pero también tendremos más falsos positivos. El incrementar -n ocasionará tambien que maq use más memoria y sea más lento. Si queremos guardar las secuencias que no pueden ser mapeadas (para posteriores análisis, por ejemplo para tratar de ensamblarlas) nos será útil la opción -u. Podemos intentar entonces:
maq map -n 2 -e 70 -u solexa_3.a.unmap solexa_3.a.map ref.bfa solexa_3.a.bfq
-- maq-0.7.1
[ma_load_reads] loading reads...
...
# tardó 2:30 minutos en palenque
Las secuencias que no pudieron ser mapeadas con los parámetros que elegimos las podemos ver directamente, pues es un archivo de texto:
more solexa_3.a.unmap
Si queremos a estas alturas ver el resultado del alineamiento, tenemos que extraerlo con el siguiente comando:
maq mapview solexa_3.a.map > solexa_3.a.map.txt
head solexa_3.a.map.txt
IL2_30_2_3_400_762 chr 121 - 0 0 85 85 85 0 0 1 0 30 cGgGCAGATACTTTAACCAATATAGGAATA 0I2IIIIIIIIIIIIIIIIIIIIIIIIIII
IL2_30_2_17_1003_180 chr 168 + 0 0 76 76 76 0 0 1 0 30 AATGACAGAGTACACAAcAtcCaTgaAcCg IIIIIIIIIIIIIIIII5I0*I0I.,I8G%
Las columnas de este archivo se describen a continuación:
1 identificador del read
2 nombre de la secuencia de referencia
3 posición en la secuencia de referencia
4 cadena de la secuencia de referencia
5-6 información para mapeo de secuencias apareadas
7-9 calidad de mapeo, idénticas por no ser secuencias apareadas
10 número de mismatches del mejor hit
11 suma de calidades de los mismatches del mejor hit
12 número de hits con 0-mismatches en las primeras 24 bases
13 número de hits con 1-mismatches en las primeras 24 bases
14 longitud del read
15 secuencia del read
16 calidades del read
El archivo .map es uno de los centrales del proceso de maq. Cuando tengamos todos los archivos .map correspondientes a los archivos fastq iniciales, podemos combinarlos usando:
maq mapmerge all.map solexa_3.a.map solexa_3.b.map solexa_3.c.map ...
Otras funciones para extraer información del archivo de mapeo:
maq mapcheck -s ref.bfa all.map
Number of reference sequences: 3
Length of reference sequences exlcuding gaps: 5133713
Length of gaps in the reference sequences: 0
Length of non-gap regions covered by reads: 4845450
Length of 24bp unique regions of the reference: 4694282
Reference nucleotide composition: A: 24.01%, C: 25.90%, G: 25.97%, T: 24.12%
Reads nucleotide composition: A: 23.58%, C: 26.32%, G: 26.37%, T: 23.73%
Average depth across all non-gap regions: 33.773
Average depth across 24bp unique regions: 33.854
A C G T : AC AG AT CA CG CT GA GC GT TA TC TG : 0? 1? 2? 3? 4? : 0? 1? 2? 3? 4?
1 25.3 25.4 25.0 24.3 : 0 0 0 3 0 1 2 1 0 2 1 1 : 4 12 20 30 933 : 248 169 0 0 0
2 24.6 24.8 25.9 24.7 : 1 1 0 2 0 0 0 0 0 1 1 2 : 5 15 25 35 919 : 200 33 7 3 1
3 24.6 25.8 25.0 24.6 : 1 1 0 4 0 1 0 0 0 1 1 2 : 6 18 29 40 906 : 191 30 8 4 1
4 23.9 25.7 26.7 23.7 : 1 1 0 3 0 0 0 0 0 1 1 3 : 6 17 28 38 909 : 204 34 9 4 1
5 25.1 24.2 26.6 24.1 : 2 2 1 5 1 1 2 1 1 2 2 4 : 10 17 28 37 907 : 405 41 8 4 1
...
25 23.6 26.6 26.4 23.5 : 15 5 2 9 5 2 1 3 3 3 7 8 : 126 206 161 120 385 : 99 8 2 1 1
26 23.6 26.0 26.9 23.6 : 13 9 3 12 13 4 1 2 3 3 6 13 : 149 209 159 117 364 : 114 10 3 2 1
27 23.2 26.7 26.8 23.4 : 23 9 3 12 10 4 1 5 5 4 10 13 : 178 238 166 114 304 : 117 10 3 2 1
28 23.4 26.0 26.9 23.7 : 21 11 6 16 14 8 2 4 8 4 8 12 : 203 248 165 109 273 : 121 12 3 2 1
29 22.7 26.7 27.1 23.5 : 46 19 9 20 20 10 3 10 10 7 18 20 : 245 265 161 100 228 : 169 18 4 2 1
30 22.9 26.1 27.0 24.0 : 39 20 16 27 25 20 4 11 16 8 16 20 : 278 274 156 93 199 : 173 25 0 0 0
<-contenido de bases-> <-frecuencia de cambios observados-> <-calidades reads -> <-cambios/calidad->
maq mapstat all.map
-- Total number of reads: 5791504
-- Sum of read length: 173745120
-- Error rate: 0.016155
-- Average read length: 30.00
-- Mismatch statistics:
-- MM 0 4051603
-- MM 1 1036434
-- MM 2 423348
-- MM 3 206417
-- MM 4 64118
-- MM 5 9584
-- Mapping quality statistics:
-- MQ 00-09 146191 146191
-- MQ 10-19 150298 150298
-- MQ 20-29 18452 18452
-- MQ 30-39 240481 240481
-- MQ 40-49 389043 389043
-- MQ 50-59 638169 638169
-- MQ 60-69 475203 475203
-- MQ 70-79 1216844 1216844
-- MQ 80-89 2312666 2312666
-- MQ 90-99 204157 204157
Para obtener un resúmen por cada base de nuestra referencia, sobre la cobertura y las secuencias que ahí maapean, podemos generar un archivo llamado pileup.
maq pileup -v ref.bfa all.map | head -n 20
chr 1 A 0 @ @ @
chr 2 G 0 @ @ @
chr 3 A 0 @ @ @
chr 4 G 0 @ @ @
chr 5 A 0 @ @ @
chr 6 T 1 @, @I @v
chr 7 T 1 @, @I @v
chr 8 A 2 @,g @I$ @vy
chr 9 C 3 @,., @III @vyh
chr 10 G 4 @,.,. @IIII @vyhz
chr 11 T 4 @,.,. @IIII @vyhz
chr 12 C 5 @,.,.. @II8II @vyhz{
chr 13 T 6 @,.,.., @I4II1I @vyhz{{
chr 14 G 9 @,.,..,,,. @IIIIIIII9 @vyhz{{uv{
chr 15 G 9 @,.,..,,,. @IIIIIIII< @vyhz{{uv{
chr 16 T 9 @,.,..,,,. @HIIIIIIII @vyhz{{uv{
chr 17 T 10 @,.,..,,,., @I?IIIIIIII @vyhz{{uv{v
chr 18 G 11 @,.,..,,,.,, @IIIIIIIIIII @vyhz{{uv{vw
chr 19 C 11 @,.,..,,,.,, @IIHIIIIIIII @vyhz{{uv{vw
chr 20 A 11 @,.,..,,,.,, @IIIIIIIIIII @vyhz{{uv{vw
El cual incluye el cromosoma, posición, base de referencia, cobertura en esa posición y 3 columnas codificadas: secuencias de los reads, calidad de los reads, calidad de mapeo. Estas últimas columnas siempre empiezan con "@". La secuencia de los reads se muestra en "." o "," cuando la base es igual a la de referencia si no, se muestra la base. Para indicar la cadena "+" se usa la "," o una letra mayúscula, para la cadena "-" se usa el "." o una letra minúscula.
Generación del consenso de mapeo y predicción de SNPs
El siguiente archivo central del proceso maq es el consenso (.cns). Lo obtenemos así:
maq assemble -N 1 -s consensus.cns ref.bfa all.map
A partir del consenso, existen varias funciones para extraer archivos con distinto tipo de información. Podemos por ejemplo extraer el consenso completo, es decir la secuencia del nuevo genoma de acuerdo al consenso de mapeo, junto con sus valores de calidad de cada base, en formato fastq:
maq cns2fq consensus.cns > cns.fq
head(cns.fq)
@chr
nnnnnttaCGTCTGGTTGCAAGAGATCATAACAGGGGAAATTGATTGAAAATAAATATAT
CGCCAGCAGCACATGAACAAGTTTCGGAATGTGATCAATTTAAAAATTTATTGACTTAGG
CGGGCAGATACTTTAACCAATATAGGAATACAAGACAGACAAATAAAAATGACAGAGTAC
ACAACATCCATGAACCGCATCAGCACCACCACCATTACCACCATCACCATTACCACAGGT
AACGGTGCGGGCTGACGCGTACAGGAAACACAGAAAAAAGCCCGCACCTGAACAGTGCGG
GCTTTTTTTTCGACCAGAGATCACGAGGTAACAACCATGCGAGTGTTGAAGTTCGGCGGT
ACATCAGTGGCAAATGCAGAACGTTTTCTGCGTGTTGCCGATATTCTGGAAAGCAATTCC
AGGCAAGGGCAGGTAGCGACCGTACTTTCCGCCCCCGCGAAAATTACCAACCATCTGGTG
GCGATGATTGAAAAAACTATCGGCGGCCAGGATGCTTTGCCGAATATCAGCGATGCCGAA
El objetivo de muchos estudios de resecuenciación es encontrar los SNPs. Para obtener todas las posiciones que varían entre el consenso y el genoma de referencia usamos cns2snp:
maq cns2snp consensus.cns > cns.snp
head cns.snp
chr 12478 C T 255 39 1.00 63 62 C 255 Y
chr 17187 C T 255 33 1.00 63 62 N 255 N
chr 22623 T C 255 36 1.00 63 62 A 255 M
los campos de esta salida se describen a continuación:
1 nombre del cromosoma de referencia
2 posición en el cromosoma
3 base de referencia
4 base consenso
5 calidad del consenso
6 profundidad de la cobertura
7 número promedio de hits de los reads que cubren esta posición
8 calidad de mapeo más alta de los reads que cubren esta posición
9 calidad mínima del consenso de las 3 bases de cada lado del sitio
10 segunda base consenso
11 log likelihood ratio de la segunda y tercera base
12 tercera base consenso
En particular el campo 5 nos define la calidad del consenso y por lo tanto la confianza que tenemos de que es un SNP verdadero. También es bueno considerar la calidad de las bases vecinas, definidas en el campo 9. Una región de baja calidad puede implicar un problema de mapeo, no variabilidad puntual. La columna 7 no da información sobre la redundancia de esta posición. Si el número es cercano a 1, esta región es única, pero si es mayor puede implicar una duplicación o una región repetitiva.
Debido que estos SNP potenciales pueden contener muchos falsos positivos, conviene filtrarlos para elegir los mejores candidatos. Para esto maq incluye un script de perl que nos ayuda:
maq.pl SNPfilter
Usage: maq.pl SNPfilter [options]
Options: -d INT minimum depth to call a SNP [3]
-D INT maximum depth (<=254), otherwise ignored [256]
-n INT minimum neighbouring quality [20]
-Q INT required max mapping quality of the reads covering the SNP [40]
-q INT minimum consensus quality [20]
-w INT size of the window in which SNPs should be filtered out [5] (cerca de un indel)
-S FILE splitread output [null]
-F FILE indelpe output [null]
-f FILE indelsoa output [null]
-s INT indelsoa score (= left_clip + right_clip - across) [3]
-m INT indelsoa: max number of reads mapped across the indel [1]
-W INT window size for filtering dense SNPs [10]
-N INT maximum number of SNPs in a window [2]
-a alternative filter for single end reads
Recuerden que no podemos usar información de indeles porque no contamos con secuencias apareadas. Podemos usar algunas de las reglas que sugieren en el artículo de maq:
maq.pl SNPfilter -d 4 -n 20 -Q 40 -q 40 -W 10 -N 3 cns.snp > filtered.snp
wc -l cns.snp filtered.snp
7489 cns.snp
256 filtered.snp
Como podrán ver, se puede eliminar un gran numero de SNPs de esta manera.
Visualización del mapeo con maqview
En un momento dado vamos a querer explorar manualmente ciertas regiones de mapeo, para visualizar mejor la cobertura, o la calidad del mapeo al rededor de un SNP potencial. Para esto nos sirven programas como maqview.
Primero debemos indexar los archivos de mapeo y consenso:
maqindex -i -c consensus.cns all.map
Maqview es un programa que utliza OpenGL para darnos una interfaz gráfica al mapeo. Nos permite navegar rápidamente por nuestro genoma de referencia, observando las secuencias que han sido mapeadas, posiciones de error, calidad de secuencia, calidad de mapeo, SNPs, etc.
maqview -c consensus.cns all.map
Un ejemplo mostrando un SNP
Una de las modalidades gráficas, los colores indican la cadena de mapeo
Ejercicios:
Escriban un script para poder trabajar fácilmente con maq y múltiples archivos fastq de entrada.
Procesen 10 de los archivos fastq, hasta el punto de predecir SNPs.
Comparen sus predicciones con las de este archivo: filtered.snp
Si asumen que esos son los verdaderos SNPs, qué pueden decir sobre falsos/verdaderos positivos/negativos?
Prueben algunas combinaciones de filtrado, para ver cómo cambia esto.
Subscribe to:
Posts (Atom)