gene_expression_analysis.py

#!/usr/bin/env python3
# -*- coding: utf-8 -*-
"""
Differential Gene Expression and Pathway Analysis Pipeline
---------------------------------------------------------
This script provides a comprehensive framework for analyzing gene expression data,
identifying differentially expressed genes, and performing pathway enrichment analysis
with AI-powered interpretation of results.

Author: Cline
Date: December 4, 2025
"""

import os
import sys
import argparse
import logging
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
import seaborn as sns
from scipy import stats
import statsmodels.stats.multitest as smm
from sklearn.decomposition import PCA
from sklearn.preprocessing import StandardScaler
from Bio import Entrez
import gseapy as gp
import mygene
import json
import openai
from dotenv import load_dotenv

# Set up logging
logging.basicConfig(
    level=logging.INFO,
    format='%(asctime)s - %(name)s - %(levelname)s - %(message)s',
    handlers=[logging.StreamHandler()]
)
logger = logging.getLogger('gene_expression_analysis')

# Load environment variables for API keys
load_dotenv()


class GeneExpressionAnalyzer:
    """
    Class for analyzing gene expression data, performing differential expression analysis,
    and conducting pathway enrichment analysis.
    """
    
    def __init__(self, config=None):
        """
        Initialize the analyzer with configuration.
        
        Args:
            config (dict, optional): Configuration parameters for analysis
        """
        self.config = config or {}
        self.data = None
        self.metadata = None
        self.normalized_data = None
        self.diff_expr_results = None
        self.pathway_results = None
        self.openai_client = None
        
        # Set Entrez email for NCBI API
        if 'entrez_email' in self.config:
            Entrez.email = self.config['entrez_email']
        
        # Initialize OpenAI client if API key is available
        if os.environ.get('OPENAI_API_KEY'):
            try:
                from openai import OpenAI
                self.openai_client = OpenAI(api_key=os.environ.get('OPENAI_API_KEY'))
            except ImportError:
                logger.warning("OpenAI package not available. AI interpretation disabled.")
        else:
            logger.warning("OPENAI_API_KEY not found in environment variables. AI interpretation disabled.")
            
        # Initialize MyGene client for gene annotations
        self.mg = mygene.MyGeneInfo()
        
    def load_data(self, expression_file, metadata_file=None, format='csv', sample_col='sample_id', gene_id_col='gene_id'):
        """
        Load gene expression data and optional metadata.
        
        Args:
            expression_file (str): Path to expression data file
            metadata_file (str, optional): Path to metadata file
            format (str): File format ('csv', 'tsv', or 'excel')
            sample_col (str): Name of the column containing sample IDs
            gene_id_col (str): Name of the column containing gene IDs
            
        Returns:
            self: For method chaining
        """
        logger.info(f"Loading expression data from {expression_file}")
        
        try:
            if format.lower() == 'csv':
                self.data = pd.read_csv(expression_file)
            elif format.lower() == 'tsv':
                self.data = pd.read_csv(expression_file, sep='\t')
            elif format.lower() == 'excel':
                self.data = pd.read_excel(expression_file)
            else:
                raise ValueError(f"Unsupported format: {format}")
                
            # Set gene IDs as index if specified
            if gene_id_col in self.data.columns:
                self.data = self.data.set_index(gene_id_col)
            
            logger.info(f"Loaded expression data with {self.data.shape[0]} genes and {self.data.shape[1]} samples")
            
            # Load metadata if provided
            if metadata_file:
                logger.info(f"Loading metadata from {metadata_file}")
                if format.lower() == 'csv':
                    self.metadata = pd.read_csv(metadata_file)
                elif format.lower() == 'tsv':
                    self.metadata = pd.read_csv(metadata_file, sep='\t')
                elif format.lower() == 'excel':
                    self.metadata = pd.read_excel(metadata_file)
                
                logger.info(f"Loaded metadata for {self.metadata.shape[0]} samples")
                
                # Validate that sample IDs in metadata match expression data
                if sample_col in self.metadata.columns:
                    metadata_samples = set(self.metadata[sample_col])
                    data_samples = set(self.data.columns)
                    
                    if not metadata_samples.issubset(data_samples):
                        missing = metadata_samples - data_samples
                        logger.warning(f"Found {len(missing)} samples in metadata not present in expression data")
        
        except Exception as e:
            logger.error(f"Error loading data: {str(e)}")
            raise
            
        return self
    
    def preprocess_data(self, normalize_method='log2', filter_low_expr=True, min_samples=3):
        """
        Preprocess gene expression data.
        
        Args:
            normalize_method (str): Normalization method ('log2', 'z-score', None)
            filter_low_expr (bool): Whether to filter low-expression genes
            min_samples (int): Minimum number of samples with expression for gene filtering
            
        Returns:
            self: For method chaining
        """
        if self.data is None:
            logger.error("No data loaded. Please load data first.")
            return self
            
        logger.info("Preprocessing expression data")
        
        # Create a copy to avoid modifying original data
        self.normalized_data = self.data.copy()
        
        # Filter low-expression genes
        if filter_low_expr:
            logger.info(f"Filtering genes with low expression in less than {min_samples} samples")
            before_genes = self.normalized_data.shape[0]
            
            # Define a gene as expressed if value > 0
            expressed = (self.normalized_data > 0).sum(axis=1)
            self.normalized_data = self.normalized_data[expressed >= min_samples]
            
            after_genes = self.normalized_data.shape[0]
            logger.info(f"Filtered out {before_genes - after_genes} low-expression genes")
        
        # Apply normalization if specified
        if normalize_method:
            logger.info(f"Applying {normalize_method} normalization")
            
            if normalize_method == 'log2':
                # Add small constant to avoid log(0)
                self.normalized_data = np.log2(self.normalized_data + 1)
                
            elif normalize_method == 'z-score':
                # Z-score normalization (standardization) gene-wise
                self.normalized_data = StandardScaler().fit_transform(self.normalized_data.T).T
                self.normalized_data = pd.DataFrame(
                    self.normalized_data,
                    index=self.data.index,
                    columns=self.data.columns
                )
        
        logger.info(f"Preprocessing complete. Data shape: {self.normalized_data.shape}")
        return self
    
    def explore_data(self, output_dir='.', prefix=''):
        """
        Explore and visualize the data.
        
        Args:
            output_dir (str): Directory to save plots
            prefix (str): Prefix for output file names
            
        Returns:
            self: For method chaining
        """
        if self.normalized_data is None:
            logger.warning("No normalized data available. Using raw data.")
            data_to_explore = self.data
        else:
            data_to_explore = self.normalized_data
        
        if data_to_explore is None:
            logger.error("No data available for exploration.")
            return self
            
        logger.info("Exploring gene expression data")
        
        # Create output directory if it doesn't exist
        os.makedirs(output_dir, exist_ok=True)
        
        # 1. Distribution of expression values
        plt.figure(figsize=(10, 6))
        sns.histplot(data_to_explore.values.flatten(), kde=True)
        plt.title('Distribution of Gene Expression Values')
        plt.xlabel('Expression Value')
        plt.ylabel('Frequency')
        plt.tight_layout()
        plt.savefig(os.path.join(output_dir, f"{prefix}expression_distribution.png"))
        plt.close()
        
        # 2. Principal Component Analysis
        if data_to_explore.shape[1] >= 3:  # Need at least 3 samples for meaningful PCA
            logger.info("Performing PCA")
            
            # Transpose to have samples as rows
            pca = PCA(n_components=2)
            pca_result = pca.fit_transform(data_to_explore.T)
            variance_explained = pca.explained_variance_ratio_ * 100
            
            plt.figure(figsize=(10, 8))
            
            if self.metadata is not None and 'condition' in self.metadata.columns:
                # Match sample order between PCA and metadata
                samples = data_to_explore.columns
                meta_subset = self.metadata[self.metadata['sample_id'].isin(samples)]
                meta_subset = meta_subset.set_index('sample_id').loc[samples].reset_index()
                
                # Create scatter plot with condition colors
                sns.scatterplot(
                    x=pca_result[:, 0],
                    y=pca_result[:, 1],
                    hue=meta_subset['condition'],
                    s=100
                )
                plt.legend(title='Condition')
            else:
                # Basic PCA plot without condition colors
                plt.scatter(pca_result[:, 0], pca_result[:, 1], s=100)
                
                # Annotate sample names
                for i, sample in enumerate(data_to_explore.columns):
                    plt.annotate(sample, (pca_result[i, 0], pca_result[i, 1]))
            
            plt.title('PCA of Gene Expression Data')
            plt.xlabel(f'PC1 ({variance_explained[0]:.2f}%)')
            plt.ylabel(f'PC2 ({variance_explained[1]:.2f}%)')
            plt.tight_layout()
            plt.savefig(os.path.join(output_dir, f"{prefix}pca_plot.png"))
            plt.close()
            
            # Save PC loadings for potential future analysis
            pc_df = pd.DataFrame(pca.components_.T, index=data_to_explore.index, columns=['PC1', 'PC2'])
            pc_df.to_csv(os.path.join(output_dir, f"{prefix}pca_loadings.csv"))
            
        # 3. Sample correlation heatmap
        plt.figure(figsize=(12, 10))
        corr_matrix = data_to_explore.corr()
        sns.heatmap(corr_matrix, annot=False, cmap='coolwarm', vmin=-1, vmax=1)
        plt.title('Sample Correlation Heatmap')
        plt.tight_layout()
        plt.savefig(os.path.join(output_dir, f"{prefix}sample_correlation.png"))
        plt.close()
        
        logger.info(f"Data exploration plots saved to {output_dir}")
        return self
    
    def differential_expression(self, group1, group2, method='t-test', fdr=0.05, lfc=1.0):
        """
        Perform differential expression analysis between two groups.
        
        Args:
            group1 (list): List of sample IDs for group 1
            group2 (list): List of sample IDs for group 2
            method (str): Statistical method ('t-test', 'wilcoxon')
            fdr (float): False discovery rate threshold
            lfc (float): Log fold-change threshold
            
        Returns:
            self: For method chaining
        """
        if self.normalized_data is None:
            logger.warning("No normalized data available. Please preprocess the data first.")
            return self
        
        logger.info(f"Performing differential expression analysis: {group1[0]}... vs {group2[0]}...")
        
        # Extract data for each group
        group1_data = self.normalized_data[group1]
        group2_data = self.normalized_data[group2]
        
        results = []
        
        # Calculate statistics for each gene
        for gene in self.normalized_data.index:
            g1_expr = group1_data.loc[gene].values
            g2_expr = group2_data.loc[gene].values
            
            # Calculate mean expression
            mean_g1 = np.mean(g1_expr)
            mean_g2 = np.mean(g2_expr)
            
            # Calculate log fold change
            log_fold_change = mean_g2 - mean_g1  # Already log2 transformed data
            
            # Apply statistical test
            if method == 't-test':
                stat, pvalue = stats.ttest_ind(g1_expr, g2_expr, equal_var=False)
            elif method == 'wilcoxon':
                stat, pvalue = stats.ranksums(g1_expr, g2_expr)
            else:
                raise ValueError(f"Unsupported statistical method: {method}")
            
            results.append({
                'gene_id': gene,
                'mean_group1': mean_g1,
                'mean_group2': mean_g2,
                'log2_fold_change': log_fold_change,
                'statistic': stat,
                'p_value': pvalue
            })
        
        # Convert to DataFrame
        self.diff_expr_results = pd.DataFrame(results)
        
        # Calculate adjusted p-values (FDR)
        self.diff_expr_results['padj'] = smm.multipletests(
            self.diff_expr_results['p_value'], method='fdr_bh'
        )[1]
        
        # Annotate results with gene symbols
        if len(self.diff_expr_results) > 0:
            try:
                gene_ids = self.diff_expr_results['gene_id'].tolist()
                annotations = self.mg.querymany(gene_ids, scopes='entrezgene,ensembl.gene', 
                                              fields='symbol,name', species='human')
                
                # Create a dictionary for mapping gene IDs to symbols and names
                symbol_map = {}
                name_map = {}
                for ann in annotations:
                    if 'query' in ann:
                        gene_id = ann['query']
                        symbol = ann.get('symbol', 'Unknown')
                        name = ann.get('name', 'Unknown')
                        symbol_map[gene_id] = symbol
                        name_map[gene_id] = name
                
                # Add the annotations to the results
                self.diff_expr_results['gene_symbol'] = self.diff_expr_results['gene_id'].map(symbol_map)
                self.diff_expr_results['gene_name'] = self.diff_expr_results['gene_id'].map(name_map)
            except Exception as e:
                logger.warning(f"Error annotating genes: {str(e)}")
        
        # Identify significant genes
        sig_genes = self.diff_expr_results[
            (self.diff_expr_results['padj'] < fdr) & 
            (abs(self.diff_expr_results['log2_fold_change']) > lfc)
        ]
        
        logger.info(f"Found {len(sig_genes)} significantly differentially expressed genes "
                   f"(FDR < {fdr}, |log2FC| > {lfc})")
        
        return self
    
    def visualize_diff_expr(self, output_dir='.', prefix=''):
        """
        Visualize differential expression results.
        
        Args:
            output_dir (str): Directory to save plots
            prefix (str): Prefix for output file names
            
        Returns:
            self: For method chaining
        """
        if self.diff_expr_results is None:
            logger.error("No differential expression results available.")
            return self
            
        logger.info("Visualizing differential expression results")
        
        # Create output directory if it doesn't exist
        os.makedirs(output_dir, exist_ok=True)
        
        # 1. Volcano plot
        plt.figure(figsize=(10, 8))
        
        # Plot all genes in gray
        plt.scatter(
            self.diff_expr_results['log2_fold_change'],
            -np.log10(self.diff_expr_results['padj']),
            alpha=0.5,
            s=30,
            color='gray'
        )
        
        # Highlight significant genes
        sig_genes = self.diff_expr_results[
            (self.diff_expr_results['padj'] < 0.05) & 
            (abs(self.diff_expr_results['log2_fold_change']) > 1)
        ]
        
        if len(sig_genes) > 0:
            plt.scatter(
                sig_genes['log2_fold_change'],
                -np.log10(sig_genes['padj']),
                alpha=0.8,
                s=50,
                color='red'
            )
            
            # Label top significant genes
            top_genes = sig_genes.nsmallest(10, 'padj')
            for _, gene in top_genes.iterrows():
                label = gene.get('gene_symbol', gene['gene_id'])
                plt.annotate(
                    label,
                    (gene['log2_fold_change'], -np.log10(gene['padj'])),
                    fontsize=9
                )
        
        plt.axhline(-np.log10(0.05), linestyle='--', color='black', alpha=0.3)
        plt.axvline(-1, linestyle='--', color='black', alpha=0.3)
        plt.axvline(1, linestyle='--', color='black', alpha=0.3)
        
        plt.xlabel('Log2 Fold Change')
        plt.ylabel('-log10(adjusted p-value)')
        plt.title('Volcano Plot of Differential Expression')
        plt.tight_layout()
        plt.savefig(os.path.join(output_dir, f"{prefix}volcano_plot.png"))
        plt.close()
        
        # 2. Heatmap of top differentially expressed genes
        if self.normalized_data is not None and len(sig_genes) > 0:
            # Get top genes by significance
            top_n = min(50, len(sig_genes))
            top_genes = sig_genes.nsmallest(top_n, 'padj')
            top_gene_ids = top_genes['gene_id'].tolist()
            
            # Extract expression data for these genes
            if len(top_gene_ids) > 0:
                heatmap_data = self.normalized_data.loc[top_gene_ids]
                
                # Generate labels for genes
                if 'gene_symbol' in top_genes.columns:
                    labels = top_genes.set_index('gene_id')['gene_symbol'].to_dict()
                    heatmap_data = heatmap_data.rename(index=labels)
                
                plt.figure(figsize=(12, min(10, 0.25 * top_n + 2)))
                sns.clustermap(
                    heatmap_data,
                    cmap='viridis',
                    z_score=0,  # Z-score rows (genes)
                    figsize=(12, min(14, 0.25 * top_n + 4)),
                    dendrogram_ratio=0.1,
                    yticklabels=True
                )
                plt.savefig(os.path.join(output_dir, f"{prefix}heatmap_top_deg.png"))
                plt.close()
        
        # 3. MA plot
        plt.figure(figsize=(10, 8))
        
        # Calculate mean expression
        self.diff_expr_results['mean_expr'] = (
            self.diff_expr_results['mean_group1'] + 
            self.diff_expr_results['mean_group2']
        ) / 2
        
        # Plot all genes in gray
        plt.scatter(
            self.diff_expr_results['mean_expr'],
            self.diff_expr_results['log2_fold_change'],
            alpha=0.5,
            s=30,
            color='gray'
        )
        
        # Highlight significant genes
        if len(sig_genes) > 0:
            plt.scatter(
                sig_genes['mean_expr'],
                sig_genes['log2_fold_change'],
                alpha=0.8,
                s=50,
                color='red'
            )
            
            # Label top significant genes
            top_genes = sig_genes.nsmallest(10, 'padj')
            for _, gene in top_genes.iterrows():
                label = gene.get('gene_symbol', gene['gene_id'])
                plt.annotate(
                    label,
                    (gene['mean_expr'], gene['log2_fold_change']),
                    fontsize=9
                )
        
        plt.axhline(0, linestyle='--', color='black', alpha=0.3)
        plt.axhline(1, linestyle='--', color='black', alpha=0.3)
        plt.axhline(-1, linestyle='--', color='black', alpha=0.3)
        
        plt.xlabel('Mean Expression')
        plt.ylabel('Log2 Fold Change')
        plt.title('MA Plot of Differential Expression')
        plt.tight_layout()
        plt.savefig(os.path.join(output_dir, f"{prefix}ma_plot.png"))
        plt.close()
        
        # Save the results to CSV
        self.diff_expr_results.to_csv(os.path.join(output_dir, f"{prefix}differential_expression_results.csv"), index=False)
        
        logger.info(f"Differential expression visualizations saved to {output_dir}")
        return self
    
    def pathway_analysis(self, gene_set='KEGG_2021_Human', organism='human', output_dir='.', prefix=''):
        """
        Perform pathway enrichment analysis on differentially expressed genes.
        
        Args:
            gene_set (str): Gene set collection to use
            organism (str): Organism name
            output_dir (str): Directory to save results
            prefix (str): Prefix for output file names
            
        Returns:
            self: For method chaining
        """
        if self.diff_expr_results is None:
            logger.error("No differential expression results available.")
            return self
        
        logger.info(f"Performing pathway analysis using {gene_set}")
        
        # Create output directory if it doesn't exist
        os.makedirs(output_dir, exist_ok=True)
        
        # Find significant genes for enrichment
        sig_genes = self.diff_expr_results[
            (self.diff_expr_results['padj'] < 0.05) & 
            (abs(self.diff_expr_results['log2_fold_change']) > 1)
        ]
        
        if len(sig_genes) == 0:
            logger.warning("No significant genes found for pathway analysis.")
            return self
        
        # Extract gene symbols for enrichment
        if 'gene_symbol' in sig_genes.columns:
            gene_list = sig_genes['gene_symbol'].dropna().tolist()
        else:
            gene_list = sig_genes['gene_id'].tolist()
        
        if len(gene_list) == 0:
            logger.warning("No valid gene symbols available for enrichment.")
            return self
            
        try:
            # Perform enrichment analysis
            enr = gp.enrichr(
                gene_list=gene_list,
                gene_sets=gene_set,
                organism=organism,
                outdir=None,  # Don't save output files yet
                no_plot=True
            )
            
            # Get the results as dataframe
            self.pathway_results = enr.results
            
            logger.info(f"Found {len(self.pathway_results)} enriched pathways")
            
            # Filter for significant pathways
            sig_pathways = self.pathway_results[self.pathway_results['Adjusted P-value'] < 0.05]
            
            if len(sig_pathways) > 0:
                logger.info(f"Found {len(sig_pathways)} significantly enriched pathways (FDR < 0.05)")
                
                # Plot top enriched pathways
                if len(sig_pathways) > 0:
                    top_pathways = sig_pathways.head(20)  # Top 20 pathways
                    
                    plt.figure(figsize=(12, 8))
                    sns.barplot(
                        x='-log10(Adjusted P-value)',
                        y='Term',
                        data=top_pathways.assign(**{
                            '-log10(Adjusted P-value)': -np.log10(top_pathways['Adjusted P-value'])
                        }).sort_values('-log10(Adjusted P-value)', ascending=True)
                    )
                    plt.title(f'Top Enriched Pathways ({gene_set})')
                    plt.tight_layout()
                    plt.savefig(os.path.join(output_dir, f"{prefix}enriched_pathways.png"))
                    plt.close()
                    
                    # Network visualization of enriched pathways
                    try:
                        gp.plot.dotplot(
                            enr.results.head(20),
                            title=f'Enriched Pathways ({gene_set})',
                            figsize=(12, 10),
                            cmap='viridis',
                            ofname=os.path.join(output_dir, f"{prefix}pathway_dotplot.png")
                        )
                    except Exception as e:
                        logger.warning(f"Error creating pathway dotplot: {str(e)}")
                        
                    # Save the results to CSV
                    self.pathway_results.to_csv(os.path.join(output_dir, f"{prefix}pathway_results.csv"), index=False)
            
        except Exception as e:
            logger.error(f"Error in pathway analysis: {str(e)}")
        
        return self
    
    def ai_interpretation(self, output_dir='.', prefix=''):
        """
        Use OpenAI to interpret differential expression and pathway results.
        
        Args:
            output_dir (str): Directory to save interpretation
            prefix (str): Prefix for output file names
            
        Returns:
            self: For method chaining
        """
        if self.openai_client is None:
            logger.warning("OpenAI client not available. Skipping AI interpretation.")
            return self
        
        if self.diff_expr_results is None:
            logger.warning("No differential expression results available for interpretation.")
            return self
            
        logger.info("Generating AI interpretation of analysis results")
        
        # Create output directory if it doesn't exist
        os.makedirs(output_dir, exist_ok=True)
        
        try:
            # Find significant genes for interpretation
            sig_genes = self.diff_expr_results[
                (self.diff_expr_results['padj'] < 0.05) & 
                (abs(self.diff_expr_results['log2_fold_change']) > 1)
            ]
            
            # Extract the top 50 significant genes for interpretation
            if len(sig_genes) > 0:
                top_genes = sig_genes.nsmallest(50, 'padj')
                
                # Format the gene data for the AI
                if 'gene_symbol' in top_genes.columns and 'gene_name' in top_genes.columns:
                    gene_info = []
                    for _, gene in top_genes.iterrows():
                        gene_info.append({
                            "gene_symbol": gene.get('gene_symbol', gene['gene_id']),
                            "gene_name": gene.get('gene_name', "Unknown"),
                            "log2_fold_change": float(gene['log2_fold_change']),
                            "padj": float(gene['padj'])
                        })
                else:
                    gene_info = top_genes[['gene_id', 'log2_fold_change', 'padj']].to_dict('records')
            
                # Format pathway results if available
                pathway_info = []
                if self.pathway_results is not None:
                    sig_pathways = self.pathway_results[self.pathway_results['Adjusted P-value'] < 0.05]
                    if len(sig_pathways) > 0:
                        top_pathways = sig_pathways.head(20)
                        pathway_info = top_pathways[['Term', 'Adjusted P-value', 'Genes']].to_dict('records')
                
                # Create the prompt for the AI
                prompt = f"""
                I need you to interpret the results of a differential gene expression analysis. 
                
                Here are the top differentially expressed genes (adjusted p-value < 0.05, |log2 fold change| > 1):
                
                {json.dumps(gene_info, indent=2)}
                
                """
                
                if pathway_info:
                    prompt += f"""
                    Here are the top enriched pathways:
                    
                    {json.dumps(pathway_info, indent=2)}
                    
                    """
                
                prompt += """
                Please provide a comprehensive interpretation including:
                1. A summary of the key differentially expressed genes and their biological functions
                2. The main biological pathways and processes affected
                3. Potential biological implications of these changes
                4. Suggestions for further experimental validation
                
                Keep the interpretation scientifically accurate but accessible.
                """
                
                # Generate the AI interpretation
                response = self.openai_client.chat.completions.create(
                    model="gpt-4",
                    messages=[
                        {"role": "system", "content": "You are a bioinformatics expert specialized in gene expression analysis. Provide scientifically accurate interpretations of genomic data."},
                        {"role": "user", "content": prompt}
                    ],
                    temperature=0.2,
                    max_tokens=1500
                )
                
                # Extract the interpretation
                interpretation = response.choices[0].message.content
                
                # Save the interpretation to a file
                with open(os.path.join(output_dir, f"{prefix}ai_interpretation.txt"), "w") as f:
                    f.write(interpretation)
                
                logger.info(f"AI interpretation saved to {os.path.join(output_dir, f'{prefix}ai_interpretation.txt')}")
                
        except Exception as e:
            logger.error(f"Error generating AI interpretation: {str(e)}")
        
        return self
    
    def run_full_analysis(self, expression_file, metadata_file=None, group1=None, group2=None, 
                         output_dir='results', prefix=''):
        """
        Run the full analysis pipeline.
        
        Args:
            expression_file (str): Path to expression data file
            metadata_file (str, optional): Path to metadata file
            group1 (list): List of sample IDs for group 1
            group2 (list): List of sample IDs for group 2
            output_dir (str): Directory to save results
            prefix (str): Prefix for output file names
            
        Returns:
            self: For method chaining
        """
        logger.info("Starting full analysis pipeline")
        
        # Create output directory if it doesn't exist
        os.makedirs(output_dir, exist_ok=True)
        
        # 1. Load data
        self.load_data(expression_file, metadata_file)
        
        # 2. Preprocess data
        self.preprocess_data()
        
        # 3. Explore data
        self.explore_data(output_dir, prefix)
        
        # 4. Differential expression analysis
        if group1 is not None and group2 is not None:
            self.differential_expression(group1, group2)
            self.visualize_diff_expr(output_dir, prefix)
            
            # 5. Pathway analysis
            self.pathway_analysis(output_dir=output_dir, prefix=prefix)
            
            # 6. AI interpretation
            self.ai_interpretation(output_dir, prefix)
        
        logger.info(f"Full analysis complete. Results saved to {output_dir}")
        return self


