Python for Bioinformatics

Bioinformatics Tutorials & Code Walkthroughs

Practical, hands-on tutorials crafted for researchers, students, and computational biologists to build real-world computational workflows.

1. Parsing FASTA and FASTQ Sequences with Biopython

from Bio import SeqIO

# Calculate average GC content and length distribution across FASTA
fasta_file = "sequences.fasta"
total_length = 0
gc_count = 0
sequence_count = 0

for record in SeqIO.parse(fasta_file, "fasta"):
    seq = record.seq.upper()
    seq_len = len(seq)
    total_length += seq_len
    gc_count += seq.count("G") + seq.count("C")
    sequence_count += 1
    print(f"ID: {record.id} | Length: {seq_len} bp | GC: {(seq.count('G')+seq.count('C'))/seq_len*100:.2f}%")

avg_gc = (gc_count / total_length) * 100 if total_length > 0 else 0
print(f"\nTotal Sequences: {sequence_count}")
print(f"Mean GC Content: {avg_gc:.2f}%")

2. Filtering VCF Variants by Depth & Quality with Python

import gzip

def filter_vcf(input_vcf, output_vcf, min_qual=30, min_dp=10):
    open_func = gzip.open if input_vcf.endswith(".gz") else open
    with open_func(input_vcf, "rt") as infile, open(output_vcf, "w") as outfile:
        for line in infile:
            if line.startswith("#"):
                outfile.write(line)
                continue
            fields = line.strip().split("\t")
            qual = float(fields[5]) if fields[5] != "." else 0.0
            info = fields[7]
            
            # Extract DP from INFO
            dp = 0
            for tag in info.split(";"):
                if tag.startswith("DP="):
                    dp = int(tag.split("=")[1])
                    break
                    
            if qual >= min_qual and dp >= min_dp:
                outfile.write(line)

filter_vcf("raw_variants.vcf.gz", "filtered_high_confidence.vcf")