PRÁCTICA DE ALINEAMIENTO o MAPEO DE SECUENCIAS RNA-Seq con TopHat

Prededentes:

En la carpeta de la UCO /usr/local/uco/bms/practica3, hay 3 carpetas:

Objetivos de la práctica

  1. Buscar con Google el programa bowtie2, e inspeccionar la página para encontrar cómo instalarlo y obtener el manual que te indica cómo se usa.
  2. Lo mismo con el programa Tophat.
  3. Reconocer el contenido del archivo chr1.fa (que contiene el cromosoma 1 humano, versión de anotación 18 que usaremos con genoma de referencia)
  4. Reconocer el contenido de su correspondiente archivo de anotaciones tipo gtf. Nota: usaremos solo la secuencia del cromosoma 1, y no del genoma humano completo para no generar archivos de gran tamaño y extrema exigencias de computación)
  5. Reconocer el contenido de los archivos fastq correspondientes a un experimento de RNA-Seq. Ver su contenido y la información que nos proporciona.
    1. hay que identificar si son secuencias pareadas o no
    2. Si son archivos sincronizados, esto es, si tienen el mismo número de lecturas cada uno, y si los nombres de las lecturas coinciden en orden. Para ello usad cat en combinación con wc -l, y head y tail
  6. Usar el programa adecuado de bowtie2 para indexar el archivo fasta del cromosoma 1 humano que usaremos como genoma de referencia, algo que es imprescindible para mapear rápidamente las secuencias obtenidas de un experimento RNA-Seq
  7. Una vez creado el archivo índice con bowtie2, realizaremos el mapeo con TopHat, uno de los programas más conocidos de mapeo transcriptómico. Señalar que hay otros mapeadores, como STAR, HISAT, o incluso otras versiones como TopHat2, HISAT2.
  8. Veremos por primera vez lo que es un archivo BAM/SAM, que es el archivo que TopHat y otros mapeadores dan como resultado, y lo que son los archivos BED, ya que TopHat proporciona información de las delecciones y las inserciones con archivos de este formato BED
  9. Analizaremos la información que TopHat da sobre el mapeo, como el porcentaje de lecturas que mapean, las delecciones, las inserciones, los junctions...
  10. Haremos un control de calidad del mapeo obtenido analizando el archivo BAM con el programa samstat.
  11. Aprenderemos cómo podemos manipular y extraer información de los archivos SAM/BAM usando el conjunto de programas incluido en samtools. Aprenderemos a ver el contenido interno de dichos archivos, la estadística que aporta, como extraer información de los archivos SAM, y un largo etcétera
  12. Aprenderemos las bases de cómo interpretar la información del apartado FLAGS con ayuda de la información proporcionada por esta página WEB
  13. Aprenderemos a obtener información del contenido del apartado CIGAR y de otros campos
  14. Aprenderemos a extraer alguna información útil de los archivos SAM/BAM con ayuda de cut, grep, sort, uniq, etc..
  15. Aprenderemos a ordenar primero y luego a indexar el archivo de mapeo BAM/SAM con ayuda de samtools generando el correspondiente archivo de índices .bai. Muchos programas exigen que los archivos BAM estén ordenados e indexados, como el programa IGV (Integrative Genome Viewer) que veremos en el siguiente punto. El indexado u organización de estos archivos son necesarios para que la información que contienen se muevan fluidos por los programas
  16. Cargaremos el genoma de referencia, su correspondiente archivo gtf y el archivo BAM de resultados indexado en el programa IGV, y veremos qué información nos porporciona dicho programa. Veremos la cobertura del mapeo, inserciones, delecciones y variantes en comparación con el genoma de referencia. Añadiremos los archivos junctions.bed para ver qué información nos da. Para ello visitaremos en parte la ayuda del programa IGV, como ésta que te da información de los junctions
  17. Y si nos da tiempo, haremos un contaje de secuencias RNA-Seq mapeadas, usando el archivo bam/sam y el programa htseq-count
  18. No podremos hacer un análisis de expresión diferencial porque solo disponemos de una muestra, pero podemos sentar las bases de lo que hay que hacer

