miércoles, 31 de diciembre de 2014

Ensamblado de un genoma de novo

En esta entrada llevaremos a cabo el ensamblado de un genoma de novo y comparemos posteriormente el resultado con un genoma de referencia. Para realizar el ensamblado haremos uso de la herramienta Velvet, que se basa en las gráficas de Brujin y en una serie de algoritmos para realizar el ensamblado. Posteriormente, la comparación con el genoma de referencia se realizará con Mauve, una herramienta de alineamiento múltiple de genomas desarrollada por la Universidad de Wisconsin-Madison.

Los pasos que hemos seguido para llevar a cabo esta práctica quedan detallados en la página web de la asignatura, al que podéis acceder en este enlace.


Primera parte: Ensamblado con Velvet

Para realizar el ensamblado nos hemos basado en unas secuencias pareadas, obtenidas con un equipo Illumina HiSeq2000 a partir de una estirpe de E. Coli TY-2482.

Tras descargar las secuencias, el primer paso es realizar un análisis preliminar de las secuencias para comprobar la calidad de las mismas. Como ya se ha explicado en entradas anteriores, el control de la calidad de la secuenciación podrá realizarse con el programa FastQC.

Hemos comprobado que la calidad de las secuencias es óptima. Podemos empezar entonces el ensamblado. Velvet es una de tantas herramientas del ámbito de las Biociencias que sólo es accesible en Linux, por lo que tendremos que trabajar en este sistema operativo.
Tras descargar el software y compilarlo (los detalles de este proceso los podéis encontrar en el enlace a la práctica de la página web de la asignatura) tenemos los dos comandos que necesitaremos para el ensamblado: Velveth y Velvetg. El primero de ellos se encarga de la construcción del perfil de datos que empleará posteriormente Velvetg, incluyendo archivos con información sobre cada una de las secuencias. Velvetg se encarga de construir las gráficas de Brujin a partir del k-mer especificado anteriormente a Velveth. 

Como acabamos de mencionar, el ensamblador Velvet hace uso de las gráficas de Brujin. Este método surge para superar los inconvenientes que emplea el método basado en los solapamientos, ya que realiza una fragmentación más amplia del genoma y tiene en cuenta todas las posibilidades de ensamblaje de los nodos.  Es decir, a diferencia del método de solapamiento, las gráficas de Brujin nos proporcionan muchas más opciones de posibles secuencias y de todas las propuestas selecciona aquella que es más larga y coherente con la hipótesis. 

Sin embargo, a pesar del claro avance que supone, no está libre de ambigüedades, puesto que el ensamblaje de nodos puede resolverse por distintos caminos. Por ello, se aconseja realizar varios ensamblajes con distintos valores de k. En nuestro caso particular, hemos partido de los siguientes valores de k: 27, 29, 33, 35 y 37. 
En este punto, es importante destacar que programas como Velvet únicamente trabajan bien con valores impares de k mer, puesto que de esta forma se evita que un fragmento sea el reverso complementario de sí mismo y porque entre secuencias palindrómicas, el nucleótido central es distinto.

Actualmente se han descrito alrededor de diez parámetros que evalúan la calidad del ensamblado. Uno de los más sonados es el valor de N50, que se define como "el tamaño del último contig que, sumado a las respectivas longitudes de los contigs de tamaño superior, se obtiene una longitud igual o superior a la mitad del genoma en estudio". Así pues, podemos suponer que cuanto mayor es el N50 de nuestro ensamblado mayor es la calidad del mismo, puesto que mayor es la longitud de los contigs y menor número de huecos quedará. 

Sin embargo, esta correlación no es siempre acertada ya que puede tratarse de que estemos realizando un ensamblado incorrecto y los fragmentos que se están uniendo no se corresponda con el orden real del genoma que estamos analizando. Esto lo corroboraremos comparando nuestro ensamblado con un genoma de confianza (para lo que usaremos Mauve). 

Así pues, hemos obtenido los siguientes datos tras lanzar Velvet con nuestra secuenciación: 





Segunda parte: Ordenando los Scaffolds con Mauve

Mauve es un sistema ideado para construir eficientemente múltiples alineamientos de un genoma, teniendo en cuenta eventos evolutivos a gran escala como reordenaciones o inversiones. En esta práctica nos centraremos en la interfaz básica del programa. 

En primer lugar y siguiendo los pasos descritos en la página web de la asignatura, introducimos el genoma de referencia y seguidamente el resultado de el ensamblado realizado por Velvet, ambos en extensión ".fa". El programa irá realizando sucesivos alineamientos, de los cuales el último es el resultado final y el alineamiento válido. 

