Home » bondad predicción de genes

bondad predicción de genes

Contenidos

El objetivo de esta práctica será el de cuantificar la bondad de la predicción del programa GeneMark.

Para guiarnos, nos fijaremos en la figura 1 de la siguiente publicación: Burset M1, Guigó R. (1996) Evaluation of gene structure prediction programs.

guigo_1996

Protocolo:

1) Descargar la secuencia del genoma completo de Rickettsia massiliae que se usará para la predicción de genes. Cambiar el nombre de la secuencia a ‘CP000683’ y eliminar la descripción

2) usar el ORFfinder para detectar ORFs sobre una subsecuencia

3) Descargar las secuencias de las CDS de los genes anotados de Rickettsia massiliae. Estos datos se usarán para obtener la anotación de la secuencia con la que tenemos que comparar la predicción.

4) Predecir los genes mediante el GeneMark y guardar el resultado en formato ‘GFF’. Eliminar la última línea en blanco.

5) Eliminar el encabezado (las primeras 9 líneas) del fichero en GFF (resultado).

6) Obtener la anotación

El fichero con las secuencias de las CDS contiene para cada secuencia una anotación en formato fasta.
>lcl|CP000683.1_cds_RMA_0001_1 [gene=RMA_0001] [protein=hypothetical protein] [protein_id=ABV84350.1] [location=1..822]

El objetivo será extraer las coordenadas dadas por el tag location=1..822. Mediante una línea de comando:

cat Rickettsia_cds.fa | awk ‘{ if($1 ~ />/){ if($0 ~ /location=complement/){ split($0,f,”location=complement\\(“); split(f[2],f2,”\\.\\.”); gsub(“]”, “”,f2[2]); gsub(“\\)”, “”,f2[2]); print “CP000683\t”f2[1]-1″\t”f2[2]”\t-\t0\t-“; } else{ split($0,f,”location=”); split(f[2],f2,”\\.\\.”); gsub(“]”, “”,f2[2]); print “CP000683\t”f2[1]”\t”f2[2]”\t-\t0\t+”;} }} ‘ > rickettsia_cds_coord.txt

obtendremos una salida que consta de las coordenadas en formato bed

–desglosando el comando 
escribir los encabezados
cat sequence-2.txt | awk ‘{ if($0 ~ />/) { print $0 } }’

escribir los encabezados de las CDS de genes localizados en la hebra (-)
cat sequence-2.txt | awk ‘{ if($0 ~ />/) { if($0 ~ /location=complement/ ){ print $0 } }}’

cortar la cadena de caracteres en ‘location=complement(‘; Las coordenadas estaran en el indice 2 del array f (f[2])
cat sequence-2.txt | awk ‘{ if($0 ~ />/) { if($0 ~ /location=complement/ ){ split($0,f,”location=complement\\(“); print f[2] } }}’

cortar las coordenadas en ‘..’;
cat sequence-2.txt | awk ‘{ if($0 ~ />/) { if($0 ~ /location=complement/ ){ split($0,f,”location=complement\\(“); split(f[2],f1,”\\.\\.”) ; print f1[1]” “f1[2] ;} }}’

eliminar ‘([‘ y escribir las coordenadas en formato BED
cat sequence-2.txt | awk ‘{ if($0 ~ />/) { if($0 ~ /location=complement/ ){ split($0,f,”location=complement\\(“); split(f[2],f1,”\\.\\.”) ; gsub(“\\)]”,””,f1[2]); print “CP000683\t”f1[1]-1″\t”f1[2]”\t-\t0\t-” ;} }}’

redireccionar la salida a un fichero
cat sequence-2.txt | awk ‘{ if($0 ~ />/) { if($0 ~ /location=complement/ ){ split($0,f,”location=complement\\(“); split(f[2],f1,”\\.\\.”) ; gsub(“\\)]”,””,f1[2]); print “CP000683\t”f1[1]-1″\t”f1[2]”\t-\t0\t-” ;} }}’ > annotation_-.bed
ahora para la hebra (+)
cat sequence-2.txt | awk ‘{ if($0 ~ />/) { if($0 ~ /location=[0-9]/ ){ split($0,f,”location=”); print f[2] } }}’

cat sequence-2.txt | awk ‘{ if($0 ~ />/) { if($0 ~ /location=[0-9]/ ){ split($0,f,”location=”); split(f[2],f1,”\\.\\.”) ; gsub(“]”,””,f1[2]); print “CP000683\t”f1[1]-1″\t”f1[2]”\t-\t0\t+” ; } }}’ > annotation_+.txt

 

unir los dos ficheros

cat annotation_+.txt annotation_-.txt > annotation.txt

 

  • Expresiones condicionales en awk

4 Awk If Statement Examples ( if, if else, if else if, 😕 )

  • funciones tipo ‘String’ (cadena de caracteres)

https://www.gnu.org/software/gawk/manual/html_node/String-Functions.html

  • arrays en awk

http://www.delorie.com/gnu/docs/gawk/gawk_119.html

 

7) Usar Galaxy para comparar la predicción con la anotación (mirrors)

  • Subir los datos a Galaxy: Get Data –> Upload File, especificando directamente los formatos BED (para la anotación) y GFF para la predicción.
  • Convertir GFF en BED: Herramienta lápiz –> Convert Format
  • Filtrar por hebra: Filter and Sort –> Filter (c6 == ‘+’ and c6 ==’-‘)
  • Comparar las listas: Operate on Genomic Intervals –> Coverage

 

¿Cómo podemos calcular la sensibilidad y especificidad?

Anuncios

Fecha de entrega del trabajo personal extendido al 5 de marzo

Comienzo de las clases: 4 de febrero, 12h

Aula-Inf O01 (Decanato)

—————————

Para un mayor aprovechamiento de la asignatura, se recomienda acudir a clase con ordenador portatil o tableta.