Na Prática: Fundamentos do Biopython

pip install biopython

from Bio import SeqIO          # Parsing de arquivos e E/S de sequências
from Bio.Seq import Seq        # O objeto Seq central
from Bio import Entrez         # Acesso ao banco de dados do NCBI
from Bio import Align          # Alinhamento par a par e múltiplo
from Bio import pairwise2      # Alinhamento par a par legado
from Bio.PDB import PDBParser  # Parsing de estruturas

from Bio.Seq import Seq

# Criar uma sequência
dna = Seq("ATGCGATCGATCGATCG")

# Operações semelhantes a string funcionam como esperado
print(len(dna))           # 18
print(dna[0:6])           # ATGCGA
print(dna.count("ATG"))   # 1

# Operações específicas de biologia
print(dna.complement())           # TACGCTAGCTAGCTAGC
print(dna.reverse_complement())   # CGATCGATCGATCGCAT

# Transcrição: DNA → RNA (T → U)
rna = dna.transcribe()
print(rna)  # AUGCGAUCGAUCGAUCG

# Tradução: RNA → proteína
protein = rna.translate()
print(protein)  # MRSS* (onde * é o códon de parada)

# DNA → proteína direto
protein = dna.translate()
print(protein)       # MRSS*
print(protein[:-1])  # Remover códon de parada: MRSS

# Por padrão, Seq é sem tipo — funciona para DNA, RNA ou proteína
# Para sequências de RNA, use U em vez de T
rna = Seq("AUGCGAUCG")

# MutableSeq permite modificação in-place
from Bio.Seq import MutableSeq
mseq = MutableSeq("ATGCGATCG")
mseq[0] = "G"     # substituição pontual
print(mseq)       # GTGCGATCG

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

record = SeqRecord(
    Seq("ATGCGATCGATCG"),
    id="NM_007294",
    name="BRCA1",
    description="Homo sapiens BRCA1 mRNA",
    annotations={"molecule_type": "mRNA"}
)

print(record.id)           # NM_007294
print(record.seq)          # ATGCGATCGATCG
print(record.description)  # Homo sapiens BRCA1 mRNA

from Bio import SeqIO

# Analisar um arquivo FASTA
records = list(SeqIO.parse("sequences.fasta", "fasta"))
for record in records:
    print(f"{record.id}: {len(record.seq)} pb")

# Analisar um arquivo GenBank
for record in SeqIO.parse("brca1.gb", "genbank"):
    print(record.id)
    print(record.description)
    print(f"Features: {len(record.features)}")

# Ler um único registro
record = SeqIO.read("single.fasta", "fasta")  # levanta ValueError se >1 registro

# Escrever em FASTA
SeqIO.write(records, "output.fasta", "fasta")

# Converter entre formatos
with open("output.gb", "w") as f:
    SeqIO.write(records, f, "genbank")

# Converter FASTQ → FASTA
records = SeqIO.parse("reads.fastq", "fastq")
SeqIO.write(records, "reads.fasta", "fasta")

record = SeqIO.read("brca1.gb", "genbank")

for feature in record.features:
    print(f"Tipo: {feature.type}")
    print(f"Localização: {feature.location}")
    print(f"Fita: {feature.location.strand}")  # 1 (forward) ou -1 (reverse)

# Filtrar por features CDS
cds_features = [f for f in record.features if f.type == "CDS"]

for cds in cds_features:
    gene = cds.qualifiers.get("gene", ["unknown"])[0]
    protein_id = cds.qualifiers.get("protein_id", ["?"])[0]

    # Extrair a sequência CDS (trata junções e fita)
    cds_seq = cds.extract(record.seq)
    protein = cds_seq.translate(to_stop=True)

    print(f"{gene} ({protein_id}): {len(protein)} aa")

from Bio.SeqFeature import SimpleLocation

feature = record.features[0]
loc = feature.location

print(loc.start)   # início baseado em 0 (ExactPosition)
print(loc.end)     # fim baseado em 0 (ExactPosition)
print(loc.strand)  # 1 ou -1

# Verificar se é composto (com splicing)
from Bio.SeqFeature import CompoundLocation
if isinstance(loc, CompoundLocation):
    for part in loc.parts:
        print(f"  Éxon: {part.start}–{part.end}")

from Bio.SeqUtils import gc_fraction

seq = Seq("ATGCGATCGATCGATCG")
gc = gc_fraction(seq)
print(f"Conteúdo GC: {gc:.1%}")  # 47.1%

# Composição manual
from collections import Counter
counts = Counter(str(seq))
print(counts)

