Manejo de Secuencias con python

Una vez instalados python y biopython, vamos a utilizarlos para manejar las secuencias utilizando el archivo de covid19 descargado previamente. En python, un texto es una secuencia de cadena literal (string) y se puede manipular fácilmente, entonces para comenzar, vamos utilizar algunos de los ejemplos de la introducción a linux para entender cómo hacerlos con python. Los primeros ejemplos simples no requieren de librerías adicionales, funcionan con los comando básicos de python. Para esto entonces habra python

python

Cada vez que entramos a la terminal debemos activa el entorno de conda que creamos source activate python3.9

# El comando print imprime el texto entre comillas
print("Hello world!")

# El mismo resultado se obtiene sumando las palabras y concatenando el texto
helloworld = 'Hello, ' + 'world!'
print(helloworld)

Ya que un genoma es una secuencia de texto, se puede manejar de forma similar. Para facilitar el código utilizamos la librería de Biopython que ya tiene implementadas funciones específicas para secuencias.

Una vez dentro de python (y cada vez que ingresamos) debemos volver a llamar las librería y definir las variables

import Bio
from Bio.SeqRecord import SeqRecord
from Bio import SeqIO

# define una variable para abrir el archivo fasta como input
dna_file = SeqIO.read('Sars_cov.dna.fa', "fasta")

# define una variable para la secuencia completa del covid
covid19_seq = dna_file.seq

# cuenta la longitud del genoma de covid
print(f'El genoma del virus que causa Covid-19 (SARS-CoV-2) consiste de {len(covid19_seq)} bases.')

El comando print imprime en la pantalla el texto literal (entre comillas ‘’) y la función len() con el parámetro covid19_seq que definimos anteriormente. La f dentro de print es un prefijo que formatea la cadena literal (para más información).

En este caso nos da la longitud del genoma de covid que es de 29903 (el genoma humano está en la escala de miles de millones ó GB).

La variable covid19_seq es una lista ordenada de bases (string), es decir que la letra de cada posición se puede saber por su índice

# La primera posición empieza con el índice 0
covid19_seq[0]
# Segunda posición
covid19_seq[1]
# Las primeras 50 posiciones
covid19_seq[0:50]
# ó
covid19_seq[:50]
# Y de atrás para delante
covid19_seq[-50]

Asímismo, los condicionales y bucles funcionan similar que en la terminal. Por ejemplo

## Condicionales lógicos
5 == 5
5 == 4
4 < 5

# El mismo ejemplo de la terminal de imprimir los número de 1-10 con for loop
for i in range(1,10):
    print(i)

# El mismo ejemplo de la terminal de imprimir los número de 1-10 con while loop
i = 1 
while i<10:
    print(i)
    i += 1

Ahora vamos a utilizar los bucles y condicionales para buscar el codón de inicio de transcripción para el genoma de covid19

# Define el codón de inicio
start_codon = 'AUG'

# Escanea la secuencia hasta encontrar el patrón
for i in range(len(covid19_seq)):
    if covid19_seq[i:i+3] == start_codon:
        print('The start codon starts at index', i)
        break
else:
    print('Codon not found in sequence.')

Por qué este script no encontró un codón de inicio? Porque el codón está en términos del RNA mensajero y el input que le dimos fue el del genoma. Para transcribirlo, vamos a definir otra variable

covid_mRNA = covid19_seq.transcribe()

Ejercicio Vuelva a correr el script anterior para el RNA mensajero y verifique si la posición de indicada si contiene el códon de inicio

Escribir secuencias en un archivo

