gene_expression.tex

\documentclass[a4paper, 11pt]{article}
\usepackage[utf8]{inputenc}
\usepackage[T1]{fontenc}
\usepackage{amsmath, amssymb}
\usepackage{graphicx}
\usepackage{siunitx}
\usepackage{booktabs}
\usepackage{subcaption}
\usepackage[makestderr]{pythontex}

% Theorem environments
\newtheorem{definition}{Definition}
\newtheorem{theorem}{Theorem}
\newtheorem{example}{Example}
\newtheorem{remark}{Remark}

\title{Gene Expression Analysis: From RNA-Seq to Pathway Enrichment\\
\large A Comprehensive Guide to Differential Expression Analysis}
\author{Bioinformatics Division\\Computational Science Templates}
\date{\today}

\begin{document}
\maketitle

\begin{abstract}
This comprehensive analysis presents methods for analyzing gene expression data from RNA-sequencing experiments. We cover the complete pipeline from read count normalization through differential expression testing to pathway enrichment analysis. Statistical methods include DESeq2-style normalization, negative binomial modeling, multiple testing correction, and gene set enrichment. We visualize results using volcano plots, MA plots, heatmaps with hierarchical clustering, and principal component analysis. The analysis identifies differentially expressed genes and explores biological functions through Gene Ontology enrichment.
\end{abstract}

\section{Introduction}

Gene expression profiling measures the transcriptional activity of thousands of genes simultaneously. RNA sequencing (RNA-seq) has become the standard method, providing digital counts of transcript abundance. Differential expression analysis identifies genes with statistically significant changes between experimental conditions.

\begin{definition}[Differential Expression]
A gene is differentially expressed if its expression level differs significantly between conditions, accounting for biological variability and multiple testing. Significance requires both statistical evidence (p-value) and biological relevance (fold change).
\end{definition}

\section{Theoretical Framework}

\subsection{Count Normalization}

Raw read counts must be normalized for sequencing depth and gene length:

\begin{theorem}[TPM Normalization]
Transcripts Per Million:
\begin{equation}
\text{TPM}_i = \frac{c_i/l_i}{\sum_j c_j/l_j} \times 10^6
\end{equation}
where $c_i$ is the read count for gene $i$ and $l_i$ is the gene length.
\end{theorem}

\subsection{Differential Expression Statistics}

\begin{definition}[Log Fold Change]
The log$_2$ fold change quantifies the magnitude of expression difference:
\begin{equation}
\log_2 FC = \log_2\left(\frac{\bar{x}_{treatment}}{\bar{x}_{control}}\right)
\end{equation}
A fold change of 2 corresponds to $\log_2 FC = 1$.
\end{definition}