def find_orfs(seq: Seq, min_length: int = 100) -> list:
    orfs = []
    for strand, nuc in [(+1, seq), (-1, seq.reverse_complement())]:
        for frame in range(3):
            i = frame
            while i < len(nuc) - 2:
                codon = nuc[i:i+3]
                if codon == "ATG":
                    for j in range(i, len(nuc) - 2, 3):
                        if nuc[j:j+3] in ("TAA", "TAG", "TGA"):
                            orf_len = j + 3 - i
                            if orf_len >= min_length:
                                orfs.append({
                                    "start": i, "end": j + 3,
                                    "strand": strand, "frame": frame,
                                    "length": orf_len, "seq": nuc[i:j+3]
                                })
                            break
                i += 3
    return sorted(orfs, key=lambda x: -x["length"])

from Bio import Align

aligner = Align.PairwiseAligner()
aligner.mode = "global"        # Needleman-Wunsch
# aligner.mode = "local"       # Smith-Waterman

aligner.match_score = 2
aligner.mismatch_score = -1
aligner.open_gap_score = -2
aligner.extend_gap_score = -0.5

seq1 = Seq("ATGCGATCG")
seq2 = Seq("ATGCAATCG")

alignments = aligner.align(seq1, seq2)
best = alignments[0]

print(f"Pontuação: {best.score}")
print(best)

from Bio.Blast import NCBIWWW, NCBIXML

result_handle = NCBIWWW.qblast(
    "blastn",          # programa: blastn, blastp, blastx, tblastn, tblastx
    "nt",              # banco de dados: nt, nr, refseq_rna, etc.
    "ATGCGATCGATCG",   # sequência query
    hitlist_size=10
)

blast_records = NCBIXML.parse(result_handle)
blast_record = next(blast_records)

for alignment in blast_record.alignments:
    for hsp in alignment.hsps:
        print(f"Match: {alignment.title[:60]}")
        print(f"Pontuação: {hsp.score}, E-valor: {hsp.expect:.2e}")
        print(f"Identidade: {hsp.identities}/{hsp.align_length}")
        break

from Bio import AlignIO
from collections import Counter

# Analisar um arquivo de alinhamento Clustal
alignment = AlignIO.read("sequences.aln", "clustal")

print(f"Comprimento do alinhamento: {alignment.get_alignment_length()}")
print(f"Número de sequências: {len(alignment)}")

# Obter uma coluna do alinhamento
col_10 = alignment[:, 10]   # coluna 10, todas as sequências
print(f"Coluna 10: {col_10}")

# Calcular conservação por posição
def column_identity(alignment, position: int) -> float:
    col = alignment[:, position]
    most_common = Counter(col).most_common(1)[0][1]
    return most_common / len(col)

from Bio.PDB import PDBParser
import urllib.request

# Baixar e analisar
urllib.request.urlretrieve("https://files.rcsb.org/download/1JM7.pdb", "1JM7.pdb")
parser = PDBParser(QUIET=True)
structure = parser.get_structure("1JM7", "1JM7.pdb")

# Navegar pela hierarquia
model = structure[0]              # primeiro MODEL
chain = model["A"]                # cadeia A
residue = chain[100]              # resíduo 100
atom = residue["CA"]              # carbono alfa
print(f"Coordenadas CA: {atom.get_vector()}")

# Todos os resíduos na cadeia A
residues = [r for r in chain.get_residues() if r.get_id()[0] == " "]
print(f"Resíduos padrão na cadeia A: {len(residues)}")

def get_ca_distance(res1, res2) -> float:
    ca1 = res1["CA"].get_vector()
    ca2 = res2["CA"].get_vector()
    return (ca1 - ca2).norm()

res10 = chain[10]
res50 = chain[50]
dist = get_ca_distance(res10, res50)
print(f"Distância Cα: {dist:.2f} Å")

from Bio import Entrez, SeqIO
from Bio.SeqUtils import gc_fraction
from Bio.Seq import Seq
import time

Entrez.email = "seu@email.com"

def analyze_gene(gene_name: str, organism: str = "Homo sapiens") -> dict:
    handle = Entrez.esearch(
        db="nucleotide",
        term=f"{gene_name}[Gene Name] AND {organism}[Organism] AND RefSeq[Filter] AND mRNA[Filter]",
        retmax=1
    )
    result = Entrez.read(handle)
    handle.close()

    if not result["IdList"]:
        return {"error": "Gene não encontrado"}

    time.sleep(0.4)
    handle = Entrez.efetch(db="nucleotide", id=result["IdList"][0], rettype="gb", retmode="text")
    record = SeqIO.read(handle, "genbank")
    handle.close()

    cds_features = [f for f in record.features if f.type == "CDS"]
    analysis = {
        "accession": record.id,
        "total_length_bp": len(record.seq),
        "gc_content": gc_fraction(record.seq),
        "cds_count": len(cds_features),
    }

    if cds_features:
        cds_seq = cds_features[0].extract(record.seq)
        protein = cds_seq.translate(to_stop=True)
        analysis["protein_length_aa"] = len(protein)
        analysis["protein_gc"] = gc_fraction(cds_seq)

    return analysis

result = analyze_gene("BRCA1")
for key, value in result.items():
    print(f"{key}: {value}")