8: Data Analysis and Working with Data in Python

Chapter 8: Data Analysis and Working with Data in Python

Next-Generation Sequencing (NGS) produces massive volumes of structured tables: sample manifests, gene quantification matrices, variant call sets, and clinical outcome registries. Mastering data cleaning, regex pattern extraction, and robust tabular file manipulation is fundamental to bioinformatics.

1. Working with Complex Biological Datasets

Biological datasets frequently feature multi-delimiter structures, headers with metadata (such as VCF header lines starting with ##), and missing values. Processing them requires clean data pipelines.

2. Importing and Exporting Data

Streaming large genomic tables line-by-line prevents out-of-memory errors on large sequencing datasets:

# Memory-efficient streaming of large gene expression tables
def process_expression_data(input_file, output_file, min_cpm=1.0):
    with open(input_file, "r") as infile, open(output_file, "w") as outfile:
        header = infile.readline()
        outfile.write(header)
        
        passed_genes = 0
        for line in infile:
            parts = line.strip().split("\t")
            gene_id = parts[0]
            expression_values = [float(x) for x in parts[1:]]
            
            # Keep genes expressed above threshold in at least 50% of samples
            if sum(x >= min_cpm for x in expression_values) >= len(expression_values) / 2:
                outfile.write(line)
                passed_genes += 1
                
    print(f"Filtering Complete: {passed_genes} genes retained.")

process_expression_data("raw_counts.tsv", "filtered_expression.tsv")

3. Cleaning and Preprocessing Data

Genomic datasets frequently contain incomplete entries, duplicate identifiers, or inconsistent chromosome names (e.g. chr1 vs 1). Normalizing these ensures downstream compatibility with tools like GATK and IGV.

4. Regular Expressions in Biology (re module)

Regular expressions allow precise pattern matching across nucleotide motifs, restriction enzyme sites, and identifier parsing.

import re

# 1. Matching degenerate IUPAC motifs (EcoRI: GAATTC)
ecori_pattern = re.compile(r"GAATTC", re.IGNORECASE)
sequence = "ATGCGAATTCAGCTAGAATTCGG"
for match in ecori_pattern.finditer(sequence):
    print(f"EcoRI Cleavage Site at position {match.start()}..{match.end()}")

# 2. Extracting Depth (DP) and Allele Frequency (AF) from VCF INFO field
info_field = "AC=2;AF=0.048;AN=42;DP=182;FS=1.2;MQ=59.4"
dp_match = re.search(r"DP=(\d+)", info_field)
af_match = re.search(r"AF=([0-9.]+)", info_field)

depth = int(dp_match.group(1)) if dp_match else 0
af = float(af_match.group(1)) if af_match else 0.0
print(f"Variant Metrics -> Read Depth: {depth}x | Allele Frequency: {af*100:.2f}%")

5. Manipulating Data

Filtering, sorting, and combining patient cohorts by somatic mutation status and sequencing depth thresholds.

6. Plotting Data Basics

Previewing sequencing coverage distributions, quality curves, and insert size distributions using standard Python scripting.