El alineamiento entre las secuencias seleccionadas se nos muestra como un panel horizontal junto con una escala que coordenadas del genoma. Cada bloque de colores representa una región de la secuencia del nuestro genoma alineado con parte del otro genoma, siendo presumiblemente homóloga al mismo y libre de reordenamiento. 

Cuando uno de los bloques de nuestro genoma problema se localiza por encima de la línea central del esquema, indica que  se orienta hacia delante con respecto a la secuencia del genoma de referencia. En contraposición, los bloques que se localizan por debajo de la línea central representan regiones que se localizan en orientación inversa con respecto a la secuencia de referencia. Por último, las regiones que no forman parte de ningún bloque constituyen secuencias que no presentan homología con el genoma de referencia. 

Además, dentro de cada bloque de color, Mauve dibuja un perfil de similitud de la secuencia del genoma, donde cuanto mayor es la altura del perfil mayor es la conservación en esa región. 

A continuación se expondrán los cinco alineamientos realizados en función de los contigs obtenidos con las k mer indicadas en la primera parte, por orden desde k= 27 hasta k=37: 











De todos los alineamientos realizados podemos observar como claramente el que más homología comparte con el genoma de referencia es el alineamiento k = 27, en contraposición con lo que auguraba el N50 de los respectivos ensamblados realizados con velvet. Esto nos demuestra cómo dicho parámetro debe ser secundario a la hora de establecer el alineamiento más óptimo pues, como en este caso, se puede tratar de un alineamiento "ficticio" en el que los distintos segmentos se han ensamblado erróneamente. 




jueves, 4 de diciembre de 2014

Control de calidad de secuenciación con SOLiD (II): tras el cribado

En este caso llevamos a cabo tres cribados en el control de la secuenciación y un filtrado, para la que utilizamos las herramientas de trimming y filtering de FastX-Toolkit, respectivamente: 
  1. Selección de secuencias con una calidad igual o superior a Q= 20.
  2. Selección de secuencias con una calidad igual o superior a Q= 28.
  3. Selección de secuencias con una calidad igual o superior a Q= 28 pero eliminando las lecturas de menos de 47 bases. 
  4. Filtrado de lecturas que tengan menos del 90% de bases con una calidad superior a Q= 20. 
A continuación se muestra la evolución en la aceptabilidad en la relación de módulos proporcionada por FastQC y se comentarán los aspectos más interesantes. 




1.  Selección de secuencias con una calidad igual o superior a Q = 20

En la primera pestaña de este análisis comprobamos que se han tenido que eliminar 9.233 para conseguir las condiciones de calidad establecidas. Las secuencias que se han mantenido presentan distintas longitudes (1-50) como veremos más adelante. 

El módulo Per base sequence content observamos notables diferencias con respecto al respectivo análisis en la secuencia original, de forma que al eliminar las secuencias "erróneas" se ha conseguido disminuir la diferencia en el contenido de A y T o G y C (ahora las diferencias son superiores al 10% pero inferiores al 20%). 

Otro aspecto que se ha visto modificado significativamente es el contenido en GC de la secuencia, que se ve claramente disminuido, alejándose más de lo que marca la secuencia de referencia. Una posible explicación puede ser a la aparición en este análisis de secuencias de distintos tamaños lo que haga que no coincida el contenido GC de estas nuevas lecturas con el de referencia, en las mismas posiciones. Cabe destacar además que el porcentaje del contenido en GC con respecto al total de secuencias no se ve alterado con respecto al análisis original (en ambos casos es de 55%). 




Por último, merece la pena destacar el módulo Sequence length distribution, en cuyo gráfico podemos observar (como ya comentamos anteriormente) la distribución diferencial de la longitud de las lecturas entre la secuencia original (izquierda) y tras el cribado (derecha):




2.  Selección de secuencias con una calidad igual o superior a Q = 28

En este análisis se han visto significativamente modificados, con respecto al cribado anterior, los parámetros: per sequence quality scores, overepresented sequences  y Kmer content

En cuando a los valores de calidad de secuencia, éstos se ven muy superados con respecto tanto al análisis anterior como a la secuencia original. Con estas nuevas condiciones de Q, se consigue que el valor de calidad media más frecuente sea superior a 27.



Curiosamente, las secuencias sobrerepresentadas se ven aumentadas. Una posible explicación es que la eliminación de las secuencias que no alcanzaban una Q= 28 haya dado lugar a una modificación en la relación secuencias sobrerepresentadas/secuencias totales de tal forma que se haya alcanzado el límite para el que el programa lanza el aviso (warning) en este parámetro. 

