From 03af32b9b57f87abaa5762b1b89e79f9a4b53e5b Mon Sep 17 00:00:00 2001 From: "copilot-swe-agent[bot]" <198982749+Copilot@users.noreply.github.com> Date: Wed, 25 Mar 2026 20:05:06 +0000 Subject: [PATCH 1/3] Initial plan From 56f3fea7daeee4a1824d1c9d8ab5875a457c66d4 Mon Sep 17 00:00:00 2001 From: "copilot-swe-agent[bot]" <198982749+Copilot@users.noreply.github.com> Date: Wed, 25 Mar 2026 20:14:41 +0000 Subject: [PATCH 2/3] Add DGE analysis scripts and notebook for RNA-seq and microarray datasets Co-authored-by: mejian1 <55107518+mejian1@users.noreply.github.com> Agent-Logs-Url: https://github.com/mejian1/ExopherGeneExpressionProfiling/sessions/b9803320-3ce5-489f-8205-4260d2b43455 --- RNAseq_DGE_analysis.ipynb | 630 ++++++++++++++++++++++++++++++++++++++ microarray_DGE_GSE27677.R | 175 +++++++++++ microarray_DGE_GSE9720.R | 153 +++++++++ 3 files changed, 958 insertions(+) create mode 100644 RNAseq_DGE_analysis.ipynb create mode 100644 microarray_DGE_GSE27677.R create mode 100644 microarray_DGE_GSE9720.R diff --git a/RNAseq_DGE_analysis.ipynb b/RNAseq_DGE_analysis.ipynb new file mode 100644 index 0000000..2775f0f --- /dev/null +++ b/RNAseq_DGE_analysis.ipynb @@ -0,0 +1,630 @@ +{ + "nbformat": 4, + "nbformat_minor": 0, + "metadata": { + "colab": { + "provenance": [], + "authorship_tag": "ABX9TyN/TKwX0Wa6g8aUZvub0HyK", + "include_colab_link": true + }, + "kernelspec": { + "name": "python3", + "display_name": "Python 3" + }, + "language_info": { + "name": "python" + } + }, + "cells": [ + { + "cell_type": "markdown", + "metadata": { + "id": "view-in-github", + "colab_type": "text" + }, + "source": [ + "\"Open" + ] + }, + { + "cell_type": "markdown", + "metadata": { + "id": "intro_md" + }, + "source": [ + "# RNA-seq Differential Gene Expression Analysis\n", + "\n", + "This notebook performs a complete differential gene expression (DGE) analysis pipeline\n", + "for RNA-seq data. The workflow covers:\n", + "\n", + "1. **Install dependencies** – PyDESeq2 and supporting libraries\n", + "2. **Load & inspect raw count data** – from TSV/CSV files or simulated data\n", + "3. **Quality control** – library size distribution, PCA, sample correlation\n", + "4. **Normalization** – size-factor estimation (DESeq2 median-of-ratios)\n", + "5. **Differential expression** – Wald test via PyDESeq2\n", + "6. **Results filtering** – apply adjusted p-value and fold-change thresholds\n", + "7. **Visualization** – Volcano plot, MA plot, heatmap of top DEGs, PCA\n", + "8. **Export results** – write DEG tables to TSV files\n", + "\n", + "**Dataset context**: *C. elegans* sbp-1 RNAi knockdown RNA-seq experiment\n", + "(GEO accession GSE70692). Upload the following files before running:\n", + "- `GSE70692_sbp1RNAi_raw_counts.tsv` – integer read-count matrix\n", + " (genes × samples, with a `sample_info.tsv` metadata file)\n", + "\n", + "If raw count files are not available the notebook generates a synthetic\n", + "dataset that mirrors the structure of the experiment so all code cells\n", + "can still be executed end-to-end.\n" + ] + }, + { + "cell_type": "markdown", + "metadata": { + "id": "install_md" + }, + "source": [ + "## 1. Install Dependencies\n" + ] + }, + { + "cell_type": "code", + "metadata": { + "id": "install_cell" + }, + "source": [ + "# Install PyDESeq2 and additional plotting / analysis libraries\n", + "!pip install pydeseq2 adjustText --quiet\n", + "!pip install matplotlib seaborn scikit-learn scipy statsmodels --quiet" + ], + "execution_count": null, + "outputs": [] + }, + { + "cell_type": "markdown", + "metadata": { + "id": "load_md" + }, + "source": [ + "## 2. Load and Inspect Count Data\n" + ] + }, + { + "cell_type": "code", + "metadata": { + "id": "load_cell" + }, + "source": [ + "import os\n", + "import numpy as np\n", + "import pandas as pd\n", + "\n", + "COUNTS_FILE = '/content/GSE70692_sbp1RNAi_raw_counts.tsv'\n", + "METADATA_FILE = '/content/sample_info.tsv'\n", + "\n", + "# ---------------------------------------------------------------\n", + "# If real data files are not uploaded, create a synthetic dataset\n", + "# that mirrors the sbp-1 RNAi experiment structure:\n", + "# 3 control samples vs 3 sbp-1 RNAi samples, ~20 000 genes\n", + "# ---------------------------------------------------------------\n", + "if not os.path.exists(COUNTS_FILE):\n", + " print(\"[INFO] Count file not found – generating synthetic dataset.\")\n", + " rng = np.random.default_rng(42)\n", + " n_genes, n_samples = 20_000, 6\n", + " sample_ids = [f\"ctrl_{i}\" for i in range(1, 4)] + [f\"sbp1_{i}\" for i in range(1, 4)]\n", + " gene_ids = [f\"WBGene{str(i).zfill(8)}\" for i in range(1, n_genes + 1)]\n", + "\n", + " # Base counts drawn from negative binomial\n", + " base_mu = rng.gamma(shape=1, scale=100, size=n_genes)\n", + " counts = rng.negative_binomial(n=10, p=10 / (10 + base_mu[:, None]),\n", + " size=(n_genes, n_samples))\n", + "\n", + " # Introduce ~1 000 truly differentially expressed genes\n", + " n_de = 1_000\n", + " de_idx = rng.choice(n_genes, n_de, replace=False)\n", + " fold_changes = rng.choice([-1, 1], n_de) * rng.uniform(1, 4, n_de)\n", + " counts[de_idx, 3:] = np.clip(\n", + " (counts[de_idx, :3] * (2 ** fold_changes[:, None])).astype(int), 0, None\n", + " )\n", + "\n", + " counts_df = pd.DataFrame(counts, index=gene_ids, columns=sample_ids)\n", + " counts_df.to_csv(COUNTS_FILE, sep='\\t')\n", + "\n", + " metadata = pd.DataFrame({'condition': ['control'] * 3 + ['sbp1_RNAi'] * 3},\n", + " index=sample_ids)\n", + " metadata.to_csv(METADATA_FILE, sep='\\t')\n", + " print(f\"[INFO] Synthetic counts: {counts_df.shape[0]} genes × {counts_df.shape[1]} samples\")\n", + "\n", + "# Load count matrix (genes as rows, samples as columns)\n", + "counts_df = pd.read_csv(COUNTS_FILE, sep='\\t', index_col=0)\n", + "metadata = pd.read_csv(METADATA_FILE, sep='\\t', index_col=0)\n", + "\n", + "print(\"Count matrix shape:\", counts_df.shape)\n", + "print(\"\\nSample metadata:\")\n", + "display(metadata)\n", + "print(\"\\nFirst 5 genes:\")\n", + "display(counts_df.head())" + ], + "execution_count": null, + "outputs": [] + }, + { + "cell_type": "markdown", + "metadata": { + "id": "qc_md" + }, + "source": [ + "## 3. Quality Control\n", + "\n", + "Inspect library sizes, per-gene count distributions, and sample-to-sample\n", + "correlations before normalization.\n" + ] + }, + { + "cell_type": "code", + "metadata": { + "id": "qc_cell" + }, + "source": [ + "import matplotlib.pyplot as plt\n", + "import seaborn as sns\n", + "from sklearn.preprocessing import StandardScaler\n", + "from sklearn.decomposition import PCA\n", + "\n", + "conditions = metadata['condition']\n", + "palette = {'control': '#1B9E77', 'sbp1_RNAi': '#D95F02'}\n", + "colors = [palette[c] for c in conditions]\n", + "\n", + "fig, axes = plt.subplots(1, 3, figsize=(18, 5))\n", + "\n", + "# --- Library sizes ---\n", + "lib_sizes = counts_df.sum(axis=0) / 1e6\n", + "axes[0].bar(range(len(lib_sizes)), lib_sizes, color=colors)\n", + "axes[0].set_xticks(range(len(lib_sizes)))\n", + "axes[0].set_xticklabels(lib_sizes.index, rotation=45, ha='right')\n", + "axes[0].set_xlabel('Sample')\n", + "axes[0].set_ylabel('Library size (M reads)')\n", + "axes[0].set_title('Library Sizes')\n", + "from matplotlib.patches import Patch\n", + "axes[0].legend(handles=[Patch(color=v, label=k) for k, v in palette.items()])\n", + "\n", + "# --- Log2 count distribution (pseudo-count + 1) ---\n", + "log_counts = np.log2(counts_df + 1)\n", + "axes[1].boxplot(log_counts.values, labels=log_counts.columns,\n", + " patch_artist=True,\n", + " boxprops=dict(facecolor='lightblue'))\n", + "axes[1].set_xlabel('Sample')\n", + "axes[1].set_ylabel('log2(count + 1)')\n", + "axes[1].set_title('Count Distribution per Sample')\n", + "plt.setp(axes[1].get_xticklabels(), rotation=45, ha='right')\n", + "\n", + "# --- Sample correlation heatmap ---\n", + "corr_mat = log_counts.corr()\n", + "sns.heatmap(corr_mat, ax=axes[2], annot=True, fmt='.3f',\n", + " cmap='coolwarm', vmin=0.8, vmax=1.0,\n", + " linewidths=0.5)\n", + "axes[2].set_title('Sample Correlation (log2 counts)')\n", + "\n", + "plt.tight_layout()\n", + "plt.show()\n", + "print(\"Minimum library size:\", lib_sizes.min().round(2), \"M reads\")\n", + "print(\"Maximum library size:\", lib_sizes.max().round(2), \"M reads\")" + ], + "execution_count": null, + "outputs": [] + }, + { + "cell_type": "markdown", + "metadata": { + "id": "deseq_md" + }, + "source": [ + "## 4 & 5. Normalization and Differential Expression with PyDESeq2\n", + "\n", + "PyDESeq2 implements the DESeq2 method (Love et al. 2014) in Python.\n", + "It estimates size factors using the median-of-ratios method and fits\n", + "a negative binomial GLM per gene, followed by a Wald test.\n" + ] + }, + { + "cell_type": "code", + "metadata": { + "id": "deseq_cell" + }, + "source": [ + "from pydeseq2.dds import DeseqDataSet\n", + "from pydeseq2.ds import DeseqStats\n", + "\n", + "# PyDESeq2 expects samples as rows and genes as columns\n", + "counts_T = counts_df.T\n", + "\n", + "# Remove genes with zero counts in all samples\n", + "counts_T = counts_T.loc[:, counts_T.sum(axis=0) > 0]\n", + "print(f\"Genes after filtering zero-count genes: {counts_T.shape[1]}\")\n", + "\n", + "# Build DESeq2 dataset\n", + "dds = DeseqDataSet(\n", + " counts = counts_T,\n", + " metadata = metadata,\n", + " design_factors = 'condition',\n", + " refit_cooks = True,\n", + " n_cpus = 4,\n", + ")\n", + "\n", + "# Run DESeq2 (size-factor estimation + dispersion + GLM fitting)\n", + "dds.deseq2()\n", + "\n", + "# Extract normalized counts for downstream visualization\n", + "norm_counts = dds.layers['normed_counts']\n", + "norm_counts_df = pd.DataFrame(norm_counts,\n", + " index = counts_T.index,\n", + " columns= counts_T.columns)\n", + "print(\"\\nSize factors:\")\n", + "print(dds.size_factors)\n" + ], + "execution_count": null, + "outputs": [] + }, + { + "cell_type": "code", + "metadata": { + "id": "wald_cell" + }, + "source": [ + "# Wald test: sbp1_RNAi vs control\n", + "stat_res = DeseqStats(\n", + " dds,\n", + " contrast = ('condition', 'sbp1_RNAi', 'control'),\n", + " alpha = 0.05,\n", + " cooks_filter = True,\n", + " independent_filter = True,\n", + ")\n", + "stat_res.run_wald_test()\n", + "stat_res.p_value_adjustment() # Benjamini-Hochberg FDR\n", + "\n", + "# Collect results into a DataFrame\n", + "results = stat_res.results_df.copy()\n", + "results.index.name = 'gene_id'\n", + "results = results.sort_values('padj', na_position='last')\n", + "\n", + "print(f\"Total genes tested: {len(results)}\")\n", + "print(f\"Significant DEGs (padj < 0.05): {(results['padj'] < 0.05).sum()}\")\n", + "display(results.head(10))" + ], + "execution_count": null, + "outputs": [] + }, + { + "cell_type": "markdown", + "metadata": { + "id": "filter_md" + }, + "source": [ + "## 6. Filter and Summarize DEG Results\n" + ] + }, + { + "cell_type": "code", + "metadata": { + "id": "filter_cell" + }, + "source": [ + "PADJ_THRESHOLD = 0.05\n", + "LFC_THRESHOLD = 1.0 # |log2 fold change| >= 1 (i.e., >=2-fold change)\n", + "\n", + "sig = results.dropna(subset=['padj'])\n", + "sig = sig[sig['padj'] < PADJ_THRESHOLD]\n", + "\n", + "deg_up = sig[sig['log2FoldChange'] > LFC_THRESHOLD]\n", + "deg_down = sig[sig['log2FoldChange'] < -LFC_THRESHOLD]\n", + "deg_all = sig[sig['log2FoldChange'].abs() >= LFC_THRESHOLD]\n", + "\n", + "print(f\"Up-regulated DEGs (padj<{PADJ_THRESHOLD}, logFC>{LFC_THRESHOLD}): {len(deg_up)}\")\n", + "print(f\"Down-regulated DEGs (padj<{PADJ_THRESHOLD}, logFC<-{LFC_THRESHOLD}): {len(deg_down)}\")\n", + "print(f\"Total DEGs (padj<{PADJ_THRESHOLD}, |logFC|>={LFC_THRESHOLD}): {len(deg_all)}\")\n", + "\n", + "print(\"\\nTop up-regulated genes:\")\n", + "display(deg_up.sort_values('log2FoldChange', ascending=False).head(10))\n", + "\n", + "print(\"\\nTop down-regulated genes:\")\n", + "display(deg_down.sort_values('log2FoldChange').head(10))" + ], + "execution_count": null, + "outputs": [] + }, + { + "cell_type": "markdown", + "metadata": { + "id": "viz_md" + }, + "source": [ + "## 7. Visualization\n", + "\n", + "### 7a. Volcano Plot\n" + ] + }, + { + "cell_type": "code", + "metadata": { + "id": "volcano_cell" + }, + "source": [ + "import matplotlib.pyplot as plt\n", + "import numpy as np\n", + "\n", + "plot_df = results.dropna(subset=['padj', 'log2FoldChange']).copy()\n", + "# Clip p-values to avoid log10(0); 1e-100 is well within float64 range\n", + "plot_df['-log10_padj'] = -np.log10(plot_df['padj'].clip(lower=1e-100))\n", + "\n", + "# Assign color categories\n", + "def assign_color(row):\n", + " if row['padj'] < PADJ_THRESHOLD and row['log2FoldChange'] >= LFC_THRESHOLD:\n", + " return 'Up-regulated'\n", + " elif row['padj'] < PADJ_THRESHOLD and row['log2FoldChange'] <= -LFC_THRESHOLD:\n", + " return 'Down-regulated'\n", + " else:\n", + " return 'Not significant'\n", + "\n", + "plot_df['category'] = plot_df.apply(assign_color, axis=1)\n", + "color_map = {'Up-regulated': '#D95F02', 'Down-regulated': '#1B9E77', 'Not significant': '#AAAAAA'}\n", + "\n", + "fig, ax = plt.subplots(figsize=(10, 7))\n", + "for cat, grp in plot_df.groupby('category'):\n", + " ax.scatter(grp['log2FoldChange'], grp['-log10_padj'],\n", + " c=color_map[cat], alpha=0.5, s=8, label=cat, rasterized=True)\n", + "\n", + "# Threshold lines\n", + "ax.axhline(-np.log10(PADJ_THRESHOLD), color='black', linestyle='--', linewidth=0.8, alpha=0.7)\n", + "ax.axvline( LFC_THRESHOLD, color='black', linestyle='--', linewidth=0.8, alpha=0.7)\n", + "ax.axvline(-LFC_THRESHOLD, color='black', linestyle='--', linewidth=0.8, alpha=0.7)\n", + "\n", + "# Label top 10 genes by -log10(padj)\n", + "top_genes = plot_df.nlargest(10, '-log10_padj')\n", + "for _, row in top_genes.iterrows():\n", + " ax.annotate(row.name, (row['log2FoldChange'], row['-log10_padj']),\n", + " fontsize=7, ha='center', va='bottom',\n", + " xytext=(0, 4), textcoords='offset points')\n", + "\n", + "ax.set_xlabel('log2 Fold Change (sbp-1 RNAi / control)', fontsize=12)\n", + "ax.set_ylabel('-log10(adjusted p-value)', fontsize=12)\n", + "ax.set_title('Volcano Plot – sbp-1 RNAi vs Control (RNA-seq)', fontsize=13)\n", + "ax.legend(markerscale=2)\n", + "ax.grid(True, linestyle='--', alpha=0.3)\n", + "plt.tight_layout()\n", + "plt.savefig('volcano_plot_sbp1_RNAi.png', dpi=150, bbox_inches='tight')\n", + "plt.show()\n", + "print(\"Volcano plot saved to volcano_plot_sbp1_RNAi.png\")" + ], + "execution_count": null, + "outputs": [] + }, + { + "cell_type": "markdown", + "metadata": { + "id": "ma_md" + }, + "source": [ + "### 7b. MA Plot (log fold change vs mean expression)\n" + ] + }, + { + "cell_type": "code", + "metadata": { + "id": "ma_cell" + }, + "source": [ + "fig, ax = plt.subplots(figsize=(9, 6))\n", + "\n", + "for cat, grp in plot_df.groupby('category'):\n", + " ax.scatter(np.log2(grp['baseMean'] + 1), grp['log2FoldChange'],\n", + " c=color_map[cat], alpha=0.4, s=6, label=cat, rasterized=True)\n", + "\n", + "ax.axhline(0, color='darkblue', linewidth=1)\n", + "ax.set_xlabel('log2 Mean Expression (baseMean)', fontsize=12)\n", + "ax.set_ylabel('log2 Fold Change', fontsize=12)\n", + "ax.set_title('MA Plot – sbp-1 RNAi vs Control', fontsize=13)\n", + "ax.legend(markerscale=3)\n", + "ax.grid(True, linestyle='--', alpha=0.3)\n", + "plt.tight_layout()\n", + "plt.savefig('MA_plot_sbp1_RNAi.png', dpi=150, bbox_inches='tight')\n", + "plt.show()\n", + "print(\"MA plot saved to MA_plot_sbp1_RNAi.png\")" + ], + "execution_count": null, + "outputs": [] + }, + { + "cell_type": "markdown", + "metadata": { + "id": "heatmap_md" + }, + "source": [ + "### 7c. Heatmap of Top 50 DEGs\n" + ] + }, + { + "cell_type": "code", + "metadata": { + "id": "heatmap_cell" + }, + "source": [ + "import seaborn as sns\n", + "\n", + "top50 = deg_all.nsmallest(50, 'padj').index\n", + "\n", + "# Z-score normalize across samples\n", + "heat_data = norm_counts_df.loc[:, top50].T\n", + "mean_expr = heat_data.mean(axis=1).values[:, None]\n", + "std_expr = heat_data.std(axis=1).values[:, None]\n", + "heat_z = (heat_data - mean_expr) / std_expr\n", + "\n", + "# Column colour bar for conditions\n", + "col_colors = pd.Series(\n", + " [palette[conditions[s]] for s in heat_z.columns],\n", + " index=heat_z.columns\n", + ")\n", + "\n", + "g = sns.clustermap(\n", + " heat_z,\n", + " col_colors = col_colors,\n", + " cmap = 'vlag',\n", + " center = 0,\n", + " method = 'average',\n", + " figsize = (10, 14),\n", + " dendrogram_ratio = (0.1, 0.2),\n", + " cbar_pos = (0, 0.2, 0.03, 0.4),\n", + ")\n", + "g.fig.suptitle('Heatmap – Top 50 DEGs (Z-score)', y=1.01, fontsize=14)\n", + "plt.savefig('heatmap_top50_DEGs.png', dpi=150, bbox_inches='tight')\n", + "plt.show()\n", + "print(\"Heatmap saved to heatmap_top50_DEGs.png\")" + ], + "execution_count": null, + "outputs": [] + }, + { + "cell_type": "markdown", + "metadata": { + "id": "pca_md" + }, + "source": [ + "### 7d. PCA of Normalized Counts\n" + ] + }, + { + "cell_type": "code", + "metadata": { + "id": "pca_cell" + }, + "source": [ + "from sklearn.decomposition import PCA\n", + "from sklearn.preprocessing import StandardScaler\n", + "\n", + "# Use top 500 most variable genes\n", + "top_var_genes = norm_counts_df.var(axis=0).nlargest(500).index\n", + "pca_input = norm_counts_df[top_var_genes]\n", + "\n", + "scaler = StandardScaler()\n", + "pca = PCA(n_components=2)\n", + "pca_coords = pca.fit_transform(scaler.fit_transform(pca_input))\n", + "\n", + "fig, ax = plt.subplots(figsize=(7, 6))\n", + "for cond, grp_idx in metadata.groupby('condition').groups.items():\n", + " mask = [s in grp_idx for s in norm_counts_df.index]\n", + " ax.scatter(pca_coords[mask, 0], pca_coords[mask, 1],\n", + " label=cond, color=palette[cond], s=100, edgecolors='k')\n", + "\n", + "for i, sample in enumerate(norm_counts_df.index):\n", + " ax.annotate(sample, (pca_coords[i, 0], pca_coords[i, 1]),\n", + " fontsize=8, ha='center', va='bottom',\n", + " xytext=(0, 6), textcoords='offset points')\n", + "\n", + "ax.set_xlabel(f'PC1 ({pca.explained_variance_ratio_[0]:.1%} variance)', fontsize=12)\n", + "ax.set_ylabel(f'PC2 ({pca.explained_variance_ratio_[1]:.1%} variance)', fontsize=12)\n", + "ax.set_title('PCA of Normalized Counts (top 500 variable genes)', fontsize=13)\n", + "ax.legend()\n", + "ax.grid(True, linestyle='--', alpha=0.3)\n", + "plt.tight_layout()\n", + "plt.savefig('PCA_normalized_counts.png', dpi=150, bbox_inches='tight')\n", + "plt.show()\n", + "print(\"PCA plot saved to PCA_normalized_counts.png\")" + ], + "execution_count": null, + "outputs": [] + }, + { + "cell_type": "markdown", + "metadata": { + "id": "bar_md" + }, + "source": [ + "### 7e. Bar Chart – Up vs Down Regulated Gene Counts\n" + ] + }, + { + "cell_type": "code", + "metadata": { + "id": "bar_cell" + }, + "source": [ + "fig, ax = plt.subplots(figsize=(6, 5))\n", + "categories_bar = ['Up-regulated', 'Down-regulated']\n", + "counts_bar = [len(deg_up), len(deg_down)]\n", + "bar_colors_bar = ['#D95F02', '#1B9E77']\n", + "\n", + "bars = ax.bar(categories_bar, counts_bar, color=bar_colors_bar, width=0.5)\n", + "for bar in bars:\n", + " ax.text(bar.get_x() + bar.get_width() / 2,\n", + " bar.get_height() + max(counts_bar) * 0.02,\n", + " str(int(bar.get_height())),\n", + " ha='center', va='bottom', fontweight='bold', fontsize=12)\n", + "\n", + "ax.set_ylabel('Number of Genes', fontsize=12)\n", + "ax.set_title(f'DEGs (padj < {PADJ_THRESHOLD}, |logFC| ≥ {LFC_THRESHOLD})\\nsbp-1 RNAi vs Control',\n", + " fontsize=12)\n", + "ax.set_ylim(0, max(counts_bar) * 1.15)\n", + "ax.grid(axis='y', linestyle='--', alpha=0.4)\n", + "plt.tight_layout()\n", + "plt.savefig('bar_chart_DEG_counts.png', dpi=150, bbox_inches='tight')\n", + "plt.show()\n", + "print(\"Bar chart saved to bar_chart_DEG_counts.png\")" + ], + "execution_count": null, + "outputs": [] + }, + { + "cell_type": "markdown", + "metadata": { + "id": "export_md" + }, + "source": [ + "## 8. Export Results\n" + ] + }, + { + "cell_type": "code", + "metadata": { + "id": "export_cell" + }, + "source": [ + "# All results (sorted by adjusted p-value)\n", + "results.to_csv('GSE70692_sbp1RNAi_DESeq2_all_genes.tsv', sep='\\t')\n", + "\n", + "# Significant DEG subsets\n", + "deg_up.to_csv('GSE70692_sbp1RNAi_DESeq2_DEGs_up.tsv', sep='\\t')\n", + "deg_down.to_csv('GSE70692_sbp1RNAi_DESeq2_DEGs_down.tsv', sep='\\t')\n", + "deg_all.to_csv('GSE70692_sbp1RNAi_DESeq2_DEGs_P0.05_FC2.tsv', sep='\\t')\n", + "\n", + "print(\"Results exported:\")\n", + "print(f\" All genes: GSE70692_sbp1RNAi_DESeq2_all_genes.tsv ({len(results)} rows)\")\n", + "print(f\" Up-regulated DEGs: GSE70692_sbp1RNAi_DESeq2_DEGs_up.tsv ({len(deg_up)} rows)\")\n", + "print(f\" Down-regulated DEGs:GSE70692_sbp1RNAi_DESeq2_DEGs_down.tsv ({len(deg_down)} rows)\")\n", + "print(f\" DEGs (|logFC|>=1): GSE70692_sbp1RNAi_DESeq2_DEGs_P0.05_FC2.tsv ({len(deg_all)} rows)\")" + ], + "execution_count": null, + "outputs": [] + }, + { + "cell_type": "markdown", + "metadata": { + "id": "summary_md" + }, + "source": [ + "## Summary\n", + "\n", + "| Step | Method | Output |\n", + "|------|--------|--------|\n", + "| Normalization | DESeq2 median-of-ratios | Size factors, normalized count matrix |\n", + "| Differential expression | Negative binomial GLM + Wald test | log2 fold change, p-value, adjusted p-value |\n", + "| Multiple testing correction | Benjamini-Hochberg FDR | `padj` column |\n", + "| Thresholds | padj < 0.05 AND \\|log2FC\\| ≥ 1 | Up/down DEG lists |\n", + "| Visualizations | Volcano, MA, Heatmap, PCA, Bar chart | PNG files |\n", + "\n", + "### Key References\n", + "- Love MI, Huber W, Anders S (2014). *Moderated estimation of fold change and\n", + " dispersion for RNA-seq data with DESeq2*. Genome Biology, 15:550.\n", + "- Muzellec B, Telenczuk M, Cabeli V, Andreux M (2023). *PyDESeq2: a python\n", + " package for bulk RNA-seq differential expression analysis*. Bioinformatics.\n", + "- GSE70692: sbp-1 RNAi C. elegans RNA-seq dataset. NCBI GEO.\n" + ] + } + ] +} \ No newline at end of file diff --git a/microarray_DGE_GSE27677.R b/microarray_DGE_GSE27677.R new file mode 100644 index 0000000..e7b4a56 --- /dev/null +++ b/microarray_DGE_GSE27677.R @@ -0,0 +1,175 @@ +# ============================================================ +# Differential Gene Expression Analysis - Microarray Dataset +# GSE27677: "A Fasting-Responsive Signaling Pathway that +# Extends Life Span in C. elegans" +# DOI: 10.1016/j.celrep.2012.12.018 +# Platform: GPL200 (Affymetrix C. elegans Genome Array) +# ============================================================ +# Requirements: BiocManager::install(c("GEOquery","limma","Biobase")) + +library(Biobase) +library(GEOquery) +library(limma) + +# ---------------------------------------------------------------- +# 1. Load series and platform data from GEO +# ---------------------------------------------------------------- +gset <- getGEO("GSE27677", GSEMatrix = TRUE, AnnotGPL = TRUE) +if (length(gset) > 1) idx <- grep("GPL200", attr(gset, "names")) else idx <- 1 +gset <- gset[[idx]] + +# Sanitize feature variable labels +fvarLabels(gset) <- make.names(fvarLabels(gset)) + +# ---------------------------------------------------------------- +# 2. Assign sample groups +# gsms string: one character per sample, X = exclude +# Groups: +# 0 = control (fed/untreated) +# 1 = fasting 3 hr +# 2 = fasting 6 hr +# 3 = fed (refeeding control) +# 4 = fasting (long-term) +# ---------------------------------------------------------------- +gsms <- "0000001122XXXXXXXXXXXX333444XXXXXXXXXXXXXXXXXX" +sml <- strsplit(gsms, split = "")[[1]] + +# Exclude samples marked as "X" +sel <- which(sml != "X") +sml <- sml[sel] +gset <- gset[, sel] + +# ---------------------------------------------------------------- +# 3. Log2-transform expression values (if not already on log scale) +# ---------------------------------------------------------------- +ex <- exprs(gset) +qx <- as.numeric(quantile(ex, c(0, 0.25, 0.5, 0.75, 0.99, 1.0), na.rm = TRUE)) +LogC <- (qx[5] > 100) || + (qx[6] - qx[1] > 50 && qx[2] > 0) || + (qx[2] > 0 && qx[2] < 1 && qx[4] > 1 && qx[4] < 2) +if (LogC) { + ex[ex <= 0] <- NaN + exprs(gset) <- log2(ex) +} + +# ---------------------------------------------------------------- +# 4. Build design matrix and fit linear model with limma +# ---------------------------------------------------------------- +sml <- paste0("G", sml) +fl <- as.factor(sml) +gset$description <- fl +design <- model.matrix(~ description + 0, gset) +colnames(design) <- levels(fl) + +fit <- lmFit(gset, design) + +# Contrasts: +# G4-G0: long-term fasting vs control (primary comparison) +# G1-G0: 3 hr fasting vs control +# G2-G1: 6 hr fasting vs 3 hr fasting (sequential time-point) +# G3-G2: refed vs 6 hr fasting (sequential time-point) +# G4-G3: long-term fasting vs refed (sequential time-point) +cont.matrix <- makeContrasts( + G4 - G0, # long-term fasting vs control + G1 - G0, # 3 hr fasting vs control + G2 - G1, # 6 hr vs 3 hr fasting + G3 - G2, # refed vs 6 hr fasting + G4 - G3, # long-term fasting vs refed + levels = design +) +fit2 <- contrasts.fit(fit, cont.matrix) +fit2 <- eBayes(fit2, proportion = 0.01) + +# ---------------------------------------------------------------- +# 5. Extract top differentially expressed genes (primary contrast: +# long-term fasting vs control, G4-G0) +# ---------------------------------------------------------------- +tT <- topTable(fit2, adjust = "fdr", sort.by = "B", number = 250, + coef = "G4 - G0") + +# Select informative columns +available_cols <- intersect( + c("ID", "adj.P.Val", "P.Value", "F", "t", "B", "logFC", + "Gene.symbol", "Gene.title", "Gene.Symbol"), + colnames(tT) +) +tT <- subset(tT, select = available_cols) +write.table(tT, file = "GSE27677_fasting_DGE_top250.tsv", row.names = FALSE, sep = "\t") +cat("Top 250 DEGs written to GSE27677_fasting_DGE_top250.tsv\n") + +# Full table (all genes, primary contrast) +tT_all <- topTable(fit2, adjust = "fdr", sort.by = "B", number = Inf, + coef = "G4 - G0") +write.table(tT_all, file = "GSE27677_fasting_DGE_all_genes.tsv", row.names = FALSE, sep = "\t") +cat("All-gene table written to GSE27677_fasting_DGE_all_genes.tsv\n") + +# ---------------------------------------------------------------- +# 6. Significant DEG subsets (adj.P.Val < 0.05, |logFC| >= 1) +# ---------------------------------------------------------------- +sig_up <- tT_all[!is.na(tT_all$adj.P.Val) & tT_all$adj.P.Val < 0.05 & tT_all$logFC > 0, ] +sig_dn <- tT_all[!is.na(tT_all$adj.P.Val) & tT_all$adj.P.Val < 0.05 & tT_all$logFC < 0, ] +sig_fc2 <- tT_all[!is.na(tT_all$adj.P.Val) & tT_all$adj.P.Val < 0.05 & abs(tT_all$logFC) >= 1, ] + +write.table(sig_up, file = "GSE27677_fasting_DEGs_up.tsv", row.names = FALSE, sep = "\t") +write.table(sig_dn, file = "GSE27677_fasting_DEGs_down.tsv", row.names = FALSE, sep = "\t") +write.table(sig_fc2, file = "GSE27677_fasting_DEGs_P0.05_FC2.tsv", row.names = FALSE, sep = "\t") +cat(sprintf("Up-regulated DEGs (fasting vs control): %d\n", nrow(sig_up))) +cat(sprintf("Down-regulated DEGs (fasting vs control): %d\n", nrow(sig_dn))) +cat(sprintf("DEGs |logFC|>=1 (fasting vs control): %d\n", nrow(sig_fc2))) + +# ---------------------------------------------------------------- +# 7. Diagnostic Plots +# ---------------------------------------------------------------- +dT <- decideTests(fit2, adjust.method = "fdr", p.value = 0.05, lfc = 0) + +# -- P-value histogram (primary contrast) -- +hist(tT_all$adj.P.Val, + col = "grey", border = "white", + xlab = "Adjusted P-value (FDR)", + ylab = "Number of genes", + main = "GSE27677 – Adjusted P-value Distribution (Fasting vs Control)") + +# -- Venn diagram across contrasts -- +vennDiagram(dT, circle.col = palette()) + +# -- QQ plot of moderated t-statistics -- +t_good <- which(!is.na(fit2$F)) +qqt(fit2$t[t_good], fit2$df.total[t_good], + main = "Moderated t-statistic QQ Plot (GSE27677)") + +# -- Volcano plot (primary contrast: G4-G0) -- +ct <- 1 # first contrast: G4-G0 +volcanoplot(fit2, coef = ct, + main = colnames(fit2)[ct], + pch = 20, + highlight = length(which(dT[, ct] != 0)), + names = rep("+", nrow(fit2))) + +# -- MA (mean-difference) plot -- +plotMD(fit2, column = ct, status = dT[, ct], + legend = FALSE, pch = 20, cex = 1, + main = "MA Plot – Fasting vs Control (GSE27677)") +abline(h = 0) + +# ---------------------------------------------------------------- +# 8. Box-and-whisker and density plots across samples +# ---------------------------------------------------------------- +labels_map <- c(G0 = "control", G1 = "fasting 3 hr", + G2 = "fasting 6 hr", G3 = "refeeding control", + G4 = "fasting (long-term)") + +ex_ord <- exprs(gset)[, order(fl)] +fl_ord <- fl[order(fl)] +palette(c("#dfeaf4", "#f4dfdf", "#f2cb98", "#b7d9b7", "#d7a3d0")) +par(mar = c(7, 4, 2, 1)) +boxplot(ex_ord, + boxwex = 0.6, notch = TRUE, + main = paste0("GSE27677 / ", annotation(gset), " selected samples"), + outline = FALSE, las = 2, + col = fl_ord) +legend("topleft", labels_map[levels(fl)], fill = palette(), bty = "n") + +par(mar = c(4, 4, 2, 1)) +plotDensities(exprs(gset), group = fl, + main = paste0("GSE27677 / ", annotation(gset), " Value Distribution"), + legend = "topright") diff --git a/microarray_DGE_GSE9720.R b/microarray_DGE_GSE9720.R new file mode 100644 index 0000000..c88d898 --- /dev/null +++ b/microarray_DGE_GSE9720.R @@ -0,0 +1,153 @@ +# ============================================================ +# Differential Gene Expression Analysis - Microarray Dataset +# GSE9720: "The Mediator Subunit MDT-15 Confers Metabolic +# Adaptation to Ingested Material" +# DOI: 10.1371/journal.pgen.1000021 +# Platform: GPL5859 (Affymetrix C. elegans Genome Array) +# ============================================================ +# Requirements: BiocManager::install(c("GEOquery","limma","Biobase")) + +library(Biobase) +library(GEOquery) +library(limma) + +# ---------------------------------------------------------------- +# 1. Load series and platform data from GEO +# ---------------------------------------------------------------- +gset <- getGEO("GSE9720", GSEMatrix = TRUE, AnnotGPL = FALSE) +if (length(gset) > 1) idx <- grep("GPL5859", attr(gset, "names")) else idx <- 1 +gset <- gset[[idx]] + +# Sanitize feature variable labels so they are valid R names +fvarLabels(gset) <- make.names(fvarLabels(gset)) + +# ---------------------------------------------------------------- +# 2. Assign sample groups +# GSE9720 has 5 samples: GSM246680 (N2 control), GSM246681 (fat-6 +# RNAi), GSM246682 (fat-7 RNAi), GSM246683 (fat-6/fat-7 RNAi), +# GSM246684 (mdt-15 RNAi). +# gsms encodes one character per sample: +# 0 = control (G0), 1 = mdt-15 RNAi (G1), X = exclude +# Update this string to reflect the sample order in your download. +# ---------------------------------------------------------------- +gsms <- "0XXX1" # 1 control (GSM246680), exclude fat-6/7 RNAi (GSM246681-3), 1 mdt-15 RNAi (GSM246684) +sml <- strsplit(gsms, split = "")[[1]] + +# Exclude samples marked as "X" +sel <- which(sml != "X") +sml <- sml[sel] +gset <- gset[, sel] + +# ---------------------------------------------------------------- +# 3. Log2-transform expression values (if not already on log scale) +# ---------------------------------------------------------------- +ex <- exprs(gset) +qx <- as.numeric(quantile(ex, c(0, 0.25, 0.5, 0.75, 0.99, 1.0), na.rm = TRUE)) +LogC <- (qx[5] > 100) || + (qx[6] - qx[1] > 50 && qx[2] > 0) || + (qx[2] > 0 && qx[2] < 1 && qx[4] > 1 && qx[4] < 2) +if (LogC) { + ex[ex <= 0] <- NaN + exprs(gset) <- log2(ex) +} + +# ---------------------------------------------------------------- +# 4. Build design matrix and fit linear model with limma +# ---------------------------------------------------------------- +sml <- paste0("G", sml) # e.g. "G0", "G1" +fl <- as.factor(sml) +gset$description <- fl +design <- model.matrix(~ description + 0, gset) +colnames(design) <- levels(fl) + +fit <- lmFit(gset, design) + +# Contrast: mdt-15 RNAi (G1) vs control (G0) +cont.matrix <- makeContrasts(G1 - G0, levels = design) +fit2 <- contrasts.fit(fit, cont.matrix) +fit2 <- eBayes(fit2, proportion = 0.01) + +# ---------------------------------------------------------------- +# 5. Extract top differentially expressed genes +# ---------------------------------------------------------------- +tT <- topTable(fit2, adjust = "fdr", sort.by = "B", number = 250) + +# Select informative columns (adjust column names to what is available) +available_cols <- intersect(c("ID", "adj.P.Val", "P.Value", "t", "B", "logFC", "ORF", "Gene.Symbol"), + colnames(tT)) +tT <- subset(tT, select = available_cols) +write.table(tT, file = "GSE9720_mdt15_DGE_top250.tsv", row.names = FALSE, sep = "\t") +cat("Top 250 DEGs written to GSE9720_mdt15_DGE_top250.tsv\n") + +# Full table (all genes) +tT_all <- topTable(fit2, adjust = "fdr", sort.by = "B", number = Inf) +write.table(tT_all, file = "GSE9720_mdt15_DGE_all_genes.tsv", row.names = FALSE, sep = "\t") +cat("All-gene table written to GSE9720_mdt15_DGE_all_genes.tsv\n") + +# ---------------------------------------------------------------- +# 6. Significant DEG subsets (adj.P.Val < 0.05, |logFC| >= 1) +# ---------------------------------------------------------------- +dT <- decideTests(fit2, adjust.method = "fdr", p.value = 0.05, lfc = 0) +sig_up <- tT_all[tT_all$adj.P.Val < 0.05 & tT_all$logFC > 0, ] +sig_dn <- tT_all[tT_all$adj.P.Val < 0.05 & tT_all$logFC < 0, ] +sig_fc2 <- tT_all[tT_all$adj.P.Val < 0.05 & abs(tT_all$logFC) >= 1, ] + +write.table(sig_up, file = "GSE9720_mdt15_DEGs_up.tsv", row.names = FALSE, sep = "\t") +write.table(sig_dn, file = "GSE9720_mdt15_DEGs_down.tsv", row.names = FALSE, sep = "\t") +write.table(sig_fc2, file = "GSE9720_mdt15_DEGs_P0.05_FC2.tsv", row.names = FALSE, sep = "\t") +cat(sprintf("Up-regulated DEGs: %d\n", nrow(sig_up))) +cat(sprintf("Down-regulated DEGs: %d\n", nrow(sig_dn))) +cat(sprintf("DEGs |logFC|>=1: %d\n", nrow(sig_fc2))) + +# ---------------------------------------------------------------- +# 7. Diagnostic Plots +# ---------------------------------------------------------------- + +# -- P-value histogram -- +hist(tT_all$adj.P.Val, + col = "grey", border = "white", + xlab = "Adjusted P-value (FDR)", + ylab = "Number of genes", + main = "GSE9720 – Adjusted P-value Distribution") + +# -- Venn diagram of up/down regulated genes -- +vennDiagram(dT, circle.col = palette()) + +# -- QQ plot of moderated t-statistics -- +t_good <- which(!is.na(fit2$F)) +qqt(fit2$t[t_good], fit2$df.total[t_good], + main = "Moderated t-statistic QQ Plot (GSE9720)") + +# -- Volcano plot -- +ct <- 1 # first (and only) contrast +volcanoplot(fit2, coef = ct, + main = colnames(fit2)[ct], + pch = 20, + highlight = length(which(dT[, ct] != 0)), + names = rep("+", nrow(fit2))) + +# -- MA (mean-difference) plot -- +plotMD(fit2, column = ct, status = dT[, ct], + legend = FALSE, pch = 20, cex = 1, + main = "MA Plot – MDT-15 RNAi vs Control (GSE9720)") +abline(h = 0) + +# ---------------------------------------------------------------- +# 8. Box-and-whisker plot of expression across samples +# ---------------------------------------------------------------- +ex_ord <- exprs(gset)[, order(fl)] +fl_ord <- fl[order(fl)] +palette(c("#1B9E77", "#7570B3")) +par(mar = c(7, 4, 2, 1)) +boxplot(ex_ord, + boxwex = 0.6, notch = TRUE, + main = paste0("GSE9720 / ", annotation(gset)), + outline = FALSE, las = 2, + col = fl_ord) +legend("topleft", levels(fl), fill = palette(), bty = "n") + +# -- Expression value density plot -- +par(mar = c(4, 4, 2, 1)) +plotDensities(exprs(gset), group = fl, + main = paste0("GSE9720 / ", annotation(gset), " Value Distribution"), + legend = "topright") From 6e9b32c63973152db1f83fd136aa0a3ca83631d1 Mon Sep 17 00:00:00 2001 From: Driscoll Lab Date: Fri, 10 Apr 2026 19:27:55 -0400 Subject: [PATCH 3/3] Enhance README with updates and structure improvements Expanded README with additional contributors, dataset details, and high-impact improvements for clarity and usability. --- README.md | 54 +++++++++++++++++++++++++++++++++++++++++++++++++++++- 1 file changed, 53 insertions(+), 1 deletion(-) diff --git a/README.md b/README.md index 1d550b3..1b4083e 100644 --- a/README.md +++ b/README.md @@ -1,9 +1,61 @@ # ExopherGeneExpressionProfiling Driscoll Lab at Rutgers University Molecular Biology, Biochemistry, Genetics, functional genomics and bioinformatics -For additional information and consultation regarding the code, please email me, Nelson Mejia, or Mark Saba if you'd like. +For additional information and consultation regarding the code, please email me, Nelson Mejia, Mark Saba, Durga, Julissa, and James Saba who is our computer scientist. This repository contains the R scripts and Python code to analyze two gene expression microarray profiles. The first Dataset is titled "A Fasting-Responsive Signaling Pathway that Extends Life Span in C. elegans" DOI: 10.1016/j.celrep.2012.12.018 NCBI GEO Accession Number GSE27677 +Harvald E, Sprenger R, Dall K..Multi-omics Analyses of Starvation Responses Reveal a Central Role for Lipoprotein Metabolism in Acute Starvation Survival in C. elegans +Cell Systems, 2017; 5, 38-52.e4 GEO accession number; GSE98919 The second Dataset is titled "The Mediator Subunit MDT-15 Confers Metabolic Adaptation to Ingested Material" DOI: 10.1371/journal.pgen.1000021 MNCBI GEO Accession Number GSE9720 Each dataset was analyzed slightly differently given the microarrays used in each study, so for each script, the process of normalization, transformation, and final processing varies, so please pay close attention to which script is being used and execute the line of code with caution. This repository is still under construction, so please check back for updates! + + + +of note, +Please follow up with these crucial updates! + +High-impact README improvements (bioinformatics / differential expression) +Even without seeing it yet, these are the most common upgrades that materially help: + +Clear purpose + scope (top 10 lines) + +1–2 sentence summary of what the folder/pipeline does (inputs → method → outputs). +Explicitly state organism, tissue/cell type, and comparison being tested (if known). +Reproducibility section + +Exact toolchain: R/Python version, key packages + versions, OS. +A renv.lock / environment.yml / requirements.txt mention and how to restore. +Inputs/Outputs contract + +Bullet list of required input files with examples (counts matrix format, metadata columns, reference genome annotation). +A table of outputs (DE table columns, normalized counts, plots) + where they are written. +One “Quickstart” command path + +The minimal commands to run end-to-end, ideally copy/paste: +install deps +run +expected output location +Method details (brief but precise) + +Which DE method (DESeq2 / edgeR / limma-voom / Seurat, etc.) +Filtering thresholds, normalization, multiple testing correction (FDR), and model formula. +Troubleshooting + FAQ + +Common errors (missing sample IDs, non-integer counts for DESeq2, batch effect confusion). +How to validate the design matrix / contrasts. +Project structure + +Short tree of directories and what each contains. +Citations + +Cite the main method paper(s) and any relevant dataset references. +If you have a DOI / preprint / manuscript, link it. +Make it skimmable + +Use consistent headings, add a TOC, keep paragraphs short, add code blocks. +Add a “Results sanity checks” section + +MA plot / volcano expectations +sample clustering/PCA +how to interpret log2FC and adjusted p-values