Treballarem amb els formats textuals de seqüenciació d'ADN, ARN i proteïnes més habituals en bioinformàtica, el FASTA i el Genbank; i també repassarem d'altres de coneguts.

Format FASTA

Probablement és el format de fitxer més utilitzat per a seqüències i un dels tipus de formats de fitxer més comuns en bioinformàtica.

El format de fitxer FASTA té els seus orígens en el programa FAST, utilitzat per a l'alineació de seqüències.

El format de fitxer es defineix simplement com un fitxer de text pla amb una o més entrades que consisteixen en una línia amb un símbol > seguit d'una línia de definició identificativa única, o defline, i una o més línies de dades de seqüència.

Crear un fitxer de text fasta és molt fàcil, tant en un editor de text pla com notepad o VSCode.

Podem crear un fitxer unifasta (una sola seqüencia) que es digui uniseq.fasta amb l'editor VSCode, simplement copia el text i guarda'l.

>Seqüència aminoàcids de prova
MTHCP*MTI*

O bé crear un fitxer multifasta (més d'una seqüència) anomenat sequences.fa des del terminal Linux:

echo ">a
ACGCGTACGTGACGACGATCG
>b
ATTTCGCGACTCTGCCTACGCTAC
>c
GGGAAACCTTTTTTT" > sequences.fa

Bingo! Ja tenim un fitxer multifasta :) amb les seqüències a, b i c.

El requisit fonamental és que el fitxer sigui plain text per tal que es pugui tractar amb qualsevol aplicació de processament de textos o llenguatge de programació.

Per tant, aquests fitxers es tracten millor en editors de text com nano, sublime o VSCode.

Per veure un fitxer FASTA des de la línia d'ordres sense editar-lo, pots fer servir la aplicació cat.

cat uniseq.fasta
>Seqüència aminoàcids de prova
MTHCP*MTI*

També pots llegir-les linia a linia amb more i less; que això és ideal per a fitxers grans:

less sequences.fa

Per sortir de less només has d’apretar la tecla q (“quit”).

També pots filtrar les primeres i últimes files amb les comandes: head i tail.

Característiques.

Realment no hi ha cap restricció sobre la seqüència de la definició en el format estàndard, excepte que la definició ha de seguir immediatament el > sense cap espai intermedi.

Per tant, dins de d’seqüència, hi ha dos coses:

  • Un comentari/descripció que comença per > i es una sola linea
  • La seqüència de bases dividides per línies de ...[70] caràcters, pot variar.