Finalmente se produce un importante aumento de las secuencias Kmer, lo que se corresponde también con el aumento del porcentaje de secuencias sobrerepresentadas y que podrían también deberse a un aumento en la relación de posible ADN contaminante o de preferencia en la PCR frente al conjunto de secuencias totales. Es decir, que secuencias ajenas a nuestro ADN de interés presentaran la suficiente calidad en la secuenciación para no ser eliminadas del procesamiento, por lo que su representación respecto al total de lecturas se haya visto aumentada. 


3.  Selección de secuencias con una calidad igual o superior a Q = 28 y de longitud superior a 47 bases 

En este cribado cabe mencionar la notable disminución en el total de secuencias (ahora un total de 44.119) que ha obligado las nuevas condiciones de cribado establecidas. 
Con respecto al análisis anterior, los apartados que se han visto modificados -ambos favorablemente- ha sido per sequence GC content y overepresented sequences, por lo que podemos suponer que se debe a la eliminación de secuencias de longitud inferior a 47 bases, de ADN contaminante o preferentes para la PCR, que provocaban las alteraciones comentadas.




4.  Filtrado de lecturas que tengan menos del 90% de bases con una calidad superior a Q= 20. 

Podemos comprobar como este análisis es el que mejores parámetros (más acercados a una correcta secuenciación) proporciona (menor número de módulos marcados como warning o failure). Es especialmente llamativo el gráfico de contenido en GC, que se asemeja increíblemente con el propio de la secuencia de referencia: 




Sin embargo, es necesario remarcar que el total de secuencias procesadas ha disminuido hasta 25.686 (casi una quinta parte del total de lecturas en la secuencia original), debido a las estrictas condiciones de calidad que hemos establecido. 

Esto nos lleva a una conclusión general en cuanto a los controles de calidad de la secuenciación: a la hora de realizar un control de calidad de secuenciación tenemos que llegar a un compromiso entre las condiciones de calidad establecidas y el número de lecturas a ser procesadas, puesto que una notable disminución del total de secuencias analizadas conducirá a un control de calidad poco significativo o representativo de la secuenciación problema origina. 

miércoles, 3 de diciembre de 2014

Control de calidad de secuenciación con SOLiD

Como ya hicimos anteriormente para Illumina, comenzaremos con la descripción de los parámetros que proporciona FastQC a partir del análisis de una secuenciación SOLiD (a la que podéis acceder clickando en este enlace). El análisis de la secuencia original (previa a la aplicación de los cribados) nos dejó la siguiente relación de parámetros.




1.  Basic Statistics

De este apartado, los aspectos que nos interesan destacar son el tipo de archivo y la codificación (al igual que con la secuencia de Illumina, es una archivo de base calls y de código Illumina 1.5), el total de secuencias y su longitud (125.000 secuencias de 50 bases) y el contenido en GC (55%). 

2.  Per base sequence quality

Este módulo viene marcado como un failure, lógico ya que si atendemos al gráfico adjunto, una gran cantidad de lecturas presentan niveles de calidad realmente bajos (en contadas ocasiones superan el umbral de Q= 20). 




3.  Per sequence quality scores

Podemos observar como, a diferencia de lo que encontramos en el análisis inicial de Illumina, en este caso existe un subconjunto de alrededor de 4.500 lecturas que presentan bajos niveles de calidad (valor de 5 en una escala de 1-33).

En nuestro caso, el análisis está marcado como warning, debido a que, como podemos observar en el gráfico, el valor de calidad media más frecuente se encuentra por debajo de 27.




4.  Per base sequence content

Este módulo viene acompañado del símbolo de failure, lo que como ya hemos visto en post anteriores, indica que las diferencias en contenido entre A y T o entre G y C son superiores al 20%. Los principales motivos suelen ser secuencias sobrerepresentadas (adaptadores o secuencias de ARN), fragmentación y/o composición de la librería sesgada, etc. 




5. Per sequence GC content

Podemos comprobar como el contenido en GC de nuestra secuenciación se corresponde casi perfectamente con la distribución de referencia. 




6. Per base N content

En este caso, el análisis viene marcado como failure, lo que se debe a que en alguna posición el contenido de Ns proporcionadas por el método de secuenciación es superior al 10%. Podemos comprobar, observando el gráfico, como únicamente se cumple para la posición 48 de las lecturas, por lo que podemos atribuirlo a un error del método de secuenciación (error físico).




