Python for Bioinformatics: Essential Biopython Recipes for FASTA, GenBank, and BLAST Parsing

>NC_000913.3 Escherichia coli str. K-12 substr. MG1655, complete genome
AGCTTTTCATTCTGACTGCAACGGGCAATATGTCTCTGTGTGGATTAAAAAAAGAGTGTCTGATAGCAGC
TTCTGAACTGGTTACCTGCCGTGAGTAAATTAAAATTTTATTGACTTAGGTCACTAAATACTTTAACCAA

LOCUS       NC_000913            4641652 bp    DNA     circular BCT 24-MAY-2024
DEFINITION  Escherichia coli str. K-12 substr. MG1655, complete genome.
ACCESSION   NC_000913
VERSION     NC_000913.3
FEATURES             Location/Qualifiers
     source          1..4641652
                     /organism="Escherichia coli str. K-12 substr. MG1655"
     gene            190..255
                     /gene="thrL"
                     /locus_tag="b0001"
     CDS             190..255
                     /gene="thrL"
                     /locus_tag="b0001"
                     /codon_start=1
                     /transl_table=11
                     /product="thr operon leader peptide"
                     /translation="MKRISTTITTTITITTGTLARK"

[FASTA File] --------> [SeqIO.parse()] --------> [GC Skew / Filtering]
                                                        |
[GenBank File] ------> [SeqFeature Extraction] -> [CDS Translation Table 11]
                                                        |
[NCBI Entrez] -------> [EFetch / ESearch] -----> [SearchIO BLAST Parsing]

from Bio import SeqIO
from Bio.SeqUtils import gc_fraction

def stream_and_filter_fasta(input_fasta, output_fasta, min_len=500, min_gc=0.45):
    """
    Streams a FASTA file record-by-record, filtering by length and GC content.
    Memory footprint remains constant (under 50 MB) regardless of file size.
    """
    retained_count = 0
    total_count = 0

    with open(output_fasta, "w") as out_handle:
        # SeqIO.parse returns a generator
        for record in SeqIO.parse(input_fasta, "fasta"):
            total_count += 1
            seq_len = len(record.seq)
            gc_val = gc_fraction(record.seq) # Returns float between 0.0 and 1.0

            if seq_len >= min_len and gc_val >= min_gc:
                # Modify header to record QC metrics
                record.description = f"{record.description} | len={seq_len} gc={gc_val:.3f}"
                SeqIO.write(record, out_handle, "fasta")
                retained_count += 1

    print(f"Processed {total_count} records. Retained {retained_count} records ({retained_count/total_count*100:.1f}%).")

# Example execution
stream_and_filter_fasta("unfiltered_contigs.fasta", "qc_filtered_contigs.fasta")

# Create an on-disk index file (zero RAM consumption)
record_dict = SeqIO.index_db("contigs.idx", ["unfiltered_contigs.fasta"], "fasta")

# Retrieve any sequence instantaneously in O(1) time
target_record = record_dict["contig_10492"]
print(f"Accession: {target_record.id}, Length: {len(target_record.seq)}")
record_dict.close()

from Bio.Seq import Seq

dna = Seq("ATGCGTACCGGTTTCGACGACTAG")

# 1. Reverse Complement
rev_comp = dna.reverse_complement()
print("Reverse Complement:", rev_comp) # CTAGTCGTCGAAACCGGTACGCAT

# 2. In Silico Transcription (DNA -> mRNA)
mrna = dna.transcribe()
print("mRNA:", mrna) # AUGCGUACCGGUUUCGACGACUAG

# 3. Translation with Genetic Code Tables
# Table 1: Standard Eukaryotic Nuclear | Table 11: Bacterial, Archaeal, Plant Plastid
protein = dna.translate(table=11, to_stop=True)
print("Translated Peptide (Table 11):", protein) # MRTVFDD

