11: Analyzing Biological Data with Python

Chapter 11: Analyzing Biological Data with Python

Translating raw multi-omics measurements into biological insight requires rigorous quantitative pipelines. This chapter details transcriptomic normalization (TPM, RPKM), dimensionality reduction (PCA, t-SNE), protein-protein interaction network modeling, and false-discovery-rate statistical corrections.

1. Gene Expression Analysis & Normalization (TPM Calculation)

Transcripts Per Million (TPM) normalizes read counts for both gene length and sequencing depth:

import pandas as pd
import numpy as np

def calculate_tpm(counts_df, gene_lengths_kb):
    """Calculates Transcripts Per Million (TPM) from raw count matrix."""
    # 1. Normalize for gene length (Reads Per Kilobase)
    rpk = counts_df.div(gene_lengths_kb, axis=0)
    
    # 2. Normalize for library size (per million scaling factor)
    scaling_factors = rpk.sum(axis=0) / 1e6
    tpm = rpk.div(scaling_factors, axis=1)
    return tpm

# Example usage
counts = pd.DataFrame({"SampleA": [150, 4200, 80], "SampleB": [180, 5100, 95]}, index=["Gene1", "Gene2", "Gene3"])
lengths = pd.Series([1.2, 3.5, 0.8], index=["Gene1", "Gene2", "Gene3"])
tpm_matrix = calculate_tpm(counts, lengths)
print("TPM Expression Matrix:\n", tpm_matrix)

2. Clustering & PCA of Gene Expression Data

Principal Component Analysis (PCA) reveals cohort variance, batches, and biological groupings:

from sklearn.decomposition import PCA
from sklearn.preprocessing import StandardScaler

# Log-transform and standardize expression
log_tpm = np.log2(tpm_matrix + 1).T
scaled_data = StandardScaler().fit_transform(log_tpm)

pca = PCA(n_components=2)
principal_components = pca.fit_transform(scaled_data)
print(f"Variance Explained: PC1={pca.explained_variance_ratio_[0]*100:.1f}%, PC2={pca.explained_variance_ratio_[1]*100:.1f}%")

3. Network Analysis with Python (NetworkX)

Constructing Protein-Protein Interaction (PPI) networks and computing topological hub centralities:

import networkx as nx

# Create Protein-Protein Interaction network
ppi = nx.Graph()
ppi.add_edge("TP53", "MDM2", weight=0.98)
ppi.add_edge("TP53", "CDKN1A", weight=0.95)
ppi.add_edge("TP53", "BAX", weight=0.91)
ppi.add_edge("EGFR", "GRB2", weight=0.88)
ppi.add_edge("GRB2", "SOS1", weight=0.92)

# Degree centrality and betweenness
degrees = dict(ppi.degree())
betweenness = nx.betweenness_centrality(ppi)
print("Top Hub Proteins by Degree:", sorted(degrees.items(), key=lambda x: x[1], reverse=True))

4. Statistical Testing & FDR Corrections

Applying Student’s t-test and Benjamini-Hochberg False Discovery Rate (FDR) corrections across thousands of hypothesis tests:

from scipy.stats import ttest_ind
from statsmodels.stats.multitest import multipletests

# Compute p-values across 1,000 genes
p_values = [0.001, 0.04, 0.12, 0.0003, 0.05, 0.89]  # Sample p-values
reject, corrected_padj, _, _ = multipletests(p_values, alpha=0.05, method="fdr_bh")

for p, padj, is_sig in zip(p_values, corrected_padj, reject):
    print(f"Raw p={p:.4f} -> FDR padj={padj:.4f} | Significant: {is_sig}")