7. Kmer content

Este módulo viene acompañado del símbolo warning. Cabe destacar en este caso cómo la mayor parte de secuencias kmer se localizan al comienzo de las lecturas, lo que se debe probablemente -como ya se explicó en pos anteriores- al empleo de cebadores aleatorios. 




El resto de parámetros (distribución de la longitud de las secuencias, nivel de duplicaciones, secuencias sobrerepresentadas y contenido de adaptadores) se corresponde con una correcta secuenciación por lo que no se incidirá sobre ellos en este post. 

martes, 2 de diciembre de 2014

Control de calidad de secuenciación con Illumina (II): tras el cribado

Un vez hemos entendido y descrito los diferentes parámetros que nos proporciona FastQC en  relación a una secuenciación problema, podemos utilizar FastX-Toolkit que, como avanzamos en el anterior post, nos permitirá procesar el archivo de secuenciación antes de lanzar un FastQC. En este caso nos centraremos en la herramienta de trimming. En total hemos realizado cuatro cribados, a saber: 
  1. Selección de secuencias con una calidad igual o superior a Q= 20. 
  2. Selección de secuencias con una calidad igual o superior a Q= 28.
  3. Selección de secuencias con una calidad igual o superior a Q= 28 pero eliminando las lecturas de menos de 30 bases. 
  4. Selección de secuencias con una calidad igual o superior a Q= 28 pero eliminando las lecturas de menos de 35 bases. 
A continuación, se muestra un esquema con la evolución en la "aceptabilidad" de los distintos análisis en función del cribado realizado, y se comentará para cada uno de ellos los aspectos más significativos. 





1.  Selección de secuencias con una calidad igual o superior a Q = 20

Tal como mencionamos en el anterior post, se produce una progresiva caída en la calidad de la secuenciación a medida que se avanza en las posiciones de las lecturas. Por ello, con este cribado se pretende eliminar aquellas lecturas que presentan una calidad menor de Q= 20. 

Así, si atendemos al apartado basic statistics comprobaremos que el total de secuencias procesadas se ha reducido de las 25.000 iniciales a 24.890, por lo que se han eliminado un total de 110 secuencias que no cumplían con los requisitos establecidos de calidad. Podemos advertir este aumento de la calidad global (superior a Q= 20 en todas las posiciones) en el apartado Per base sequency quality:




También se ha visto modificado el módulo Sequence length distribution que aparece ahora con el símbolo warning, debido a que no todas las secuencias son de la misma longitud. Esto se debe precisamente a lo que ya comentamos en el post anterior: el programa ha acortado la longitud de las secuencias con baja calidad, eliminando las últimas posiciones y ha conservado aquellas que tras el procesamiento superaban el valor mínimo de Q requerido. 




2.  Selección de secuencias con una calidad igual o superior a Q = 28.

Tras este procesamiento (eliminación de secuencias con una calidad inferior a Q= 20), se registraron un total de 24.865 (se eliminaron 135 lecturas). En cuanto al resto de parámetros no se producen cambios significativos respecto al primer cribado. 


3.  Selección de secuencias con una calidad igual o superior a Q = 28 y una longitud igual o superior a 30 pb. 

Estableciendo estas condiciones sí encontramos algunas diferencias importantes con respecto al análisis de la secuencia original y a los dos cribados anteriores, en concreto, un aumento en la aceptabilidad de los módulos contenido en GC  y contenido de Kmer.

En cuanto al contenido en GC observamos que las diferencias entre nuestra secuenciación y la de referencia del organismos están menos marcadas que en los anteriores análisis. 




Para el caso del contenido de Kmer, el símbolo de warning nos informa de que ha tenido lugar un aumento del p-value del test binomial respecto a los anteriores análisis, pero que éste es aún menor del 0.01.

Una de las causas principales de la aparición de Kmer al inicio de la biblioteca se debe al empleo de primers aleatorios, lo que deriva de un muestreo previo incompleto de los posibles cebadores aleatorios. 




4.  Selección de secuencias con una calidad igual o superior a Q = 28 y una longitud igual o superior a 35 pb. 

Por último, tras el procesamiento realizado bajo estas condiciones se produjo una pérdida de la mejora obtenida para los parámetros contenido en GC  y contenido de Kmer que se obtuvo en el análisis anterior, donde se permitían secuencias de una longitud a partir de 30 pb. 

viernes, 28 de noviembre de 2014

Control de calidad de secuenciación con Illumina

¡Bienvenidos de nuevo!