Si es un fitxer MULTI-FASTA (més d'una línia), és separen les seqüències amb una linea de capçalera ( > ).

Qualsevol tipus de seqüència o conjunt de seqüències es pot posar en un fitxer FASTA: els 3.000 milions de parells de bases del genoma humà complet, la seqüència d'ARN completa d'un virus o la seqüència d'ADN del promotor del gen del teu ratolí favorit són exemples vàlids.

On trobar fitxers FASTA.

Ja sabeu que un dels millors llocs on aconseguir fitxers .fasta d’organismes reals de forma legal és el repositori públic de l’NCBI, que ens facilita un cercador molt complet.

Exposem algusns enllaços directes com a exemples.

Gènoma complet del SarsCov2:

Gen de la orquídea:

Gen insulina a l’ésser humà.

Gen de l'insulina al gos:

Com descarregar-los ?

Per automatitzar el nostre programari, pots descarregar-te fitxers fasta com qualsevol fitxer.

Des del terminal (bash):

wget -O ls_orchid.fasta https://raw.githubusercontent.com/biopython/biopython/refs/heads/master/Doc/examples/ls_orchid.fasta

També pots veure com fer-ho amb Python i la llibreria urllib3; per tal de descarregar-lo únicament si no hi és.

https://xtec.dev/python/http/#cache

⚠️ Però compte! Hi ha restriccions si els volem baixar d'NCBI. ⚠️

NCBI té obert un repositori FTP amb fitxers de les seves bases de dades (com els articles PMC o PUBMED) que es troba a:

Però ni des d'aquest repositori ni des de programari fet amb Python o bash se'ns permet baixar directament qualsevol fitxer fasta.

Més endavant veurem com fer-ho amb un programari específic que distribueix NCBI/Entrez:

Lectura i tractament de fitxers fasta amb GNU/Linux (Bash).

Mostrar contingut fitxers

Obre un terminal en Linux (la icona és una finestra negra) o bé prem la drecera Ctrl + Alt + T, que funciona en distribucions basades en Debian (Ubuntu, Mint, PopOS …)

Introdueix aquesta comanda, que ens permet veure la capçalera i el contingut de les seqüències del fitxer sequences.fa (que hem creat abans):

grep -E '^\S+$' sequences.fa | awk '{printf "%s ", $0; getline; print $0}'

Resultat:

>a ACGCGTACGTGACGACGATCG
>b ATTTCGCGACTCTGCCTACGCTAC
>c GGGAAACCTTTTTTT

Gràcies al llenguatge del shell awk podem aconseguir-ho.

Una altra manera és usar un script de bash; que anomenarem readfasta.sh:

#!/bin/bash
echo 'Llegim tots els fitxers fasta de la carpeta actual.'
for i in *.fasta; do
    echo "Núm Seqüències"
    grep -c "^>" $i
    echo "Info Seqüències"
    # Capçalera seqüència.
    grep "^>" $i
    # Contingut seqüència:
    grep -v "^>" $i
done

Recorda com executar-lo:

chmod u+x readfasta.sh
./readfasta.sh

Recompte dades de seqüències

Una altra ordre útil de GNU/Linux és wc (abreviatura de "recompte de paraules") per comptar el nombre de línies d'un fitxer; juntament amb el programa de filtratge d'expressions grep:

 $ grep '>' sequences.fasta | wc
  3 3 9

El símbol pipa | (pipe) envia la sortida d'un programa (grep '>' sequences.fasta) a l'entrada d'un altre programa (wc), ja que amb aquest exemple estem enviant la sortida de grep al wc per comptar.

Important! Aneu amb compte per assegurar-vos que el símbol > estigui inclòs entre cometes en aquesta ordre. Si no s'inclouen les cometes, es sobreescriurà el fitxer!

Guardar resultats en nous fitxers

A més a més; un símbol > a la línia d'ordres redirigeix ​​la sortida d'un programa a un fitxer i sobreescriu el fitxer.

Per exemple, podem redirigir la sortida de l'exemple anterior així:

 $ grep ">" sequences.fasta | wc > numSequences.txt

D'aquesta manera no imprimeix res a la pantalla, sinó emmagatzema la informació en un fitxer.

Transcripció i edició de seqüències

També podem realitzar algunes operacions, com la transcripció; amb la comanda sed (stream editor).

 cat seq_3.fasta | sed ‘s/T/U/g’

Lectura de fitxers amb Python

Podem llegir fàcilment fitxers textuals amb Python, i agafar tota la informació de les cadenes.

# Llegir el fitxer amb 'with open' per assegurar que es tanca correctament.
with open('sequences.fa', 'r') as f:
    lines = f.read().splitlines()

# Diccionari per emmagatzemar les seqüències
genes = {}

# Variables per mantenir el context mentre iterem
gene_id = None
sequence = []

for line in lines:
    if line.startswith('>'):  # Les capçaleres comencen amb '>'
        if gene_id:
            # Guarda la seqüència acumulada abans de processar una nova capçalera
            genes[gene_id] = ''.join(sequence)
        gene_id = line[1:].strip()  # Guarda el nou identificador
        sequence = []  # Reinicia la seqüència
    else:
        sequence.append(line.strip())  # Afegeix la línia de seqüència

# Guarda l'última seqüència
if gene_id:
    genes[gene_id] = ''.join(sequence)

# Mostra els resultats
for gene_id, sequence in genes.items():
    print(f'ID: {gene_id}, Sequence: {sequence}, Seq.Long: {len(sequence)})')

Lectura de fitxers amb (Bio)Python

Si només necessitem llegir els fitxers i obtenir informació molt bàsica aquestes solucions que aporta GNU/Linux són bones i eficients.

Però en el nostre cas volem fer operacions complexes, com obtenir estadístiques, alinear seqüències, crear arbres… i per això és preferible usar Python.

Pots usar mètodes de Python per obrir i obtenir informació bàsica d’un fitxer fasta i obtenir estadístiques bàsiques (longitud, percentatge gc…) amb el suport de llibreries com Biopython.

Per exemple, aquest codi llegirà el contingut del fitxer sequence.fa, guardarà el contingut en l’objecte sequences i el mostrarà per pantalla.

from Bio import SeqIO

file = "sequences.fa"
sequences = SeqIO.parse(open(file), "fasta")

for record in sequences:
   print(record.id, record.seq)

El resultat és

a ACGCGTACGTGACGACGATCG
b ATTTCGCGACTCTGCCTACGCTAC
c GGGAAACCTTTTTTT

Pots observar que també hem introduït un altre objecte anomenat record, i cada seqüència del nostre fitxer FASTA s'emmagatzema en aquest objecte.

La crida als mètodes id i seq de l’objecte record retorna el text de defline i un objecte Seq.

Estructura dels SeqRecord de BioPython.

La classe SeqIO, ens servirà per llegir fitxers fasta, genbank i altres i fins i tot per a crear-ne de nous.

Tenim 2 mètodes de lectura:

read() --> Llegeix fitxers (.fasta i .genbank) que tenen només UNA (seqüència o fitxa) Resultat = SeqRecord

parse() --> Llegeix fitxers (.fasta i .genbank) que tenen només UNA O MOLTES (seqüències o fitxes de genbank)

Resultat = Llista SeqRecord's

Cada objecte SeqRecord conté, en el cas dels fitxers fasta:

  • Objecte Seq, amb la seqüència completa.
  • id. El que surt just després del símbol >
  • name. Metadada primera linia amb el nom.
  • description. Metadada primera linia amb descripció addicional.

Tenim el mètode d’escritura write, si volem guardar ràpidament les nostres investigacions (tot i que habitualment s'escriu als Genbank, no els fasta)

write(seq,info,format)
  • seq → La seqüència o llista de seqüències en format SeqRecord o format Seq.
  • info → Aquí hi ha la informació de la capçalera, es pot passar en un diccionari que contingui id, name, … Com abans, si volem un multifasta li hem de passar una llista.
  • format → ‘fasta’, ‘genbank’ ...

Script lectura de fitxer fasta millorat.

Si necessites un programa que en comptes de definir-li el fitxer .fasta dins del codi li puguis passar un fitxer .fasta que t’has descarregat prèviament pots usar aquest codi; que un cop llegit el fitxer transcriu a ARN totes les seqüències (tant si el fasta té 1 o diverses):

import sys
import os
from Bio import SeqIO, Seq

# Validacions
if len(sys.argv) != 2:
    print("Instruccions ús: " + sys.argv[0] + " <fitxer fasta>")
    sys.exit()
fastaFile = sys.argv[1]
if not os.path.isfile(fastaFile):
    print(f"Error: El fitxer {fastaFile} no existeix!")
    sys.exit()

print(f"Fitxer {fastaFile} llegit correctament. Transcrivint seqüències.")
sequences = SeqIO.parse(fastaFile, 'fasta')
for record in sequences:
    print(record.id, record.seq.transcribe())

print("Seqüències transcrites! Fins aviat!")

Aquí sota teniu la prova del seu funcionament:

>> $ python3 read_fasta.py sequences.fa
Fitxer sequences.fa llegit correctament. Transcrivint seqüències.
a ACGCGUACGUGACGACGAUCG
b AUUUCGCGACUCUGCCUACGCUAC
c GGGAAACCUUUUUUU
Seqüències transcrites! Fins aviat!
>> $ python3 read_fasta.py
Instruccions ús: read_fasta.py <fitxer fasta>
>> $ python3 read_fasta.py xd
Error: El fitxer xd no existeix!

Lectura de fitxers fasta d'aminoàcids.

Copia't aquesta seqüència en un fitxer fasta, disponible a https://www.ncbi.nlm.nih.gov/protein/AAD44166.1?report=fasta

echo ">gi|5524211|gb|AAD44166.1| cytochrome b [Elephas maximus maximus]
LCLYTHIGRNIYYGSYLYSETWNTGIMLLLITMATAFMGYVLPWGQMSFWGATVITNLFSAIPYIGTNLV
EWIWGGFSVDKATLNRFFAFHFILPFTMVALAGVHLTFLHETGSNNPLGLTSDSDKIPFHPYYTIKDFLG
LLILILLLLLLALLSPDMLGDPDNHMPADPLNTPLHIKPEWYFLFAYAILRSVPNKLGGVLALFLSIVIL
GLMPFLHTSKHRSMMLRPLSQALFWTLTMDLLTLTWIGSQPVEYPYTIIGQMASILYFSIILAFLPIAGX
IENY" > elephas_prog.fasta

Ara, prova el mateix programa que has creat per llegir el fitxer i veuràs que, tret de la transcrició que la ignorarà, funciona exactament igual.

Fitxer de prova:

>> $python3 read_fasta.py elephas_prog.fasta 
Fitxer elephas_prog.fasta llegit correctament. Transcrivint seqüències.
gi|5524211|gb|AAD44166.1| LCLYUHIGRNIYYGSYLYSEUWNUGIMLLLIUMAUAFMGYVLPWGQMSFWGAUVIUNLFSAIPYIGUNLVEWIWGGFSVDKAULNRFFAFHFILPFUMVALAGVHLUFLHEUGSNNPLGLUSDSDKIPFHPYYUIKDFLGLLILILLLLLLALLSPDMLGDPDNHMPADPLNUPLHIKPEWYFLFAYAILRSVPNKLGGVLALFLSIVILGLMPFLHUSKHRSMMLRPLSQALFWULUMDLLULUWIGSQPVEYPYUIIGQMASILYFSIILAFLPIAGXIENY
Seqüències transcrites! Fins aviat!

Execució de l'script:

python3 read_fasta.py sequences.fa
Fitxer sequences.fa llegit correctament. Transcrivint seqüències.
a 0.6190476190476191 TRT*RRS
b 0.5416666666666666 ISRLCLRY
c 0.3333333333333333 GKPFF
Seqüències transcrites! Fins aviat!

Activitat.

1.- Crea un programa com a script que llegeixi un fitxer fasta per paràmetre, calculi el percentatge GC i tradueixi les 30 primeres bases.

Ha de contolar que l'usuari hagi enviat un fitxer fasta ben formatat (amb les capçaleres que comencin per >).

import sys
import os
from Bio import SeqIO
from Bio.SeqUtils import gc_fraction

# Validacions
if len(sys.argv) != 2:
    print("Instruccions ús: " + sys.argv[0] + " <fitxer fasta>")
    sys.exit()
fastaFile = sys.argv[1]
if not os.path.isfile(fastaFile):
    print(f"Error: El fitxer {fastaFile} no existeix!")
    sys.exit()

print(f"Fitxer {fastaFile} llegit correctament. Transcrivint seqüències.")
sequences = SeqIO.parse(fastaFile, 'fasta')
for record in sequences:
    print(record.id, gc_fraction(record.seq,"weighted"), record.seq[0:30].translate())

print("Seqüències transcrites! Fins aviat!")

2.- Crea un nou programa que tradueixi un seqüència d’ADN a una de proteïnes. Al final, ha de guardar la proteïna en un altre fitxer que tingui el mateix nom juntament amb el text "_newprot". La capçalera ha de ser la mateixa que el fasta de nucleòtids, afegint el text ", translated by Biopython".

import sys
from Bio.Seq import Seq
from Bio import SeqIO
from Bio.SeqRecord import SeqRecord
import os

def tradueix_adn_a_proteines(input_file):
    """Tradueix seqüències d'ADN a proteïnes i les guarda en un nou fitxer FASTA."""
    try:
        # Llegeix el fitxer FASTA
        sequences = SeqIO.parse(open(input_file), "fasta")
        print(f"Fitxer llegit: {input_file}")
        
        # Inicialitza una llista per les seqüències traduïdes
        seq_records_proteina = []
        
        # Traducció de seqüències
        for seq_record in sequences:
            try:
                # Traduir a proteïna
                seq_proteina = seq_record.seq.translate()
                # Actualitzar la capçalera
                nova_capcalera = f"{seq_record.description}, translated by Biopython"
                # Crear un nou SeqRecord
                seq_record_proteina = SeqRecord(seq_proteina, id=seq_record.id, description=nova_capcalera)
                seq_records_proteina.append(seq_record_proteina)
            except Exception as e:
                print(f"Error en traduir {seq_record.id}: {e}")
        
        # Guarda les seqüències traduïdes
        output_file = os.path.splitext(input_file)[0] + "_newprot.fasta"
        SeqIO.write(seq_records_proteina, output_file, "fasta")
        print(f"Les seqüències traduïdes s'han guardat a {output_file}")
    except Exception as e:
        print(f"Error en processar el fitxer {input_file}: {e}")

if __name__ == "__main__":
    # Comprova si s'ha passat un argument d'entrada
    if len(sys.argv) != 2:
        print(f"Ús: {sys.argv[0]} <fitxer fasta d'ADN o ARN>")
        sys.exit(1)

    # Nom del fitxer d'entrada
    input_file = sys.argv[1]
    tradueix_adn_a_proteines(input_file)

Ho provem amb un fitxer fasta creat per nosaltres:

echo ">Demofa by Miquel
AGACTAGCAGACTATCGGTACGTCACATAG" > act2-nuc.fasta
python3 bioprot.py act2-nuc.fasta

Veiem que efectivament ha guardat un nou fitxer amb la traducció:

cat demo_newprot.fasta 
>Demofa by Miquel, translated by Biopython
RLADYRYVT*

Fastaiterator

Per a llegir fitxers multifasta grans, ens convé usar la classe FastaIterator per millorar el rendiment.

És més eficient en termes de memòria RAM quan tractes amb fitxers grans, ja que només es llegeix una seqüència a la vegada.

Anem a provar com funciona amb el multifasta ls_orchid.fasta i calcularem el temps d'execució.

from Bio import SeqIO
import time

MULTIFASTA_FILE = "ls_orchid.fasta"

def processa_multifasta(file_path):
    """
    Llegeix un fitxer multifasta utilitzant SeqIO.parse i mostra la longitud de cada seqüència.
    """
    try:
        print("Processant el fitxer multifasta...\n")
        
        start_time = time.time()
        
        # Iterador de SeqRecords des del fitxer FASTA
        for record in SeqIO.parse(file_path, "fasta"):
            print(f"ID: {record.id}")
            print(f"Descripció: {record.description}")
            print(f"Longitud: {len(record.seq)}")
            print("--------")
        
        elapsed_time = time.time() - start_time
        print(f"\nTemps total per processar el fitxer: {elapsed_time:.2f} segons")
    
    except Exception as e:
        print(f"Error en processar el fitxer: {e}")

if __name__ == "__main__":
    processa_multifasta(MULTIFASTA_FILE)

Instruccions per executar el programa amb mesurament del temps:

  1. Guarda el programa en un fitxer, per exemple, process_multifasta.py.
  2. Assegura't que tens un fitxer multifasta a la mateixa carpeta, com ara ls_orchid.fasta.
  3. Executa el programa amb la comanda time des de la línia de comandes:
time python process_multifasta.py

La sortida del programa, llegint el fitxer ls_orchid.fasta, seria:

Longitud: 744
--------
ID: gi|2765564|emb|Z78439.1|PBZ78439
Descripció: gi|2765564|emb|Z78439.1|PBZ78439 P.barbatum 5.8S rRNA gene and ITS1 and ITS2 DNA
Longitud: 592
--------

Temps total per processar el fitxer: 0.06 segons

Format FastQ

Aquí, treballarem amb fitxers FASTQ, la sortida en format estàndard que utilitzen els seqüenciadors moderns. Aprendràs a treballar amb puntuacions de qualitat per base.

Què és el format FASTQ?**

El format FASTQ s’utilitza àmpliament per informar de dades de seqüenciació d’alt rendiment. Cada entrada consta de quatre línies:

  1. Identificador: Comença amb @ i conté informació sobre la seqüència.
  2. Seqüència: La seqüència de nucleòtids corresponent.
  3. Separador: Una línia amb + (i un text amb comentaris, opcionals).
  4. Qualitat: Una cadena que codifica les puntuacions de qualitat per a cada nucleòtid.

Exemple d’entrada FASTQ:

@SRR001666.1 071112_SLXA-EAS1_s_7:5:1:817:345 length=36
GGGTGATGGCCGCTGCCGATGGCGTCAAATCCCACC
+SRR001666.1 071112_SLXA-EAS1_s_7:5:1:817:345 length=36
IIIIIIIIIIIIIIIIIIIIIIIIIIIIII9IG9IC

Càlcul de puntuacions PHRED

Les puntuacions de qualitat en un fitxer FASTQ segueixen l’estàndard PHRED, que es defineix com:

S_PHRED = -10 * log10(p_err)

  • p_err: Probabilitat d'error en la identificació d'una base.
  • Una puntuació PHRED alta significa una probabilitat baixa d'error.

Conversió de puntuacions PHRED a caràcters ASCII:

  1. Cada caràcter de la cadena de qualitat representa una puntuació PHRED.
  2. La puntuació es calcula amb: S_PHRED = ord(Q_char) - 33
    • ord(Q_char): Retorna el valor ASCII (en decimal) del caràcter.
    • El -33 ajusta per assegurar que el primer caràcter imprimible correspon a un valor de puntuació 0.

Probabilitat d'error derivada:

p_err = 10^(-S_PHRED / 10)

Per exemple, el primer caràcter imprimible és el !, que dins de la taula ASCII és el 33 llavors aquest és el que indica una qualitat més baixa.

I si tenim els caràcters B i E; quin té qualitat més alta ? Doncs E perquè té un número més alt de la taula ASCII (el 69) que la B (el 66).

Codi exemple per calcular el % probabilitat d'error d'un caràcter de fastq

def calcular_probabilitats_error(puntuacions_qualitat):
    probabilitats = [10 ** (-(ord(caracter) - 33) / 10.0) for caracter in puntuacions_qualitat]
    return probabilitats

# Exemple d'ús:
puntuacions_qualitat = ['!', 'A', 'a', '%', '-']
probabilitats_error = calcular_probabilitats_error(puntuacions_qualitat)

for i in range(len(puntuacions_qualitat)):
    print(f"Puntuació: {puntuacions_qualitat[i]}, Probabilitat d'error: {probabilitats_error[i]:.4f}")

Resultats:

Puntuació: !, Probabilitat d'error: 1.0000
Puntuació: A, Probabilitat d'error: 0.0006
Puntuació: a, Probabilitat d'error: 0.0000
Puntuació: %, Probabilitat d'error: 0.3981
Puntuació: -, Probabilitat d'error: 0.0631

Lectura i processament de fitxers FASTQ:

Biopython ens facilita la lectura del fitxer i també el càlcul de les probabilitats d'error i el valor s_phred de la qualitat dels caràcters.

Això és possible gràcies a que la classe Seq te un atribut que es diu letter_annotations i s'omple quan llegim un fitxer fastq.

from Bio import SeqIO

file = "sequences.fq"
data = SeqIO.parse(file, "fastq")

for record in data:
    print(f"ID: {record.id}")
    print(f"Seqüència: {record.seq}")
    print(f"Puntuacions PHRED: {record.letter_annotations['phred_quality']}")

# Càlcul de probabilitats d'error amb list comprehension
p_err_array = [round(10 ** (-float(Q) / 10.0),6) for Q in record.letter_annotations["phred_quality"]]
print(f"% p_err: {p_err_array}")

Sortida

ID: SRR993731.910
Seqüència: TGTCCTCTGGCTCCAGGTCTCATGATGAAAAAATTTATGGAGTCCTGGACA
Puntuacions PHRED: [26, 30, 27, 35, 26, 11, 17, 33, 32, 28, 18, 32, 10, 17, 32, 36, 26, 27, 27, 32, 27, 36, 36, 37, 29, 36, 31, 37, 27, 37, 37, 40, 32, 30, 35, 35, 34, 35, 27, 35, 26, 35, 24, 33, 24, 30, 2, 2, 35, 22, 36]

Com descarregar i llegir fitxers fastq comprimits.

Aquest és un fitxer força petit (27 MB) i representa part de les dades seqüenciades d'una dona yoruba (NA18489). Si consultes el 1.000 Genomes Project, veuràs que la gran majoria dels fitxers FASTQ són molt més grans (fins a dues ordres de magnitud més grans).

Per obtenir més informació sobre el 1.000 Genomes Project, pots visitar la seva pàgina oficial:

import os
import subprocess
import gzip
from Bio import SeqIO

url = "ftp://ftp.1000genomes.ebi.ac.uk/vol1/ftp/phase3/data/NA18489/sequence_read/SRR003265.filt.fastq.gz"
filename = "SRR003265.filt.fastq.gz"

if not os.path.exists(filename):
    print(f"Descarregant {filename}...")
    # Utilitzar el comandament 'wget' mitjançant subprocess, ja que sinó NCBI no ens deix.
    subprocess.run(["wget", "-nd", url])
    print(f"{filename} descarregat correctament.")
else:
    print(f"El fitxer {filename} ja existeix, no es descarregarà.")

# Obrir i llegir el fitxer gzip un cop descarregat
with gzip.open(filename, 'rt', encoding='utf-8') as f:
    recs = SeqIO.parse(f, 'fastq')
    rec = next(recs)  # Agafar el primer registre (read)
    print(f"ID: {rec.id}")
    print(f"Descripció: {rec.description}")
    print(f"Seqüència: {rec.seq}")
    print(f"Puntuacions de qualitat: {rec.letter_annotations['phred_quality']}")

Indexació d'un fitxer FASTQ - com llegir fitxers FASTQ grans més ràpidament?

Els fitxers FASTQ solen ser molt grans, amb milions de lectures. A causa de la gran quantitat de dades, no es poden carregar tots els registres a la memòria alhora. Per això, quan es realitzen tasques de filtratge i retallat, podem iterar sobre el fitxer observant només un SeqRecord a la vegada.

No obstant això, de vegades no es pot utilitzar un gran bucle o un iterador - potser necessitaràs un accés aleatori a les lectures. En aquest cas, la funció Bio.SeqIO.index() pot resultar molt útil, ja que et permet accedir a qualsevol lectura en el fitxer FASTQ pel seu nom.

Així que això és útil quan saps l'ID de la seqüència a la qual vols accedir, ja que et permet llegir-la de manera instantània, en comptes de fer un bucle sobre el fitxer fins trobar la seqüència desitjada.

Descarregar un fitxer FASTQ de ~1GB - triga uns 2 minuts

wget ftp://ftp.sra.ebi.ac.uk/vol1/fastq/SRR494/SRR494102/SRR494102.fastq.gz

Activitat.

3.- Agafa un fitxer fastq d'NCBI o una altra font confiable i mostra el seu contingut i la seva qualitat.

Recursos útils

Format Genbank

Genbank conté la majoria de les seqüències d'ADN, ARN i proteïnes conegudes del món, i també emmagatzema informació bibliogràfica per a les seqüències.

Va ser creat l’any 1979 i finalitzat l’any 1982. 10 anys més tard es va integrar totalment dins l’NCBI.

La versió actual de Genbank, a partir d'octubre de 2015, conté 188,372,017 seqüències, que inclouen 202,237,081,559 nucleòtids.

Pots accedir a les entrades de Genbank des de la seva pàgina principal:

També pots cercar les entrades .genbank des de diverses bases de dades d’NCBI: com Nucleotide, Protein… de la mateixa manera que ho feiem amb els fitxers .fasta

Genbank és de naturalesa arxivística, de manera que pot contenir entrades redundants.

Genbank és una base de dades, però també és un format de fitxer molt detallat.

Com descarregar manualment fitxer Genbank i llegir-lo amb Biopython.

Ens descarregarem un exemple manualment per provar, del gen complet mitocondrial de l’elefant:

Si ens el descarreguem amb l’wget ens surt en format XML i per ara no ens interessa; per tant, millor ho fem del portal web de l’NCBI (opcions Complete Record, File, Genbank)

Hem renombrat el fitxer per a distigir-lo.

mv sequences.gb NC_005129.2_elephas_mitoc.genbank

Ara sí, podem usar el següent codi:

from Bio.Seq import Seq
from Bio import SeqIO

# Llegeix multifasta o multigenbank, retorna un Iterador.
demo_genbank = list(SeqIO.parse("NC_005129.2_elephas_mitoc.genbank", "genbank"))
print(demo_genbank)

Que com a sortida mostra un resum:

[SeqRecord(seq=Seq('GTTAATGTAGCTTAAAACAAAAGCAAGGTACTGAAAATACCTAGACGAGTATAT...CTC'), id='NC_005129.2', name='NC_005129', description='Elephas maximus mitochondrion, complete genome', dbxrefs=['BioProject:PRJNA927338'])]

Com descarregar fitxer Genbank llegir-lo amb Biopython i utilitat NCBI-Entrez.

Ara, aprendrem a utilitzar el mòdul de Biopython Entrez per recuperar automàticament un fitxer en format Genbank.

Si ja coneixes un número d'accés gi per al gen, pots recuperar les dades directament de NCBI.

És molt interessant usar-lo si volem descarregar o revisar moltes seqüències alhora.

⚠️ Heu de saber que l’únic requeriment és assegurar-se de no fer més de 3 operacions per minut, i que cal indicar l’email si no tens una API Key. ⚠️

Aquí tens un exemple, amb una part de la sortida (el fitxer Genbank complet és massa gran per mostrar-lo aquí).

from Bio import Entrez
Entrez.email = "mamoro10@xtec.cat"
handle = Entrez.efetch(db="protein", id='NP_003173', rettype="gb", retmode="text")
print(handle.read())

Sortida (parcial):

 LOCUS       NP_003173    129 aa    DNA_input linear  PRI 28-NOV-2015
 DEFINITION  protachykinin-1 isoform beta precursor [Homo sapiens].
 ACCESSION   NP_003173
 VERSION     NP_003173.1 GI:4507341
 DBSOURCE    REFSEQ: accession NM_003182.2
 KEYWORDS    RefSeq.
 SOURCE      Homo sapiens (human)
   ORGANISM  Homo sapiens
             Eukaryota; Metazoa; Chordata; Craniata; Vertebrata; Euteleostomi;
             Mammalia; Eutheria; Euarchontoglires; Primates; Haplorrhini;
             Catarrhini; Hominidae; Homo.
   REFERENCE 1 (residues 1 to 129)
   AUTHOR    Agaeva,G.A., Agaeva,U.T. and Godjaev,N.M.
   TITLE     [Particularities of Spatial Organization of Human Hemokinin-1
             and Mouse/Rat Hemokinin-1 Molecules]
   JOURNAL   Biofizika 60 (3), 457-470 (2015)
   PUBMED    26349209
   REMARK    GeneRIF: The spatial structures of human, mouse, and rat
             hemokinin-1 protein isoforms have been presented.
 ...
 FEATURES    Location/Qualifiers
   source    1..129
             /organism="Homo sapiens"
             /db_xref="taxon:9606"
             /chromosome="7"
             /map="7q21-q22"
   Protein   1..129
             /product="protachykinin-1 isoform beta precursor"
             /note="neuropeptide gamma; neuropeptide K; tachykinin,
             precursor 1 (substance K, substance P, neurokinin 1,
             neurokinin 2, neuromedin L, neurokinin alpha, neuropeptide
             K, neuropeptide gamma); tachykinin 2; protachykinin;
             preprotachykinin; neurokinin A; protachykinin-1; PPT"
 ...

   CDS       1..129
             /gene="TAC1"
             /gene_synonym="Hs.2563; NK2; NKNA; NPK; TAC2"
             /coded_by="NM_003182.2:247..636"
             /note="isoform beta precursor is encoded by transcript variant beta"
             /db_xref="CCDS:CCDS5649.1"
             /db_xref="GeneID:6863"
             /db_xref="HGNC:HGNC:11517"
             /db_xref="HPRD:08876"
             /db_xref="MIM:162320"
   ORIGIN
     1 mkilvalavf flvstqlfae eiganddlny wsdwydsdqi keelpepfeh llqriarrpk
    61 pqqffglmgk rdadssiekq vallkalygh gqishkrhkt dsfvglmgkr alnsvayers
   121 amqnyerrr

Podem millorar-lo una mica; assegurant-nos que només es descarrega si no el tenim:

from Bio import Entrez, SeqIO
import os

GENBANK_NAME = "SarsCov2.gb"
DATABASE = "nuccore"
ACCESSION = "NC_045512"
Entrez.email = "mamoro10@xtec.cat"

filepath = os.path.join(os.getcwd(), GENBANK_NAME)

if not os.path.exists(filepath):
    with Entrez.efetch(db=DATABASE, id=ACCESSION, rettype="gb", retmode="text") as handle:
        with open(filepath, 'w') as genbank_file:
            genbank_file.write(handle.read())
    print(f"Genbank descarregat i guardat com a {GENBANK_NAME}")
else:
    print(f"Ja existeix el fitxer {GENBANK_NAME}.")
    demo_genbank = list(SeqIO.parse(GENBANK_NAME, "genbank"))
    print(demo_genbank)

Estructura d’un fitxer Genbank.

El format Genbank està pensat perquè nosaltres el poguem llegir sense eines especials fins a cert punt, però els fitxers tenen un format que els permet analitzar-los amb eines de programari, com el de Biopython.

En resum, l’esquema del format Genbank és:

Anem a centrar-nos en explicar els 2 primers blocs del Genbank; acompanyats d’un exemple que des del 2020 ha preocupat a tot el món i el coneixem força bé. Aquí el tenim:

Bloc 1 del Genbank: Header

Els GenBank sempre comencen amb la paraula LOCUS i contenen els següents camps al Header o annotations (com es coneixen en Biopython):

  1. Locus: Identificació del lloc.
  2. Definition: Descripció textual.
  3. ACCESSION: Identificador únic de la fitxa.
  4. Version: Notació general formada per l’ID + el número de versió.
  5. DBLINK: Projecte d’on prové la seqüència.
  6. KEYWORDS: Paraules clau per als cercadors.
  7. SOURCE: Descripció científica que pot incloure taxonomies (en format de llista).
  8. REFERENCE: Informació d’altres persones i publicacions relacionades. Inclou camps com:
    • AUTHORS
    • TITLE
    • JOURNAL
    • PUBMED (si escau).
  9. COMMENTS: Text lliure o informacions diverses.

Bloc 2 del Genbank: Gene Annotation

Aquest bloc, conegut en Biopython com a features i també com Gene Annotation, conté informació sobre els gens de la seqüència. Inclou dades sobre cadenes codificants (CDS), promotors, pèptids o regions no traduïdes. A continuació, els camps més importants:

  1. source: Informació general sobre la seqüència (longitud, organisme, tipus, etc.).
  2. gene: Detalls específics del gen.
  3. CDS:
    • Cadena codificant (indica la regió codificant).
    • translation: Seqüència proteica derivada.
    • Sovint inclou diversos CDS o combinacions (joins).
    • Dins de cada CDS:
      • Gen associat.
      • Codó d’inici (si és l’1 es tradueix amb la lletra M).
      • ID de la proteïna funcional.
  4. Altres features:
    • 5'UTR: Regió inicial no codificant (per exemple: 1..265 no codificat com a gene).
    • mat_peptide: Regions concretes de la proteïna que tenen una funció específica (per exemple, unir-se a altres proteïnes).

Per més informació sobre els conceptes i la taula de codons, pots consultar:
Gene Expression

Bloc 3 del Genbanki i final: Seqüència

El tercer és la cadena d'ADN, ARN o aminoàcids. A Genbank surt en blocs de 10 en 10 per facilitat la lectura en paper, però BioPython ajunta la cadena automàticament.

La línia // marca el final d'un fitxer genbank.

RefSeq

RefSeq (Reference Sequence) és una base de dades completa i no redundant de seqüències biològiques, mantinguda pel NCBI, que inclou ADN, ARN i proteïnes. Aquí tens els seus principals detalls i característiques:

  1. RefSeq vs. GenBank:

    • RefSeq: Seqüències seleccionades i anotades amb gran cura i precisió.
    • GenBank: Format de fitxer i base de dades que inclou seqüències, revisades o no.

  2. Connexió amb altres bases de dades:

    • RefSeq enllaça proteïnes i seqüències d'àcids nucleics corresponents a un gen específic.

  3. Estadístiques:

    • RefSeq manté aproximadament 100,000,000 entrades, cadascuna amb un identificador únic.

  4. Identificadors RefSeq:

    • XR_: Seqüències d'ARN que no són ARN missatgers.
    • XM_: ARN missatgers.
    • XP_: Proteïnes.
    • NP_: Proteïnes conegudes i confirmades.
    • NC_: Cadenes de nucleòtids (normalment ADN).
    • NR_: Seqüències d'ARN que no són ARN missatgers.
    • NM_: ARN missatgers.

Exemple de lectura d'un registre RefSeq

from Bio import Entrez
Entrez.email = "mamoro10@xtec.cat"

record = Entrez.read(Entrez.esearch(db="protein", term="NP_003173"))
print(record['IdList'])

handle = Entrez.efetch(db="protein", id=record["IdList"][0], rettype="fasta")
print(handle.read())

Sortida del codi:

['4507341']
>gi|4507341|ref|NP_003173.1| protachykinin-1 isoform beta precursor [Homo sapiens]
MKILVALAVFFLVSTQLFAEEIGANDDLNYWSDWYDSDQIKEELPEPFEHLLQRIARRPKPQQFFGLMGK
RDADSSIEKQVALLKALYGHGQISHKRHKTDSFVGLMGKRALNSVAYERSAMQNYERRR

Altres bases de dades del NCBI

El NCBI és un recurs immens amb diverses bases de dades útils. A continuació, algunes destacades:

  1. Entrez: Motor de cerca global entre bases de dades NCBI (detall a la pròxima activitat).
  2. Gene: Base de dades de gens amb informació sobre nomenclatura, RefSeq, mapes, vies, variacions i fenotips.
  3. GEO: Gene Expression Omnibus. Base de dades d'expressió gènica i dades de seqüències d'alt rendiment.
  4. OMIM: Informació sobre trets heretables, fenotips genètics i variants al·lèliques associades a malalties.
  5. Taxonomy: Base de dades de classificació i nomenclatura (representa el 10% de les espècies del planeta).
  6. PubMed: Motor de cerca per accedir a MEDLINE (citacions i resums de recerca biomèdica).
  7. PMC: Arxiu electrònic d'articles científics de text complet, amb accés gratuït al contingut.

Exemple de lectura d'un article de PMC

Resum: A Systematic Review: Do the Use of Machine Learning, Deep Learning, and Artificial Intelligence Improve Patient Outcomes in Acute Myocardial Ischemia Compared to Clinician-Only Approaches?

Llegeix l’article complet:
PubMed

Altres Bases de Dades Biològiques

Ensembl

Ensembl és una base de dades europea especialitzada en seqüències biològiques, desenvolupada per l'European Bioinformatics Institute (EBI), part de l'EMBL.

Accessible a través de http://ensembl.org, proporciona una gran col·lecció de dades biològiques, similar a la base de dades NCBI/GenBank. A més, ofereix l'eina BioMart, que facilita la recuperació d'informació com ara gens, transcrits i seqüències de proteïnes.

UCSC Genome Bioinformatics

La University of California, Santa Cruz (UCSC) manté la base de dades UCSC Genome Bioinformatics, que proporciona recursos i eines per a la recerca genòmica, principalment en models animals.

Disponible a http://genome.ucsc.edu, inclou un navegador genòmic, pàgines de descàrregues i programari com BLAT per a alineaments ràpids de seqüències, amb integració directa als navegadors genòmics per a visualitzar els resultats.

UniProt

La Universal Protein Resource (UniProt) és una base de dades de seqüències de proteïnes creada per col·laboració entre l'EMBL, l'Swiss Institute of Bioinformatics, i el Protein Information Resource (PIR).

Inclou múltiples bases de dades, com UniRef, que organitza seqüències en clústers segons el percentatge d'identitat. Per exemple, UniRef100 agrupa seqüències idèntiques de com a mínim 11 residus, i UniRef90 fa clústers amb un 90% o més d'identitat.

DNA Data Bank of Japan (DDBJ)

És una base de dades pública que recopila, gestiona i distribueix dades de seqüències genòmiques. Com a membre del International Nucleotide Sequence Database Collaboration (INSDC), treballa conjuntament amb l'NCBI/GenBank als EUA i l'EMBL-EBI/ENA a Europa.

Accessible a través de https://www.ddbj.nig.ac.jp/index-e.html, és un recurs crític per a la recerca genòmica global.

També proporciona eines avançades com BLAST i sistemes per explorar les dades genòmiques. És especialment valuós per a dades genòmiques d'espècies asiàtiques, oferint un punt de vista complementari a altres bases de dades globals.

Nota final
Aquestes bases de dades són essencials per a la biologia computacional i la genòmica, ja que proporcionen accés a seqüències i eines per analitzar-les. Els investigadors poden triar entre Ensembl, UCSC, UniProt o DDBJ segons el tipus de dades o eines que necessitin.

Llegir tots els blocs del fitxer Genbank

Si volem tractar parts importants del Genbank (l’id Accession, l’organisme, taxonomia, la seqüència d’ADN, algunes de les cadenes codificants de proteïnes (CDS), la llista d’autors, comentaris, entre d’altres), el codi anterior ens és insuficient.

Per això us proporcionem una millora del codi anterior.

Per simplificar-lo, seguim amb el Genbank de l’exemple 1 (usant el mètode SeqIO.read, també pot funcionar amb parse).

from Bio                import SeqIO          # SeqIO is a module (several classes)
from Bio.SeqRecord      import SeqRecord      # SeqRecord is a class
#from Bio.SeqFeature     import SeqFeature, FeatureLocation # Inner classes
from pprint             import pformat, pp    # Així veiem millor els diccionaris.

# Llegim un fitxer, que és unigenbank.
GENBANK_NAME: str = "NC_005129.2_elephas_mitoc.genbank"
with open(GENBANK_NAME, "r") as file:
    demo_genbank: SeqRecord = SeqIO.read(GENBANK_NAME, "genbank")
    # També funciona amb:
    # demo_genbank: list[SeqRecord] = list(SeqIO.read(GENBANK_NAME, "genbank"))

    print('SeqRecord from genbank file:')
    # Com hem fet amb el fasta, podem llegir dades del SeqRecord del genbank.
    # Annotacions, del bloc 1.
    # Ens retorna un diccionari. Els camps compostos (com autors) són una llista.
    pp(demo_genbank.annotations)

    # Si només ens interessen les referències i els seus autors.
    authors = []
    for reference in demo_genbank.annotations.get('references', []):
        authors.extend(reference.authors.split(", "))
    # Si volem eliminar duplicats
    authors = set(authors)
    print("Llista d'autors:", authors)

    # Solució més completa:
    references = demo_genbank.annotations.get('references', [])
    for i, reference in enumerate(references, 1):
        authors = reference.authors.split(", ")
        title = reference.title
        journal = reference.journal
        pubmed = reference.pubmed_id

        print(f"Referència {i}:")
        print(f"  Autors: {authors}")
        print(f"  Títol: {title}")
        print(f"  Journal: {journal}")
        print(f"  PubMed ID: {pubmed}")

    # Si volem obtenir totes les features en general.
    # Per no allargar la sortida, posem les 3 primeres.
    features = demo_genbank.features[0:3]
    for feature in features:
        print("Tipo de característica:", feature.type)
        print("Ubicació:", feature.location)

        # Accedeix als qualifiers (etiquetes) de la característica
        qualifiers = feature.qualifiers
        for key, value in qualifiers.items():
            print(f"{key}: {value}")

        print()
       
    # Si només ens interessen els comentaris.
    comments = demo_genbank.annotations.get('comment', None)
    if comments:
        print("Llista de comentaris:", comments)

    # Per llegir seqüència de nucleòtids, del bloc 3, és igual que el fasta.
    # Fixeu-vos que surt tota junta, a diferència del genbank.
    print('Primers 100 caràcters Seqüència ')
    print(demo_genbank.seq[0:100])

Editar un fitxer Genbank amb Biopython.

De la mateixa manera que passa amb els fasta i fastq, podem editar la informació del genbank fàcilment amb el mètode write del mòdul SeqIO.

Això és per a facilitar sinèrgies amb la comunitat científica.

Nosaltres no ho farem habitualment, però si us fa falta veiem un exemple. El codi:

from Bio import SeqIO
from Bio.SeqFeature import Reference
from pprint import pp

file_path = "./data/example.genbank"
with open(file_path, "r") as file:
    # Lee el archivo GenBank
    record = SeqIO.read(file, "genbank")

    # Obtiene la lista de referencias existentes o crea una nueva lista
    references = record.annotations.get('references', [])

    # Crea una nueva referencia
    new_reference = Reference()
    new_reference.authors = 'Amoros,M.A.'
    new_reference.title = 'Example Omnics'
    new_reference.location = f"bases 1 to {len(record.seq)}"

    # Agrega la nueva referencia a la lista de referencias
    references.insert(0, new_reference)
    record.annotations['references'] = references

# Ruta al archivo GenBank de salida
output_file_path = "./data/example_updated.genbank"
with open(output_file_path, "w") as output_file:
    SeqIO.write(record, output_file, "genbank")

print(f"Nueva referencia y versión agregadas a {output_file path}")

Genbank original:

LOCUS       123456                      9 bp    DNA              UNK 01-JAN-2024
DEFINITION  GATA seq.
ACCESSION   123456
VERSION     123456
KEYWORDS    .
SOURCE      Random bacteria
  ORGANISM  Random bacteria
            Bacteria; Random bacteria
REFERENCE   1  (bases 1 to 9)
  AUTHORS   Bacterio,P.
  TITLE     El sulfat atòmic
FEATURES             Location/Qualifiers
ORIGIN
        1 atggattga
//

El nou genbank ha de tenir les següents referències:

REFERENCE   1
  AUTHORS   Amoros,M.A.
  TITLE     Example Omnics
REFERENCE   2  (bases 1 to 9)
  AUTHORS   Bacterio,P.
  TITLE     El sulfat atòmic

Activitat

4.- Utilitzant el genbank del gènoma de referència del virus SarsCov2; respon als següents punts usant Biopython:

  • Quina longitud té la secuencia ?
  • Mostra’n els 30 primers caràcters.
  • Mostrar quantes annotations (bloc 1) hi ha, l’accession i la referència més recent (la primera del fitxer genbank).
  • Mostra el type i la location de totes les features del genbank.

Hauries de provar si funciona amb altres fitxers genbank.

Solució

5.- Crea una API per tal de tenir uns mètodes GET que serviran per consultar informació fasta i genbank que ja tinguem guardats al servidor. Retornaran la resposta en format JSON.

6.- Crea un portal web que proporcioni un formular per a permeti mostrar informació d'un genbank que ja tinguem al servidor.

Per exemple, pots tenir els fitxers de seqüències de la Orquídea i els del SarsCov2. Assegura’t que el seu nom coincideixi amb l’accession o codi; o que almenys el contingui.

Informació que podem aconseguir:

Dels fasta:

  • La descripció o descripcions (si és multifasta).
  • Mostrar la seqüència o les seqüències (si és multifasta).
  • Tractament errors: Si no existeix un fasta amb el codi indicat s’ha d’informar a l’usuari i proposar-li un accession que existeixi.

Dels genbank:

  • Títol
  • Accesion id
  • Enllaç directe a l’NCBI (no ho fa Biopython però és fàcil d’aconseguir)
  • Referència més recent.
  • Número de features
  • Type i location de tots els CDS.
  • 30 primers caràcters de la cadena d’origen.
  • Tractament errors: Si no existeix un genbank amb l’accession indicat s’ha d’informar a l’usuari i proposar-li un accession que existeixi.

Solució (parcial)

Eines per a resoldre la activitat

En l'anterior projecte tens la manera d'aconseguir enviar per formulari dades al servidor:

I fins i tot tens un visualitzador de seqüències que et permet visualitzar totes les features d'una seqüència o d'un Genbank sencer:

Altres formats:

TODO Incorporar algún exemple més.