def main():
    """Main function for command line execution."""
    parser = argparse.ArgumentParser(description="Differential Gene Expression and Pathway Analysis Pipeline")
    
    parser.add_argument("--expr", required=True, help="Path to expression data file (CSV, TSV, or Excel)")
    parser.add_argument("--meta", help="Path to metadata file (optional)")
    parser.add_argument("--format", default="csv", choices=["csv", "tsv", "excel"], help="File format")
    parser.add_argument("--output", default="results", help="Output directory for results")
    parser.add_argument("--prefix", default="", help="Prefix for output files")
    parser.add_argument("--group1", nargs="+", help="Sample IDs for group 1")
    parser.add_argument("--group2", nargs="+", help="Sample IDs for group 2")
    parser.add_argument("--entrez-email", help="Email for NCBI Entrez API")
    parser.add_argument("--gene-set", default="KEGG_2021_Human", help="Gene set to use for pathway analysis")
    parser.add_argument("--organism", default="human", help="Organism")
    
    args = parser.parse_args()
    
    # Configure and run analysis
    config = {
        'entrez_email': args.entrez_email
    }
    
    analyzer = GeneExpressionAnalyzer(config)
    
    try:
        analyzer.run_full_analysis(
            expression_file=args.expr,
            metadata_file=args.meta,
            group1=args.group1,
            group2=args.group2,
            output_dir=args.output,
            prefix=args.prefix
        )
    except Exception as e:
        logger.error(f"Error in analysis: {str(e)}")
        sys.exit(1)


if __name__ == "__main__":
    main()