Volvemos después de este largo parón para iniciarnos en  la asignatura de Biología de Sistemas.

El primer problema con el que nos encontramos es que muchos de los software que vamos a utilizar únicamente están disponibles para Linux. La opción más recomendable para los usuarios de Mac o Windows, es instalar VirtualBox, un software que me permitirá instalar un sistema operativo Linux dentro de vuestro propio sistema operativo, y así trabajar con estas herramientas. Si tenéis alguna dificultad con la instalación, aquí os dejo un tutorial que me ha sido muy útil para instalarlo. 

Dicho esto, la primera tarea que vamos a realizar será un Control de calidad de la secuenciación con FastQC, una de las principales herramientas que tenemos actualmente para realizar controles de la calidad de secuenciación. Podéis descargar dicho programa en el siguiente enlace. También descargaremos FastX-Toolkit, un software que incluye diversas herramientas para el pre-procesamiento de los archivos fasta y fastq  y que podéis descargar en este enlace.

En las próximas entradas del blog valoraremos la calidad de la secuenciación obtenida por dos plataformas distintas: Illumina (que trataremos en este y el siguiente post) y SOLiD (que trataremos en los dos siguientes).

Comenzamos pues, por Illumina. Una vez hemos descargado nuestra secuencia problema en formato '.fastq' (un formato que incluye tanto la secuencia como sus correspondientes puntuaciones de calidad, en código ASCII) nos disponemos a lanzar un FastQC, que nos proporciona la siguiente información relativa a la calidad de secuenciación de nuestra secuencia problema:


Junto a la denominación de la pestaña podemos encontrar hasta tres símbolos: el tic verde nos indica que el análisis se corresponde con el de una correcta secuenciación. Mientras que el warning (símbolo ámbar, no mostrado en la imagen anterior) es una "alarma" que indica que en ese aspecto en concreto la secuenciación no es correcta. En el extremo final tenemos el símbolo rojo de failure (mala secuenciación). hsd

1.  Basic Statistics:

En este apartado se presentan simplemente las características básicas de nuestra secuencia problema: shd

  • Nombre del archivo: nombre original de la secuencia problema. 
  • Tipo del archivo: base calls o colorspace. En este caso -base calls- indica que el código empleado para correlacionar los picos del cromatograma se compone de las distintas bases nucelotídicas del ADN. 
  • Codificación: indica el rango de caracteres ASCII empleado para especificar la calidad de la secuenciación. En este caso, se corresponde con Illumina 1.5 (desde "C=3" hasta "h=40").
  • Secuencias totales: hace referencia al total de secuencias procesadas en el análisis (25.000).
  • Secuencias filtradas: incluye el total de secuencias que han sido eliminadas atendiendo a las condiciones de cribado ("0", puesto que aún no se han establecido condiciones de trimming). 
  • Longitud de la secuencia: proporciona el rango de longitud de las secuencias (desde las más cortas a las más largas). En este caso, aparece únicamente el valor "38" lo que indica que todas las secuencias son de esa misma longitud. 
  • Contenido de GC (%): indica el contenido total de citosina y guanina en las secuencias. En este caso es de un 45%. 


2.  Per base sequence quality:

Este apartado proporciona una visión general de la de calidad de secuenciación para las distintas bases localizadas en cada una de las 38 posiciones (eje X), teniendo en cuenta todas las secuencias procesadas. El gráfico informa sobre muchos parámetros (mediana, rango intercuartil para cada posición, etc.) pero lo que más nos interesa en este caso es el valor de quality score (eje Y) relativo a cada posición. Observamos como dicho valor desciende a medida que avanzamos en las posiciones de la secuencia, algo lógico dada la naturaleza de los sistemas de secuenciación que se emplean actualmente, caracterizados por una alta eficiencia en las bases iniciales y una progresiva pérdida de ésta a medida que avanza en la secuencia. 

En nuestro caso particular, el programa nos proporciona un failure, que aparece cuando el cuartil inferior para cualquier base es menor de 5 o si el valor de la mediana para cualquier base es menor de 20. En definitiva, a partir de la posición 34 podemos afirmar que la calidad de la secuenciación es muy baja, por lo que deberíamos o suprimir las secuencias cuyas últimas bases presenten un valor de Q inferior a 20 o, más apropiado, acortar la longitud de las secuencias con baja calidad en las últimas posiciones y conservar las que, tras este procesamiento, superen un valor de Q de 20.


3.  Per tile sequence quality:

En este caso, un alto contenido en colores fríos indica una alta calidad de la secuenciación, tal como se cumple en nuestro caso particular. En caso de que encontráramos teselas de distintas tonalidades, supondríamos que esos "descensos" en la calidad de secuenciación se deben a incidencias de carácter físico en dicho proceso (burbujas en el equipo, etc.). 


4.  Per sequence quality scores:

En este apartado se incluye un gráfico que representan las puntuaciones medias de calidad (eje X) frente al total de lecturas. En definitiva, nos permite visualizar rápidamente si un conjunto de nuestras lecturas presenta de forma significativa un bajo valor de calidad. En nuestro caso, la mayor parte de las secuencias presenta un aceptable valor de calidad (alrededor de Q=30) y muy pocas secuencias se encuentran en un bajo rango de calidad, lo que se corresponde a lo deducido en el apartado anterior. 

En este análisis, dado que no existen diferencias significativas en la secuenciación en función del origen de los distintos genomas, bajos niveles de calidad se correlacionan con problemas físicos (sistemáticos) del método de secuenciación (burbujas en el equipo, final de la célula de flujo, etc.). Las secuencias con baja calidad pueden eliminarse, siempre que constituyen un bajo porcentaje del total. En caso de que dicho porcentaje sea representativo, debería repetirse de nuevo la secuenciación completa. 

5.  Per base sequence content:


Proporciona un gráfico que representa la proporción de cada una de las bases de DNA en las distintas posiciones (1-38) de todas las secuencias. Esperaríamos encontrar un porcentaje aproximadamente similar entre bases complementarias (A con T y C con G) puesto que los cebadores deberían unirse a ambas cadenas con la misma afinidad. En este análisis, marcado con el símbolo de failure la diferencia de correspondencia entre las bases complementarias en ambas cadenas es superior al 20%. 



Entre los posibles motivos más probables se encuentra la preferencia diferencial de amplificación de las dos cadenas. Asimismo, si este gráfico nos muestra un porcentaje de CG distinto al característico de nuestra muestra en particular, supondremos que se trata de una interferencia por ADN contaminante. 



6.  Per base sequence GC content:

Se trata de un análisis complementario al anterior, ya que presenta la distribución de contenido en GC a lo largo de todas las secuencias, obtenida a partir de la secuenciación experimental, frente a la misma distribución teórica.

En nuestro caso, encontramos importantes diferencias que, como se ha expuesto anteriormente, atribuimos a una contaminación de ADN con distinto contenido en GC respecto a nuestra muestra problema. De hecho, el programa FastQC nos proporciona el símbolo de failure junto al análisis, que aparece cuando la suma de las desviaciones de la distribución normal representan más del 30% del total de las lecturas. 

7.  Per base N content:

En esta pestaña se proporciona un gráfico que indica el porcentaje de "bases N" para cada posición, ya que el método de secuenciación empleado es incapaz de correlacionar el pico del cromatograma con una de las cuatro bases nucleotídicas con la suficiente confianza (elevado ruido de fondo). En nuestro caso, el secuenciador no ha incluido Ns en ningún caso, lo que indica que ha atribuido con la suficiente seguridad una base para cada posición de cada lectura de la secuencia.



8.  Sequence Length Distribution:

Esta pestaña presenta una distribución que relaciona el total de lecturas de la secuenciación con su respectiva longitud. Así pues, observamos que en nuestro caso se obtiene un único pico, ya que todas las lecturas originadas (alrededor de 25.000) son de la misma longitud (38 pb). 


9.  Sequence Duplication Levels:

En este caso, comparamos el contenido de secuencias repetidas en relación con el total de secuencias procesadas. Teóricamente, dado que en el proceso de secuenciación la rotura de la secuencia es aleatoria, no esperaríamos encontrar duplicidad en las secuencias. Es decir, aún centrándonos en una región determinada, no obtendríamos fragmentos idénticos, sino solapantes entre sí. En la práctica, un pequeño porcentaje de duplicidad suele correlacionarse con una alta cobertura en la secuenciación, mientras que un elevado porcentaje de esta duplicidad hace referencia a problemas en la secuenciación.

En el análisis de nuestra secuencia hemos obtenido de nuevo un failure, debido a que las secuencias repetidas representan más del 50% del total de lecturas. La causa más probable para este error es la inclusión de ADN contaminante tras el proceso de rotura de nuestra secuencia original, de forma que dicha secuencia contaminante es preferente para la amplificación por PCR y de esta forma da como resultado un significativo porcentaje de secuencias repetidas.



10.  Overrepresented sequences:

En este módulo se presentan un listado de todas las secuencias que representan más del 0.1% del total. En general, estas librerías de secuencias no deberían contener un gran número de secuencias sobrerepresentadas, siendo el motivo más probable de aparición la contaminación con ADN exógeno aunque también puede tratarse de una secuencia biológicamente significativa. 


En el análisis de nuestra secuencia (marcado con el símbolo de error) estas secuencias sobrerepresentadas suponen más del 1% del total. El programa hace además un alineamiento de estas secuencias en distintas bases de datos y nos proporciona el origen más probable de dichas secuencias. 

11.  Kmer content

El módulo Kmer content surge bajo la asunción de que ningún pequeño fragmento se encuentre posicionado en las distintas lecturas de forma sesgada. Por ello, este análisis mide el total de subsecuencias de k bases (7 en este caso) dentro de cada lectura y aplica un test binomial para determinar desviaciones significativas de lo que sería una cobertura uniforme para todas las posiciones. Como podemos observar en la correspondiente gráfica, estos kmer aparecen como picos discretos de enriquecimiento en determinadas posiciones.


lunes, 17 de junio de 2013

Análisis de secuencia

En esta entrada seleccionaremos y analizaremos una secuencia nucleotídica, teniendo en cuenta los parámetros y herramientas que hemos estudiado a lo largo de toda la asignatura práctica de Biosíntesis de Macromoléculas. La secuencia problema así como las cuestiones que nos proponemos determinar se encuentran accesibles en el siguiente enlace.

En primer lugar, utilizaremos el programa Chromas Lite® que nos proporciona un cromatograma de la secuencia problema, del que se adjunta un fragmento a continuación:




Posteriormente, exportamos la secuencia en formato FASTA. De esta forma, podremos trabajar con la secuencia en cualquiera de las herramientas bioinformáticas para obtener información de las macromoléculas, en general, y los ácidos nucleicos, en particular.

1. Análisis de secuencias con EMBOSS:
A continuación vamos a adentrarnos en la información que podemos obtener de la secuencia problema haciendo uso de JEMBOSS, una interfaz gráfica del paquete de programas EMBOSS (European Molecular Biology Open Software Suite). 

Valoración de CDS con Plotorf
En primer lugar, accedemos al paquete informático JEMBOSS, que está disponible online.
Buscamos en el menú de la izquierda la opción Plotorf, que nos permitirá identificar los CDS de nuestra secuencia problema. Los resultados obtenidos se muestran a continuación: 

Se nos proporcionan las seis posibles fases de lectura para la secuencia nucleotídica y, resaltado en color, los fragmentos de dicha secuencia que codifican para una proteína, de donde podemos deducir que, con mayor probabilidad, la fase de lectura para la secuencia problema es la número cuatro. 

Obtención de marcos abiertos de lectura con Sixpack
A continuación, nos disponemos a traducir la secuencia a sus seis posibles fases de lectura. Para ello, haremos uso de una herramienta también de EMBOSS: EMBOSS Sixpack. Encontraréis el correspondiente enlace para acceder a ella en el menú de la izquierda de la página principal de JEMBOSS online.

Tras lanzar el Sixpack nos aparecerá una interfaz que esquematiza las secuencias peptídicas (en código de una letra) de las seis posibles fases de lectura de nuestra secuencia nucleotídica (en la siguiente imagen se muestra un fragmento de la misma). Se adjunta además, un recuento del total de ORF's (marcos de lectura abiertos) para cada una de dichas fases. 


Puntuación de posibles regiones codificantes: Tcode
Para finalizar con EMBOSS vamos a emplear la herramienta Tcode. Haciendo uso de una serie de algoritmos, Tcode es capaz de identificar regiones codificantes de proteínas en una secuencia dada. Para acceder a esta herramienta basta con seguir el procedimiento llevado a cabo para el resto de programas de EMBOSS

Tcode hace una lista de los posibles CDS de nuestra secuencia problema (que muestra finalmente en un gráfico), dándoles una puntuación (Testcode value) que le permite clasificarlos en coding (codificante, por encima de la línea verde), non- coding (no codificante, por debajo de la línea roja) o no opinion (puede tratarse de un fragmento codificante o no codificante, entre las líneas verde y roja).


2. Comparación de secuencias: 
Anotaciones: 
Como hemos podido comprobar en entradas anteriores, una de las herramientas más útiles para la comparación de secuencias es BLAST (Basic Local Alignment Search Tool), un programa bioinformático de alineamiento de secuencias de tipo local (ADN, ARN, proteínas).