EJECUCIÓN DE LA PRÁCTICA

  1. Como se supone que ya hemos aprendido a instalar y compilar programas como velvet bajo linux (que no es más que leerse las instrucciones y usar en su caso la orden make, o la configure/make), le he pedido al administrador del sistema de la UCO que intale por nosotros los siguientes programas. Todos los programas están definidos en el PATH y por tanto se pueden invocar sin necesidad de usar el "./"
  1. Si tienes un ordenador propio que tiene instalado Linux, haz una copia de las secuencias localmente. Te hará falta el genoma de referencia, el GTF, y las dos secuencias stranded
  2. Si no tienes un ordenador con una instalación local de Linux habrás de trabajar como es acostumbrado o bien en los ordenadores de las aulas de la UCO o mediante una conexión con ts.uco.es
  3. Antes que nada, deberíamos hacer un control de calidad de las secuencias con ayuda de FastQC. Lo haremos sobre las secuencias incluidas en la carpeta /usr/local/uco/bms/practica3/stranded. Debeis saber:
    1. Si las secuencias son pareadas o no.
    2. Y si son secuencias pareadas, si los archivos están sincronizados, esto es, que cada lectura tenga su compañero (mate) en el otro archivo. O dicho de otro modo, que no haya lecturas huérfanas. Lo podéis hacer con un wc -l (para ver el número total de líneas), con un head y un tail, para comprobar que el encabezamiento o la cola del archivo son los mismos. Si veis que algo difiere de la normalidad, usad un programa que lo sincronize, como el programa reformat.sh del paquete BBMap
  4. Si todo es correcto, podemos hacer el mapeo de las lecturas de RNA-Seq.
  5. Como TopHat mapea usando bowtie2, tenemos que indexar con bowtie2 el genoma de referencia
    1. Haz un "ls -l" de la carpeta /usr/local/uco/bms/practica3/reference para ver el contenido. Analiza qué es cada archivo
    2. Como no tenemos permiso para escribir en la carpeta pública /usr/local/uco/bms, y tenemos los problemas de espacio en nuestras carpetas locales le he pedido al administrador de la UCO que nos grabe en la carpeta reference los archivos de índice creado por bowtie2-build. Los encontramos en la carpeta reference.
    3. Lanza la orden bowtie2-build path_al_archivo/reference/chr1.fa chr1 (tendrás que hacerlo desde un sitio donde tengas permiso de escritura)
    4. Cuando veas que está funcionando, puedes optar por dejarlo andar o hacer un Control-C para detener el proceso. Y es que puedes o no tener espacio en tu cuenta personal para que se puedan crear esos archivos de índice. Por eso, y por prevenir, le he pedido al administrador que cree los archivos por nosotros que se ha guardado en /usr/local/uco/bms/practica3/reference
    5. Mira en la carpeta que archivos nuevos se han creado
  6. Mapea las secuencias contenidas en /usr/local/uco/bms/practica3/stranded con tophat usando la siguiente información
    1. que el insert size de los fragmentos shotgun es de 450, con una desviación estándar de 350 (siempre referido a bases)
    2. Ejecuta Tophat y mira todos los parámetros con los que se puede configurar. Anótalos. Puedes lanzar la orden $tophat >& ayuda.tophat para crear un archivo que luego con $less ayuda.tophat puedas ver cómodamente
    3. Vamos a crear un script para hacerlo todo más cómodo. Si tienes fallos o quieres añadir o quitar algo, lo puedes hacer con el editor vi o pico
      1. Crea una carpeta en el directorio temporal con mkdir /tmp/tu_login (donde tu_login es el nombre de tu cuenta, en mi caso, tu_login es bb1rofra)
      2. Escribe cat > tophat.sh (enter) o con pico tophat.sh
      • #! usr/bin/bash (enter)
      • $tophat -o /tmp/tu_login/TOPHAT \(enter)
      • --no-coverage-search \(enter)
      • -r 450 # La distancia entre las lecturas en el fragmento
      • --mate-std-dev 350 \(enter) # La desviación estándar de la distancia
      • -p 4 \(enter) # Usará cuatro núcleos. Si tienes más, aumenta el número
      • -G /usr/local/uco/bms/practica3/reference/gencode_v18_chr1.gtf \(enter)
      • --library-type fr-firststrand \(enter) # Esto es porque es stranded con el método dUTP
      • /usr/local/uco/bms/practica3/reference/chr1 \(enter) # Este es el índice
      • /usr/local/uco/bms/practica3/stranded/H1hesc_R1_chr1.fastq \(enter)
      • /usr/local/uco/bms/practica3/stranded/H1hesc_R2_chr1.fastq(enter)
      • (acaba el archivo y terminalo dándole a un Control + D)
      • Recuerda darle permisos de ejecución con chmod. En ese caso has de ejecutarlo como siempre, con el ./tophat.sh. O eso, o lo ejecutas con la orden bash tophat.sh aunque no tenga permiso de ejecución

       

      • Explicación
        • - o /tmp/tu_login/TOPHAT creará esa carpeta donde los resultados que TopHat genere serán guardados. Recuerda que antes de ejecutar este script deberás crear la carpeta /tmp/tu_login
        • --no coverage-search es en un principio prescindible. El cálculo de la cobertura exige tiempo y por eso no lo haremos en esta práctica (fíjate que hay un doble guión al principio)
        • -r 450 indica el tamano interno desconocido de los fragmentos. La separación entre los finales de las lecturas pareadas en un fragmento. Explicaremos en clase exactamente qué parte es la que hay que poner
        • --mate-std-dev 350 indica a tophat como de variable son las secuencias pareadas en longitud (fíjate el doble guión del principio)
        • Al introducir el archivo GTF en tophat se consiguen dos cosas. Por un lado introducir el nombre del gen al que mapea dentro del archivo resultado. Otro proposito es tener en cuenta la información de los exones e intrones que hay en la secuencia GTF como ayuda en el mapeo y para hacer el splicing. De esa forma se hace un mapeo asistido porque permite al programa generar un transcriptoma canónico, el especificado en el archivo GTF
        • --library-type fr-firststrand indica qué librería se está usando. En este caso se usa una librería stranded hecha con dUTP. Eso lo sabéis vosotros cuando hacéis la librería al pactarlo con la empresa de secuenciación (fíjate el doble guión del principio)
        • /usr/local/uco/bms/practica3/reference/chr1 senala al programa el lugar donde están los archivos índice del genoma de referencia. Como el índice tiene el mismo nombre que la secuencia fasta, se reconocen mutuamente
        • y en las dos ultimas partes debes indicar el path y el nombre de las dos secuencias pareadas
  7. Cuando ejecutes tophat, ve mirando lo que aparece en la pantalla y describe más o menos lo que está haciendo. Descubriras que
    1. hace algún filtrado de las secuencias por calidad.
    2. cómo invoca a bowtie2 para hacer los alineamientos para comprobar que el programa está disponible. Recuerda que TopHat no hace los mapeos, solo gestiona el uso del programa mapeador (bowtie2) para hacerlo sensible a la presencia de splicing o ayustamiento
    3. como usa la información del gtf para generar un transcriptoma, y como lo indexa con bowtie2-build. La generación de este transcriptoma permite mapear las lecturas "canónicas", las ya conocidas y reconocidas, con lo que ahorra bastante tiempo porque no tiene necesidad de "romper" las lecturas en las "junctions"
    4. que empieza a mapear las lecturas izquierda y derecha con ese transcriptoma
    5. que guarda por un lado las lecturas que han mapeado con ese transcriptoma y las que no han mapeado
    6. que luego "rompe" las lecturas que no han mapeado en 2-3 trozos (depende de la longitud de las lecturas) y las vuelve a mapear
  8. Cuanto termine, haz como mínio una copia del archivo accepted_hits.bam en tu cuenta personal. Ahora trabajaremos con ella