Con frecuencia estamos interesados en un subconjunto de los genes o proteínas para realizar un análisis a profundidad de ellos. En este caso vamos a seleccionar las glicoproteínas que son aquellas que codifican la espícula. Esta espícula (“spike”) es la encargada de mediar la entrada del virus hacia el hospedero y tiene una tasa de mutación mas alta que los otros genes, lo que presumiblemente dió origen a cepas mas virulentas (https://en.wikipedia.org/wiki/Coronavirus_spike_protein). Para esto vamos a llamar al archivo de proteínas dentro de una variable. Recordemos que el archivo se puede llamar sin añadir la ruta completa si abrimos python dentro de la carpeta donde se encuentra el archivo. De lo contrario debemos darle la ruta

from Bio.Seq import Seq
from Bio.SeqRecord import SeqRecord

prot_file = "/Users/laura/docencia/eia/bioinfo/clase_python/Sars_cov.prot.fa"

for record in SeqIO.parse(prot_file, "fasta"):
    if (record.description.split()[9] == "glycoprotein"):
        print(record)

Cuantas secuencias contiene el archivo? Vamos a imprimir su id, su tamaño y su descripción:

for record in SeqIO.parse(prot_file, "fasta"):
    print(record.id,len(record.seq),record.description)

La descripción tiene mucha información importante pero esta no es muy clara, entonces vamos a extraer el fragmento que contiene la función. Para esto dividimos la descripción en items, utilizando la función split.(). Esta función separa una lista en items. Por ejemplo, yo puedo descomponer mi nombre en las letras que lo componen:

for l in "laura":
    print(l.split())

Asímismo, se puede separar la descripción en las palabras que la componen.

for record in SeqIO.parse(prot_file, "fasta"):
    print(record.description.split())

y luego le damos el número del item de interés [] (es el item 9 pero python comienza a contar desde 0, entonces es el 8)

from Bio.Seq import Seq
from Bio.SeqRecord import SeqRecord

for record in SeqIO.parse(prot_file, "fasta"):
    if (record.description.split()[9] == "glycoprotein"):
        print(record)

Hasta aquí hemos seleccionado las secuencias que contienen el término “glycoprotein”. Para escribirlas en un archivo en un formato adecuado de secuencias vamos a utilizar la función de biopython Bio.SeqIO.write() y a cambiar un poco el estilo.

glycoproteins = (record for record in SeqIO.parse(prot_file, "fasta") if (record.description.split()[9] == "glycoprotein"))
SeqIO.write(glycoproteins, "sars_covid19_glycoproteins.fa", "fasta")

Donde eperas que esté el archivo? Chequéalo mirando el directorio actual (“get current working directory”). Verifica que si se haya generado con el nombre y formatos correctos, y la información que deseamos que contenga.

import os

os.getcwd()
os.listdir()

Por último, vamos a escribir los comando que utilizamos en un script de python (con extensión “.py”) para dejar registro de lo que hicimos y pueda ser repetible

#!/usr/local/bin/

# Script por Laura Salazar Jaramillo
# Genera un subconjunto de glicoproteinas de SARS-covid 19

# Importamos las librerías necesarias
import os
import sys

import Bio
from Bio import SeqIO
from Bio.Seq import Seq

# Archivos de entrada (inputs)
dna_file = "/home/lsalazarj/bicomp2024/Sars_cov.dna.fa"
prot_file = "Sars_cov.prot.fa"


# Listamos el tamano del genoma
for record in SeqIO.parse(dna_file, "fasta"):
    print("La longitud en nc del genoma es:",len(record.seq))

# Listamos los id y longitudes de las proteínas
for record in SeqIO.parse(prot_file, "fasta"):
    print(record.id,len(record.seq),record.description.split()[9])

# Seleccionamos las glicoproteínas y las guardamos en un archivo aparte
glycoproteins = (record for record in SeqIO.parse(prot_file, "fasta") if (record.description.split()[9] == "glycoprotein"))
SeqIO.write(glycoproteins, "sars_covid19_glycoproteins.fa", "fasta")

Este script se puede ejecutar en la terminal escribiendo los comandos

python glycoproteins.py > glycoproteins.out

Alternativamente se pueden cambiar los permisos del archivo (chmod +x glycoproteins.py) y ejecutarlo

./glycoproteins.py

Ejercicio con secuencias de Citocromo B Vamos a utilizar las secuencias de Citocromo B descargadas previamente para repasar el manejo de secuencias utilizando python.

• 1. Cree una carpeta que se llame CytB_seqs

• 2. Entre a la carpeta

• 3. Copie desde allí las secuencias de Citocromo (o descárguelas nuevamente)

• 4. Active los entornos de conda

• 5. Comience Python

• 6. Llame las librerías

• 7. Defina los archivos input

• 8. Realice un loop sobre las secuencias de DNA para imprimir los nombres de los ID y los tamaños del gen: qué tanto difieren?

• 9. Repita los putnos 7 y 8 pero con el archivo de proteínas: qué tanto difieren éstas?

• 10. Ahora viene el reto de la clase: Utilice las secuencias de DNA para traducir las secuencias y guarde este archvivo.

• 11. Compare las secuencias obtenidas con las del archivo de las proteínas: son iguales?

• 12. Utilice la sección 3.8 del tutorial de Biopython y repita el punto 10

mkdir CytB_seqs
cd CytB_seqs
cp /Users/laura/docencia/eia/bioinfo/Online_BioInf/talleres/CytBDNA.txt .
cp /Users/laura/docencia/eia/bioinfo/Online_BioInf/talleres/CytBProt.txt .

import Bio
from Bio.SeqRecord import SeqRecord
from Bio import SeqIO

cytdna_file = "CytBDNA.txt"
cytprot_file = "CytBProt.txt"

for record in SeqIO.parse(cytdna_file, "fasta"):
        print(record.id,len(record.seq))

for record in SeqIO.parse(cytprot_file, "fasta"):
        print(record.id,len(record.seq))

for record in SeqIO.parse(cytdna_file, "fasta"):
        mRNA=record.seq.transcribe()
        print(record.id,mRNA)
        print(record.id,mRNA.translate("Vertebrate Mitochondrial"))