Práctica 1

Recap

Obetivos:

  • Escribir un script que lea una secuencia de nucleótidos en memoria para calcular distintas frecuencias
  • Entender el script base propuesto que calcula la frecuencia para un nucleótido
  • Resolver las cuestiones planteadas haciendo uso de las pistas

A continuación tenéis el script base (también disponible en /home/biocomp/practica1/script_base.py).

# asignar la ruta del fichero de entrada
# (es decir: DONDE esta el fichero que queremos leer)

input_path = "/home/biocomp/practica1/secuencia.fa"

# abrimos el fichero correspondiente en modo lectura

with open(input_path,"r") as ifile:
    # guardamos las líneas en la lista lines
    lines = ifile.readlines()

# inicializamos una variable vacía
sequence = ""
# iteramos todas las líneas del fichero
for l in lines:
    if not l.startswith(">"): # comprobamos si la línea comienza con > (formato fasta)
        sequence = sequence + l.rstrip() # rstrip elimina el salto de línea al final (si lo hay)

# ahora ya tenemos la secuencia completa cargada en sequence
# usamos el método count para contar cuantas veces aparece una letra concreta

C_occurrences = sequence.count("C")
G_occurrences = sequence.count("G")

print("Number of G: ",G_occurrences)
print("Number of C: ",C_occurrences)

Preguntas/Problemas

  1. Modificar el script para obtener la ruta del fichero de entrada del usuario (Dificultad: Fácil)
  2. Determinar la longitud de la secuencia y sacar el resultado en pantalla (Dificultad: Fácil)
  3. Calcular el %G+C de la secuencia (el operador para la multiplicación es *, mientras el de la división es /) y escribir el resultado en pantalla (Dificultad: Fácil)
  4. Generalizar el contenido en bases (poder calcular %A+T, o %C+T, etc) (Dificultad: Media)
  5. Calcular el número y porcentaje de todos los n-meros en la secuencia (para un n definido). Por ejemplo, para n=2 calcula dímeros para n=3 calcula trímeros, etc. (Dificultad: Difícil)
  6. Calcular la frecuencia esperada para cada n-mero y la ratio entre las frecuencias observadas y esperadas
  7. Incluir un parametro de paso (step) para poder calcular los n-meros de forma solapantes o no-solapantes
  8. Usar este fichero (codon_AA) para generar una estatística del uso de codones. Escribir un fichero de salida con las siguientes columnas:
    • codon
    • aminoácido
    • número de codones (las veces que se usa el codon en la(s) secuencia(s) de entrada
    • frecuencia relativa dentro de los codones degenerados.
    • frecuencia del codon (relativa al número total de codones)
  9. Plotear la ratio obs/esp para los n-meros

 

Funciones y pistas que pueden ayudar para resolver la cuestión 4 y 5

Cuestión 4:

## definir los nucleótidos que queremos analizar
nucleotides = ["A", "T"]
# de esta forma se puede calcular la frecuencia de cada uno dentro de un bucle for 
# luego se suma todo y se divide por el número de bases
n_counts = 0

for nt in nucleotides:

# haz los cálculos aquí

n_counts = n_counts  + numero_calculado

# % = n_counts/total de nts


Cuestión 5


# para resolver este problema hay dos claves: usar diccionarios e iterar todas 


# las subsecuencias presentes para un n determinado 


# por ejemplo para trinucleótidos podría empezar así


subseq_dict = {} #inicializamos un diccionario vacio, las keys son los trinuc y el valor su frecuencia


for i,s in enumerate(sequence):


current_sub = sequence[i:i+3]


#seleccionamos los nucleotidos entre la posicion que iteramos y 3 mas adelante


if curren_sub in subseq_dict.keys():


# suma uno al valor que tuviera esa key


else:


subseq_dict[current_sub] = 1


# inicializamos esta key con valor 1 # para completar el ejercicio debeis generalizarlo (i.e. que se pueda usar cualquier N di-,tri-nucleotidos etc # en el bucle for (tal cual esta), os dara un error porque las ultimas iteraciones superan la longitud de la secuencia # si una secuencia tiene longitud=4, intentar


acceder a la posicion 5 o mas provoca un fallo


cuestión 6


la frecuencia esperada la podemos obtener a partir de las frecuencias observadas de los mononucleótidos


cuestión 7


import matplotlib.pyplot as plt

plt.bar(xAxis,yAxis)

plt.title(‘OBS/EXP ratio of ‘+str(nmerLength)+”-mers”)

plt.xlabel(‘nmer’)

plt.ylabel(‘obs/exp’)

plt.show()