Podemos distinguir entre varias modalidades de BLAST. En este caso realizaremos BlastN (Nucleotide Blast), que compara una secuencia nucleotídica con bases de datos que contengan también secuencias nucleotídicas.


Se nos muestra, en un diagrama de barras con código de colores (representativo de la puntuación de alineamiento), las secuencias más similares como una distribución de alineamientos: cada barra representa una secuencia, ordenadas por porcentaje decreciente de similitud.


La primera de las barras (máxima similitud) se corresponde con el CDS completo de la proteína ROS1 de Arabidopsis thaliana, por lo que deducimos que nuestra secuencia problema se corresponde con el gen que codifica dicha proteína.


3. Información sobre la proteína:
Clickamos en el número de accesión correspondiente a dicha proteína  (NM_129207.4) y accedemos a una interfaz que nos proporciona mayor información (locus del gen, clasificación del organismo, artículo y publicación...), además se muestra la secuencia proteínica (en código de una letra) deducible a partir de nuestra secuencia nucleotídica problema, y que se adjunta a continuación:

Entre los aspectos más relevantes que podemos encontrar en las anotaciones, se encuentra la información relativa a la funcionalidad de la proteína. En este caso en particular, encontramos que nuestra proteína se trata de un represor transcripcional que actúa en procesos de silenciameiento génico, mediante la desmetilación de regiones del promotor diana. Interacciona físicamente con otras proteínas, como RPA2/ROR1. Se ha encontrado en mutantes de Ros1 un aumento de la metilación en varios promotores. 
Entre los loci afectados por ros1, algunos se ven afectados en la metilación de citosinas aunque la mayor parte de ellos se afectan en otros nucleótidos. 

En la interfaz de las anotaciones aparece un menú a la derecha con enlaces a varias bases de datos donde podremos completar la búsqueda y ampliar la información. Una de las opciones es Map viewer, que de una forma gráfica sitúa nuestro gen en su correspondiente cromosoma (en este caso, cromosoma 2), junto con el resto de loci de  dicho cromosoma, pudiendo saber qué genes se encuentran más cerca de ros1 y, por tanto, con una alta probabilidad segregarán juntos. 

 InterPro:
A continuación, vamos a hacer uso de Interpro (en concreto, Interproscan) para identificar dominios de interés en la proteína ROS1. 
Para ello accedemos en primer lugar a la página web de Interpro, haciendo click en el enlace, y buscamos la herramienta InterproScan y copiamos la secuencia peptídica que obtuvimos al realizar el Blast en el recuadro de inserción, clickamos en submit  y ya sólo queda esperar los resultados. 



Encontramos esquematizadas las distintas regiones con dominios, mostradas bajo un código de colores (según la base de datos de procedencia). Observamos algunos dominios HHH (varias hélices alfa) y HTH (hélice-giro-hélice) que se relacionan con una DNA glicosilasa (acción de reparación del daño en el ADN) que podemos inferir que pertenece a la familia de endonucleasas III. 
Haciendo click en los recuadros al comienzo de cada representación, obtendremos información más detallada, así como los números de accesión GO.

BlastX y BlastP
Para finalizar esta práctica resumen vamos a realizar dos útlimos Blast

- Blastx: este programa usa como entrada una secuencia nucleotídica que traduce en sus seis posibles marcos de lectura, y compara estas secuencias con proteínas recogidas en bases de de datos (suele utilizarse cuando se sospecha que la secuencia de entrada codifica para una proteína, pero no se sabe exactamente cuál es su producto). 

Los resultados que obtenemos tras lanzar el Blastx se esquematizan a continuación: 
Obtenemos que la secuencia con mayor porcentaje de similitud se corresponde con la proteína Ros1 de Arabidpsis thaliana, como era de esperar. El resto de líneas que aparacen en el Blast hacen referencia a proteínas hipotéticas con secuencia recogida en las bases de datos. 

- BlastP: este programa compara una secuencia de aminoácidos con otras secuencias, también de aminoácidos, recogidas en bases de datos, para realizar alineamientos introduciendo gaps (huecos) según matrices de sustitución BLOSUM o PAM. 

Tras realizar el BlastP obtenemos los siguientes resultados: 

En este caso, obtenemos un mayor número de coincidencias. Las dos primeras secuencias se corresponden con secuencias de la proteína ROS1 de Arabidopsis, almacenadas en distintas bases de datos. Del mismo modo que cuando llevamos a cabo el BlastX, el resto de coincidencias (con menor porcentaje de similitud) se corresponde con secuencias catalogadas como proteínas hipotéticas.