def find_orfs(sequence, min_protein_len=50, table_id=11):
    """
    Discovers all candidate open reading frames (ORFs) across
    all 6 reading frames (3 forward, 3 reverse).
    """
    orfs = []
    seq_len = len(sequence)

    # Iterate forward (strand +1) and reverse complement (strand -1)
    for strand, nuc_seq in [(+1, sequence), (-1, sequence.reverse_complement())]:
        for frame in range(3):
            trans_seq = str(nuc_seq[frame:].translate(table=table_id))
            trans_len = len(trans_seq)
            start_pos = 0

            while start_pos < trans_len:
                start_codon = trans_seq.find("M", start_pos)
                if start_codon == -1:
                    break
                stop_codon = trans_seq.find("*", start_codon)
                if stop_codon == -1:
                    break

                peptide_len = stop_codon - start_codon
                if peptide_len >= min_protein_len:
                    peptide = trans_seq[start_codon:stop_codon]
                    orfs.append({
                        "strand": strand,
                        "frame": frame + 1,
                        "length_aa": peptide_len,
                        "peptide": peptide
                    })
                start_pos = stop_codon + 1

    return sorted(orfs, key=lambda x: x["length_aa"], reverse=True)

# Test ORF Finder
test_seq = Seq("ATGAAGCGTATTTCTACCACCATTACTACCACCATCACCATTACCACCGGTACTCTGGCTCGTAACTGA")
print("Identified ORFs:", find_orfs(test_seq, min_protein_len=10))

from Bio import SeqIO
import pandas as pd

def extract_genbank_features(gbk_path):
    """
    Parses a GenBank flat file and extracts all annotated CDS features,
    automatically resolving spliced exons, locus tags, and protein products.
    """
    extracted_data = []

    for record in SeqIO.parse(gbk_path, "genbank"):
        chromosome_id = record.id

        for feature in record.features:
            if feature.type == "CDS":
                # Extract qualifiers safely using .get()
                locus_tag = feature.qualifiers.get("locus_tag", ["N/A"])[0]
                gene_name = feature.qualifiers.get("gene", ["N/A"])[0]
                product = feature.qualifiers.get("product", ["Hypothetical Protein"])[0]

                # feature.extract handles single intervals and spliced CompoundLocations
                cds_nucleotides = feature.extract(record.seq)

                # Check for documented translation or translate directly
                translation = feature.qualifiers.get("translation", [None])[0]
                if not translation:
                    translation = str(cds_nucleotides.translate(table=11, to_stop=True))

                extracted_data.append({
                    "Chromosome": chromosome_id,
                    "Locus_Tag": locus_tag,
                    "Gene": gene_name,
                    "Start": int(feature.location.start) + 1, # Convert 0-based to 1-based
                    "End": int(feature.location.end),
                    "Strand": "+" if feature.location.strand == 1 else "-",
                    "CDS_Length_bp": len(cds_nucleotides),
                    "Protein_Length_aa": len(translation),
                    "Product": product,
                    "Sequence": translation
                })

    df = pd.DataFrame(extracted_data)
    print(f"Extracted {len(df)} annotated coding sequences.")
    return df

# Example usage
# df_genes = extract_genbank_features("ecoli_k12.gbk")
# df_genes.to_csv("annotated_cds_table.csv", index=False)

from Bio import Entrez
from Bio import SeqIO
import time

# Configure NCBI Entrez API credentials
Entrez.email = "researcher@bioinformaticsdaily.com" # Required by NCBI
Entrez.api_key = "YOUR_NCBI_API_KEY_OPTIONAL"        # Raises rate limit from 3 to 10 req/sec

def download_refseq_genomes(organism_query, max_records=5):
    """
    Searches the NCBI Nucleotide database for complete bacterial genomes
    and downloads corresponding annotated GenBank files.
    """
    print(f"Searching NCBI for: {organism_query}")

    # 1. ESearch: Execute text search query
    search_handle = Entrez.esearch(
        db="nucleotide",
        term=f"{organism_query}[Organism] AND complete genome[Title] AND refseq[filter]",
        retmax=max_records,
        sort="relevance"
    )
    search_results = Entrez.read(search_handle)
    search_handle.close()

    id_list = search_results["IdList"]
    print(f"Found {len(id_list)} matching records. IDs: {id_list}")

    records = []
    # 2. EFetch: Retrieve complete records in GenBank format
    for uid in id_list:
        print(f"Fetching record UID: {uid}...")
        fetch_handle = Entrez.efetch(
            db="nucleotide",
            id=uid,
            rettype="gbwithparts", # Complete GenBank record including full sequence
            retmode="text"
        )
        record = SeqIO.read(fetch_handle, "genbank")
        fetch_handle.close()
        records.append(record)

        # Respect NCBI rate limiting
        time.sleep(0.35)

    return records

# Example: Fetch Salmonella enterica reference genomes
# genomes = download_refseq_genomes("Salmonella enterica", max_records=2)
# print("Downloaded genome size:", len(genomes[0].seq), "bp")