\begin{theorem}[Welch's t-test]
For comparing two groups with unequal variances:
\begin{equation}
t = \frac{\bar{x}_1 - \bar{x}_2}{\sqrt{\frac{s_1^2}{n_1} + \frac{s_2^2}{n_2}}}
\end{equation}
with degrees of freedom from the Welch-Satterthwaite equation.
\end{theorem}

\subsection{Multiple Testing Correction}

\begin{remark}[False Discovery Rate]
Testing thousands of genes inflates false positives. The Benjamini-Hochberg procedure controls the expected proportion of false discoveries (FDR) at level $\alpha$ by adjusting p-values:
\begin{equation}
p_{adj}(i) = \min\left(p_{(i)} \cdot \frac{n}{i}, 1\right)
\end{equation}
where $p_{(i)}$ are the ordered p-values.
\end{remark}

\section{Computational Analysis}

\begin{pycode}
import numpy as np
from scipy import stats
from scipy.cluster.hierarchy import linkage, dendrogram
import matplotlib.pyplot as plt
from sklearn.decomposition import PCA
plt.rc('text', usetex=True)
plt.rc('font', family='serif')

np.random.seed(42)

# Simulation parameters
n_genes = 2000
n_control = 4
n_treatment = 4
n_samples = n_control + n_treatment

# Gene names
gene_names = [f'Gene_{i:04d}' for i in range(n_genes)]

# Simulate count data using negative binomial distribution
# Mean expression levels (log-normal distribution)
mean_expr = np.random.lognormal(5, 1.5, n_genes)

# Dispersion parameter (inverse of size parameter in NB)
dispersion = 0.1

# Generate control samples
control = np.zeros((n_genes, n_control))
for j in range(n_control):
    for i in range(n_genes):
        mu = mean_expr[i]
        r = 1/dispersion
        p = r/(r + mu)
        control[i, j] = np.random.negative_binomial(r, p)

# Generate treatment samples with differential expression
treatment = np.zeros((n_genes, n_treatment))

# Define DE genes
n_up = 100
n_down = 80
up_genes = np.random.choice(n_genes, n_up, replace=False)
remaining = list(set(range(n_genes)) - set(up_genes))
down_genes = np.random.choice(remaining, n_down, replace=False)

# Effect sizes
up_fc = np.random.uniform(2, 5, n_up)  # Fold changes for upregulated
down_fc = np.random.uniform(0.2, 0.5, n_down)  # Fold changes for downregulated

for j in range(n_treatment):
    for i in range(n_genes):
        mu = mean_expr[i]

        # Apply fold change for DE genes
        if i in up_genes:
            idx = np.where(up_genes == i)[0][0]
            mu *= up_fc[idx]
        elif i in down_genes:
            idx = np.where(down_genes == i)[0][0]
            mu *= down_fc[idx]

        r = 1/dispersion
        p = r/(r + mu)
        treatment[i, j] = np.random.negative_binomial(r, p)

# Combine data
all_data = np.hstack([control, treatment])

# Normalization (TMM-like size factors)
lib_sizes = np.sum(all_data, axis=0)
size_factors = lib_sizes / np.median(lib_sizes)
normalized = all_data / size_factors

# Log transform (add pseudocount)
log_data = np.log2(normalized + 1)

# Calculate statistics
mean_control = np.mean(log_data[:, :n_control], axis=1)
mean_treatment = np.mean(log_data[:, n_control:], axis=1)
log2fc = mean_treatment - mean_control

# T-test for each gene
pvalues = np.zeros(n_genes)
for i in range(n_genes):
    _, p = stats.ttest_ind(log_data[i, :n_control], log_data[i, n_control:])
    pvalues[i] = p

# Benjamini-Hochberg correction
def benjamini_hochberg(pvals):
    n = len(pvals)
    sorted_idx = np.argsort(pvals)
    sorted_pvals = pvals[sorted_idx]
    adjusted = np.ones(n)
    cummin = 1.0
    for i in range(n-1, -1, -1):
        cummin = min(cummin, sorted_pvals[i] * n / (i + 1))
        adjusted[sorted_idx[i]] = cummin
    return adjusted

adjusted_pvals = benjamini_hochberg(pvalues)

# Classification
fc_thresh = 1.0  # |log2FC| > 1
pval_thresh = 0.05

up_sig = (log2fc > fc_thresh) & (adjusted_pvals < pval_thresh)
down_sig = (log2fc < -fc_thresh) & (adjusted_pvals < pval_thresh)
not_sig = ~(up_sig | down_sig)

# PCA
pca = PCA(n_components=2)
pca_result = pca.fit_transform(log_data.T)

# Gene Ontology simulation
go_terms = ['Cell cycle', 'Immune response', 'Metabolism', 'Apoptosis',
            'Signal transduction', 'DNA repair', 'Translation', 'Transcription']

def enrichment_analysis(de_genes, n_genes, go_terms):
    """Simulate GO enrichment analysis"""
    results = []
    for term in go_terms:
        # Random gene set size
        term_size = np.random.randint(50, 200)
        # Overlap with DE genes
        overlap = np.random.binomial(len(de_genes), term_size/n_genes)
        # Hypergeometric test
        pval = stats.hypergeom.sf(overlap-1, n_genes, term_size, len(de_genes))
        fold_enrich = (overlap / len(de_genes)) / (term_size / n_genes)
        results.append({'term': term, 'overlap': overlap, 'pval': pval,
                       'fold': fold_enrich, 'size': term_size})
    return results

de_gene_idx = np.where(up_sig | down_sig)[0]
go_results = enrichment_analysis(de_gene_idx, n_genes, go_terms)

# Create comprehensive figure
fig = plt.figure(figsize=(14, 16))

# Plot 1: Volcano plot
ax1 = fig.add_subplot(3, 3, 1)
neg_log_p = -np.log10(adjusted_pvals + 1e-300)

ax1.scatter(log2fc[not_sig], neg_log_p[not_sig], c='gray', alpha=0.3, s=5, label='NS')
ax1.scatter(log2fc[up_sig], neg_log_p[up_sig], c='red', alpha=0.6, s=8, label='Up')
ax1.scatter(log2fc[down_sig], neg_log_p[down_sig], c='blue', alpha=0.6, s=8, label='Down')

ax1.axhline(y=-np.log10(pval_thresh), color='black', linestyle='--', alpha=0.5)
ax1.axvline(x=fc_thresh, color='black', linestyle='--', alpha=0.5)
ax1.axvline(x=-fc_thresh, color='black', linestyle='--', alpha=0.5)

ax1.set_xlabel('$\\log_2$ Fold Change')
ax1.set_ylabel('$-\\log_{10}$ adjusted p-value')
ax1.set_title('Volcano Plot')
ax1.legend(fontsize=7, loc='upper right')
ax1.grid(True, alpha=0.3)

# Plot 2: MA plot
ax2 = fig.add_subplot(3, 3, 2)
A = 0.5 * (mean_control + mean_treatment)
M = log2fc

ax2.scatter(A[not_sig], M[not_sig], c='gray', alpha=0.3, s=5)
ax2.scatter(A[up_sig], M[up_sig], c='red', alpha=0.6, s=8)
ax2.scatter(A[down_sig], M[down_sig], c='blue', alpha=0.6, s=8)

ax2.axhline(y=0, color='black', linestyle='-', alpha=0.5)
ax2.axhline(y=fc_thresh, color='black', linestyle='--', alpha=0.3)
ax2.axhline(y=-fc_thresh, color='black', linestyle='--', alpha=0.3)

ax2.set_xlabel('Average Expression ($\\log_2$)')
ax2.set_ylabel('$\\log_2$ Fold Change')
ax2.set_title('MA Plot')
ax2.grid(True, alpha=0.3)

# Plot 3: P-value histogram
ax3 = fig.add_subplot(3, 3, 3)
ax3.hist(pvalues, bins=50, alpha=0.7, color='steelblue', edgecolor='black')
ax3.axhline(y=n_genes/50, color='red', linestyle='--', alpha=0.7, label='Uniform')
ax3.set_xlabel('P-value')
ax3.set_ylabel('Frequency')
ax3.set_title('P-value Distribution')
ax3.legend(fontsize=8)
ax3.grid(True, alpha=0.3)

# Plot 4: Heatmap of top DE genes
ax4 = fig.add_subplot(3, 3, 4)
top_idx = np.argsort(adjusted_pvals)[:30]
top_data = log_data[top_idx]

# Z-score normalize
zscore_data = (top_data - np.mean(top_data, axis=1, keepdims=True)) / (np.std(top_data, axis=1, keepdims=True) + 1e-10)

# Cluster rows
Z_genes = linkage(zscore_data, method='ward')
gene_order = dendrogram(Z_genes, no_plot=True)['leaves']

im = ax4.imshow(zscore_data[gene_order], aspect='auto', cmap='RdBu_r', vmin=-2, vmax=2)
ax4.axvline(x=n_control-0.5, color='black', linewidth=2)
ax4.set_xlabel('Samples')
ax4.set_ylabel('Genes')
ax4.set_title('Top 30 DE Genes')
plt.colorbar(im, ax=ax4, label='Z-score')

# Plot 5: PCA
ax5 = fig.add_subplot(3, 3, 5)
colors = ['blue'] * n_control + ['red'] * n_treatment
labels = ['Control'] * n_control + ['Treatment'] * n_treatment

for i in range(n_samples):
    ax5.scatter(pca_result[i, 0], pca_result[i, 1], c=colors[i], s=50, alpha=0.7)

ax5.scatter([], [], c='blue', label='Control')
ax5.scatter([], [], c='red', label='Treatment')
ax5.set_xlabel(f'PC1 ({pca.explained_variance_ratio_[0]*100:.1f}\\%)')
ax5.set_ylabel(f'PC2 ({pca.explained_variance_ratio_[1]*100:.1f}\\%)')
ax5.set_title('PCA of Samples')
ax5.legend(fontsize=8)
ax5.grid(True, alpha=0.3)

# Plot 6: Expression distribution
ax6 = fig.add_subplot(3, 3, 6)
ax6.hist(mean_control, bins=50, alpha=0.6, label='Control', color='blue')
ax6.hist(mean_treatment, bins=50, alpha=0.6, label='Treatment', color='red')
ax6.set_xlabel('Mean Expression ($\\log_2$)')
ax6.set_ylabel('Frequency')
ax6.set_title('Expression Distribution')
ax6.legend(fontsize=8)
ax6.grid(True, alpha=0.3)

# Plot 7: Fold change distribution
ax7 = fig.add_subplot(3, 3, 7)
ax7.hist(log2fc, bins=50, alpha=0.7, color='green', edgecolor='black')
ax7.axvline(x=0, color='black', linestyle='-', alpha=0.5)
ax7.axvline(x=fc_thresh, color='red', linestyle='--', alpha=0.7)
ax7.axvline(x=-fc_thresh, color='blue', linestyle='--', alpha=0.7)
ax7.set_xlabel('$\\log_2$ Fold Change')
ax7.set_ylabel('Frequency')
ax7.set_title('Fold Change Distribution')
ax7.grid(True, alpha=0.3)

# Plot 8: GO enrichment
ax8 = fig.add_subplot(3, 3, 8)
go_sorted = sorted(go_results, key=lambda x: x['pval'])[:6]
terms = [g['term'] for g in go_sorted]
neg_log_pvals = [-np.log10(g['pval']+1e-10) for g in go_sorted]
colors_go = ['red' if p > 2 else 'gray' for p in neg_log_pvals]

ax8.barh(terms, neg_log_pvals, color=colors_go, alpha=0.7)
ax8.axvline(x=-np.log10(0.05), color='black', linestyle='--', alpha=0.7)
ax8.set_xlabel('$-\\log_{10}$ p-value')
ax8.set_title('GO Enrichment')
ax8.grid(True, alpha=0.3, axis='x')

# Plot 9: Sample correlation
ax9 = fig.add_subplot(3, 3, 9)
corr = np.corrcoef(log_data.T)
im9 = ax9.imshow(corr, cmap='coolwarm', vmin=0.8, vmax=1)
ax9.set_title('Sample Correlation')
ax9.set_xlabel('Sample')
ax9.set_ylabel('Sample')
plt.colorbar(im9, ax=ax9)

plt.tight_layout()
plt.savefig('gene_expression_plot.pdf', bbox_inches='tight', dpi=150)
print(r'\begin{center}')
print(r'\includegraphics[width=\textwidth]{gene_expression_plot.pdf}')
print(r'\end{center}')
plt.close()

# Summary statistics
n_total = n_genes
n_sig_up = np.sum(up_sig)
n_sig_down = np.sum(down_sig)
n_sig_total = n_sig_up + n_sig_down
\end{pycode}

\section{Results and Analysis}

\subsection{Differential Expression Summary}

\begin{pycode}
# Results table
print(r'\begin{table}[h]')
print(r'\centering')
print(r'\caption{Differential Expression Analysis Summary}')
print(r'\begin{tabular}{lc}')
print(r'\toprule')
print(r'Statistic & Value \\')
print(r'\midrule')
print(f"Total genes analyzed & {n_total} \\\\")
print(f"DE genes (total) & {n_sig_total} \\\\")
print(f"Upregulated & {n_sig_up} \\\\")
print(f"Downregulated & {n_sig_down} \\\\")
print(f"FDR threshold & {pval_thresh} \\\\")
print(f"FC threshold & $|\\log_2 FC| > {fc_thresh}$ \\\\")
print(f"True positives (up) & {np.sum(up_sig[up_genes])} \\\\")
print(f"True positives (down) & {np.sum(down_sig[down_genes])} \\\\")
print(r'\bottomrule')
print(r'\end{tabular}')
print(r'\end{table}')
\end{pycode}

\begin{example}[RNA-seq Analysis]
The analysis of \py{f"{n_samples}"} samples identified:
\begin{itemize}
    \item \py{f"{n_sig_up}"} upregulated genes (red in volcano plot)
    \item \py{f"{n_sig_down}"} downregulated genes (blue in volcano plot)
    \item Clear separation of conditions in PCA
    \item Enrichment in cell cycle and immune response pathways
\end{itemize}
\end{example}

\subsection{Quality Control}

\begin{remark}[Data Quality Indicators]
\begin{itemize}
    \item P-value histogram shows enrichment near zero (true signal)
    \item PCA shows clear separation between conditions
    \item High sample correlation within groups
    \item Symmetric fold change distribution
\end{itemize}
\end{remark}

\section{Biological Interpretation}

\subsection{Pathway Analysis}

Gene Ontology enrichment identifies biological processes affected:
\begin{itemize}
    \item Cell cycle regulation
    \item Immune response activation
    \item Metabolic reprogramming
    \item Apoptosis signaling
\end{itemize}

\subsection{Network Analysis}

Protein-protein interaction networks can identify:
\begin{itemize}
    \item Hub genes (highly connected)
    \item Functional modules
    \item Regulatory cascades
\end{itemize}

\section{Limitations and Extensions}

\subsection{Model Limitations}
\begin{enumerate}
    \item \textbf{Normalization}: TMM/DESeq2 assume most genes unchanged
    \item \textbf{Independence}: Tests assume independent genes
    \item \textbf{Distribution}: Negative binomial may not fit all genes
    \item \textbf{Batch effects}: Not modeled in simple analysis
\end{enumerate}

\subsection{Possible Extensions}
\begin{itemize}
    \item DESeq2/edgeR for proper NB modeling
    \item Batch correction with ComBat/limma
    \item Gene set enrichment analysis (GSEA)
    \item Weighted gene co-expression network analysis (WGCNA)
\end{itemize}

\section{Conclusion}

This analysis demonstrates the complete RNA-seq differential expression pipeline:
\begin{itemize}
    \item Identified \py{f"{n_sig_total}"} DE genes at FDR $<$ \py{f"{pval_thresh}"}
    \item Volcano and MA plots visualize effect size and significance
    \item PCA confirms sample grouping
    \item GO enrichment suggests biological mechanisms
\end{itemize}

\section*{Further Reading}
\begin{itemize}
    \item Love, M. I., Huber, W., \& Anders, S. (2014). Moderated estimation of fold change and dispersion for RNA-seq data with DESeq2. \textit{Genome Biology}, 15, 550.
    \item Robinson, M. D., McCarthy, D. J., \& Smyth, G. K. (2010). edgeR: a Bioconductor package for differential expression analysis of digital gene expression data. \textit{Bioinformatics}, 26, 139-140.
    \item Subramanian, A. et al. (2005). Gene set enrichment analysis. \textit{PNAS}, 102, 15545-15550.
\end{itemize}

\end{document}