VIENDO LOS ALINEAMIENTOS CON IGV

  1. Abre IGV invocándolo directamente desde Linux. IGV es un programa creado en JAVA que funciona con cualquier sistema operativo. Dispones de un tutorial y de una guía de usuario en este enlace. El programa se ha instalado en la UCO, pero puede descargarse localmente e incluso ejecutarse desde la propia página IGV
  2. Abre el programa y ejecuta File / New Session para empezar desde cero
  3. Si ves que hay cargado datos, pulsa con el botón derecho del ratón sobre la columna izquierda que contiene esos datos y dale a la opción Remove Track. Esto hay que hacerlo porque IGV se carga de forma predeterminada con el genoma humano completo
  4. Con la opción Genomes, carga el archivo chr1.fa. Antes, deberás hacer una copia de él en tu carpeta principal, porque para usarlo, tendrás que tener permisos de escritura en dicha carpeta.
  5. Con la opción File, ve incorporando dentro de IGV el correspondiente archivo gtf de anotaciones. El programa te avisará de que tiene que indexarlo, pero lo hará él automaticamente siempre que tengas una copia de este archivo en tus carpetas con permisos de escritura. En este momento tienes un archivo fasta con la secuencia cargada en el programa, y además, su correspondiente archivo de anotaciones que describe dónde están los genes, sus nombres y propiedades, y los distintos elementos que contiene, en particular los intrones.
    1. Juega con el nivel de zoom (hay una barra con un menos y un más en la parte superior derecha) para ver la información de los genes. Es importante que tengas en cuenta el nivel de zoom, porque a elevados valores, no serás capaz de ver nada.
    2. Juega con el botón derecho del menú para ver cuántas isoformas hay descritos en estos genes (opción Expanded o Squished)
    3. Juega con las otras opciones para ver lo que hace el programa
    4. Acude a la ayuda o Guía de Usuario del programa
    5. Acude al canal de Youtube que explica muchas de las propiedades que tiene este programa
  6. Ahora ejecuta esta orden samtools index accepted_hits.bam. Esto sirve para indexar (ordenar) el archivo bam para que funcione correcta y más rapidamente dentro del programa. Se va a generar un archivo de índice con la terminación *bai. Ten en cuenta que los archivos BAM contienen una gran cantidad de información. La indexación permite su uso más rápido y racional
  7. Ahora carga el archivo accepted-hits.bam. Se cargarán a la vez la información de los alineamientos, y el archivo índice. SI no has generado un archivo índice con el mismo nombre que el bam, el programa te lo va a indicar
  8. Verás como lo acepta sin problemas, siempre que mantengas el archivo bai de índice en el mismo directorio que su correspondiente archivo bam del que procede
  9. Maneja las diferentes partes del programa y familiarizate sobre todo con los niveles de zoom hasta ser capaz de ver las lecturas
  10. Selecciona en el campo que está a la izquierda de Go la zona chr1:226040038-226045165 (aunque el propio programa le ponga comas, no hace falta ponerlas al escribirlas)
  11. Mira la información de la cobertura de la secuenciación y algunos colores en algunas columnas. Debes tener en cuenta que son muy pocas las lecturas que hemos incluido en el mapeo. Puedes optar por mapear las lecturas de las secuencias pareadas que usamos para ensamblar el genoma de E.coli y un genoma de Ecoli para ver una cobertura mayor
    1. Mira la correspondencia de las lecturas con la de los intrones indicados por el programa
    2. Fíjate en los límites de las lecturas que definen los exones/intrones
    3. Centra la página donde ves los colorines en las columnas de la cobertura ¿A qué crees que se corresponden?
    4. Juega con el nivel de zoom hasta llegar a ver las bases
    5. Juega con el menú del botón derecho del raton y los niveles de zoom para ver las variantes
    6. Mira en la parte de abajo donde se señala la posición de la secuencia indicada por el ratón. Busca la base 226040705 e indica qué es lo que ves
    7. Mira en esta región chr1:226,037,699-226,037,738 e indica qué significa los colores que aparecen en la cobertura del alineamiento
    8. Si re preguntas por qué ves lecturas en color rojo, mira esta página sobre el color dependiendo del tamaño del inserto. Si de veras quieres sacarle partido a este programa, no tendrás mas remedio que leerte la guía de usuario del programa
    9. Mira en ESTA PÁGINA para ver como se interpretan las orientaciones de las lecturas y cómo se descubren las reorganizaciones génicas, duplicaciones en tandem, translocaciones en el mismo cromosoma, etc
    10. Mira esta página para ver cómo se pueden ver las "splice junctions" a través de las Sashimi plots. Más información sobre las Sashimi plots en esta página
    11. Ya que TopHat te da los archivos BED tanto de las inserciones como de las delecciones,
  12. Busca información sobre cómo se representan estas lecturas mapeadas desde ESTA PÁGINA
  13. Busca con IGV estas zonas
    1. chr1:150226611-150247125
    2. chr1:18745 (para que veas que también se puede buscar una sola base)
    3. Busca mutaciones o fallos en la secuenciación seleccionado esta región en IGV chr1:150233921-150234000
    4. Explora los límites de las lecturas jugando con los niveles de zoom y mira si corresponden a los límites de los exones señalados por el archivo gtf (por ejemplo chr1:150233801-150234120)
    5. Carga el archivo gtf para ver cómo las lecturas se ajustan a los perfiles de los exones, y cómo muchas de las lecturas se han tenido que "romper" para alinearlas al genoma de referencia dejando un intrón en su interior
    6. Carga el archivo de junctions (junctions.bed) que TopHat te ha proporcionado, y juega con el botón derecho del ratón para verlo
  14. Busca información sobre el programa y por qué dibuja con diferentes colores las lecturas (http://software.broadinstitute.org/software/igv/DefaultDisplay)
  15. (opcional) Lanza el programa tophat con las mismas secuencias, pero esta vez senala --library-type fr-unstranded. Cárgalo todo en IGV y dime si ves algunas diferencias. Se generará un nuevo archivo bam que deberás indexar previamente con samtools para conseguirlo representar dentro de IGV

 

TRABAJANDO CON LOS ARCHIVOS DE ALINEAMIENTO O MAPEO EN FORMATO SAM/BAM

  1. Ejecuta por primera vez la orden $samtools
  2. Mira los diferentes subprogramas u opciones que tiene. En un principio te agrupa los subprogramas que tiene por funciones
  3. Entre en las diferentes opciones de las opciones invocando samtools y la correspondiente opción (por ejemplo, usando samtools view)
  4. Busca una opción con samtools que te de la estadística del archivo accepted_hits.bam que ha creado tophat (samtools flagstat accepted_hits.bam)
  5. Mira el contenido de este archivo bam con la orden samtools view (usando o no un pipe con las órdenes head, tail, more o less)
  6. Busque la orden que le permita ver la cabecera del archivo accepted_hits.bam con la orden samtools view (nota: esto es por que observe que si trabaja con un BAM, la orden samtools view no aporta información de la cabecera del archivo)
  8. Busque el modo de ver sólo una región del archivo accepted_hits.bam
  9. Mira ESTA PÁGINA y vamos a jugar un poco con la información que podemos obtener del archivo accepted_hits.bam. Por ejemplo:
    1. Estudie el formato del archivo accepted_hits.bam con la orden samtools view | less -S
    2. Haga corresponder cada campo (columna) con la información que proporciona
    3. Vamos a usar la información proporcionada en esta página o la suminsitrada por samtools flags N (donde N es el valor del FLAG) para interpretar lo que tenemos
      1. Interprete lo que significa el FLAG 99
      2. Interprete lo que significa el FLAG 147
      3. Interprete lo que significa el FLAG 577
      4. Como alternativa a la página anterior, lance la orden samtools flags para ver qué ordenes puede dar
        1. Ejemplo: samtools flags UNMAP indicará que el código hexadecimal es el 0X4 ó el decimal 4 para indicar un par concordante
        2. Mire lo que hace la opción samtools view -f y samtools view -F
          1. ejemplo: samtools view -f 4 | head (o wc -l) dará información o el número de lecturas que no mapean
          2. samtools view -F 4 dará información de todas las lecturas que mapean, porque -F (en mayúsculas) da la información contraria que -f
            1. Investigue lo que significa buscar con el FLAG 66. ?¿Cuantas lecturas hay? (samtools view -f 66 accepted_hits.bam).
            2. Extraiga estas lecturas en un archivo para poder verlas representadas luego en IGV
            3. Mire lo que significa el flag 6 y el flag 7.
            4. Busque con -F 6 las lecturas que son mapeadas y pareadas
            5. Busque con -F 7 las lecturas que son mapeadas, pareadas y concordantes o correctamente mapeadas
            6. Busque las lecturas que son mapeadas, pareadas y no concordantes con la combinación -F 6 y -f 1
            7. Haga sumas y restas con las opciones 4, 5 y 6 a ver si salen las cuentas
    4. Vamos a usar la información del campo RNAME para ver qué lecturas no han sido mapeadas. Buscaremos en el archivo de especificaciones y aplicaremos una orden samtools view adecuada
    5. Vamos a usar la información del campo CIGAR
      1. Para ver cuantas lecturas contienen intrones (col 6, CIGAR, buscando la letra N). Use para ello la orden cut para aislar la sexta columna, y la orden grep para encontrar lineas con N
      2. Para ver la longitud de los intrones
      3. para ver cuantas lecturas tienen inserciones relativas al genoma de referencia
      4. para ver cuál es la longitud de esas inserciones
    6. Mire el significado del valor RNAME, y busque cuantas lecturas han quedado sin mapear
  10. Ordena por nombres o por coordenadas el archivo accepted_hits.bam con la orden adecuada de samtools.
    1. Si lo ordenas por coordenadas, lo puedes usar directamente en un programa explorador de genomas como IGV o IGB
    2. Si lo ordenas por nombres, lo puedes usar con programas como HTSeq-count para obtener una tabla que contiene los nombres de los genes y las cuentas, listo para ser usado en RNA-Seq
  11. Vamos a preparar el archivo BAM con samtools par poder usarlo con otros programas o para ver información del mapeo (una ayuda la puedes obrtener desde los propios comandos, o desde esta página)
    1. Indexa el archivo accepted_hits.bam con la orden samtools index accepted_hits.bam. Verás que se crea el archivo de índice .bai. Como toda indexación, esto permite acelarar y optimizar el uso de estos archivos
    2. Indexa el archivo chr1.fasta con la orden samtools faidx chr1.fa. Como veis, también se puede indexar los archivos fasta
    3. Mira la cabecera del archivo con esta orden samtools view -H accepted_hits.bam. Deduzca qué información obtiene de dicha cabecera: Observará los cromosomas que contiene, así como información de si el archivo está ordenado y el criterio de ordenación, que puede ser por coordenadas o por alfabeticamente por el nombre del gen
    4. Mira una parte del archivo sam/bam con la orden samtools view accepted_hits.bam chr1:150226611-150247125
    5. Genera un nuevo archivo que contenga solo la parte anterior usando un ">" y el un nuevo archivo con un nuevo nombre de archivo
    6. Mira la ayuda de samtools tview (explora el modo de uso)
      1. Visualiza los alineamientos en formato texto en una posición determinada con samtools tview -d T -p chr1:150233900-150234034 accepted_hits.bam chr1.fa (corresponde a una región donde hay lecturas)
      2. Visualiza los alineamientos en un modo interactivo con samtools tview -d C -p chr1:150233900-150234034 accepted_hits.bam chr1.fa (corresponde a una región donde hay lecturas)
    7. Mira los datos en formato mpileup (apilados) con la orden samtools mpileup -f chr1.fa accepted_hits.bam (requiere que el archivo fasta esté indexado con samtools faidx)
    8. Genera un archivo de variantes en formato VCF no comprimido con la orden samtools mpileup -v -u -f chr1.fa accepted_hits.bam > variantes.vcf
      1. Aprovecha que has generado el archivo de variantes, para cargarlo en IGV siguiendo las instrucciones indicadas en esta página
      2. Podrias investigar cómo se puede hacer un diagnóstico con los archivos vcf. Busca y mira páginas que usen arhivos vcf de variantes
    9. Genera el mismo archivo de variantes en formato BCF comprimido con la orden samtools mpileup -g -f chr1.fa accepted_hits.bam > variantes.bcf
    10. Elimina los duplicados de PCR con la orden samtools rmdup accepted_hits.bam > accepted_hits_noPCR.sam
    11. Convierte el archivo archivo accepted_hits.bam al formato de texto normal SAM para ser capaz de ver y leer su contenido. Busca en samtools como hacerlo.
    12. Regenere desde un archivo BAM todos las lecturas que han sido mapeadas en un archivo fastq. Búsca cómo se hace.
    13. Regenere desde un archivo BAM todas las lecturas a un archivo fasta. Búsca cómo se hace.
    14. Mira que orden hay que dar para concatenar o unir dos o más archivos BAM en un único archivo
  12. Obtén solo la cabecera del archivo BAM/SAM e indica qué información proporciona. Deberás buscar en la ayuda de samtools view cómo hacerlo.
  13. Ahora haremos un análisis control de calidad del archivo bam de mapeo con samstat. Busca este programa con Google e invocalo directamente con samstat y el nombre del archivo bam. Verás que se genera un archivo html. Abrelo y analizalo.
  14. Si quieres sacarle más partido a los archivos BAM/SAM mira las utilidades del conjunto de programas Piccard
  15. Y ahora jugaremos con Rstudio (que está instalado en la UCO con ciertas limitaciones) y esta página WEB que te da información sobre R
  16. Y hágamos ahora una extracción de cuentas con el programa htseq-count htseq-count -f bam -s yes -m intersection-strict accepted_hits.bam gencode_v18_chr1.gtf > Contaje_genes.txt. Deberás ver la ayuda del programa para entender que hacen estos cualificadores (-f, -s, -m...)

 

Tareas y prácticas opcionales y más información

  1. Las especificaciones oficiales de un archivo SAM (en formato pdf)
  2. Puedes ver qué significa los valores de FLAG desde todos estos enlaces
    1. desde esta página (solo funciona en decimal)
    2. Desde esta otra página en la que ves los flags en binario, decimal y hexadecimal
    3. Desde esta otra (muy recomendada) para ver algunos valores de FLAGS más característicos (en formato pdf)
    4. O en esta slideshare (a partir de la diapositiva 32). En esta diapositiva se explica también lo que es un SAM/BAM
  3. Puedes descargar una publicación PDF que te indica todo lo relacionado con archivos SAM/BAM
  4. Aquí puedes ver información sobre el pileup incluido en SAM
  5. Mira esta página para obtener un pequeno tutorial de cómo usar samtools
  6. Mira esta página para ver cómo se puede utilizar samtools para extraer información sobre el mapeo
  7. Otro tutorial (magnifico) de cómo usar samtools
  8. Cómo puedes obtener información de las variantes usando samtools para generar un archivo VCF
  9. Mira esta página si quieres ver un tutorial sobre IGV
  10. Un PDF que da indicaciones de cómo hacer un experimento de expresión diferencial con R y el paquete DESeq2
  11. La ayuda del propio paquete DESeq2 ya da las pautas para aprender a realizar un estudio de expresión diferencial
  12. Este otro tutorial ya hace algo más que un análisis de expresión diferencial: realiza enriquecimientos funcionales para interpretar lo que se está expresando diferencialmente
  13. Otro magnifico tutorial que obtiene datos de lecturas desde la base de datos NCBI GEO
  14. y este último como se integra RStudio para hacer otro tutorial

Significado de los flags

BitHx BitDec Info  ------ --------------------------------------------------------------------------  
0x0001	1 the read is paired in sequencing, no matter whether it is mapped in a pair  
0x0002	2 the read is mapped in a proper pair (depends on the protocol, normally inferred during alignment) 1  
0x0004	4 the query sequence itself is unmapped  
0x0008	8 the mate is unmapped 1  
0x0010	16 strand of the query (0 for forward; 1 for reverse strand)  
0x0020	32 strand of the mate 1  
0x0040	64 the read is the first read in a pair (see below)  
0x0080	128 the read is the second read in a pair (see below)  
0x0100	256 the alignment is not primary (a read having split hits may have multiple primary alignment records)  
0x0200	512 the read fails platform/vendor quality checks  
0x0400	1024 the read is either a PCR duplicate or an optical duplicate