from Bio.Blast.Applications import NcbiblastpCommandline
from Bio import SearchIO
import pandas as pd

def run_and_parse_blast(query_fasta, blast_db, evalue_threshold=1e-5):
    """
    Executes local blastp and parses results into a structured pandas DataFrame.
    """
    output_xml = "blast_results.xml"

    # 1. Construct and execute the blastp command
    blastp_cline = NcbiblastpCommandline(
        cmd="blastp",
        query=query_fasta,
        db=blast_db,
        evalue=evalue_threshold,
        outfmt=5, # XML output
        out=output_xml
    )
    print("Executing command:", blastp_cline)
    stdout, stderr = blastp_cline()

    # 2. Parse XML with Bio.SearchIO
    parsed_hits = []
    for query_result in SearchIO.parse(output_xml, "blast-xml"):
        for hit in query_result:
            for hsp in hit.hsps: # High-scoring Segment Pairs
                # Calculate percent identity
                pct_identity = (hsp.ident_num / hsp.aln_span) * 100

                parsed_hits.append({
                    "Query_ID": query_result.id,
                    "Query_Length": query_result.seq_len,
                    "Hit_ID": hit.id,
                    "Hit_Description": hit.description,
                    "E_Value": hsp.evalue,
                    "Bit_Score": hsp.bitscore,
                    "Identity_Pct": pct_identity,
                    "Alignment_Span": hsp.aln_span,
                    "Query_Start": hsp.query_start,
                    "Query_End": hsp.query_end,
                    "Hit_Start": hsp.hit_start,
                    "Hit_End": hsp.hit_end
                })

    df = pd.DataFrame(parsed_hits)
    return df

# Example schema demonstration
print("SearchIO parser ready for automated comparative pipelines.")

GC Skew = (G - C) / (G + C)

import numpy as np
import matplotlib.pyplot as plt
import pandas as pd

def calculate_and_plot_gc_skew(sequence_str, window_size=10000, step_size=2000):
    """
    Computes sliding-window GC Skew and cumulative GC skew across a chromosome,
    plotting replication origin (minimum cumulative skew) and terminus.
    """
    positions = []
    skew_values = []
    seq_len = len(sequence_str)

    for i in range(0, seq_len - window_size, step_size):
        subseq = sequence_str[i:i + window_size].upper()
        g = subseq.count("G")
        c = subseq.count("C")
        skew = (g - c) / (g + c) if (g + c) > 0 else 0

        positions.append(i + window_size // 2)
        skew_values.append(skew)

    positions = np.array(positions) / 1e6 # Convert to Megabases (Mb)
    skew_values = np.array(skew_values)
    cum_skew = np.cumsum(skew_values)

    # Generate publication multi-panel plot
    fig, (ax1, ax2) = plt.subplots(2, 1, figsize=(12, 7), sharex=True)

    # Panel 1: Windowed GC Skew
    ax1.plot(positions, skew_values, color="#2563eb", linewidth=0.8, alpha=0.8)
    ax1.axhline(0, color="gray", linestyle="--", linewidth=0.8)
    ax1.set_ylabel("Windowed GC Skew", fontsize=11)
    ax1.set_title("Chromosomal GC Skew Analysis (Window = 10 kb)", fontsize=13)
    ax1.grid(True, linestyle=":", alpha=0.6)

    # Panel 2: Cumulative GC Skew (Inflection points indicate origin/terminus)
    ax2.plot(positions, cum_skew, color="#dc2626", linewidth=1.5)
    ori_index = np.argmin(cum_skew)
    ax2.axvline(positions[ori_index], color="black", linestyle="--", label=f"Predicted oriC (~{positions[ori_index]:.2f} Mb)")
    ax2.set_xlabel("Genomic Position (Mb)", fontsize=11)
    ax2.set_ylabel("Cumulative GC Skew", fontsize=11)
    ax2.legend(loc="upper right", frameon=True)
    ax2.grid(True, linestyle=":", alpha=0.6)

    plt.tight_layout()
    plt.savefig("GC_Skew_Analysis.pdf", dpi=300)
    plt.show()

# Example: Generate simulated chromosome visualization
simulated_genome = "ATGCGTACCGGT" * 100000
calculate_and_plot_gc_skew(simulated_genome, window_size=5000, step_size=1000)