{ "nbformat": 4, "nbformat_minor": 0, "metadata": { "colab": { "provenance": [] }, "kernelspec": { "name": "python3", "display_name": "Python 3" }, "language_info": { "name": "python" } }, "cells": [ { "cell_type": "markdown", "source": [ "The colorectal cancer (CRC) cell line RNA-seq data was obtained from the expanded Cancer Cell Line Encyclopedia (CCLE) project, submitted by Ghandi et al. (Broad Institute/Novartis). The raw paired-end transcriptomic sequencing files for 56 distinct colorectal cancer cell lines were accessed via the NCBI BioProject database under the accession number PRJNA523380 (Sequence Read Archive: SRP186687)." ], "metadata": { "id": "9NZ-yvqybbY7" } }, { "cell_type": "markdown", "source": [ "### CODE CELL -1\n", "## 1. Data Acquisition and Cohort Standardization\n", "\n", "To ensure the highest statistical confidence in our `rMATS` alternative splicing analysis, transcriptomic data was carefully selected to eliminate technical batch effects between the tumor and normal cohorts.\n", "\n", "* **The Malignant Cohort:** Colorectal cancer cell lines from the Cancer Cell Line Encyclopedia (CCLE, BioProject Accession: **PRJNA523380**).\n", "* **The Normal Cohort:** 18 samples of fresh-frozen human normal colorectal mucosa (GEO Accession: **GSE50760**).\n", "\n", "### Rationale for Standardization:\n", "1. **Library Preparation Parity (Poly-A Selection):** Both the CCLE cohort and GSE50760 utilized oligo-dT magnetic beads for poly-A selection. This guarantees we are not introducing a \"Total RNA vs. Poly-A\" batch effect, which would artificially skew the Percent Spliced In (PSI) ratios of our target lncRNA scaffolds.\n", "2. **Exclusion of FFPE Tissue:** Formalin fixation chemically cross-links and degrades RNA, creating artificial dropouts that algorithms like `rMATS` misinterpret as skipped exons. By strictly utilizing fresh-frozen mucosal tissue (GSE50760), we preserved the integrity of the exon-intron boundaries.\n", "3. **Epithelial/Mucosal Matching:** Colorectal carcinoma originates in the epithelial lining. Selecting GSE50760 ensures the normal baseline consists of colonic mucosa rather than full-thickness bowel resections, mathematically eliminating smooth muscle and stromal RNA noise from the baseline comparison." ], "metadata": { "id": "PPyHW1TZkbQQ" } }, { "cell_type": "markdown", "source": [ "2.1. Data Acquisition and Cohort Standardization\n", "Raw paired-end RNA-sequencing data for the colorectal cancer cell lines were accessed from the Sequence Read Archive (SRA) via the NCBI Entrez Direct API. To ensure statistical comparability across the cohort and mitigate sequencing depth bias during differential alternative splicing analysis, a strict library size threshold was applied. SRA accession files exceeding 10 GB in compressed format were excluded from the primary analysis. Files of this magnitude typically represent extreme ultra-deep sequencing runs (>150 million reads) or Total RNA libraries heavily burdened with ribosomal RNA (rRNA) retention.\n", "\n", "Evidence/Justification: It is well-documented in transcriptomics that mixing extreme read-depth outliers with standard-depth libraries (20-50 million reads) artificially inflates the false discovery rate (FDR) for novel splice junctions. Algorithms that calculate Percent Spliced In (PSI) are highly sensitive to library size confounding; un-normalized excessive read depth in a subset of samples can result in the spurious detection of alternative splicing events. Therefore, enforcing this maximum size threshold standardizes the statistical power across the cohort and reduces technical batch effects while preserving robust coverage of polyadenylated transcripts, including long non-coding RNAs such as GAS5 and PVT1. Data retrieval was performed utilizing the SRA Toolkit (v3.0.0) prefetch utility to ensure complete data transfer, followed by fasterq-dump for local FASTQ extraction.\n", "\n", "2.2. Quality Control and Pre-processing\n", "Raw sequencing reads underwent rigorous quality assessment using FastQC. Adapter sequences, poly-G tails, and low-quality bases were trimmed utilizing fastp (v0.23.0) with default heuristic sliding-window quality filtering. Reads falling below the quality threshold or minimum length constraints were discarded to prevent aberrant splice junction mapping.\n", "\n", "2.3. Spliced-Aware Genome Alignment\n", "Cleaned reads were mapped to the human reference genome (GRCh38 primary assembly) with annotations provided by GENCODE (Release 40). Alignment was executed using the Spliced Transcripts Alignment to a Reference (STAR) software (v2.7.10a). The STAR index was generated with a splice junction overhang (--sjdbOverhang) of 50 to optimize mapping across annotated splice sites.\n", "\n", "To maximize the stringency of novel junction discovery—a critical requirement for non-coding RNA splicing analysis—STAR was executed with highly customized parameters. Spurious splice junctions were filtered out using the --outFilterType BySJout command. The tolerance for mismatches was strictly capped at 5% of the read length (--outFilterMismatchNoverLmax 0.05). To capture accurate exon-intron boundaries, the minimum overhang for unannotated junctions was set to 4 bases (--alignSJoverhangMin 4), and multimapping was restricted to a maximum of 10 loci (--outFilterMultimapNmax 10).\n", "\n", "2.4. Post-processing and Cloud Storage\n", "To ensure structural integrity of the output, STAR was configured to output unsorted BAM files, which were subsequently coordinate-sorted and indexed using Samtools (samtools sort and samtools index) allocated with 24 CPU threads to prevent memory fragmentation. The finalized, coordinate-sorted BAM files and their respective .bai indices were securely transferred to a Google Cloud Storage bucket for downstream differential splicing evaluation.\n", "\n", "Below is the code of the run_pipeline.sh file which we used to download the tumor cell line files and perform star alignment on them and store the .bam and .bai files in cloud bucket fo further analysis\n", "```\n", "%%writefile run_pipeline.sh\n", "#!/bin/bash\n", "set -euo pipefail\n", "\n", "# Define your Cloud Bucket paths\n", "GCS_BUCKET=\"gs://crcvsnormal\"\n", "GCS_BAMDIR=\"${GCS_BUCKET}/BAMs\"\n", "FAILED_LOG=\"/content/failed_downloads.txt\"\n", "\n", "echo \"==========================================================\"\n", "echo \"[STEP 1] Installing Dependencies & Setting up...\"\n", "echo \"==========================================================\"\n", "apt-get update -qq && apt-get install -y -qq sra-toolkit fastqc rna-star samtools ncbi-entrez-direct\n", "\n", "if [ ! -x /usr/local/bin/fastp ]; then\n", " wget -q http://opengene.org/fastp/fastp && chmod a+x ./fastp && mv ./fastp /usr/local/bin/\n", "fi\n", "\n", "mkdir -p ~/.ncbi\n", "echo '/LIBS/GUID = \"ebe12217-19ad-45ce-8c3b-7f154fc1d4d3\"' > ~/.ncbi/user-settings.mkfg\n", "echo '/repository/user/main/public/root = \"/content/ncbi_cache\"' >> ~/.ncbi/user-settings.mkfg\n", "\n", "BASEDIR=\"/content/colorectal_ccle\"\n", "BAMDIR=\"${BASEDIR}/BAMs\"\n", "QC_DIR=\"${BASEDIR}/QC\"\n", "FASTQ_DIR=\"${BASEDIR}/FASTQ\"\n", "REF_DIR=\"${BASEDIR}/reference_genome\"\n", "STAR_IDX_DIR=\"${BASEDIR}/STAR_GRCh38_Index\"\n", "DISK_TMP_BASE=\"${BASEDIR}/star_tmp\"\n", "\n", "mkdir -p \"$BAMDIR\" \"$QC_DIR\" \"$FASTQ_DIR\" \"$DISK_TMP_BASE\" \"$REF_DIR\" \"$STAR_IDX_DIR\"\n", "\n", "echo \"==========================================================\"\n", "echo \"[STEP 2] Genome Index Verification...\"\n", "echo \"==========================================================\"\n", "if [ ! -f \"${STAR_IDX_DIR}/SAindex\" ]; then\n", " echo \"Genome index not found! Building index...\"\n", " cd \"$REF_DIR\"\n", " wget -qO GRCh38.primary_assembly.genome.fa.gz https://ftp.ebi.ac.uk/pub/databases/gencode/Gencode_human/release_40/GRCh38.primary_assembly.genome.fa.gz\n", " wget -qO gencode.v40.annotation.gtf.gz https://ftp.ebi.ac.uk/pub/databases/gencode/Gencode_human/release_40/gencode.v40.annotation.gtf.gz\n", " gunzip -f GRCh38.primary_assembly.genome.fa.gz gencode.v40.annotation.gtf.gz\n", " STAR --runThreadN 30 --runMode genomeGenerate --genomeDir \"$STAR_IDX_DIR\" --genomeFastaFiles \"${REF_DIR}/GRCh38.primary_assembly.genome.fa\" --sjdbGTFfile \"${REF_DIR}/gencode.v40.annotation.gtf\" --sjdbOverhang 50\n", "fi\n", "\n", "STAR --genomeLoad Remove --genomeDir \"$STAR_IDX_DIR\" || true\n", "\n", "echo \"==========================================================\"\n", "echo \"[STEP 3] Processing Samples...\"\n", "echo \"==========================================================\"\n", "# NEW: Replaced curly braces with 'tail' and 'cut' to prevent bash EOF errors\n", "tail -n +2 /content/56rnaseqcclecolorectal.txt | cut -d',' -f1 | tr -d '\\r' > /content/rnaseq_srr_only.txt\n", "\n", "while read -u 3 srr; do\n", " [ -z \"$srr\" ] && continue\n", " echo \"----------------------------------------------------------\"\n", " echo \"Processing Sample: $srr\"\n", " echo \"----------------------------------------------------------\"\n", "\n", " if grep -Fxq \"$srr\" \"$FAILED_LOG\" 2>/dev/null; then\n", " echo \"[BLACKLIST] $srr previously failed to download. Skipping permanently...\"\n", " continue\n", " fi\n", "\n", " if gsutil ls \"${GCS_BAMDIR}/${srr}_Aligned.sortedByCoord.out.bam.bai\" >/dev/null 2>&1; then\n", " echo \"[CLOUD RESUME] BAM and Index for $srr already exist in your bucket. Skipping...\"\n", " continue\n", " fi\n", "\n", " # NEW: Replaced curly braces with 'sed' and 'cut'\n", " SIZE_MB=$(esearch -db sra -query \"$srr\" | efetch -format runinfo 2>/dev/null | sed -n '2p' | cut -d',' -f8)\n", " if [[ -n \"$SIZE_MB\" && \"$SIZE_MB\" -gt 10000 ]]; then\n", " echo \"[BLACKLIST] Skipping $srr. File size is ${SIZE_MB} MB.\"\n", " continue\n", " fi\n", "\n", " C1=\"${FASTQ_DIR}/${srr}_1.cleaned.fastq\"\n", " C2=\"${FASTQ_DIR}/${srr}_2.cleaned.fastq\"\n", "\n", " if [ -f \"$C1\" ] && [ -f \"$C2\" ]; then\n", " echo \"[SMART RESUME] Cleaned FASTQ files found! Skipping download...\"\n", " else\n", " echo \"[CLEANUP] Removing partial files for fresh download...\"\n", " rm -f \"${FASTQ_DIR}/${srr}\"*.fastq || true\n", " rm -rf \"${DISK_TMP_BASE}/${srr}_STARtmp\" || true\n", " rm -rf /content/ncbi_cache/* ~/ncbi/public/sra/* || true\n", "\n", " MAX_RETRIES=5\n", " DOWNLOAD_SUCCESS=0\n", " for i in $(seq 1 $MAX_RETRIES); do\n", " echo \"[PREFETCH] Attempt $i to fetch $srr...\"\n", " if prefetch \"$srr\" --max-size 10G; then\n", " DOWNLOAD_SUCCESS=1\n", " break\n", " else\n", " echo \"[WARNING] Prefetch connection dropped. Retrying in 10 seconds...\"\n", " sleep 10\n", " fi\n", " done\n", "\n", " if [ \"$DOWNLOAD_SUCCESS\" -eq 0 ]; then\n", " echo \"[ERROR] Failed to prefetch $srr after 5 attempts. Permanently blacklisting.\"\n", " echo \"$srr\" >> \"$FAILED_LOG\"\n", " continue\n", " fi\n", "\n", " echo \"[EXTRACT] Prefetch complete! Extracting FASTQs locally...\"\n", " if ! fasterq-dump --split-files \"$srr\" --threads 30 --mem 40G --outdir \"$FASTQ_DIR\"; then\n", " echo \"[ERROR] Failed to extract FASTQs for $srr. Permanently blacklisting.\"\n", " echo \"$srr\" >> \"$FAILED_LOG\"\n", " continue\n", " fi\n", "\n", " R1=\"${FASTQ_DIR}/${srr}_1.fastq\"\n", " R2=\"${FASTQ_DIR}/${srr}_2.fastq\"\n", "\n", " if [ ! -f \"$R1\" ] || [ ! -f \"$R2\" ]; then\n", " continue\n", " fi\n", "\n", " echo \"[CLEANUP] FASTQs extracted. Wiping SRA cache to save disk space...\"\n", " rm -rf /content/ncbi_cache/* ~/ncbi/public/sra/* ./$srr || true\n", "\n", " echo \"[FASTQC] Running FastQC...\"\n", " fastqc -t 8 -q \"$R1\" \"$R2\" -o \"$QC_DIR\"\n", "\n", " echo \"[FASTP] Running fastp...\"\n", " fastp -i \"$R1\" -I \"$R2\" -o \"$C1\" -O \"$C2\" -h \"${QC_DIR}/${srr}_fastp.html\" -j \"${QC_DIR}/${srr}_fastp.json\" --thread 16 2>/dev/null\n", " rm -f \"$R1\" \"$R2\" || true\n", " fi\n", "\n", " TMPDIR=\"${DISK_TMP_BASE}/${srr}_STARtmp\"\n", " rm -rf \"$TMPDIR\"\n", "\n", " echo \"[STAR] Starting alignment for $srr (Unsorted Output)...\"\n", " STAR --genomeLoad LoadAndKeep --runThreadN 30 --genomeDir \"$STAR_IDX_DIR\" --readFilesIn \"$C1\" \"$C2\" --outSAMtype BAM Unsorted --outSAMstrandField intronMotif --outFilterMatchNminOverLread 0.0 --outFilterScoreMinOverLread 0.0 --outFilterMatchNmin 15 --seedSearchStartLmax 50 --outFilterMultimapNmax 10 --outFilterMismatchNoverLmax 0.05 --winAnchorMultimapNmax 50 --alignSJoverhangMin 4 --alignSJDBoverhangMin 1 --outFilterMultimapScoreRange 1 --outFilterType BySJout --outFileNamePrefix \"${BAMDIR}/${srr}_\" --outTmpDir \"$TMPDIR\" > /dev/null\n", "\n", " echo \"[SAMTOOLS] Sorting BAM file...\"\n", " samtools sort -@ 24 -m 2G -o \"${BAMDIR}/${srr}_Aligned.sortedByCoord.out.bam\" \"${BAMDIR}/${srr}_Aligned.out.bam\"\n", "\n", " echo \"[SAMTOOLS] Indexing BAM...\"\n", " samtools index -@ 24 \"${BAMDIR}/${srr}_Aligned.sortedByCoord.out.bam\"\n", "\n", " echo \"[GCS UPLOAD] Transferring finished BAM and Index to Cloud Storage bucket...\"\n", " gsutil cp \"${BAMDIR}/${srr}_Aligned.sortedByCoord.out.bam\" \"${BAMDIR}/${srr}_Aligned.sortedByCoord.out.bam.bai\" \"${GCS_BAMDIR}/\"\n", "\n", " echo \"[CLEANUP] Wiping all local FASTQs, BAMs, and Caches for $srr...\"\n", " rm -f \"${BAMDIR}/${srr}_\"* || true\n", " rm -f \"${FASTQ_DIR}/${srr}\"* || true\n", " rm -rf \"$TMPDIR\" || true\n", " rm -rf /content/ncbi_cache/* ~/ncbi/public/sra/* ./$srr || true\n", "\n", " echo \"[STATUS] $srr complete and safely stored in gs://crcvsnormal!\"\n", "done 3< /content/rnaseq_srr_only.txt\n", "\n", "STAR --genomeLoad Remove --genomeDir \"$STAR_IDX_DIR\" || true\n", "echo \"Pipeline Finished!\"\n", "```\n", "\n", "Below is the code of the run_normal_pipeline.sh file which we used to download the normal colorectal mucosa genome files and perform star alignment on them and store the .bam and .bai files in cloud bucket for further analysis\n", "\n", "```\n", "%%writefile run_normal_pipeline.sh\n", "#!/bin/bash\n", "set -euo pipefail\n", "\n", "GCS_BUCKET=\"gs://crcvsnormal\"\n", "GCS_BAMDIR=\"${GCS_BUCKET}/BAMSNORMAL\"\n", "FAILED_LOG=\"/content/normal_failed_downloads.txt\"\n", "\n", "echo \"==========================================================\"\n", "echo \"[STEP 1] Installing Dependencies & Setting up...\"\n", "echo \"==========================================================\"\n", "apt-get update -qq && apt-get install -y -qq sra-toolkit fastqc rna-star samtools ncbi-entrez-direct\n", "\n", "if [ ! -x /usr/local/bin/fastp ]; then\n", " wget -q http://opengene.org/fastp/fastp && chmod a+x ./fastp && mv ./fastp /usr/local/bin/\n", "fi\n", "\n", "mkdir -p ~/.ncbi\n", "echo '/LIBS/GUID = \"ebe12217-19ad-45ce-8c3b-7f154fc1d4d3\"' > ~/.ncbi/user-settings.mkfg\n", "echo '/repository/user/main/public/root = \"/content/ncbi_cache\"' >> ~/.ncbi/user-settings.mkfg\n", "\n", "BASEDIR=\"/content/colorectal_normal\"\n", "BAMDIR=\"${BASEDIR}/BAMs\"\n", "QC_DIR=\"${BASEDIR}/QC\"\n", "FASTQ_DIR=\"${BASEDIR}/FASTQ\"\n", "DISK_TMP_BASE=\"${BASEDIR}/star_tmp\"\n", "\n", "REF_DIR=\"/content/colorectal_ccle/reference_genome\"\n", "STAR_IDX_DIR=\"/content/colorectal_ccle/STAR_GRCh38_Index\"\n", "\n", "mkdir -p \"$BAMDIR\" \"$QC_DIR\" \"$FASTQ_DIR\" \"$DISK_TMP_BASE\"\n", "\n", "echo \"==========================================================\"\n", "echo \"[STEP 2] Genome Index Verification...\"\n", "echo \"==========================================================\"\n", "if [ ! -f \"${STAR_IDX_DIR}/SAindex\" ]; then\n", " echo \"Genome index not found! Building index...\"\n", "else\n", " echo \"Existing genome index found in ${STAR_IDX_DIR}. Skipping build!\"\n", "fi\n", "\n", "STAR --genomeLoad Remove --genomeDir \"$STAR_IDX_DIR\" || true\n", "\n", "echo \"==========================================================\"\n", "echo \"[STEP 3] Verifying Normal Tissue SRR List...\"\n", "echo \"==========================================================\"\n", "if [ -s /content/normal_rnaseq_srr_only.txt ]; then\n", " echo \"Pre-filtered SRR list found! Skipping NCBI search.\"\n", " echo \"Processing $(wc -l < /content/normal_rnaseq_srr_only.txt) genuine normal tissue samples.\"\n", "else\n", " echo \"[FATAL] /content/normal_rnaseq_srr_only.txt is missing. Run the Python snippet first!\"\n", " exit 1\n", "fi\n", "\n", "echo \"==========================================================\"\n", "echo \"[STEP 4] Processing Samples...\"\n", "echo \"==========================================================\"\n", "while read -u 3 srr; do\n", " [ -z \"$srr\" ] && continue\n", " echo \"----------------------------------------------------------\"\n", " echo \"Processing Sample: $srr\"\n", " echo \"----------------------------------------------------------\"\n", "\n", " if grep -Fxq \"$srr\" \"$FAILED_LOG\" 2>/dev/null; then\n", " echo \"[BLACKLIST] $srr previously failed to download. Skipping...\"\n", " continue\n", " fi\n", "\n", " if gsutil ls \"${GCS_BAMDIR}/${srr}_Aligned.sortedByCoord.out.bam.bai\" >/dev/null 2>&1; then\n", " echo \"[CLOUD RESUME] BAM and Index for $srr already exist. Skipping...\"\n", " continue\n", " fi\n", "\n", " SIZE_MB=$(esearch -db sra -query \"$srr\" 2>/dev/null | efetch -format runinfo 2>/dev/null | sed -n '2p' | cut -d',' -f8 || echo \"\")\n", " if [[ -n \"$SIZE_MB\" && \"$SIZE_MB\" -gt 10000 ]]; then\n", " echo \"[BLACKLIST] Skipping $srr. File size is ${SIZE_MB} MB.\"\n", " continue\n", " fi\n", "\n", " C1=\"${FASTQ_DIR}/${srr}_1.cleaned.fastq\"\n", " C2=\"${FASTQ_DIR}/${srr}_2.cleaned.fastq\"\n", "\n", " if [ -f \"$C1\" ] && [ -f \"$C2\" ]; then\n", " echo \"[SMART RESUME] Cleaned FASTQ files found! Skipping download...\"\n", " else\n", " echo \"[CLEANUP] Removing partial files for fresh download...\"\n", " rm -f \"${FASTQ_DIR}/${srr}\"*.fastq || true\n", " rm -rf \"${DISK_TMP_BASE}/${srr}_STARtmp\" || true\n", " rm -rf /content/ncbi_cache/* ~/ncbi/public/sra/* || true\n", "\n", " MAX_RETRIES=5\n", " DOWNLOAD_SUCCESS=0\n", " for i in $(seq 1 $MAX_RETRIES); do\n", " echo \"[PREFETCH] Attempt $i to fetch $srr...\"\n", " if prefetch \"$srr\" --max-size 10G; then\n", " DOWNLOAD_SUCCESS=1\n", " break\n", " else\n", " sleep 10\n", " fi\n", " done\n", "\n", " if [ \"$DOWNLOAD_SUCCESS\" -eq 0 ]; then\n", " echo \"$srr\" >> \"$FAILED_LOG\"\n", " continue\n", " fi\n", "\n", " echo \"[EXTRACT] Prefetch complete! Extracting FASTQs...\"\n", " if ! fasterq-dump --split-files \"$srr\" --threads 30 --mem 40G --outdir \"$FASTQ_DIR\"; then\n", " echo \"$srr\" >> \"$FAILED_LOG\"\n", " continue\n", " fi\n", "\n", " R1=\"${FASTQ_DIR}/${srr}_1.fastq\"\n", " R2=\"${FASTQ_DIR}/${srr}_2.fastq\"\n", "\n", " if [ ! -f \"$R1\" ] || [ ! -f \"$R2\" ]; then\n", " continue\n", " fi\n", "\n", " rm -rf /content/ncbi_cache/* ~/ncbi/public/sra/* ./$srr || true\n", "\n", " fastqc -t 8 -q \"$R1\" \"$R2\" -o \"$QC_DIR\"\n", " fastp -i \"$R1\" -I \"$R2\" -o \"$C1\" -O \"$C2\" -h \"${QC_DIR}/${srr}_fastp.html\" -j \"${QC_DIR}/${srr}_fastp.json\" --thread 16 2>/dev/null\n", " rm -f \"$R1\" \"$R2\" || true\n", " fi\n", "\n", " TMPDIR=\"${DISK_TMP_BASE}/${srr}_STARtmp\"\n", " rm -rf \"$TMPDIR\"\n", "\n", " STAR --genomeLoad LoadAndKeep --runThreadN 30 --genomeDir \"$STAR_IDX_DIR\" --readFilesIn \"$C1\" \"$C2\" --outSAMtype BAM Unsorted --outSAMstrandField intronMotif --outFilterMatchNminOverLread 0.0 --outFilterScoreMinOverLread 0.0 --outFilterMatchNmin 15 --seedSearchStartLmax 50 --outFilterMultimapNmax 10 --outFilterMismatchNoverLmax 0.05 --winAnchorMultimapNmax 50 --alignSJoverhangMin 4 --alignSJDBoverhangMin 1 --outFilterMultimapScoreRange 1 --outFilterType BySJout --outFileNamePrefix \"${BAMDIR}/${srr}_\" --outTmpDir \"$TMPDIR\" > /dev/null\n", "\n", " samtools sort -@ 24 -m 2G -o \"${BAMDIR}/${srr}_Aligned.sortedByCoord.out.bam\" \"${BAMDIR}/${srr}_Aligned.out.bam\"\n", " samtools index -@ 24 \"${BAMDIR}/${srr}_Aligned.sortedByCoord.out.bam\"\n", "\n", " gsutil cp \"${BAMDIR}/${srr}_Aligned.sortedByCoord.out.bam\" \"${BAMDIR}/${srr}_Aligned.sortedByCoord.out.bam.bai\" \"${GCS_BAMDIR}/\"\n", "\n", " rm -f \"${BAMDIR}/${srr}_\"* || true\n", " rm -f \"${FASTQ_DIR}/${srr}\"* || true\n", " rm -rf \"$TMPDIR\" || true\n", " rm -rf /content/ncbi_cache/* ~/ncbi/public/sra/* ./$srr || true\n", "\n", " echo \"[STATUS] $srr safely stored in gs://crcvsnormal/BAMSNORMAL!\"\n", "done 3< /content/normal_rnaseq_srr_only.txt\n", "\n", "STAR --genomeLoad Remove --genomeDir \"$STAR_IDX_DIR\" || true\n", "echo \"Normal Pipeline Finished!\"\n", "```\n", "\n" ], "metadata": { "id": "jxaMPMI9APzT" } }, { "cell_type": "markdown", "source": [ "# CODE CELL -2\n", "# 14. Differential Expression and Isoform Quantification (Cuffdiff)\n", "\n", "To quantify differential gene and isoform expression between the 36 colorectal cancer (CRC) cell lines and 18 normal colorectal mucosa samples, we utilized `cuffdiff` from the Cufflinks suite.\n", "\n", "### Methodological Rationale: Per-Chromosome Batching\n", "Given the massive computational footprint of aligning 54 deep-sequenced BAM files against the entire human genome simultaneously, executing `cuffdiff` in a single pass routinely exceeds standard high-performance computing (HPC) memory limits. Conversely, computationally aggregating the BAM files into a single \"Tumor\" and \"Normal\" track was actively avoided, as doing so destroys sample-level biological variance (degrees of freedom), forcing the algorithm to rely on blind dispersion estimations that yield highly conservative, statistically weak p-values.\n", "\n", "To preserve individual biological replicates for downstream Ordinary Least Squares (OLS) residualization while maintaining strict memory efficiency, we implemented a **per-chromosome batching strategy**.\n", "\n", "### Pipeline Steps:\n", "1. **Annotation Segmentation:** The GENCODE v40 primary assembly annotation (`.gtf`) was programmatically split into individual files corresponding to each primary chromosome. Alternate contigs and mitochondrial DNA were excluded to optimize processing time.\n", "2. **Sequential Quantification:** `cuffdiff` was executed sequentially on each chromosome utilizing 12 threads. This constrained the maximum RAM footprint to ~45GB (during Chromosome 1 processing), well within the limits of the computing environment, without sacrificing algorithmic rigor.\n", "3. **Data Stitching:** Upon completion, the localized outputs were computationally stitched back into unified, genome-wide tracking files. Specifically, `gene_exp.diff`, `isoform_exp.diff`, and `genes.read_group_tracking` were aggregated to feed directly into our dual-model residualization scripts.\n", "\n", "### Computational Environment\n", "* **Engine:** Cufflinks (Cuffdiff)\n", "* **Library Type:** fr-unstranded\n", "* **Threads:** 12 per chromosome\n", "\n", "```\n", "# This is formatted as code\n", "import os\n", "import subprocess\n", "import time\n", "import sys\n", "import glob\n", "\n", "def run_cmd(cmd):\n", " subprocess.run(cmd, shell=True, check=False)\n", "\n", "print(\"==========================================================\")\n", "print(\"RESUMING CUFFDIFF PIPELINE\")\n", "print(\"==========================================================\")\n", "\n", "# ==========================================\n", "# 1. INSTALL CUFFDIFF\n", "# ==========================================\n", "if subprocess.run(\"command -v cuffdiff\", shell=True, capture_output=True).returncode != 0:\n", " print(\"1. INSTALLING CUFFLINKS...\")\n", " run_cmd(\"sudo apt-get update -qq > /dev/null && sudo apt-get install -y cufflinks > /dev/null\")\n", " print(\" ✅ Cuffdiff installed.\")\n", "else:\n", " print(\" ✅ Cuffdiff is already installed. Skipping installation.\")\n", "\n", "# ==========================================\n", "# 2. LOCATE/DOWNLOAD GTF AND SPLIT\n", "# ==========================================\n", "print(\"2. LOCATING AND SPLITTING GENCODE ANNOTATION...\")\n", "GTF_ZIP = \"/content/gencode.v40.annotation.gtf.zip\"\n", "GTF = \"/content/gencode.v40.annotation.gtf\"\n", "GTF_SPLIT = \"/content/gtf_split\"\n", "\n", "if not os.path.exists(GTF):\n", " if os.path.exists(GTF_ZIP):\n", " print(f\" --> Unzipping {GTF_ZIP}...\")\n", " run_cmd(f\"unzip -q -o {GTF_ZIP} -d /content/\")\n", " extracted = glob.glob(\"/content/**/gencode.v40.annotation.gtf\", recursive=True)\n", " if extracted and extracted[0] != GTF:\n", " os.rename(extracted[0], GTF)\n", " else:\n", " print(\" --> Downloading fresh GENCODE v40 GTF directly from EBI FTP...\")\n", " run_cmd(\"wget -qO /content/gencode.v40.annotation.gtf.gz ftp://ftp.ebi.ac.uk/pub/databases/gencode/Gencode_human/release_40/gencode.v40.annotation.gtf.gz\")\n", " run_cmd(\"gunzip -f /content/gencode.v40.annotation.gtf.gz\")\n", "\n", "os.makedirs(GTF_SPLIT, exist_ok=True)\n", "\n", "if len(os.listdir(GTF_SPLIT)) < 24:\n", " print(\" --> Splitting GTF into individual chromosomes...\")\n", " awk_cmd = \"\"\"awk 'BEGIN{OFS=\"\\\\t\"} /^#/ {next} {chr=$1; gsub(/\\\\r$/,\"\",chr); out=sprintf(\"%s/%s.gtf\",\"'\"\"\" + GTF_SPLIT + \"\"\"'\", chr); print > out}' \"\"\" + GTF\n", " run_cmd(awk_cmd)\n", " print(\" ✅ Annotation split complete.\")\n", "else:\n", " print(\" ✅ Annotation already split. Skipping.\")\n", "\n", "# ==========================================\n", "# 3. RUN CUFFDIFF SEQUENTIALLY\n", "# ==========================================\n", "print(\"3. PREPARING CUFFDIFF EXECUTION...\")\n", "OUT_BASE = \"/content/cuffdiff_chromosomes\"\n", "os.makedirs(OUT_BASE, exist_ok=True)\n", "\n", "with open(\"/content/b1.txt\", \"r\") as f: normal_bams = f.read().strip()\n", "with open(\"/content/b2.txt\", \"r\") as f: tumor_bams = f.read().strip()\n", "\n", "# Gather valid chromosomes\n", "valid_chrs = []\n", "for f in sorted(glob.glob(f\"{GTF_SPLIT}/chr*.gtf\")):\n", " chr_name = os.path.basename(f).replace(\".gtf\", \"\")\n", " if \"_\" in chr_name or \"Un\" in chr_name or \"M\" in chr_name:\n", " continue\n", " valid_chrs.append(chr_name)\n", "\n", "total_chrs = len(valid_chrs)\n", "print(f\" ✅ Found {total_chrs} primary chromosomes to process.\")\n", "print(\"-\" * 58)\n", "\n", "for i, chr_name in enumerate(valid_chrs, 1):\n", " chr_gtf = f\"{GTF_SPLIT}/{chr_name}.gtf\"\n", " chr_out = f\"{OUT_BASE}/{chr_name}\"\n", " \n", " if os.path.exists(f\"{chr_out}/gene_exp.diff\"):\n", " print(f\"[{time.strftime('%H:%M:%S')}] [{i}/{total_chrs}] ⏩ Skipping {chr_name} (Already processed)\")\n", " continue\n", " \n", " print(f\"[{time.strftime('%H:%M:%S')}] [{i}/{total_chrs}] ⏳ Processing {chr_name}...\")\n", " os.makedirs(chr_out, exist_ok=True)\n", " \n", " # Run cuffdiff, piping output to log\n", " log_file = f\"{chr_out}/cuffdiff.log\"\n", " cmd = f\"cuffdiff --no-update-check -o {chr_out} -p 12 --library-type fr-unstranded -L Normal,Tumor {chr_gtf} {normal_bams} {tumor_bams} > {log_file} 2>&1\"\n", " \n", " # Start the process\n", " process = subprocess.Popen(cmd, shell=True)\n", " \n", " # Heartbeat loop for Colab\n", " while process.poll() is None:\n", " print(f\" ... still calculating {chr_name} ...\")\n", " sys.stdout.flush()\n", " time.sleep(30)\n", " \n", " if process.returncode == 0:\n", " print(f\"[{time.strftime('%H:%M:%S')}] [{i}/{total_chrs}] ✅ Finished {chr_name}\")\n", " else:\n", " print(f\" ❌ FAILED! Check log at {log_file}\")\n", " sys.exit(1)\n", "\n", "print(\"-\" * 58)\n", "\n", "# ==========================================\n", "# 4. STITCH RESULTS TOGETHER\n", "# ==========================================\n", "print(\"4. STITCHING PER-CHROMOSOME RESULTS...\")\n", "FINAL_OUT = \"/content/cuffdiff_final\"\n", "os.makedirs(FINAL_OUT, exist_ok=True)\n", "\n", "bases = [\"gene_exp.diff\", \"isoform_exp.diff\", \"genes.read_group_tracking\"]\n", "\n", "for base in bases:\n", " print(f\" --> Stitching {base}...\")\n", " out_file = f\"{FINAL_OUT}/{base}\"\n", " \n", " with open(out_file, 'w') as outfile:\n", " wrote_header = False\n", " for chr_name in valid_chrs:\n", " chr_base_file = f\"{OUT_BASE}/{chr_name}/{base}\"\n", " if os.path.exists(chr_base_file) and os.path.getsize(chr_base_file) > 0:\n", " with open(chr_base_file, 'r') as infile:\n", " lines = infile.readlines()\n", " if not lines: continue\n", " \n", " if not wrote_header:\n", " outfile.writelines(lines) # Write everything including header\n", " wrote_header = True\n", " else:\n", " outfile.writelines(lines[1:]) # Skip header\n", "\n", "print(\"==========================================================\")\n", "print(\"🎉 CUFFDIFF COMPLETED SUCCESSFULLY\")\n", "print(f\"Final stitched files are located in: {FINAL_OUT}\")\n", "print(\"==========================================================\")\n", "```\n", "\n" ], "metadata": { "id": "O4deOWRohKVk" } }, { "cell_type": "code", "source": [], "metadata": { "id": "2VjfLWyihwxg" }, "execution_count": null, "outputs": [] }, { "cell_type": "markdown", "source": [ "# CODE CELL-3\n", "### Bias Correction and Residualization of Differential Expression\n", "\n", "**Overview**\n", "To accurately identify true transcriptional upregulation and mitigate confounding factors such as isoform switching and the \"tumor dilution effect,\" a multi-step residualization workflow was applied to the standard Cuffdiff output[cite: 4]. This pipeline utilizes Ordinary Least Squares (OLS) regression to adjust the raw log2 fold changes for splicing burden and sample-level read variances[cite: 4].\n", "\n", "**Methodological Steps**\n", "\n", "1. **Raw Differential Expression Assessment:**\n", " * Initial differential expression values were extracted directly from the `gene_exp.diff` file[cite: 4].\n", " * Significant genes (FDR ≤ 0.05) were isolated, and a raw log2 fold change (`log2FC_raw`) was calculated to establish a baseline ranking of upregulated targets[cite: 4].\n", "\n", "2. **Splicing Burden Quantification:**\n", " * Isoform-level data was parsed from `isoform_exp.diff` to assess the degree of alternative splicing occurring within each gene locus[cite: 4].\n", " * An abundance filter was applied (mean FPKM ≥ 1.0) to remove low-expression noise[cite: 4].\n", " * A Splicing Burden Index (SBI) was computed for each gene by subtracting the number of significantly downregulated isoforms from the number of significantly upregulated isoforms[cite: 4].\n", "\n", "3. **Gene-Level Residualization (Correcting for Splicing Bias):**\n", " * To ensure that observed upregulation was not merely an artifact of excessive alternative splicing (where highly spliced genes appear artificially upregulated), an OLS regression model was constructed[cite: 4].\n", " * The model fit the raw log2 fold change against the total number of significant isoforms and the calculated SBI (`log2FC_raw ~ n_sig_iso + SBI`)[cite: 4].\n", " * The residuals from this model (`log2FC_resid_splicing`) represent the fold change variance independent of the gene's splicing complexity[cite: 4].\n", "\n", "4. **Sample-Level Residualization (Correcting for Replicate Variance):**\n", " * To account for inter-sample variance and technical batch effects, read group tracking data (`genes.read_group_tracking`) was integrated[cite: 4].\n", " * A secondary OLS model evaluated the log2-transformed FPKM values against the biological condition (tumor vs. normal) while incorporating dummy variables for individual replicates (`log2FPKM ~ condition + replicate`)[cite: 4].\n", " * The condition coefficient (`log2FC_resid_sample`) was extracted, providing a highly robust estimate of condition-driven expression changes independent of specific sample outliers[cite: 4].\n", "\n", "5. **Combined Bias-Corrected Scoring:**\n", " * The residuals derived from both the gene-level (splicing) and sample-level models were converted into standardized Z-scores to normalize their scales[cite: 4].\n", " * A final `CombinedScore` was calculated as the mean of these Z-scores[cite: 4].\n", " * Genes were then re-ranked based on this combined metric, allowing for the precise identification of loci that demonstrate true, biology-driven transcriptional upregulation in the malignant cohort[cite: 4].\n", "\n", "```\n", "# This is formatted as code\n", "!pip install -q statsmodels\n", "import pandas as pd, numpy as np\n", "import statsmodels.api as sm\n", "from pathlib import Path\n", "\n", "# ---------- Paths ----------\n", "GENE_DIFF = \"/content/gene_exp.diff\"\n", "ISO_DIFF = \"/content/isoform_exp.diff\"\n", "READ_GROUP = \"/content/genes.read_group_tracking\"\n", "\n", "pd.set_option(\"display.max_columns\", 200)\n", "pd.set_option(\"display.width\", 180)\n", "\n", "def safe_log2fc(v2, v1, eps=1e-9):\n", " v2 = pd.to_numeric(v2, errors=\"coerce\")\n", " v1 = pd.to_numeric(v1, errors=\"coerce\")\n", " return np.log2((v2 + eps) / (v1 + eps))\n", "\n", "# ============================================================\n", "# 1) RAW CUFFDIFF: top upregulated genes (gene_exp.diff)\n", "# ============================================================\n", "gene = pd.read_csv(GENE_DIFF, sep=\"\\t\")\n", "gene[\"significant\"] = gene[\"significant\"].astype(str).str.lower()\n", "\n", "gene_sig = gene[gene[\"significant\"].eq(\"yes\")].copy()\n", "gene_sig[\"log2FC_raw\"] = safe_log2fc(gene_sig[\"value_2\"], gene_sig[\"value_1\"])\n", "gene_sig = (\n", " gene_sig\n", " .replace([np.inf, -np.inf], np.nan)\n", " .dropna(subset=[\"log2FC_raw\"])\n", ")\n", "\n", "top_raw = (\n", " gene_sig\n", " .sort_values(\"log2FC_raw\", ascending=False)\n", " [[\"gene\",\"gene_id\",\"log2FC_raw\",\"p_value\",\"q_value\",\"value_1\",\"value_2\"]]\n", " .head(50).reset_index(drop=True)\n", ")\n", "top_raw.insert(0, \"Rank_Raw\", np.arange(1, len(top_raw)+1))\n", "top_raw.to_csv(\"/content/TopUpregulated_Raw_Cuffdiff.csv\", index=False)\n", "\n", "print(\"\\n=== Top 10 Upregulated (RAW Cuffdiff) ===\")\n", "print(top_raw.head(10).to_string(index=False))\n", "\n", "# ============================================================\n", "# 2) ISOFORM SUMMARY -> per-gene splicing burden metrics\n", "# ============================================================\n", "iso = pd.read_csv(ISO_DIFF, sep=\"\\t\")\n", "iso[\"significant\"] = iso[\"significant\"].astype(str).str.lower().eq(\"yes\")\n", "\n", "if \"log2(fold_change)\" in iso.columns:\n", " iso[\"log2fc\"] = pd.to_numeric(iso[\"log2(fold_change)\"], errors=\"coerce\")\n", "else:\n", " iso[\"log2fc\"] = safe_log2fc(iso.get(\"value_2\"), iso.get(\"value_1\"))\n", "\n", "if {\"value_1\",\"value_2\"}.issubset(iso.columns):\n", " iso[\"mean_fpkm\"] = (\n", " pd.to_numeric(iso[\"value_1\"], errors=\"coerce\").fillna(0) +\n", " pd.to_numeric(iso[\"value_2\"], errors=\"coerce\").fillna(0)\n", " ) / 2.0\n", " iso = iso[iso[\"mean_fpkm\"] >= 1.0].copy()\n", "else:\n", " iso[\"mean_fpkm\"] = np.nan\n", "\n", "grp = iso.groupby(\"gene\", dropna=False)\n", "per_gene = pd.DataFrame({\n", " \"n_isoforms\": grp[\"test_id\"].count(),\n", " \"n_sig_iso\": grp[\"significant\"].sum(),\n", " \"n_sig_up\": grp.apply(\n", " lambda d: (pd.to_numeric(d[\"log2fc\"], errors=\"coerce\") > 0).sum(), include_groups=False\n", " ),\n", " \"n_sig_down\": grp.apply(\n", " lambda d: (pd.to_numeric(d[\"log2fc\"], errors=\"coerce\") < 0).sum(), include_groups=False\n", " ),\n", "}).reset_index().rename(columns={\"gene\":\"Gene\"})\n", "\n", "per_gene[\"SBI\"] = per_gene[\"n_sig_up\"] - per_gene[\"n_sig_down\"]\n", "per_gene.to_csv(\"/content/PerGene_Isoform_Metrics.csv\", index=False)\n", "\n", "# ============================================================\n", "# 3) GENE-LEVEL RESIDUALIZATION vs splicing burden\n", "# ============================================================\n", "de = gene_sig.rename(columns={\"gene\":\"Gene\"})[\n", " [\"Gene\",\"gene_id\",\"log2FC_raw\",\"q_value\",\"p_value\"]\n", "]\n", "merged_genelevel = de.merge(per_gene, on=\"Gene\", how=\"left\")\n", "\n", "for c in [\"n_isoforms\",\"n_sig_iso\",\"n_sig_up\",\"n_sig_down\",\"SBI\"]:\n", " merged_genelevel[c] = merged_genelevel[c].fillna(0).astype(int)\n", "\n", "X = sm.add_constant(\n", " merged_genelevel[[\"n_sig_iso\",\"SBI\"]].astype(float).fillna(0.0)\n", ")\n", "y = merged_genelevel[\"log2FC_raw\"].astype(float)\n", "\n", "model_gene = sm.OLS(y, X, missing=\"drop\").fit(cov_type=\"HC3\")\n", "merged_genelevel[\"log2FC_resid_splicing\"] = model_gene.resid\n", "\n", "merged_genelevel.to_csv(\"/content/Residualized_DE_splicingBurden.csv\", index=False)\n", "\n", "top_resid_splice = (\n", " merged_genelevel\n", " .sort_values(\"log2FC_resid_splicing\", ascending=False)\n", " [[\"Gene\",\"gene_id\",\"log2FC_resid_splicing\",\"log2FC_raw\",\"n_sig_iso\",\"SBI\"]]\n", " .head(50).reset_index(drop=True)\n", ")\n", "top_resid_splice.insert(\n", " 0, \"Rank_Residual_Splicing\", np.arange(1, len(top_resid_splice)+1)\n", ")\n", "top_resid_splice.to_csv(\n", " \"/content/TopUpregulated_Residualized_Splicing.csv\", index=False\n", ")\n", "\n", "print(\"\\n=== Top 10 Upregulated (Residualized vs Splicing Burden) ===\")\n", "print(top_resid_splice.head(10).to_string(index=False))\n", "\n", "Path(\"/content/Residualization_geneLevel_OLS_Summary.txt\").write_text(\n", " str(model_gene.summary())\n", ")\n", "\n", "# ============================================================\n", "# 4) SAMPLE-LEVEL RESIDUALIZATION using genes.read_group_tracking\n", "# ============================================================\n", "rg = pd.read_csv(READ_GROUP, sep=\"\\t\", low_memory=False)\n", "\n", "rg[\"condition\"] = rg[\"condition\"].astype(str).str.strip()\n", "cond_map = {\n", " \"c\": \"normal\",\n", " \"control\": \"normal\",\n", " \"t\": \"tumor\",\n", " \"test\": \"tumor\",\n", "}\n", "rg = rg[rg[\"condition\"].isin(cond_map.keys())].copy()\n", "rg[\"condition_norm\"] = rg[\"condition\"].map(cond_map)\n", "\n", "if \"FPKM\" not in rg.columns:\n", " raise ValueError(\"Expected column 'FPKM' not found in genes.read_group_tracking\")\n", "\n", "M = rg[[\"tracking_id\",\"condition_norm\",\"replicate\",\"FPKM\"]].copy()\n", "M = M.rename(columns={\n", " \"tracking_id\": \"gene_id\",\n", " \"condition_norm\": \"condition\"\n", "})\n", "\n", "M[\"log2FPKM\"] = np.log2(pd.to_numeric(M[\"FPKM\"], errors=\"coerce\") + 1e-6)\n", "\n", "gene_symbol_map = gene[[\"gene_id\",\"gene\"]].drop_duplicates()\n", "M = M.merge(gene_symbol_map, on=\"gene_id\", how=\"left\")\n", "M = M.rename(columns={\"gene\":\"Gene\"})\n", "\n", "M[\"cond_Tumor\"] = (M[\"condition\"] == \"tumor\").astype(int)\n", "M[\"replicate\"] = M[\"replicate\"].astype(str)\n", "rep_dummies = pd.get_dummies(M[\"replicate\"], prefix=\"rep\", drop_first=True)\n", "M = pd.concat([M, rep_dummies], axis=1)\n", "\n", "design_cols = [\"cond_Tumor\"] + [c for c in M.columns if c.startswith(\"rep_\")]\n", "\n", "coef_rows = []\n", "for g, df_g in M.groupby(\"Gene\", dropna=False):\n", " if df_g[\"cond_Tumor\"].nunique() < 2:\n", " continue\n", " Z = df_g[design_cols].copy().astype(float).fillna(0.0)\n", " Z = sm.add_constant(Z)\n", " yy = df_g[\"log2FPKM\"].astype(float)\n", " if len(df_g) <= Z.shape[1]:\n", " continue\n", " try:\n", " fit = sm.OLS(yy, Z, missing=\"drop\").fit(cov_type=\"HC3\")\n", " coef_rows.append({\n", " \"Gene\": g,\n", " \"log2FC_resid_sample\": float(fit.params.get(\"cond_Tumor\", np.nan))\n", " })\n", " except Exception:\n", " continue\n", "\n", "sample_resid = pd.DataFrame(coef_rows)\n", "\n", "# ============================================================\n", "# 5) MERGE ALL + BUILD COMBINED BIAS-CORRECTED SCORE\n", "# ============================================================\n", "final = merged_genelevel.merge(sample_resid, on=\"Gene\", how=\"left\")\n", "\n", "def zseries(s):\n", " s = pd.to_numeric(s, errors=\"coerce\")\n", " m = s.mean(skipna=True)\n", " sd = s.std(skipna=True, ddof=0)\n", " if not np.isfinite(sd) or sd == 0:\n", " return pd.Series(np.nan, index=s.index)\n", " return (s - m) / sd\n", "\n", "final[\"z_splicing\"] = zseries(final[\"log2FC_resid_splicing\"])\n", "final[\"z_sample\"] = zseries(final[\"log2FC_resid_sample\"])\n", "\n", "def combined_row(row):\n", " vals = []\n", " if pd.notna(row[\"z_splicing\"]):\n", " vals.append(row[\"z_splicing\"])\n", " if pd.notna(row[\"z_sample\"]):\n", " vals.append(row[\"z_sample\"])\n", " if not vals:\n", " return np.nan\n", " return float(np.mean(vals))\n", "\n", "final[\"CombinedScore\"] = final.apply(combined_row, axis=1)\n", "\n", "final.to_csv(\"/content/Residualization_Combined_Table.csv\", index=False)\n", "\n", "top_resid_sample = (\n", " final.dropna(subset=[\"log2FC_resid_sample\"])\n", " .sort_values(\"log2FC_resid_sample\", ascending=False)\n", " [[\"Gene\",\"gene_id\",\"log2FC_resid_sample\"]]\n", " .head(50).reset_index(drop=True)\n", ")\n", "top_resid_sample.insert(\n", " 0, \"Rank_Residual_Sample\", np.arange(1, len(top_resid_sample)+1)\n", ")\n", "top_resid_sample.to_csv(\n", " \"/content/TopUpregulated_Residualized_SampleLevel.csv\", index=False\n", ")\n", "\n", "top_combined = (\n", " final.dropna(subset=[\"CombinedScore\"])\n", " .sort_values(\"CombinedScore\", ascending=False)\n", " [[\"Gene\",\"gene_id\",\"CombinedScore\",\n", " \"log2FC_resid_sample\",\"log2FC_resid_splicing\",\"log2FC_raw\"]]\n", " .head(50).reset_index(drop=True)\n", ")\n", "top_combined.insert(0, \"Rank_Combined\", np.arange(1, len(top_combined)+1))\n", "top_combined.to_csv(\n", " \"/content/TopUpregulated_BiasCorrected_Combined.csv\", index=False\n", ")\n", "\n", "t10_raw = (top_raw.head(10)[[\"Rank_Raw\",\"gene\",\"log2FC_raw\"]]\n", " .rename(columns={\"gene\":\"Gene\",\"log2FC_raw\":\"Raw_log2FC\"}))\n", "t10_spl = top_resid_splice.head(10)[\n", " [\"Rank_Residual_Splicing\",\"Gene\",\"log2FC_resid_splicing\"]\n", "]\n", "t10_smp = top_resid_sample.head(10)[\n", " [\"Rank_Residual_Sample\",\"Gene\",\"log2FC_resid_sample\"]\n", "]\n", "cmp = (t10_raw\n", " .merge(t10_spl, on=\"Gene\", how=\"outer\")\n", " .merge(t10_smp, on=\"Gene\", how=\"outer\"))\n", "cmp.to_csv(\"/content/Top10_Raw_vs_Residualized_BothModels.csv\", index=False)\n", "\n", "print(\"\\n=== Top 10 Upregulated (Residualized, Sample-level) ===\")\n", "print(top_resid_sample.head(10).to_string(index=False))\n", "\n", "print(\"\\n=== Top 10 Upregulated (Combined bias-corrected score) ===\")\n", "print(top_combined.head(10).to_string(index=False))\n", "\n", "print(\"\\nFiles written:\")\n", "for p in [\n", " \"/content/TopUpregulated_Raw_Cuffdiff.csv\",\n", " \"/content/PerGene_Isoform_Metrics.csv\",\n", " \"/content/Residualized_DE_splicingBurden.csv\",\n", " \"/content/TopUpregulated_Residualized_Splicing.csv\",\n", " \"/content/Residualization_Combined_Table.csv\",\n", " \"/content/TopUpregulated_Residualized_SampleLevel.csv\",\n", " \"/content/Top10_Raw_vs_Residualized_BothModels.csv\",\n", " \"/content/TopUpregulated_BiasCorrected_Combined.csv\",\n", " \"/content/Residualization_geneLevel_OLS_Summary.txt\",\n", "]:\n", " print(\" -\", p)\n", "```\n", "\n" ], "metadata": { "id": "tIv_57Xsmqgz" } }, { "cell_type": "markdown", "source": [ "## CODE CELL 4\n", "2. Environment Setup and High-Throughput Pipeline Initialization\n", "To perform the alternative splicing analysis, rMATS (v4.1.2) and MntJULiP were utilized. While these tools were successful in our previous esophageal study, the massive scale of this colorectal dataset (54 high-depth samples, ~315 GB) introduced unique infrastructure challenges in the Google Cloud Colab Enterprise environment.\n", "\n", "2.1 Hardware Optimization: Transition to 1,000 GB Persistent Storage\n", "Initial attempts using gcsfuse to virtually mount Google Cloud Storage buckets resulted in significant I/O latency, leading to statistical \"timeouts\" in rMATS (manifesting as P-values of 1.0) and ZeroDivisionError crashes in MntJULiP’s Bayesian model. To resolve this, the runtime was upgraded to a 1,000 GB Standard Persistent Disk. All BAM files were physically transferred from the GCS bucket to the local /content directory using the multi-threaded gsutil -m cp command.\n", "\n", "2.2 Conda Environment Setup and Dependency Resolution\n", "Deploying legacy bioinformatics tools in modern environments requires strict version isolation. We utilized Miniconda to establish isolated environments and circumvent dependency conflicts:\n", "\n", "rmats_env (Conda Installation): A dedicated Conda environment was created using Python 3.8. To ensure stability on the persistent disk, rMATS and its binary dependencies (STAR, Samtools) were installed directly into this environment from the bioconda and conda-forge channels.\n", "\n", "mntjulip_env (Virtualenv): Configured with strictly pinned legacy dependencies (numpy==1.23.5, pandas==1.3.5) to prevent mathematical variance errors.\n", "\n", "Automated Source Correction: The MntJULiP models.py was programmatically patched using sed to replace deprecated Numpy aliases (np.int, np.float) with standard Python types, preventing runtime crashes.\n", "\n", "2.3 Master Installation and Execution Block\n", "The following code represents the finalized implementation used to establish the Conda environments and execute the analysis:\n", "\n", "\n", "\n", "```\n", "# This is formatted as code\n", "# ==============================================================================\n", "# PHASE 1: CONDA INSTALLATION & ENVIRONMENT SETUP\n", "# ==============================================================================\n", "# 1. Install Miniconda and establish the rmats_env (Python 3.8)\n", "!wget -q https://repo.anaconda.com/miniconda/Miniconda3-latest-Linux-x86_64.sh\n", "!bash Miniconda3-latest-Linux-x86_64.sh -bfp /usr/local\n", "!/usr/local/bin/conda create -y -n rmats_env python=3.8 -c conda-forge -c bioconda rmats=4.1.2 star=2.7.10a samtools=1.15.1\n", "\n", "# 2. Establish mntjulip_env and apply source code patches\n", "!python3.8 -m venv /content/mntjulip_env\n", "!source /content/mntjulip_env/bin/activate && pip install -q \"numpy==1.23.5\" \"pandas==1.3.5\" \"statsmodels==0.12.2\"\n", "!sed -i 's/np\\.int\\b/int/g' /content/MntJULiP/models.py\n", "\n", "# ==============================================================================\n", "# PHASE 2: LOCAL MANIFEST GENERATION (3-COLUMN GROUPING)\n", "# ==============================================================================\n", "import os\n", "t_folder, n_folder = \"/content/BAMs\", \"/content/BAMSNORMAL\"\n", "t_files = sorted([os.path.join(t_folder, f) for f in os.listdir(t_folder) if f.endswith(\".bam\") and \"SRR8616044\" not in f])\n", "n_files = sorted([os.path.join(n_folder, f) for f in os.listdir(n_folder) if f.endswith(\".bam\")])\n", "\n", "# MntJULiP 3-column manifest\n", "with open(\"/content/all_bams_manifest_fixed.txt\", \"w\") as f:\n", " f.write(\"path\\tsample\\tcondition\\n\")\n", " for p in t_files: f.write(f\"{p}\\t{os.path.basename(p).split('_')[0]}\\tTumor\\n\")\n", " for p in n_files: f.write(f\"{p}\\t{os.path.basename(p).split('_')[0]}\\tNormal\\n\")\n", "\n", "# ==============================================================================\n", "# PHASE 3: MASTER ANALYSIS EXECUTION (rMATS -> MntJULiP)\n", "# ==============================================================================\n", "# STAGE 1: rMATS execution using absolute Conda paths\n", "!export PYTHONPATH=/usr/local/envs/rmats_env/lib/python3.8/site-packages:$PYTHONPATH && \\\n", " /usr/local/envs/rmats_env/bin/python /usr/local/envs/rmats_env/bin/rmats.py \\\n", " --b1 /content/b1.txt --b2 /content/b2.txt \\\n", " --gtf /content/colorectal_ccle/reference_genome/gencode.v40.annotation.gtf \\\n", " --od /content/rmats_out --tmp /content/rmats_tmp \\\n", " -t paired --readLength 150 --variable-read-length --nthread 30 --allow-clipping\n", "\n", "\n", "```\n", "\n" ], "metadata": { "id": "pr4Eyti47mVK" } }, { "cell_type": "markdown", "source": [ "# CODE CELL-5\n", "MntJULiP Analysis Workflow: Splicing Quantification and Differential Analysis\n", "This analysis utilizes the official MntJULiP toolkit (source: https://github.com/splicebox/MntJULiP). The installation process strictly follows the developer-provided protocol by extracting the source and executing the setup.py script. This step is critical as it compiles the internal C++ binaries (specifically the junc tool) necessary for high-throughput splice junction extraction from large-scale BAM files (315 GB).\n", "\n", "To ensure the stability of the Bayesian modeling and data merging stages for a 54-sample cohort, the environment is manually aligned to legacy versions of Dask (2023.5.0), Pandas (1.5.3), and NumPy (1.23.5). These adjustments resolve modern dependency conflicts, specifically the AssertionError in the Dask expression engine and the AttributeError from the deprecated np.int alias in the MntJULiP source code. The workflow is optimized for high-memory Google Colab instances, utilizing 16 threads and disabling the Dask Query Planner to finalize the differential splicing results in the diff_introns.txt output.\n", "\n", "```\n", "# This is formatted as code\n", "# ==============================================================================\n", "# 1. INSTALLATION & ENVIRONMENT ALIGNMENT (Official SpliceBox Protocol)\n", "# ==============================================================================\n", "# Download and extract the official MntJULiP source from SpliceBox\n", "wget https://github.com/splicebox/MntJULiP/archive/refs/tags/v1.5.2.tar.gz -O MntJULiP-1.5.2.tar.gz\n", "tar -xzf MntJULiP-1.5.2.tar.gz\n", "mkdir -p /content/mntjulip_out\n", "\n", "# Align dependencies for legacy compatibility with the MntJULiP 1.5.2 codebase\n", "pip install dask==2023.5.0 distributed==2023.5.0 pandas==1.5.3 numpy==1.23.5 \\\n", " pysam scikit-learn pystan==2.19.1.1 --quiet\n", "\n", "# Execute official compilation of C++ 'junc' binary via setup.py\n", "cd /content/MntJULiP-1.5.2\n", "python setup.py install --user --prefix=\n", "\n", "# ==============================================================================\n", "# 2. EXECUTION & DIFFERENTIAL SPLICING QUANTIFICATION\n", "# ==============================================================================\n", "cd /content\n", "export DASK_QUERY_PLANNING=False\n", "export DASK_DISTRIBUTED__QUERIES__PLANNING=False\n", "\n", "# Run main analysis pipeline; automatically reuses existing .splice.gz files\n", "# in the output directory to avoid redundant 7-hour extraction.\n", "python -u /content/MntJULiP-1.5.2/run.py \\\n", " --bam-list /content/mntjulip_manifest_fixed.tsv \\\n", " --anno-file /content/colorectal_ccle/reference_genome/gencode.v40.annotation.gtf \\\n", " --out-dir /content/mntjulip_out \\\n", " --num-threads 16\n", "```\n", "\n" ], "metadata": { "id": "liXv3C766ByK" } }, { "cell_type": "markdown", "source": [ "# CODE CELL-6\n", "3. Execution of Alternative Splicing Pipelines (rMATS and MntJULiP)\n", "Following the successful initialization of the isolated software environments, the core alternative splicing quantification was executed. Because rMATS and MntJULiP utilize fundamentally different statistical models and input requirements, the execution script was designed to toggle dynamically between their respective virtual environments while managing the heavy I/O demands of the cloud-mounted BAM files.\n", "\n", "The execution block performs the following automated steps:\n", "\n", "rMATS Execution (Conda Environment): The script activates the dedicated rmats_env (Python 3.8). To prevent memory saturation and storage overflow when processing 43 high-depth genomes, a dedicated temporary directory (--tmp /content/rmats_tmp) is explicitly defined. This ensures that the massive volume of intermediate mathematical files generated during junction mapping does not crash the local Colab instance before the final PSI values are calculated.\n", "\n", "Virtual Environment Toggling: The Conda environment is cleanly deactivated, and the strictly pinned Python virtual environment (mntjulip_env) is sourced to safely transition to the secondary pipeline.\n", "\n", "Dynamic Manifest Generation for MntJULiP: Unlike rMATS, which accepts comma-separated lists of BAM files, MntJULiP's Bayesian models strictly require a tabular clinical manifest. The script uses Bash text-replacement utilities (tr) to dynamically convert the tumor (b1.txt) and normal (b2.txt) lists into a single, combined manifest (all_bams_manifest.txt), mapping each specific BAM file path to its corresponding biological condition (\"Tumor\" or \"Normal\").\n", "\n", "MntJULiP Execution: The tool is launched utilizing the dynamically generated manifest and the standard reference genome (GENCODE v40) to apply advanced filtration to the discovered splice junctions.\n", "```%%bash\n", "# Clean up the corrupted/empty files from the previous failed run\n", "rm -rf /content/rmats_out/*\n", "rm -rf /content/rmats_tmp/*\n", "mkdir -p /content/rmats_out\n", "mkdir -p /content/rmats_tmp\n", "mkdir -p /content/mntjulip_out\n", "\n", "# ==============================================================================\n", "# 1. EXECUTE rMATS (Inside Conda Environment)\n", "# ==============================================================================\n", "echo \"[STATUS] Activating rmats_env and starting rMATS analysis...\"\n", "eval \"$(/usr/local/bin/conda shell.bash hook)\"\n", "conda activate rmats_env\n", "\n", "# Execute rMATS (ADDED --variable-read-length to save the 4.4 billion reads)\n", "rmats.py --b1 /content/b1.txt \\\n", " --b2 /content/b2.txt \\\n", " --gtf /content/colorectal_ccle/reference_genome/gencode.v40.annotation.gtf \\\n", " --od /content/rmats_out \\\n", " --tmp /content/rmats_tmp \\\n", " -t paired \\\n", " --readLength 150 \\\n", " --variable-read-length \\\n", " --cstat 0.0001 \\\n", " --nthread 30\n", "\n", "echo \"[SUCCESS] rMATS analysis complete!\"\n", "\n", "# ==============================================================================\n", "# 2. EXECUTE MntJULiP (Inside Python Virtual Environment)\n", "# ==============================================================================\n", "echo \"[STATUS] Activating mntjulip_env and starting MntJULiP analysis...\"\n", "conda deactivate\n", "source /content/mntjulip_env/bin/activate\n", "\n", "# Quick patch to install the missing Dask module\n", "echo \"[STATUS] Patching missing MntJULiP dependencies...\"\n", "pip install -q dask\n", "\n", "# Build the exact Manifest format required by MntJULiP (BAM_PATH CONDITION)\n", "echo \"[STATUS] Building MntJULiP Clinical Manifest...\"\n", "> /content/all_bams_manifest.txt\n", "for bam in $(tr ',' ' ' < /content/b1.txt); do\n", " echo -e \"$bam\\tTumor\" >> /content/all_bams_manifest.txt\n", "done\n", "for bam in $(tr ',' ' ' < /content/b2.txt); do\n", " echo -e \"$bam\\tNormal\" >> /content/all_bams_manifest.txt\n", "done\n", "\n", "# Execute MntJULiP\n", "python /content/MntJULiP/run.py \\\n", " --bam-list /content/all_bams_manifest.txt \\\n", " --anno-file /content/colorectal_ccle/reference_genome/gencode.v40.annotation.gtf \\\n", " --out-dir /content/mntjulip_out \\\n", " --num-threads 30\n", "\n", "echo \"==============================================================================\"\n", "echo \"[SUCCESS] Both rMATS and MntJULiP pipelines have finished executing!\"\n", "echo \"===============\n", "```\n", "\n" ], "metadata": { "id": "WTYBUjVkfcBM" } }, { "cell_type": "markdown", "source": [ "# CODE CELL 7\n", "# Supplementary Methods: Integrated Non-Coding RNA Splicing Pipeline\n", "\n", "### 1. Overview\n", "All splicing analyses were performed using a unified, fully reproducible Python pipeline. This workflow integrates **rMATS** (exon-centric alternative splicing) and **MntJULiP** (intron usage and differential intron regulation). By bypassing global differential expression modeling and focusing exclusively on structural transcriptomic variations, the pipeline isolates high-confidence regulatory splicing events unique to the tumor microenvironment.\n", "\n", "### 2. Computational Environment & Software Versions\n", "Analyses were performed on an enterprise-grade compute instance to ensure reproducibility and consistency across large-scale sequencing datasets.\n", "\n", "* **Operating System:** Ubuntu 22.04 LTS (x86_64)\n", "* **Python Version:** 3.10.x\n", "* **rMATS-Turbo:** v4.1.2 (for SE, RI, MXE, A5SS, and A3SS event detection)\n", "* **MntJULiP:** v1.5.2 (for Differential Splicing Analysis [DSA] and Differential Spliced Regions [DSR])\n", "* **Annotation:** GENCODE v40 (GRCh38.p13)\n", "\n", "### 3. Data Processing and Annotation\n", "To isolate the non-coding transcriptome, a master index of non-coding RNA biotypes was established using the GENCODE v40 GTF catalog. The pipeline explicitly filtered for transcripts classified under the following biotypes: `lncRNA`, `lincRNA`, `macro_lncRNA`, `bidirectional_promoter_lncRNA`, `antisense`, `sense_intronic`, `sense_overlapping`, `3prime_overlapping_ncRNA`, `processed_transcript`, `pseudogene` (including `transcribed`, `processed`, `unprocessed`, `rRNA`, `tRNA`, `snRNA`, `snoRNA` pseudogenes), and `TEC` (To Be Experimentally Confirmed).\n", "\n", "### 4. Statistical Framework & Thresholding\n", "To maintain high confidence in the candidate pool, a stringent dual-filter was applied to all detected events:\n", "* **Significance:** FDR (for rMATS) and q-value (for MntJULiP) $\\le 0.05$.\n", "* **Magnitude:** Absolute Inclusion Level Difference ($|\\Delta \\text{PSI}|$) $\\ge 0.20$.\n", "\n", "### 5. Unified Integration Logic (DSR-to-DSA Fallback)\n", "The pipeline constructs a master intron table by merging two distinct MntJULiP outputs to maximize coverage while preserving structural context:\n", "1. **DSR (`diff_spliced_introns.txt`):** Prioritized for identifying regulated intron clusters (`group_id`, `group_structure`, `group_loc`).\n", "2. **DSA (`diff_introns.txt`):** Used as the primary source for full-coverage significance statistics.\n", "\n", "Introns identified in DSR are annotated with regulatory group information (MntJULiP_source = \"DSR\"), while those unique to DSA are annotated with statistical significance metrics only (MntJULiP_source = \"DSA\").\n", "\n", "### 6. Event–Intron Anchor Matching Algorithm\n", "To identify the structural correspondence between exon-centric rMATS events and intron-centric MntJULiP dysregulation, a deterministic anchor-based matching algorithm was implemented:\n", "* **Anchor Definition:** Splice event decision points (e.g., alternative exon starts for MXE, or upstream exon starts for RI) were calculated as anchors.\n", "* **Matching Hierarchy:** Overlapping candidates were prioritized based on:\n", " 1. Anchor co-localization ($\\le 1$ bp).\n", " 2. Group structure orientation.\n", " 3. Significance (q-value).\n", " 4. Overlap length (bp).\n", "\n", "### 7. Structural and Clinical Annotation\n", "Each validated event was mapped to GENCODE v40 models. The pipeline classified events as **\"Annotated (GENCODE)\"** if splicing junctions matched existing canonical transcripts or **\"Novel (Unannotated)\"** if they deviated from known models. This allows for the differentiation between well-characterized splicing variability and previously unreported, cancer-associated structural variants.\n", "\n", "### 8. Reproducibility Statement\n", "The entire workflow was executed as a modular, deterministic Python script. It utilizes fixed external archives for all reference files (GENCODE) and software (rMATS, MntJULiP). The logic avoids external state or pseudo-randomness, ensuring that subsequent executions of this pipeline on the same raw input datasets produce identical event-level output files.\n", "\n", "```\n", "# This is formatted as code\n", "# ================================================================\n", "# Unified rMATS ↔ MntJULiP pipeline (DSR → DSA fallback)\n", "# + GENCODE intron/exon annotation\n", "#\n", "# Outputs:\n", "# /content/CRC_SplicingFrequencyRanking_byGene.csv\n", "# /content/CRC_SplicingRegulationModes_TopMostSpliced.csv\n", "# /content/CRC__With_Annotation.csv\n", "# /content/CRC_TopGenes_Combined_With_Annotation.csv\n", "# ================================================================\n", "\n", "import os, re, io, gzip, zipfile, math, json, shutil, time\n", "from pathlib import Path\n", "\n", "import pandas as pd\n", "import numpy as np\n", "\n", "pd.set_option(\"display.max_columns\", 200)\n", "pd.set_option(\"display.width\", 260)\n", "\n", "# ---------------- USER PATHS ----------------\n", "RMATS_ZIP = \"/content/rmats_results_colorectal.zip\"\n", "MNT_ZIP = \"/content/mntjulip_results_full_colorectal.zip\"\n", "GTF_ZIP = \"/content/gencode.v40.annotation.gtf.zip\"\n", "\n", "# ---------------- OUTPUTS -------------------\n", "OUT_FREQ_RANK = \"/content/CRC_SplicingFrequencyRanking_byGene.csv\"\n", "OUT_REGMODES = \"/content/CRC_SplicingRegulationModes_TopMostSpliced.csv\"\n", "OUT_COMBINED_ANNOT = \"/content/CRC_TopGenes_Combined_With_Annotation.csv\"\n", "\n", "# ---------------- PARAMS --------------------\n", "ALPHA_FDR = 0.05\n", "MIN_ABS_DPSI = 0.20\n", "MNT_Q = 0.05 \n", "NONCODING_ONLY = True \n", "TOP_PICK_N = 10 # Strictly data-driven: automatically pick the top 10\n", "COLOC_TOL = 1\n", "GENCODE_SLOP = 0 \n", "\n", "# --------------- UTILITIES ------------------\n", "def _find_first(cols, *cands):\n", " low = {c.lower(): c for c in cols}\n", " for c in cands:\n", " if c.lower() in low:\n", " return low[c.lower()]\n", " return None\n", "\n", "def _gi(v):\n", " try:\n", " if v is None or (isinstance(v, float) and np.isnan(v)):\n", " return None\n", " return int(float(str(v).strip()))\n", " except Exception:\n", " return None\n", "\n", "def _span_overlap(a, b):\n", " s = max(a[0], b[0]); e = min(a[1], b[1])\n", " return max(0, e - s)\n", "\n", "def _normalize_chr(x):\n", " return str(x).replace(\"chr\", \"\")\n", "\n", "EVENT_RE = re.compile(r\"(?:^|/|\\.)(SE|RI|A3SS|A5SS|MXE|TXI)\\.MATS\\.(?:JC|JCEC)\\.txt$\", re.IGNORECASE)\n", "\n", "def _largest_entries_by_key(zf, key_fn):\n", " bykey = {}\n", " for info in zf.infolist():\n", " key = key_fn(info)\n", " if key is None:\n", " continue\n", " if (key not in bykey) or (info.file_size > bykey[key].file_size):\n", " bykey[key] = info\n", " return bykey\n", "\n", "# ----------- GTF open & noncoding filter set -----------\n", "def _open_gtf_any_from_zip(zip_path):\n", " with zipfile.ZipFile(zip_path, 'r') as z:\n", " gtf_name = next((n for n in z.namelist() if n.endswith(\".gtf\")), None)\n", " if gtf_name is None:\n", " raise FileNotFoundError(\"No .gtf inside GTF_ZIP\")\n", " return io.TextIOWrapper(z.open(gtf_name), encoding=\"utf-8\", errors=\"ignore\")\n", "\n", "def load_noncoding_symbol_set(gtf_zip):\n", " if not NONCODING_ONLY:\n", " return None\n", " nc_types = {\n", " \"lncRNA\",\"lincRNA\",\"macro_lncRNA\",\"bidirectional_promoter_lncRNA\",\n", " \"antisense\",\"sense_intronic\",\"sense_overlapping\",\"3prime_overlapping_ncRNA\",\n", " \"processed_transcript\",\"non_coding\",\"ncRNA\",\n", " \"pseudogene\",\"transcribed_pseudogene\",\"processed_pseudogene\",\"unprocessed_pseudogene\",\n", " \"TEC\"\n", " }\n", " nc = set()\n", " with _open_gtf_any_from_zip(gtf_zip) as fh:\n", " for line in fh:\n", " if not line or line.startswith(\"#\"):\n", " continue\n", " p = line.rstrip(\"\\n\").split(\"\\t\")\n", " if len(p) < 9 or p[2] != \"gene\":\n", " continue\n", " attrs = p[8]\n", " def grab(k):\n", " m = re.search(rf'{k}\\s+\"([^\"]+)\"', attrs)\n", " return m.group(1) if m else None\n", " gname = grab(\"gene_name\")\n", " biotype = grab(\"gene_type\") or grab(\"gene_biotype\")\n", " if gname and biotype and biotype in nc_types:\n", " nc.add(gname)\n", " return nc\n", "\n", "NONCODING_SET = load_noncoding_symbol_set(GTF_ZIP) if NONCODING_ONLY else None\n", "\n", "# ----------- rMATS loaders (all + significant) -----------\n", "def _read_tsv(zf, info):\n", " with zf.open(info) as f:\n", " return pd.read_csv(f, sep=\"\\t\")\n", "\n", "def _event_key_from_row(evt, row, colnames):\n", " if evt == \"SE\":\n", " cols = [\"exonStart_0base\",\"exonEnd\",\"upstreamES\",\"upstreamEE\",\"downstreamES\",\"downstreamEE\"]\n", " elif evt == \"RI\":\n", " cols = [\"exonStart_0base\",\"exonEnd\",\"upstreamEE\",\"downstreamES\"]\n", " elif evt in (\"A3SS\",\"A5SS\"):\n", " cols = [\"longExonStart_0base\",\"longExonEnd\",\"flankingES\",\"flankingEE\",\"shortES\",\"shortEE\"]\n", " elif evt == \"MXE\":\n", " cols = [\"1stExonStart_0base\",\"1stExonEnd\",\"2ndExonStart_0base\",\"2ndExonEnd\",\n", " \"upstreamES\",\"upstreamEE\",\"downstreamES\",\"downstreamEE\"]\n", " else:\n", " cols = []\n", " cc = _find_first(colnames, \"chr\",\"chrom\",\"chromosome\") or \"chr\"\n", " ss = _find_first(colnames, \"strand\") or \"strand\"\n", " parts = [str(row.get(cc, \"\")).replace(\"chr\", \"\").strip(),\n", " str(row.get(ss, \"\")).strip()]\n", " for c in cols:\n", " v = row.get(c, None)\n", " parts.append(\"\" if v is None else (\"\" if _gi(v) is None else str(_gi(v))))\n", " parts.append(evt)\n", " return \"_\".join(parts)\n", "\n", "def load_all_rmats_events(zip_path):\n", " with zipfile.ZipFile(zip_path, 'r') as z:\n", " def k(info):\n", " m = EVENT_RE.search(info.filename)\n", " if not m:\n", " return None\n", " evt = m.group(1).upper()\n", " variant = \"JCEC\" if info.filename.endswith(\".MATS.JCEC.txt\") else \"JC\"\n", " return (evt, variant)\n", " chosen = _largest_entries_by_key(z, k)\n", " rows = []\n", " for (evt, variant), info in sorted(chosen.items()):\n", " df = _read_tsv(z, info)\n", " df.columns = [c.strip() for c in df.columns]\n", " gcol = _find_first(df.columns, \"geneSymbol\",\"GeneSymbol\",\"gene_name\",\"gene\",\"Gene\",\"GeneID\",\"gene_id\")\n", " if gcol is None:\n", " continue\n", " if NONCODING_SET is not None:\n", " df = df[df[gcol].astype(str).isin(NONCODING_SET)]\n", " if df.empty:\n", " continue\n", " colnames = list(df.columns)\n", " df[\"__event_key\"] = df.apply(lambda r: _event_key_from_row(evt, r, colnames), axis=1)\n", " df[\"gene\"] = df[gcol].astype(str)\n", " df[\"__evt_type\"] = evt\n", " df[\"__variant\"] = variant\n", " keep = [\"gene\",\"__event_key\",\"__evt_type\",\"__variant\",\n", " \"chr\",\"Chromosome\",\"chrom\",\"strand\",\n", " \"exonStart_0base\",\"exonEnd\",\"upstreamES\",\"upstreamEE\",\"downstreamES\",\"downstreamEE\",\n", " \"longExonStart_0base\",\"longExonEnd\",\"flankingES\",\"flankingEE\",\"shortES\",\"shortEE\",\n", " \"1stExonStart_0base\",\"1stExonEnd\",\"2ndExonStart_0base\",\"2ndExonEnd\"]\n", " rows.append(df[[c for c in keep if c in df.columns]].drop_duplicates())\n", " if not rows:\n", " return pd.DataFrame(columns=[\"gene\",\"__event_key\",\"__evt_type\",\"__variant\"])\n", " return pd.concat(rows, ignore_index=True).drop_duplicates([\"__event_key\"])\n", "\n", "def load_sig_rmats_events(zip_path, alpha=0.05, min_abs_dpsi=0.20):\n", " with zipfile.ZipFile(zip_path, 'r') as z:\n", " def k(info):\n", " m = EVENT_RE.search(info.filename)\n", " if not m:\n", " return None\n", " evt = m.group(1).upper()\n", " variant = \"JCEC\" if info.filename.endswith(\".MATS.JCEC.txt\") else \"JC\"\n", " return (evt, variant)\n", " chosen = _largest_entries_by_key(z, k)\n", " parts = []\n", " for (evt, variant), info in sorted(chosen.items()):\n", " df = _read_tsv(z, info)\n", " df.columns = [c.strip() for c in df.columns]\n", " gcol = _find_first(df.columns, \"geneSymbol\",\"GeneSymbol\",\"gene_name\",\"gene\",\"Gene\",\"GeneID\",\"gene_id\")\n", " fdr_col = _find_first(df.columns, \"FDR\",\"fdr\",\"qvalue\",\"q_value\",\"adjp\",\"adj_p\")\n", " dpsi_col= _find_first(df.columns, \"IncLevelDifference\",\"incLevelDifference\",\"deltaPSI\",\"dPSI\",\"delta_psi\")\n", " if gcol is None or fdr_col is None or dpsi_col is None:\n", " continue\n", " if NONCODING_SET is not None:\n", " df = df[df[gcol].astype(str).isin(NONCODING_SET)]\n", " if df.empty:\n", " continue\n", " df[fdr_col] = pd.to_numeric(df[fdr_col], errors=\"coerce\")\n", " df[dpsi_col]= pd.to_numeric(df[dpsi_col], errors=\"coerce\")\n", " df = df[(df[fdr_col] <= alpha) & (df[dpsi_col].abs() >= min_abs_dpsi)]\n", " if df.empty:\n", " continue\n", " colnames = list(df.columns)\n", " df[\"__event_key\"] = df.apply(lambda r: _event_key_from_row(evt, r, colnames), axis=1)\n", " df[\"gene\"] = df[gcol].astype(str)\n", " df[\"__evt_type\"] = evt\n", " df[\"__variant\"] = variant\n", " parts.append(df[[\"gene\",\"__event_key\",\"__evt_type\",\"__variant\",fdr_col,dpsi_col]].rename(\n", " columns={fdr_col:\"fdr\", dpsi_col:\"delta_psi\"}))\n", " if not parts:\n", " return pd.DataFrame()\n", " sig = pd.concat(parts, ignore_index=True)\n", " sig[\"__var_rank\"] = sig[\"__variant\"].map({\"JCEC\":0,\"JC\":1}).fillna(2)\n", " sig = (sig.sort_values([\"__event_key\",\"__var_rank\",\"fdr\"])\n", " .drop_duplicates(\"__event_key\")\n", " .drop(columns=\"__var_rank\"))\n", " return sig\n", "\n", "# ------------- MntJULiP: detect group columns -------------\n", "def detect_group_columns(df):\n", " def pick(regex, numeric=None, prefer_vals=None):\n", " cand = [c for c in df.columns if re.search(regex, c.lower())]\n", " if prefer_vals is not None:\n", " scored = []\n", " for c in cand:\n", " vals = df[c].dropna().astype(str).head(200)\n", " hit = any(v in {\"o\",\"i\"} for v in vals)\n", " scored.append((2 if hit else 1, c))\n", " cand = [c for _, c in sorted(scored, reverse=True)]\n", " if numeric is True:\n", " cand = [c for c in cand if pd.api.types.is_numeric_dtype(df[c])]\n", " if numeric is False:\n", " cand = [c for c in cand if not pd.api.types.is_numeric_dtype(df[c])]\n", " return cand[0] if cand else None\n", "\n", " gid = pick(r\"(group.*id|^gid$|cluster)\", numeric=False)\n", " gstruct = pick(r\"(group.*structure|structure|orient|^group$)\", prefer_vals={\"o\",\"i\"})\n", " gloc = pick(r\"(group.*loc|^loc$|locus|anchor|junction.*pos)\", numeric=True)\n", " gq = pick(r\"(group.*q|posterior|post_prob|prob|score)$\", numeric=True)\n", " return gid, gstruct, gloc, gq\n", "\n", "def _read_mnt(zf, info):\n", " with zf.open(info) as f:\n", " if info.filename.endswith(\".gz\"):\n", " return pd.read_csv(gzip.open(io.BytesIO(f.read()), \"rt\"), sep=\"\\t\")\n", " return pd.read_csv(f, sep=\"\\t\")\n", "\n", "def load_mnt_sig_introns_union(zip_path, alpha=0.05, min_abs_dpsi=0.20):\n", " with zipfile.ZipFile(zip_path, 'r') as z:\n", " p_dsr = None\n", " p_dsa = None\n", " for info in z.infolist():\n", " name = info.filename\n", " if name.endswith(\"diff_spliced_introns.txt\"):\n", " p_dsr = info\n", " elif name.endswith(\"diff_introns.txt\"):\n", " p_dsa = info\n", "\n", " if p_dsa is None:\n", " return pd.DataFrame()\n", "\n", " dall = _read_mnt(z, p_dsa)\n", " dall.columns = [c.strip() for c in dall.columns]\n", " chr_a = _find_first(dall.columns, \"chrom\",\"chr\",\"chromosome\")\n", " st_a = _find_first(dall.columns, \"start\",\"loc_start\",\"s\")\n", " en_a = _find_first(dall.columns, \"end\",\"loc_end\",\"e\")\n", " sd_a = _find_first(dall.columns, \"strand\")\n", " q_a = _find_first(dall.columns, \"q_value\",\"qvalue\",\"q\",\"FDR\",\"fdr\")\n", " g_a = _find_first(dall.columns, \"gene_name\",\"gene\",\"Gene\")\n", "\n", " if not all([chr_a, st_a, en_a, sd_a, q_a]):\n", " return pd.DataFrame()\n", "\n", " dall_proc = dall[[chr_a, st_a, en_a, sd_a, q_a] + ([g_a] if g_a else [])].copy()\n", " dall_proc = dall_proc.rename(columns={\n", " chr_a:\"chrom\", st_a:\"start\", en_a:\"end\", sd_a:\"strand\", q_a:\"MntJULiP_q\", (g_a or \"gene\"):\"gene_raw\"\n", " })\n", " dall_proc[\"chrom\"] = dall_proc[\"chrom\"].astype(str).str.replace(\"chr\",\"\")\n", " dall_proc[\"start\"] = pd.to_numeric(dall_proc[\"start\"], errors=\"coerce\")\n", " dall_proc[\"end\"] = pd.to_numeric(dall_proc[\"end\"], errors=\"coerce\")\n", " dall_proc[\"strand\"] = dall_proc[\"strand\"].astype(str)\n", " dall_proc[\"gene_dsa\"] = dall_proc[\"gene_raw\"].astype(str) if g_a else np.nan\n", "\n", " dpsi_dsa_col = _find_first(dall.columns, \"delta_psi\",\"dpsi\",\"deltaPSI\",\"IncLevelDifference\")\n", " dall_proc[\"delta_psi_dsa\"] = pd.to_numeric(dall[dpsi_dsa_col], errors=\"coerce\") if dpsi_dsa_col else np.nan\n", "\n", " intr_proc = None\n", " if p_dsr is not None:\n", " intr = _read_mnt(z, p_dsr)\n", " intr.columns = [c.strip() for c in intr.columns]\n", " chr_i = _find_first(intr.columns, \"chrom\",\"chr\",\"chromosome\")\n", " st_i = _find_first(intr.columns, \"start\",\"loc_start\",\"s\")\n", " en_i = _find_first(intr.columns, \"end\",\"loc_end\",\"e\")\n", " sd_i = _find_first(intr.columns, \"strand\")\n", " g_i = _find_first(intr.columns, \"gene_name\",\"gene\",\"Gene\")\n", "\n", " if all([chr_i, st_i, en_i, sd_i, g_i]):\n", " intr_proc = intr[[chr_i, st_i, en_i, sd_i, g_i]].copy()\n", " intr_proc = intr_proc.rename(columns={\n", " chr_i:\"chrom\", st_i:\"start\", en_i:\"end\", sd_i:\"strand\", g_i:\"gene_dsr\"\n", " })\n", " intr_proc[\"chrom\"] = intr_proc[\"chrom\"].astype(str).str.replace(\"chr\",\"\")\n", " gid_col, gstruct_col, gloc_col, gq_col = detect_group_columns(intr)\n", " intr_proc[\"group_id\"] = intr[gid_col].astype(str) if gid_col else np.nan\n", " intr_proc[\"group_structure\"] = intr[gstruct_col].astype(str) if gstruct_col else np.nan\n", " intr_proc[\"group_loc\"] = pd.to_numeric(intr[gloc_col], errors=\"coerce\") if gloc_col else np.nan\n", " intr_proc[\"group_q\"] = pd.to_numeric(intr[gq_col], errors=\"coerce\") if gq_col else np.nan\n", "\n", " dpsi_dsr_col = _find_first(intr.columns, \"delta_psi\",\"dpsi\",\"deltaPSI\",\"IncLevelDifference\")\n", " intr_proc[\"delta_psi_dsr\"] = pd.to_numeric(intr[dpsi_dsr_col], errors=\"coerce\") if dpsi_dsr_col else np.nan\n", "\n", " if intr_proc is not None:\n", " m = dall_proc.merge(intr_proc, on=[\"chrom\",\"start\",\"end\",\"strand\"], how=\"left\", suffixes=(\"\",\"_dsr\"))\n", " else:\n", " m = dall_proc.copy()\n", " for c in [\"gene_dsr\",\"group_id\",\"group_structure\",\"group_loc\",\"group_q\",\"delta_psi_dsr\"]:\n", " m[c] = np.nan\n", "\n", " m[\"gene\"] = m[\"gene_dsr\"].fillna(m[\"gene_dsa\"])\n", " m[\"delta_psi\"] = m[\"delta_psi_dsr\"]\n", " m.loc[m[\"delta_psi\"].isna(), \"delta_psi\"] = m.loc[m[\"delta_psi\"].isna(), \"delta_psi_dsa\"]\n", " m[\"MntJULiP_source\"] = np.where(m[\"group_id\"].notna(), \"DSR\", \"DSA\")\n", " m[\"MntJULiP_q\"] = pd.to_numeric(m[\"MntJULiP_q\"], errors=\"coerce\")\n", " m = m[m[\"MntJULiP_q\"] <= alpha]\n", "\n", " if \"delta_psi\" in m.columns:\n", " m = m[m[\"delta_psi\"].abs() >= min_abs_dpsi]\n", "\n", " if NONCODING_SET is not None:\n", " m = m[m[\"gene\"].astype(str).isin(NONCODING_SET)]\n", "\n", " m = m.dropna(subset=[\"start\",\"end\"]).copy()\n", " m[\"src\"] = m[\"MntJULiP_source\"]\n", " m[\"__key\"] = m[\"chrom\"].astype(str) + \":\" + m[\"start\"].astype(int).astype(str) + \"-\" + m[\"end\"].astype(int).astype(str) + \":\" + m[\"strand\"].astype(str)\n", " return m.drop_duplicates([\"gene\",\"__key\"])\n", "\n", "# ---------- Rank genes by frequency ----------\n", "rmats_all_df = load_all_rmats_events(RMATS_ZIP)\n", "mnt_all_df = load_mnt_sig_introns_union(MNT_ZIP, ALPHA_FDR, MIN_ABS_DPSI)\n", "\n", "master = pd.concat([\n", " rmats_all_df.assign(source=\"rMATS\", __event_key=rmats_all_df[\"__event_key\"]),\n", " mnt_all_df.assign(source=\"MntJULiP\", __event_key=mnt_all_df[\"__key\"])\n", "], ignore_index=True)\n", "master = master.dropna(subset=[\"gene\"]).drop_duplicates([\"gene\",\"__event_key\",\"source\"], keep=\"first\")\n", "\n", "freq = (master.groupby(\"gene\")\n", " .agg(unique_events=(\"__event_key\",\"nunique\"), total_rows=(\"__event_key\",\"size\"))\n", " .reset_index().sort_values([\"unique_events\",\"total_rows\",\"gene\"], ascending=[False,False,True])\n", " .reset_index(drop=True))\n", "freq[\"splicing_rank\"] = np.arange(1, len(freq)+1)\n", "freq.to_csv(OUT_FREQ_RANK, index=False)\n", "\n", "rmats_sig = load_sig_rmats_events(RMATS_ZIP, ALPHA_FDR, MIN_ABS_DPSI)\n", "rmats_sig_counts = rmats_sig.groupby(\"gene\")[\"__event_key\"].nunique().rename(\"rmats_sig_unique\").reset_index() if not rmats_sig.empty else pd.DataFrame(columns=[\"gene\",\"rmats_sig_unique\"])\n", "mnt_sig_counts = mnt_all_df.groupby(\"gene\")[\"__key\"].nunique().rename(\"mntjulip_any_sig_unique\").reset_index() if not mnt_all_df.empty else pd.DataFrame(columns=[\"gene\",\"mntjulip_any_sig_unique\"])\n", "\n", "reg = freq[[\"gene\",\"unique_events\",\"total_rows\",\"splicing_rank\"]].merge(rmats_sig_counts, on=\"gene\", how=\"left\").merge(mnt_sig_counts, on=\"gene\", how=\"left\")\n", "reg[\"rmats_sig_unique\"] = reg[\"rmats_sig_unique\"].fillna(0).astype(int)\n", "reg[\"mntjulip_any_sig_unique\"] = reg[\"mntjulip_any_sig_unique\"].fillna(0).astype(int)\n", "\n", "def _mode(row):\n", " r = row[\"rmats_sig_unique\"] > 0\n", " m = row[\"mntjulip_any_sig_unique\"] > 0\n", " if r and m: return \"both\"\n", " if r and not m: return \"exon_level_only\"\n", " if (not r) and m: return \"intron_only\"\n", " return \"neither\"\n", "\n", "reg[\"regulation_mode\"] = reg.apply(_mode, axis=1)\n", "reg.to_csv(OUT_REGMODES, index=False)\n", "\n", "# ---------- rMATS per-event junction segments ----------\n", "def rmats_junction_segments(evt, row):\n", " segs = []\n", " gi = _gi\n", " if evt == \"SE\":\n", " uEE = gi(row.get(\"upstreamEE\")); exS = gi(row.get(\"exonStart_0base\"))\n", " exE = gi(row.get(\"exonEnd\")); dES = gi(row.get(\"downstreamES\"))\n", " if uEE is not None and exS is not None: segs.append((min(uEE, exS), max(uEE, exS)))\n", " if exE is not None and dES is not None: segs.append((min(exE, dES), max(exE, dES)))\n", " elif evt == \"RI\":\n", " uEE = gi(row.get(\"upstreamEE\")); dES = gi(row.get(\"downstreamES\"))\n", " if uEE is not None and dES is not None: segs.append((min(uEE, dES), max(uEE, dES)))\n", " elif evt in (\"A5SS\",\"A3SS\"):\n", " lS = gi(row.get(\"longExonStart_0base\")); lE = gi(row.get(\"longExonEnd\"))\n", " sS = gi(row.get(\"shortES\")); sE = gi(row.get(\"shortEE\"))\n", " fS = gi(row.get(\"flankingES\")); fE = gi(row.get(\"flankingEE\"))\n", " if fE is not None and lS is not None: segs.append((min(fE,lS), max(fE,lS)))\n", " if lE is not None and fS is not None: segs.append((min(lE,fS), max(lE,fS)))\n", " if fE is not None and sS is not None: segs.append((min(fE,sS), max(fE,sS)))\n", " if sE is not None and fS is not None: segs.append((min(sE,fS), max(sE,fS)))\n", " elif evt == \"MXE\":\n", " uEE = gi(row.get(\"upstreamEE\")); e1S = gi(row.get(\"1stExonStart_0base\")); e1E = gi(row.get(\"1stExonEnd\"))\n", " e2S = gi(row.get(\"2ndExonStart_0base\")); e2E = gi(row.get(\"2ndExonEnd\")); dES = gi(row.get(\"downstreamES\"))\n", " if uEE is not None and e1S is not None: segs.append((min(uEE,e1S), max(uEE,e1S)))\n", " if e1E is not None and dES is not None: segs.append((min(e1E,dES), max(e1E,dES)))\n", " if uEE is not None and e2S is not None: segs.append((min(uEE,e2S), max(uEE,e2S)))\n", " if e2E is not None and dES is not None: segs.append((min(e2E,dES), max(e2E,dES)))\n", " return sorted({(s,e) for (s,e) in segs if s is not None and e is not None and e > s})\n", "\n", "# ---------- GENCODE catalog ----------\n", "def build_gencode_catalog_for_gene(gtf_zip, gene_name):\n", " exons = []\n", " with _open_gtf_any_from_zip(gtf_zip) as f:\n", " for line in f:\n", " if not line or line.startswith(\"#\"): continue\n", " parts = line.rstrip(\"\\n\").split(\"\\t\")\n", " if len(parts) < 9 or parts[2] != \"exon\": continue\n", " d = {m.group(1): m.group(2) for m in re.finditer(r'(\\S+)\\s+\"(.*?)\";', parts[8])}\n", " gname = d.get(\"gene_name\") or d.get(\"gene_id\",\"\")\n", " if gname != gene_name and gene_name not in gname: continue\n", " exons.append({\"chr\": (\"chr\"+parts[0] if not str(parts[0]).startswith(\"chr\") else parts[0]),\n", " \"start\": int(parts[3]), \"end\": int(parts[4]), \"strand\": parts[6], \"transcript_id\": d.get(\"transcript_id\",\"\")})\n", " exons_df = pd.DataFrame(exons)\n", " introns = []\n", " for tx, sub in exons_df.groupby(\"transcript_id\"):\n", " sub = sub.sort_values(\"start\")\n", " for (_, e1), (_, e2) in zip(sub.iloc[:-1].iterrows(), sub.iloc[1:].iterrows()):\n", " if e2[\"start\"] > e1[\"end\"]:\n", " introns.append({\"chr\": e1[\"chr\"], \"start\": e1[\"end\"], \"end\": e2[\"start\"], \"strand\": e1[\"strand\"], \"transcript_id\": tx})\n", " introns_df = pd.DataFrame(introns)\n", " intron_index = {(c,s): sub[[\"start\",\"end\",\"transcript_id\"]].to_records(index=False) for (c,s), sub in introns_df.groupby([\"chr\",\"strand\"])} if len(introns_df) else {}\n", " exon_index = {(c,s): sub[[\"start\",\"end\",\"transcript_id\"]].to_records(index=False) for (c,s), sub in exons_df.groupby([\"chr\",\"strand\"])}\n", " return exons_df, introns_df, exon_index, intron_index\n", "\n", "def gencode_intron_match(ch, st, a, b, intron_index, slop=0):\n", " return [tx for (s,e,tx) in intron_index.get((ch, st), []) if (s - slop) <= a <= (s + slop) and (e - slop) <= b <= (e + slop)]\n", "\n", "# ---------- Anchor helpers ----------\n", "def anchor_positions(evt, row):\n", " if evt == \"RI\":\n", " upES = _gi(row.get(\"upstreamES\"))\n", " return [upES] if upES is not None else []\n", " if evt == \"MXE\":\n", " return [x for x in [_gi(row.get(\"1stExonStart_0base\")), _gi(row.get(\"2ndExonStart_0base\"))] if x is not None]\n", " if evt == \"A5SS\":\n", " s = _gi(row.get(\"shortES\"))\n", " return [s] if s is not None else []\n", " return []\n", "\n", "def min_anchor_dist(group_loc, anchors):\n", " try: gl = int(group_loc)\n", " except: return np.inf\n", " return min(abs(gl - a) for a in anchors) if anchors else np.inf\n", "\n", "def coloc_to_anchors(group_loc, anchors, tol=1):\n", " try: gl = int(group_loc)\n", " except: return False\n", " return any(abs(gl - a) <= tol for a in anchors) if anchors else False\n", "\n", "# ---------- Pipeline Core ----------\n", "def per_gene_pipeline(gene, rmats_all_df, rmats_sig, mnt_all_df):\n", " r_any = rmats_all_df[rmats_all_df[\"gene\"] == gene].copy()\n", " if r_any.empty or rmats_sig.empty: return None\n", " chr_col = next((c for c in [\"chr\",\"chrom\",\"Chromosome\"] if c in r_any.columns), None)\n", " if chr_col is None: return None\n", " r_any[\"chrom\"] = r_any[chr_col].astype(str).map(_normalize_chr)\n", "\n", " r_coords = rmats_all_df.drop_duplicates(\"__event_key\")\n", " r_sig = rmats_sig.merge(r_coords, on=[\"__event_key\",\"gene\",\"__evt_type\"], how=\"left\")\n", " r_sig = r_sig[r_sig[\"gene\"] == gene].copy()\n", " if r_sig.empty: return None\n", "\n", " m = mnt_all_df[mnt_all_df[\"gene\"] == gene].copy()\n", " r_sig[\"chrom\"] = r_sig[chr_col].astype(str).map(_normalize_chr)\n", "\n", " pub_rows = []\n", " for _, row in r_sig.iterrows():\n", " evt = row[\"__evt_type\"]\n", " segs = rmats_junction_segments(evt, row)\n", " if not segs: continue\n", " evctr = float(np.mean([(s+e)/2.0 for (s,e) in segs]))\n", " anchors = anchor_positions(evt, row)\n", "\n", " msub = m[(m[\"chrom\"] == row[\"chrom\"]) & (m[\"strand\"] == str(row[\"strand\"]))].copy()\n", " if msub.empty:\n", " pub_rows.append({\"gene\": gene, \"event_type\": evt, \"chr\": \"chr\"+str(row[\"chrom\"]), \"strand\": row[\"strand\"], \"rMATS_event_key\": row[\"__event_key\"], \"rMATS_dPSI\": row[\"delta_psi\"], \"rMATS_FDR\": row[\"fdr\"], \"MntJULiP_intron_start\": np.nan, \"note\": \"No MntJULiP match\"})\n", " continue\n", "\n", " msub[\"Overlap_bp\"] = msub.apply(lambda r: max([_span_overlap((s,e), (int(r[\"start\"]), int(r[\"end\"]))) for s,e in segs]), axis=1)\n", " msub = msub[msub[\"Overlap_bp\"] > 0].copy()\n", " if msub.empty:\n", " pub_rows.append({\"gene\": gene, \"event_type\": evt, \"chr\": \"chr\"+str(row[\"chrom\"]), \"strand\": row[\"strand\"], \"rMATS_event_key\": row[\"__event_key\"], \"rMATS_dPSI\": row[\"delta_psi\"], \"rMATS_FDR\": row[\"fdr\"], \"MntJULiP_intron_start\": np.nan, \"note\": \"No MntJULiP match\"})\n", " continue\n", "\n", " if anchors and \"group_loc\" in msub.columns:\n", " msub[\"anchor_dist\"] = msub[\"group_loc\"].apply(lambda gl: min_anchor_dist(gl, anchors))\n", " msub[\"CoLoc\"] = msub[\"group_loc\"].apply(lambda gl: coloc_to_anchors(gl, anchors, tol=COLOC_TOL))\n", " else:\n", " msub[\"anchor_dist\"] = np.inf\n", " msub[\"CoLoc\"] = False\n", "\n", " msub = msub.sort_values([\"anchor_dist\", \"Overlap_bp\"], ascending=[True, False])\n", " best = msub.iloc[0]\n", "\n", " pub_rows.append({\n", " \"gene\": gene, \"event_type\": evt, \"chr\": \"chr\"+str(row[\"chrom\"]), \"strand\": row[\"strand\"],\n", " \"rMATS_event_key\": row[\"__event_key\"], \"rMATS_dPSI\": row[\"delta_psi\"], \"rMATS_FDR\": row[\"fdr\"],\n", " \"MntJULiP_intron_start\": int(best[\"start\"]), \"MntJULiP_intron_end\": int(best[\"end\"]),\n", " \"MntJULiP_q\": float(best.get(\"MntJULiP_q\", np.nan)), \"MntJULiP_source\": best.get(\"src\", \"\")\n", " })\n", "\n", " pub = pd.DataFrame(pub_rows).drop_duplicates([\"gene\",\"rMATS_event_key\"])\n", " if pub.empty: return None\n", "\n", " # GENCODE Support\n", " exons_df, introns_df, exon_index, intron_index = build_gencode_catalog_for_gene(GTF_ZIP, gene)\n", " def gcheck(r):\n", " if pd.isna(r[\"MntJULiP_intron_start\"]): return \"No\", \"no_match\"\n", " hit = gencode_intron_match(r[\"chr\"], str(r[\"strand\"]), int(r[\"MntJULiP_intron_start\"]), int(r[\"MntJULiP_intron_end\"]), intron_index, slop=GENCODE_SLOP)\n", " return (\"Yes\", \"exact_intron\") if hit else (\"No\", \"no_match\")\n", " pub[[\"annotated_in_gencode\",\"gencode_match_type\"]] = pub.apply(gcheck, axis=1, result_type=\"expand\")\n", "\n", " pub.to_csv(f\"/content/CRC_{gene}_With_Annotation.csv\", index=False)\n", " return pub\n", "\n", "# -------- Execute --------\n", "selected_genes = reg.head(TOP_PICK_N)[\"gene\"].tolist()\n", "print(f\"[INFO] Automatically selected top {TOP_PICK_N} genes for integration: {selected_genes}\")\n", "\n", "combined = []\n", "for g in selected_genes:\n", " res = per_gene_pipeline(g, rmats_all_df, rmats_sig, mnt_all_df)\n", " if res is not None and not res.empty:\n", " combined.append(res)\n", "\n", "if combined:\n", " pd.concat(combined, ignore_index=True).to_csv(OUT_COMBINED_ANNOT, index=False)\n", " print(f\"\\n✅ SUCCESS: Full integration table saved to {OUT_COMBINED_ANNOT}\")\n", "```\n", "\n" ], "metadata": { "id": "XtSohQLT8LJ6" } }, { "cell_type": "code", "source": [], "metadata": { "id": "R-S3f-H6pflY" }, "execution_count": null, "outputs": [] }, { "cell_type": "markdown", "source": [ "# CODE CELL-8\n", "### Clinical Validation Pipeline for Candidate Non-Coding RNAs\n", "\n", "**Objective:** To structurally validate the top-ranked alternative splicing events against standard genomic annotations and dynamically format the verified coordinates for downstream clinical correlation in the TCGA and GTEx databases.\n", "\n", "**Methodological Steps:**\n", "\n", "1. **Exon Annotation Database Construction:**\n", " To establish a baseline of known human genomic structures, the complete human transcriptome annotation was downloaded. A highly efficient, memory-mapped coordinate database was constructed by parsing this compressed GTF file. Only features explicitly annotated as `exon` were retained.\n", " * **Reference Version:** GENCODE Release 40 (GRCh38.p13).\n", " * **Data Structure:** Exon coordinates (chromosome, strand, start, end) were indexed into a hashed dictionary to enable rapid O(1) mathematical cross-referencing.\n", "\n", "2. **Target Isolation and Significance Filtering:**\n", " The previously generated frequency ranking (`Colorectal_ncRNA_Splicing_Ranking.csv`) was imported, and the top five most frequently spliced non-coding RNA loci were isolated. The comprehensive rMATS skipped-exon dataset (`SE.MATS.JC.txt`) was subsequently loaded and strictly filtered to retain only statistically significant events (False Discovery Rate [FDR] < 0.05) associated with these top target loci.\n", "\n", "3. **Coordinate Transformation and Novelty Cross-Referencing:**\n", " Bioinformatics algorithms utilize differing coordinate systems. To accurately cross-reference the junctions, a mathematical transformation was applied. The 0-based start coordinates produced by rMATS were programmatically converted to the 1-based coordinate system utilized by the GENCODE GTF. Each target exon boundary was then queried against the custom annotation database:\n", " * Events with exact boundary matches were flagged as **\"Annotated (GENCODE)\"**.\n", " * Events with divergent boundaries were flagged as **\"Novel (Unannotated)\"**.\n", "\n", "4. **Clinical Query Formatting and Export:**\n", " For each validated event, the Inclusion Level Difference (ΔPSI) was calculated to quantify the magnitude and direction of the splicing shift between tumor and normal states. The genomic coordinates were concatenated into a standardized web-query format (e.g., `chrX:start-end`) natively compatible with the TCGA SpliceSeq (COAD cohort) and GTEx Portal search architecture. This formatted data was exported as `Validation_Targets_TCGA_GTEx.csv` to facilitate independent clinical validation.\n", "\n", "**Software & Module Dependencies:**\n", "* **Python:** v3.10+\n", "* **Pandas:** v1.5+ / 2.0+ (Utilized for memory-efficient dataframe subsetting)\n", "* **Annotation Database:** GENCODE v40 (GRCh38.p13)\n", "* **Standard Python Libraries:** `os`, `gzip`, `collections.defaultdict`\n", "\n", "```\n", "import os\n", "import gzip\n", "import pandas as pd\n", "from collections import defaultdict\n", "\n", "# --- 1. File Paths ---\n", "rmats_se_file = \"rmats_colorectal_out/content/rmats_out/SE.MATS.JC.txt\"\n", "gtf_file = \"gencode.v40.annotation.gtf.gz\"\n", "ranking_file = \"Colorectal_ncRNA_Splicing_Ranking.csv\"\n", "\n", "# Ensure pre-requisite files exist\n", "if not os.path.exists(rmats_se_file) or not os.path.exists(ranking_file):\n", " raise FileNotFoundError(\"rMATS output or Ranking CSV not found. Please run the previous cell first.\")\n", "\n", "# --- 2. Build Exon Annotation Database from GENCODE v40 ---\n", "print(\"Building exact coordinate database from GENCODE v40...\")\n", "annotated_exons = defaultdict(set)\n", "\n", "with gzip.open(gtf_file, 'rt') as f:\n", " for line in f:\n", " if line.startswith('#'): continue\n", " fields = line.strip().split('\\t')\n", " if fields[2] == 'exon':\n", " chrom = fields[0]\n", " start = int(fields[3]) # GTF is 1-based\n", " end = int(fields[4])\n", " strand = fields[6]\n", " annotated_exons[(chrom, strand)].add((start, end))\n", "\n", "# --- 3. Extract Top Targets ---\n", "print(\"\\nExtracting significant events for top 5 non-coding RNAs...\")\n", "df_ranking = pd.read_csv(ranking_file)\n", "top_targets = df_ranking['Gene'].head(5).tolist() # Targeting top 5 (e.g., JPX, GAS5, SNHG17)\n", "\n", "df_rmats = pd.read_csv(rmats_se_file, sep='\\t', low_memory=False)\n", "fdr_col = next((c for c in df_rmats.columns if c.lower() in ['fdr', 'padj']), None)\n", "gene_col = next((c for c in df_rmats.columns if c.lower() in ['genesymbol', 'gene_name']), None)\n", "\n", "sig_rmats = df_rmats[df_rmats[fdr_col] < 0.05].copy()\n", "\n", "# Filter for our top target genes\n", "pattern = '|'.join(top_targets)\n", "target_events = sig_rmats[sig_rmats[gene_col].astype(str).str.contains(pattern, na=False, case=False)].copy()\n", "\n", "# --- 4. Cross-Reference and Format ---\n", "print(\"Validating annotation status and formatting for TCGA/GTEx...\")\n", "validation_records = []\n", "\n", "for idx, row in target_events.iterrows():\n", " gene = row[gene_col].replace('\"', '')\n", " chrom = row['chr']\n", " strand = row['strand']\n", " event_id = row['ID']\n", " \n", " # rMATS 0-based start, 1-based end. Convert start to 1-based for GTF match\n", " exon_start_1base = int(row['exonStart_0base']) + 1\n", " exon_end = int(row['exonEnd'])\n", " \n", " # Check if this exact exon boundaries exist in GENCODE\n", " is_annotated = (exon_start_1base, exon_end) in annotated_exons[(chrom, strand)]\n", " status = \"Annotated (GENCODE)\" if is_annotated else \"Novel (Unannotated)\"\n", " \n", " # Calculate ΔPSI (Tumor - Normal)\n", " inc_diff = row.get('IncLevelDifference', 0)\n", " \n", " validation_records.append({\n", " 'Gene': gene,\n", " 'Event_ID': event_id,\n", " 'Coordinate_Query': f\"{chrom}:{exon_start_1base}-{exon_end}\",\n", " 'Status': status,\n", " 'FDR': row[fdr_col],\n", " 'Delta_PSI': inc_diff\n", " })\n", "\n", "val_df = pd.DataFrame(validation_records)\n", "val_df.to_csv(\"Validation_Targets_TCGA_GTEx.csv\", index=False)\n", "\n", "print(\"\\n\" + \"=\"*70)\n", "print(\"TOP VALIDATION TARGETS (Extracted for Web Portal Query)\")\n", "print(\"=\"*70)\n", "print(val_df[['Gene', 'Coordinate_Query', 'Status', 'Delta_PSI']].head(15).to_string(index=False))\n", "print(\"\\n--> Full validation table saved to 'Validation_Targets_TCGA_GTEx.csv'\")\n", "# This is formatted as code\n", "```\n", "\n" ], "metadata": { "id": "tpFGj_gOJMwU" } }, { "cell_type": "markdown", "source": [ "# CODE CELL-9\n", "# Supplementary Note: Local Environment Execution and Snaptron API Query Architecture\n", "\n", "### Rationale for Local Execution Environment\n", "In previous analyses of esophageal adenocarcinoma, cloud-based environments such as Google Colab were utilized for computational pipelines. However, for the current colorectal splicing validation, execution was transitioned to a local Anaconda/Jupyter Notebook environment. This transition was necessitated by the strict Web Application Firewall (WAF) protocols implemented by the Johns Hopkins University (JHU) Snaptron server.\n", "\n", "High-volume API queries originating from known cloud datacenter IP ranges (such as those used by Google Colab) are frequently flagged as automated bot traffic and blocked. By executing the query engine locally on a residential broadband network, these cloud-specific WAF restrictions are bypassed, ensuring uninterrupted and complete data retrieval for the validation targets.\n", "\n", "### Software and Environment Specifications\n", "To ensure full reproducibility, the local environment and dependencies used for this analysis are detailed below:\n", "* **Operating System:** Windows 11 (64-bit)\n", "* **Environment:** Anaconda / Miniconda3\n", "* **Python Version:** 3.10.x\n", "* **Core Libraries:** `pandas` (v2.0+), `numpy` (v1.24+)\n", "* **Network Protocol:** Native Windows OS `curl.exe` (invoked via Python `subprocess`)\n", "\n", "### Addressing Query Syntax and the \"Zero-Match\" Controversy\n", "A critical methodological distinction must be made regarding how the Snaptron database indexes genomic coordinates, and what an empty response signifies.\n", "\n", "**1. The `chr` Prefix Artifact:**\n", "Genomic coordinates annotated via GENCODE standardly include the `chr` prefix (e.g., `chr20:38422046-38422241`). However, the Snaptron backend database is indexed strictly using integers for chromosomes. When a query is submitted containing the literal string \"chr\", the database engine executes the search but finds zero matching index keys, returning an empty set. This is a computational syntax artifact, not a reflection of biological absence. To demonstrate this transparency, the script below executes two separate queries for every target: one raw query (including `chr`) and one cleaned query (excluding `chr`).\n", "\n", "**2. Biological Absence vs. Network Failure:**\n", "Standard Python HTTP libraries (e.g., `requests`) frequently encounter `ChunkedEncodingError` exceptions when the Snaptron database drops connections prematurely during heavy database searches. To eliminate this confounding variable, a robust wrapper utilizing the native OS-level `curl` command was engineered. Because `curl` is resilient to premature connection drops, an empty return (0 matches) via this method definitively indicates that the server successfully searched the entire transcriptomic database and found no matching reads. Therefore, the absence of these specific splicing variations in normal GTEx tissue and standard TCGA junction catalogs represents a true biological negative in those cohorts, underscoring their highly specific, tumor-enriched nature in the current study.\n", "\n", "---\n", "\n", "### Execution Code\n", "```python\n", "import pandas as pd\n", "import numpy as np\n", "import subprocess\n", "import time\n", "import os\n", "\n", "print(\"--- INITIALIZING SNAPTRON OS-LEVEL QUERY ENGINE ---\")\n", "\n", "# 1. File Check\n", "targets_file = \"Validation_Targets_TCGA_GTEx.csv\"\n", "if not os.path.exists(targets_file):\n", " raise FileNotFoundError(f\"❌ Cannot find '{targets_file}'. Please place it in {os.getcwd()}\")\n", "\n", "# 2. The Robust 'curl' Wrapper (Bypasses Python HTTP library errors)\n", "def snaptron_curl_query(base_url, region, sqs_filter=None):\n", " url = f\"{base_url}?regions={region}&dtype=s\"\n", " if sqs_filter:\n", " encoded_filter = sqs_filter.replace('\"', '%22').replace(' ', '%20')\n", " url += f\"&sqs={encoded_filter}\"\n", " \n", " cmd = [\n", " \"curl.exe\",\n", " \"-s\", \n", " \"--max-time\", \"60\", \n", " \"-H\", \"User-Agent: Mozilla/5.0 (Windows NT 10.0; Win64; x64) Chrome/120.0.0.0 Safari/537.36\",\n", " url\n", " ]\n", " \n", " for attempt in range(3):\n", " try:\n", " time.sleep(1.5)\n", " result = subprocess.run(cmd, capture_output=True, text=True, check=False)\n", " text = result.stdout.strip()\n", " \n", " if not text:\n", " return {}\n", " \n", " per_sample = {}\n", " for line in text.splitlines():\n", " parts = line.split('\\t')\n", " if len(parts) >= 2:\n", " try:\n", " per_sample[parts[0]] = float(parts[1])\n", " except ValueError:\n", " pass\n", " \n", " if per_sample:\n", " return per_sample\n", " \n", " except Exception as e:\n", " print(f\" [Attempt {attempt+1}] OS command failed. Retrying...\")\n", " time.sleep(3)\n", " \n", " return {}\n", "\n", "def summarize_counts(counts_dict):\n", " if not counts_dict:\n", " return 0, 0.0\n", " vals = list(counts_dict.values())\n", " return len(vals), float(np.mean(vals))\n", "\n", "# 3. Execution\n", "df_targets = pd.read_csv(targets_file)\n", "total_targets = len(df_targets)\n", "print(f\"✅ Loaded {total_targets} targets. Starting dual-queries...\\n\")\n", "\n", "results_list = []\n", "\n", "for index, row in df_targets.iterrows():\n", " raw_region = str(row[\"Coordinate_Query\"]).strip()\n", " \n", " if pd.isna(raw_region) or raw_region == 'nan' or raw_region == '':\n", " results_list.append({\n", " \"Raw_Query\": \"None\", \"Clean_Query\": \"None\",\n", " \"GTEx_Transverse_n_RAW\": 0, \"GTEx_Transverse_n_CLEAN\": 0,\n", " \"GTEx_Sigmoid_n_RAW\": 0, \"GTEx_Sigmoid_n_CLEAN\": 0,\n", " \"TCGA_n_RAW\": 0, \"TCGA_n_CLEAN\": 0\n", " })\n", " continue\n", " \n", " # Generate the cleaned coordinate for biological testing\n", " clean_region = raw_region.replace('chr', '')\n", " \n", " print(f\"Querying [{index + 1}/{total_targets}]: {raw_region} vs {clean_region} ...\")\n", " \n", " # -- A. GTEx Transverse Colon (Testing both formats) --\n", " tr_raw = snaptron_curl_query(\"[http://snaptron.cs.jhu.edu/gtex/snaptron](http://snaptron.cs.jhu.edu/gtex/snaptron)\", raw_region, 'tissue:\"Colon - Transverse\"')\n", " tr_raw_n, _ = summarize_counts(tr_raw)\n", " \n", " tr_clean = snaptron_curl_query(\"[http://snaptron.cs.jhu.edu/gtex/snaptron](http://snaptron.cs.jhu.edu/gtex/snaptron)\", clean_region, 'tissue:\"Colon - Transverse\"')\n", " tr_clean_n, _ = summarize_counts(tr_clean)\n", " \n", " # -- B. GTEx Sigmoid Colon (Testing both formats) --\n", " sig_raw = snaptron_curl_query(\"[http://snaptron.cs.jhu.edu/gtex/snaptron](http://snaptron.cs.jhu.edu/gtex/snaptron)\", raw_region, 'tissue:\"Colon - Sigmoid\"')\n", " sig_raw_n, _ = summarize_counts(sig_raw)\n", " \n", " sig_clean = snaptron_curl_query(\"[http://snaptron.cs.jhu.edu/gtex/snaptron](http://snaptron.cs.jhu.edu/gtex/snaptron)\", clean_region, 'tissue:\"Colon - Sigmoid\"')\n", " sig_clean_n, _ = summarize_counts(sig_clean)\n", " \n", " # -- C. TCGA COAD (Testing both formats) --\n", " tc_raw = snaptron_curl_query(\"[http://snaptron.cs.jhu.edu/tcga/snaptron](http://snaptron.cs.jhu.edu/tcga/snaptron)\", raw_region, 'gdc_cases.project.project_id:\"TCGA-COAD\"')\n", " tc_raw_n, _ = summarize_counts(tc_raw)\n", " \n", " tc_clean = snaptron_curl_query(\"[http://snaptron.cs.jhu.edu/tcga/snaptron](http://snaptron.cs.jhu.edu/tcga/snaptron)\", clean_region, 'gdc_cases.project.project_id:\"TCGA-COAD\"')\n", " tc_clean_n, _ = summarize_counts(tc_clean)\n", " \n", " # Append the side-by-side proof to the results\n", " results_list.append({\n", " \"Raw_Query\": raw_region,\n", " \"Clean_Query\": clean_region,\n", " \"GTEx_Transverse_n_RAW\": int(tr_raw_n),\n", " \"GTEx_Transverse_n_CLEAN\": int(tr_clean_n),\n", " \"GTEx_Sigmoid_n_RAW\": int(sig_raw_n),\n", " \"GTEx_Sigmoid_n_CLEAN\": int(sig_clean_n),\n", " \"TCGA_n_RAW\": int(tc_raw_n),\n", " \"TCGA_n_CLEAN\": int(tc_clean_n)\n", " })\n", "\n", "# Recombine and save\n", "cohort_df = pd.DataFrame(results_list)\n", "out_annotated = pd.concat([df_targets.reset_index(drop=True), cohort_df], axis=1)\n", "\n", "out_file = \"Snaptron_API_Validated_Results_DualQuery.csv\"\n", "out_annotated.to_csv(out_file, index=False)\n", "\n", "print(\"\\n\" + \"=\"*80)\n", "print(\"✅ SNAPTRON API DUAL-VALIDATION COMPLETE\")\n", "print(f\"Results saved locally as: {out_file}\")\n", "print(\"=\"*80)" ], "metadata": { "id": "jVeSJ9wuRcJ_" } }, { "cell_type": "markdown", "source": [ "# CODE CELL-10\n", "# 12. Visualization of Alternative Splicing Architecture\n", "\n", "To highlight the structural divergence of transcriptomic dysregulation, the physical genomic coordinates of the most frequently spliced non-coding RNAs were mapped. By extracting the exact junction bounds from the integration pipeline, mutually exclusive exons (MXEs), alternative splice sites (A5SS/A3SS), and retained introns (RIs) are rendered on a single shifted axis[cite: 2]. This approach cleanly delineates the complex exon-switching events and structural instability defining the tumor microenvironment[cite: 2].\n", "\n", "### Methodological Steps\n", "\n", "**1. Configuration & Color Palettes:**\n", "A standardized color dictionary maps distinct hues to individual MXE events and specific alternative splicing mechanisms (e.g., A5SS, A3SS, RI). This guarantees visual continuity across all figures in the manuscript, matching the aesthetic established in previous EAC analyses[cite: 2]. The output directory is automatically provisioned within the runtime environment.\n", "\n", "**2. Key Parsing Logic:**\n", "The `parse_rmats_key` function deconstructs the unique `rMATS_event_key` string (e.g., `11_-_..._A5SS`) to systematically extract the event type and raw genomic coordinates. Valid numeric strings are separated from missing data flags (`-1`). Complex events like MXE are grouped into discrete tuple pairs to accurately reflect multiple alternative segments, ensuring precise spatial representation.\n", "\n", "**3. Plotting Engine:**\n", "The `plot_gene_architecture` module iterates over all parsed events for a specific gene. It identifies the absolute minimum genomic coordinate across the entire event set, establishing a fixed `anchor` parameter[cite: 2]. All subsequent junction bounds are mathematically subtracted from this anchor to shift the X-axis to zero, enhancing visual clarity. MXE variants are iteratively plotted on distinct upper horizontal lanes, while alternative mechanisms (A5SS, RI) are sequentially distributed on lower lanes[cite: 2]. Dynamic custom legends and annotations are generated programmatically prior to exporting the result as a high-resolution 300 DPI TIFF image.\n", "\n", "**4. Execution Runner:**\n", "The main entry point loads the multi-omic validation dataset (`Final_Validated_CRC_Targets_DualQuery.csv`), validates the presence of gene identifiers, and specifically subsets the dataframe for the top high-priority targets (*GAS5* and *SNHG1*). It then invokes the plotting engine sequentially.\n", "\n", "### Environment & Module Versions\n", "This visualization pipeline relies on the following standard data science packages:\n", "* **Python:** 3.10+\n", "* **pandas:** 1.5.3 or 2.0+ (For dataframe ingestion and subsetting)\n", "* **numpy:** 1.23.5 or 1.24+ (For mathematical transformations and coordinate aggregations)\n", "* **matplotlib:** 3.7.1 (For canvas generation, geometrical rendering, and figure export)\n", "\n", "```\n", "# This is formatted as code\n", "import os\n", "import pandas as pd\n", "import numpy as np\n", "import matplotlib.pyplot as plt\n", "import matplotlib.patches as mpatches\n", "\n", "# ===========================================================\n", "# 1. CONFIGURATION & COLOR PALETTES\n", "# ===========================================================\n", "INPUT_CSV = \"/content/Final_Validated_CRC_Targets_DualQuery.csv\"\n", "OUT_DIR = \"/content/Architectural_Plots\"\n", "os.makedirs(OUT_DIR, exist_ok=True)\n", "\n", "# Maintained visual continuity with established aesthetics\n", "COLORS = {\n", " \"MXE\": [\n", " \"#7b61ff\", \"#ff6b6b\", \"#4ecdc4\", \"#ffa600\",\n", " \"#1e90ff\", \"#ff1493\", \"#32cd32\", \"#8b4513\",\n", " \"#8a2be2\", \"#00ced1\"\n", " ],\n", " \"SE\": \"#a9a9a9\",\n", " \"A5SS_LONG\": \"#3b75ff\",\n", " \"A5SS_SHORT\": \"#38e5d2\",\n", " \"A3SS_LONG\": \"#000080\",\n", " \"A3SS_SHORT\": \"#4169e1\",\n", " \"RI\": \"#ff4a4a\",\n", " \"MntJULiP\": \"#ff8c00\"\n", "}\n", "\n", "# ===========================================================\n", "# 2. KEY PARSING LOGIC\n", "# ===========================================================\n", "def parse_rmats_key(key_string):\n", " \"\"\"\n", " Deconstructs the rMATS_event_key to extract raw genomic coordinates\n", " based on the event type structure.\n", " \"\"\"\n", " parts = str(key_string).split('_')\n", " if len(parts) < 4:\n", " return None, [], None\n", " \n", " evt_type = parts[-1].upper()\n", " chrom = parts[0].replace(\"chr\", \"\")\n", " \n", " # Extract only valid numeric coordinates, ignoring '-1' missing flags\n", " coords = []\n", " for p in parts[2:-1]:\n", " try:\n", " val = int(p)\n", " if val != -1:\n", " coords.append(val)\n", " except ValueError:\n", " pass\n", " \n", " segments = []\n", " \n", " # Identify the specific variant segments to plot based on event type\n", " if evt_type == \"MXE\" and len(coords) >= 4:\n", " segments.append((\"MXE\", coords[0], coords[1]))\n", " segments.append((\"MXE\", coords[2], coords[3]))\n", " elif evt_type == \"SE\" and len(coords) >= 2:\n", " segments.append((\"SE\", coords[0], coords[1]))\n", " elif evt_type in (\"A5SS\", \"A3SS\") and len(coords) >= 6:\n", " segments.append((f\"{evt_type}_LONG\", coords[0], coords[1]))\n", " segments.append((f\"{evt_type}_SHORT\", coords[4], coords[5]))\n", " elif evt_type == \"RI\" and len(coords) >= 2:\n", " segments.append((\"RI\", coords[0], coords[1]))\n", " elif len(coords) >= 2:\n", " segments.append((evt_type, coords[0], coords[1]))\n", " \n", " return evt_type, segments, chrom\n", "\n", "# ===========================================================\n", "# 3. PLOTTING ENGINE\n", "# ===========================================================\n", "def plot_gene_architecture(gene, group_df):\n", " parsed_events = []\n", " all_coords = []\n", " \n", " for _, row in group_df.iterrows():\n", " key = row.get(\"rMATS_event_key\", \"\")\n", " evt_type, segments, chrom = parse_rmats_key(key)\n", " \n", " # Fallback to MntJULiP coordinates if rMATS key parsing fails\n", " if not segments and pd.notna(row.get(\"MntJULiP_intron_start\")):\n", " try:\n", " segments = [(\"MntJULiP\", int(row[\"MntJULiP_intron_start\"]), int(row[\"MntJULiP_intron_end\"]))]\n", " evt_type = \"MntJULiP\"\n", " chrom = str(row.get(\"chr\", \"chr\")).replace(\"chr\", \"\")\n", " except Exception:\n", " pass\n", " \n", " if segments:\n", " parsed_events.append((evt_type, segments, chrom))\n", " for (_, s, e) in segments:\n", " all_coords.extend([s, e])\n", " \n", " if not parsed_events:\n", " print(f\" [!] Skipping {gene}: No valid coordinates parsed.\")\n", " return\n", "\n", " # Calculate absolute minimum coordinate to serve as the zero-anchor\n", " anchor = min(all_coords)\n", " primary_chrom = parsed_events[0][2]\n", " \n", " def shift(x):\n", " return x - anchor\n", " \n", " plt.figure(figsize=(18, 10))\n", " ax = plt.gca()\n", " \n", " y_offset = -1\n", " mxe_counter = 0\n", " legend_elements = {}\n", " \n", " # Sort events so MXEs are plotted first\n", " parsed_events.sort(key=lambda x: 0 if x[0] == \"MXE\" else 1)\n", " \n", " for evt_type, segments, chrom in parsed_events:\n", " if evt_type == \"MXE\":\n", " color = COLORS[\"MXE\"][mxe_counter % len(COLORS[\"MXE\"])]\n", " label = f\"MXE-{mxe_counter + 1}\"\n", " mxe_counter += 1\n", " if label not in legend_elements:\n", " legend_elements[label] = mpatches.Patch(color=color, label=label)\n", " \n", " # Plot MXE parts on the same horizontal lane\n", " for (seg_type, s_raw, e_raw) in segments:\n", " ax.hlines(y_offset, shift(s_raw), shift(e_raw), color=color, linewidth=12)\n", " y_offset -= 1\n", " \n", " else:\n", " for (seg_type, s_raw, e_raw) in segments:\n", " color = COLORS.get(seg_type, \"#000000\")\n", " if seg_type not in legend_elements:\n", " legend_elements[seg_type] = mpatches.Patch(color=color, label=seg_type)\n", " \n", " # Each non-MXE segment type gets its own lane\n", " ax.hlines(y_offset, shift(s_raw), shift(e_raw), color=color, linewidth=12)\n", " y_offset -= 1\n", "\n", " # Formatting and Labels\n", " mxe_title_part = f\"MXE-1 to MXE-{mxe_counter}\" if mxe_counter > 0 else \"Shifted axis\"\n", " plt.title(f\"{gene} alternative splicing architecture\\n(All Events, {mxe_title_part})\", fontsize=20, pad=20)\n", " plt.xlabel(f\"Shifted genomic coordinate (chr{primary_chrom}, {gene} locus)\", fontsize=14)\n", " plt.yticks([], [])\n", " \n", " plt.text(0, y_offset, f\"0 = chr{primary_chrom}:{anchor:,}\", fontsize=14)\n", " plt.axhline(0, color=\"black\", linewidth=1, linestyle=\"--\", alpha=0.5)\n", " \n", " plt.legend(handles=list(legend_elements.values()), title=\"Event legend\", loc=\"lower right\", fontsize=12, frameon=True)\n", " \n", " plt.tight_layout()\n", " out_file = os.path.join(OUT_DIR, f\"{gene}_Architecture_Map.tiff\")\n", " plt.savefig(out_file, dpi=300, format=\"tiff\")\n", " plt.close()\n", " print(f\" [+] Saved {gene} figure: {out_file}\")\n", "\n", "# ===========================================================\n", "# 4. EXECUTION RUNNER\n", "# ===========================================================\n", "if __name__ == \"__main__\":\n", " if not os.path.exists(INPUT_CSV):\n", " print(f\"ERROR: Cannot find {INPUT_CSV}. Ensure the file is uploaded to the Colab environment.\")\n", " else:\n", " df = pd.read_csv(INPUT_CSV)\n", " gene_col = next((c for c in df.columns if c.lower() in [\"gene\", \"genesymbol\", \"gene_name\"]), None)\n", " \n", " if gene_col:\n", " # Filtering specifically for the targets of interest\n", " target_genes = [\"GAS5\", \"SNHG1\"]\n", " for gene_name in target_genes:\n", " group_data = df[df[gene_col] == gene_name]\n", " if not group_data.empty:\n", " print(f\"Processing {gene_name}...\")\n", " plot_gene_architecture(gene_name, group_data)\n", " else:\n", " print(f\" [!] No data found for {gene_name}.\")\n", " \n", " print(f\"\\nPlots generated in {OUT_DIR}/. Use files.download() to export.\")\n", " else:\n", " print(\"ERROR: Could not identify a 'Gene' column in the CSV.\")\n", "```\n", "\n" ], "metadata": { "id": "F_RMX6cHnpTZ" } }, { "cell_type": "markdown", "source": [ "# CODE CELL-11\n", "# 14. Generation of Publication-Ready Group-Averaged Sashimi Plots\n", "\n", "To visually validate the quantitative splicing alterations identified in the integration pipeline, we utilized `rmats2sashimiplot` to render the transcriptomic read coverage and junction spanning reads.\n", "\n", "### Rationale for Group Averaging\n", "Plotting high-throughput RNA-seq data for large cohorts (e.g., all individual normal and tumor replicates) inherently causes severe vertical compression, rendering read coverage peaks and statistical text overlapping and illegible. To produce clean, publication-ready figures that accurately communicate biological divergence, we implemented **Group Averaging**.\n", "\n", "By utilizing a dynamically generated grouping file (`.gf`), the plotter aggregates the individual BAM tracks into two distinct, averaged tracks: **Normal** and **Tumor**. This mitigates individual sample noise, prevents visual overcrowding, and clearly highlights the cohort-level structural shifts in *GAS5* and *SNHG1* exon usage.\n", "\n", "### Computational Environment & Dependencies\n", "To ensure stability and reproducibility, the plotting engine was executed in a strictly controlled Python 3.10 environment. The legacy `rmats2sashimiplot` repository was cloned and programmatically converted to Python 3 syntax using the `2to3` utility.\n", "* **Plotting Engine:** `rmats2sashimiplot` (Xinglab, GitHub)\n", "* **Python:** 3.10.x\n", "* **Core Libraries:** `pysam==0.19.1`, `matplotlib==3.7.1`, `numpy==1.24.3`, `scipy==1.10.1`\n", "\n", "### Methodological Steps\n", "\n", "**1. BAM Preparation & Exclusion:**\n", "BAM files and indices are synced from Google Cloud Storage. Known corrupted or truncated samples (e.g., `SRR8616044`) are explicitly removed to prevent `samtools` indexing failures or plotting crashes. Cohort manifest files (`b1.txt` and `b2.txt`) are generated dynamically.\n", "\n", "**2. Event Extraction & Strict Parsing:**\n", "The rMATS output archive is extracted, and the raw `.MATS.txt` files are parsed. A custom Python script rebuilds the specific structural key (e.g., `11_-_..._A5SS`) to map our stringently validated dual-query targets against the raw rMATS data. Specific `.event.txt` files are generated exclusively for the *GAS5* and *SNHG1* splice junctions of interest.\n", "\n", "**3. Group Info Generation:**\n", "A grouping file (`grouping.gf`) is automatically compiled by calculating the exact number of biological replicates in the Normal and Tumor manifests, instructing the plotter to mathematically average the read densities across the respective cohorts.\n", "\n", "**4. Sashimi Plot Execution:**\n", "The plotting engine iterates through the generated event files. By utilizing the `--group-info` parameter, it renders and exports highly legible, two-track PDF visualizations to a dedicated output directory (`/content/Colorectal_Sashimi_Plots_Averaged`), avoiding redundant computations by checking for pre-existing plots.\n", "```\n", "# This is formatted as code\n", "%%bash\n", "set -e\n", "\n", "# ==========================================\n", "# 1. FILE PATHS & SETUP\n", "# ==========================================\n", "TABLE_PATH=\"C:\\Users\\jm_ma\\Downloads\\CRC_TopGenes_Combined_With_Annotation.csv\"\n", "RMATS_DIR=\"/content/rmats_output_full\"\n", "OUT_DIR_EVENTS=\"/content/rMATS_Events_To_Plot_V3\"\n", "\n", "# --> NEW DIRECTORY FOR AVERAGED PLOTS <--\n", "OUT_DIR_PLOTS=\"/content/Colorectal_Sashimi_Plots_Averaged\"\n", "\n", "echo \"=== 1. PREPARING BAMS ===\"\n", "# Skip the heavy rsync if the BAM folders already have contents\n", "if [ -z \"$(ls -A /content/bams/tumor 2>/dev/null)\" ]; then\n", " echo \"--> Folders empty. Syncing BAMS...\"\n", " mkdir -p /content/bams/tumor /content/bams/normal\n", " gsutil -m -o GSUtil:check_hashes=never rsync -r \"gs://colodata/BAMSNORMAL\" /content/bams/normal/ || true\n", " gsutil -m -o GSUtil:check_hashes=never rsync -r \"gs://colodata\" /content/bams/tumor/ || true\n", " find /content/bams/tumor -type f ! -name '*.bam' ! -name '*.bai' -delete || true\n", " \n", " echo \"--> Excluding corrupted sample SRR8616044...\"\n", " rm -f /content/bams/tumor/SRR8616044*.bam\n", " rm -f /content/bams/tumor/SRR8616044*.bai\n", "else\n", " echo \"--> BAMs already present on disk. Skipping cloud sync.\"\n", "fi\n", "\n", "echo \"--> Verifying BAM indexes (generating any missing .bai)...\"\n", "for bam in /content/bams/normal/*.bam /content/bams/tumor/*.bam; do\n", " if [ -f \"$bam\" ] && [ ! -f \"${bam}.bai\" ]; then\n", " samtools index \"$bam\"\n", " fi\n", "done\n", "\n", "echo \"--> Generating b1.txt (Normal) and b2.txt (Tumor)...\"\n", "ls -d /content/bams/normal/*.bam | paste -sd, - > /content/b1.txt\n", "ls -d /content/bams/tumor/*.bam | paste -sd, - > /content/b2.txt\n", "\n", "# ==========================================\n", "# 2. EXTRACT RMATS RESULTS\n", "# ==========================================\n", "echo \"=== 2. EXTRACTING RMATS RESULTS ===\"\n", "if [ -d \"$RMATS_DIR\" ] && [ \"$(ls -A $RMATS_DIR 2>/dev/null)\" ]; then\n", " echo \"--> rMATS results already extracted. Skipping unzip.\"\n", "else\n", " mkdir -p \"$RMATS_DIR\"\n", " unzip -q -o /content/rmats_results_colorectal.zip -d \"$RMATS_DIR\"\n", " echo \"--> Extraction complete.\"\n", "fi\n", "\n", "# ==========================================\n", "# 3. GENERATE TARGET EVENT FILES\n", "# ==========================================\n", "echo \"=== 3. STRICT PARSING OF VALIDATION TARGETS ===\"\n", "if [ -d \"$OUT_DIR_EVENTS\" ] && [ \"$(ls -A $OUT_DIR_EVENTS 2>/dev/null)\" ]; then\n", " echo \"--> Event files already generated. Skipping parsing.\"\n", "else\n", " mkdir -p \"$OUT_DIR_EVENTS\"\n", " cat << 'EOF' > /content/generate_events_v3.py\n", "import os\n", "from pathlib import Path\n", "import pandas as pd\n", "\n", "TABLE_PATH = r\"C:\\Users\\jm_ma\\Downloads\\CRC_TopGenes_Combined_With_Annotation.csv\"\n", "RMATS_DIR = Path(\"/content/rmats_output_full\")\n", "OUT_DIR = Path(\"/content/rMATS_Events_To_Plot_V3\")\n", "\n", "try:\n", " t3 = pd.read_csv(TABLE_PATH)\n", "except FileNotFoundError:\n", " t3 = pd.read_csv(\"/content/CRC_TopGenes_Combined_With_Annotation.csv\")\n", "\n", "t3 = t3.dropna(subset=['rMATS_event_key', 'gene'])\n", "target_keys = set(t3['rMATS_event_key'].astype(str).str.strip())\n", "\n", "req_coords = {\n", " \"SE\": [\"exonStart_0base\", \"exonEnd\", \"upstreamES\", \"upstreamEE\", \"downstreamES\", \"downstreamEE\"],\n", " \"RI\": [\"exonStart_0base\", \"exonEnd\", \"upstreamES\", \"upstreamEE\", \"downstreamES\", \"downstreamEE\"],\n", " \"A5SS\": [\"longExonStart_0base\", \"longExonEnd\", \"flankingES\", \"flankingEE\", \"shortES\", \"shortEE\"],\n", " \"A3SS\": [\"longExonStart_0base\", \"longExonEnd\", \"flankingES\", \"flankingEE\", \"shortES\", \"shortEE\"],\n", " \"MXE\": [\"1stExonStart_0base\", \"1stExonEnd\", \"2ndExonStart_0base\", \"2ndExonEnd\", \"upstreamES\", \"upstreamEE\", \"downstreamES\", \"downstreamEE\"],\n", "}\n", "\n", "def build_event_key(evt, row, cols):\n", " chr_col = next((c for c in cols if c.lower() in ['chr', 'chrom', 'chromosome']), 'chr')\n", " strand_col = next((c for c in cols if c.lower() == 'strand'), 'strand')\n", " chrom = str(row.get(chr_col, \"\")).replace(\"chr\", \"\").strip()\n", " strand = str(row.get(strand_col, \"\")).strip()\n", "\n", " def gi(c):\n", " if c not in cols: return \"\"\n", " v = row[c]\n", " try: return str(int(float(str(v).strip())))\n", " except: return \"\"\n", "\n", " parts = [chrom, strand]\n", " for c in req_coords.get(evt, []):\n", " parts.append(gi(c))\n", " parts.append(evt)\n", " return \"_\".join(parts)\n", "\n", "for evt in [\"SE\", \"RI\", \"A3SS\", \"A5SS\", \"MXE\"]:\n", " for variant in [\"JCEC\", \"JC\"]:\n", " target_filename = f\"{evt}.MATS.{variant}.txt\"\n", " matches = list(RMATS_DIR.rglob(target_filename))\n", " \n", " if not matches: continue\n", " \n", " path = matches[0]\n", " df = pd.read_csv(path, sep=\"\\t\", dtype=str)\n", " cols = [c.strip() for c in df.columns]\n", " df.columns = cols\n", " \n", " gcol = next((c for c in cols if c.lower() in ['genesymbol', 'gene_name', 'gene']), None)\n", " if 'ID' not in cols or not gcol: continue\n", " \n", " keys = []\n", " for _, r in df.iterrows():\n", " keys.append(build_event_key(evt, r, cols))\n", " df[\"__event_key\"] = keys\n", " \n", " df_target = df[df[\"__event_key\"].isin(target_keys)].copy()\n", " \n", " for _, mrow in df_target.iterrows():\n", " evkey = mrow[\"__event_key\"]\n", " gene = str(mrow[gcol]).replace('\"', '').strip()\n", " \n", " for c in req_coords.get(evt, []):\n", " if c in cols:\n", " val = str(mrow.get(c, \"\")).strip()\n", " if val == \"\" or val.lower() == \"nan\":\n", " mrow[c] = \"-1\"\n", " \n", " export_row = mrow.drop(\"__event_key\")\n", " export_cols = [c for c in cols if c != \"__event_key\"]\n", " \n", " out_path = OUT_DIR / f\"{gene}_{evkey}_{variant}_{evt}.event.txt\"\n", " if not out_path.exists():\n", " with open(out_path, \"w\") as f:\n", " f.write(\"\\t\".join(export_cols) + \"\\n\")\n", " f.write(\"\\t\".join(str(export_row.get(c, \"\")) for c in export_cols) + \"\\n\")\n", "EOF\n", " python3 /content/generate_events_v3.py\n", "fi\n", "\n", "# ==========================================\n", "# 4. ENVIRONMENT SETUP FOR RMATS2SASHIMIPLOT\n", "# ==========================================\n", "echo \"=== 4. INSTALLING PYTHON 3.10 & PLOTTING DEPENDENCIES ===\"\n", "if ! command -v python3.10 &> /dev/null; then\n", " sudo apt-get update -y > /dev/null 2>&1\n", " sudo apt-get install -y python3.10 python3.10-dev python3.10-distutils > /dev/null 2>&1\n", " curl -sS https://bootstrap.pypa.io/get-pip.py | python3.10 > /dev/null 2>&1\n", " python3.10 -m pip install pysam==0.19.1 matplotlib==3.7.1 numpy==1.24.3 scipy==1.10.1 > /dev/null 2>&1\n", "else\n", " echo \"--> Python 3.10 already installed. Skipping.\"\n", "fi\n", "\n", "# ==========================================\n", "# 5. CLONE SASHIMI PLOTTER\n", "# ==========================================\n", "echo \"=== 5. SETTING UP SASHIMI PLOTTER ===\"\n", "if [ ! -d \"/content/rmats2sashimiplot_fixed/src\" ]; then\n", " cd /content\n", " rm -rf rmats2sashimiplot_fixed\n", " git clone https://github.com/Xinglab/rmats2sashimiplot.git rmats2sashimiplot_fixed > /dev/null 2>&1\n", " cd rmats2sashimiplot_fixed\n", "\n", " cat << 'EOF' > convert.py\n", "import sys\n", "from lib2to3.main import main\n", "sys.argv = ['2to3', '-w', '-n', '--no-diffs', 'src']\n", "sys.exit(main(\"lib2to3.fixes\"))\n", "EOF\n", " python3.10 convert.py > /dev/null 2>&1\n", "else\n", " echo \"--> Sashimi plotter already cloned and configured. Skipping.\"\n", "fi\n", "\n", "# ==========================================\n", "# 6. GENERATE GROUPING FILE FOR AVERAGED PLOTS\n", "# ==========================================\n", "echo \"=== 6. CREATING GROUPING INFO FOR AVERAGING ===\"\n", "cat << 'EOF' > /content/make_grouping.py\n", "import os\n", "# Read the number of bams dynamically so we don't hardcode index ranges\n", "with open('/content/b1.txt', 'r') as f:\n", " b1_count = len([x for x in f.read().strip().split(',') if x])\n", "with open('/content/b2.txt', 'r') as f:\n", " b2_count = len([x for x in f.read().strip().split(',') if x])\n", "\n", "# Write the .gf format required by rmats2sashimiplot\n", "with open('/content/grouping.gf', 'w') as f:\n", " f.write(f\"Normal: 1-{b1_count}\\n\")\n", " f.write(f\"Tumor: {b1_count + 1}-{b1_count + b2_count}\\n\")\n", " \n", "print(f\"--> Created grouping.gf with Normal (1-{b1_count}) and Tumor ({b1_count + 1}-{b1_count + b2_count})\")\n", "EOF\n", "python3 /content/make_grouping.py\n", "\n", "# ==========================================\n", "# 7. RUN THE SASHIMI PLOT GENERATOR\n", "# ==========================================\n", "echo \"=== 7. RUNNING AVERAGED SASHIMI PLOTS ===\"\n", "mkdir -p \"$OUT_DIR_PLOTS\"\n", "export PYTHONPATH=$PYTHONPATH:/content/rmats2sashimiplot_fixed/src\n", "\n", "shopt -s nullglob\n", "event_files=(\"$OUT_DIR_EVENTS\"/*.event.txt)\n", "\n", "if [ ${#event_files[@]} -eq 0 ]; then\n", " echo \"--> ERROR: No event files were found. Check your target CSV format.\"\n", "else\n", " for event_file in \"${event_files[@]}\"; do\n", " filename=$(basename -- \"$event_file\")\n", " evt_type=\"${filename##*_}\"\n", " evt_type=\"${evt_type%.event.txt}\"\n", "\n", " # Skip plot generation if the PDF already exists in the NEW averaged directory\n", " if [ -f \"$OUT_DIR_PLOTS/${filename%.event.txt}/Sashimi_plot/sashimi_plot.pdf\" ]; then\n", " echo \"--> Averaged plot for $filename already exists. Skipping.\"\n", " continue\n", " fi\n", "\n", " echo \"--> Plotting Averaged $filename (Type: $evt_type)...\"\n", " \n", " # Notice we inject the --group-info parameter here and keep --l1/--l2 as a failsafe\n", " python3.10 /content/rmats2sashimiplot_fixed/src/rmats2sashimiplot/rmats2sashimiplot.py \\\n", " -o \"$OUT_DIR_PLOTS/${filename%.event.txt}\" \\\n", " --l1 Normal \\\n", " --l2 Tumor \\\n", " --event-type \"$evt_type\" \\\n", " -e \"$event_file\" \\\n", " --b1 /content/b1.txt \\\n", " --b2 /content/b2.txt \\\n", " --group-info /content/grouping.gf \\\n", " --keep-event-chr-prefix \\\n", " --min-counts 0 > /dev/null 2>&1 || echo \" [!] Warning: Plotting failed for $filename. Skipping.\"\n", " done\n", "fi\n", "\n", "echo \"=== ALL DONE: AVERAGED RESULTS SAVED ===\"\n", "echo \"--> You can find your clean PDF plots in $OUT_DIR_PLOTS\"\n", "```\n", "\n" ], "metadata": { "id": "e82puHlNFp5k" } }, { "cell_type": "markdown", "source": [ "# CODE CELL-12\n", "### Explanation of the Figure Creation Process\n", "\n", "This script focuses entirely on data visualization, processing the input CSV and producing two primary types of outputs for your target genes (GAS5 and SNHG1), resulting in exactly four figures.\n", "\n", "**Setup**\n", "The script removes the `python-docx` dependencies and only installs `poppler-utils` (system level) and `pdf2image` (Python level) to handle the PDF conversions required for the grid images.\n", "\n", "**Architecture Plotting (`plot_architecture`)**\n", "* **Data Parsing:** The `parse_rmats_key` function cracks open the raw rMATS event string in the CSV (e.g., extracting event types like MXE, SE, A3SS and bounding coordinates).\n", "* **Mapping:** The plotting function finds the earliest coordinate (the \"anchor\") to normalize the axis starting at 0, making the vast genomic numbers readable on a graph.\n", "* **Drawing:** It iterates downward (`y_offset -= 1`) drawing horizontal line segments using `matplotlib`. MXE events iterate through a custom color palette, while standard events pull from a predefined dictionary (`COLORS`).\n", "* **Output:** The plot is finalized with a legend and exported as both a `.pdf` and `.png` (to avoid Colab's common TIFF rendering errors) into the `Publication_Architecture_Plots` directory. It repeats this for both GAS5 and SNHG1 (yielding 2 architecture figures).\n", "\n", "**Sashimi Grid Tiling (`create_sashimi_grid`)**\n", "* **File Matching:** It scans the `SASHIMI_DIR` to match the rMATS event keys from your CSV to existing individual `.pdf` plots.\n", "* **Grid Calculation:** Using `math.ceil()`, it dynamically calculates how many rows and columns are needed based on the number of matches (defaulting to 3 columns, or 4 if there are more than 6 plots).\n", "* **Image Overlay & Formatting:** It converts the raw PDFs to images (`pdf2image`). To maintain a clean layout, it uses `patches.Rectangle` to essentially paint over the messy original headers with a white box. It then overlays a clean, standardized title block (Gene, Event Type, Coordinate) and an A-Z subplot label in the top left.\n", "* **Output:** The assembled multi-panel array is exported as a single, publication-ready PDF into the `Publication_Sashimi_Grids` directory for both genes (yielding the remaining 2 figures).\n", "```\n", "# This is formatted as code\n", "import os\n", "import sys\n", "import subprocess\n", "import glob\n", "import math\n", "import string\n", "import pandas as pd\n", "import numpy as np\n", "import matplotlib.pyplot as plt\n", "import matplotlib.patches as mpatches\n", "import matplotlib.patches as patches\n", "\n", "# ==========================================================\n", "# 1. BULLETPROOF INSTALLATION & KERNEL REFRESH\n", "# ==========================================================\n", "print(\"Installing system dependencies...\")\n", "subprocess.run(\"apt-get update -qq > /dev/null && apt-get install -y poppler-utils zip > /dev/null\", shell=True)\n", "\n", "try:\n", " from pdf2image import convert_from_path\n", "except ImportError:\n", " print(\"Installing python packages and refreshing kernel paths...\")\n", " subprocess.check_call([sys.executable, \"-m\", \"pip\", \"install\", \"-q\", \"pdf2image\"])\n", " import site\n", " site.main() \n", " from pdf2image import convert_from_path\n", "\n", "from google.colab import files\n", "\n", "# ==========================================\n", "# 2. CONFIGURATION & DIRECTORIES\n", "# ==========================================\n", "INPUT_CSV = \"/content/Final_Validated_CRC_Targets_DualQuery.csv\"\n", "SASHIMI_DIR = \"/content/Colorectal_Sashimi_Plots_Averaged\"\n", "\n", "OUT_ARCH_DIR = \"/content/Publication_Architecture_Plots\"\n", "OUT_GRID_DIR = \"/content/Publication_Sashimi_Grids\"\n", "os.makedirs(OUT_ARCH_DIR, exist_ok=True)\n", "os.makedirs(OUT_GRID_DIR, exist_ok=True)\n", "\n", "TARGET_GENES = [\"GAS5\", \"SNHG1\"]\n", "\n", "COLORS = {\n", " \"MXE\": [\"#7b61ff\", \"#ff6b6b\", \"#4ecdc4\", \"#ffa600\", \"#1e90ff\", \"#ff1493\", \"#32cd32\", \"#8b4513\", \"#8a2be2\", \"#00ced1\", \"#ff6347\", \"#4682b4\"],\n", " \"SE\": \"#a9a9a9\",\n", " \"A5SS_LONG\": \"#3b75ff\",\n", " \"A5SS_SHORT\": \"#38e5d2\",\n", " \"A3SS_LONG\": \"#000080\",\n", " \"A3SS_SHORT\": \"#4169e1\",\n", " \"RI\": \"#ff4a4a\",\n", " \"MntJULiP\": \"#ff8c00\"\n", "}\n", "\n", "# ==========================================\n", "# 3. ARCHITECTURE PLOTTER (Using PDF/PNG to fix fetch error)\n", "# ==========================================\n", "def parse_rmats_key(key_string):\n", " parts = str(key_string).split('_')\n", " if len(parts) < 4: return None, [], None\n", " evt_type = parts[-1].upper()\n", " chrom = parts[0].replace(\"chr\", \"\")\n", " \n", " coords = []\n", " for p in parts[2:-1]:\n", " try:\n", " val = int(p)\n", " if val != -1: coords.append(val)\n", " except ValueError: pass\n", " \n", " segments = []\n", " if evt_type == \"MXE\" and len(coords) >= 4:\n", " segments.append((\"MXE\", coords[0], coords[1]))\n", " segments.append((\"MXE\", coords[2], coords[3]))\n", " elif evt_type == \"SE\" and len(coords) >= 2:\n", " segments.append((\"SE\", coords[0], coords[1]))\n", " elif evt_type in (\"A5SS\", \"A3SS\") and len(coords) >= 6:\n", " segments.append((f\"{evt_type}_LONG\", coords[0], coords[1]))\n", " segments.append((f\"{evt_type}_SHORT\", coords[4], coords[5]))\n", " elif evt_type == \"RI\" and len(coords) >= 2:\n", " segments.append((\"RI\", coords[0], coords[1]))\n", " elif len(coords) >= 2:\n", " segments.append((evt_type, coords[0], coords[1]))\n", " \n", " return evt_type, segments, chrom\n", "\n", "def plot_architecture(gene, group_df):\n", " parsed_events = []\n", " all_coords = []\n", " \n", " for _, row in group_df.iterrows():\n", " key = row.get(\"rMATS_event_key\", \"\")\n", " evt_type, segments, chrom = parse_rmats_key(key)\n", " if segments:\n", " parsed_events.append((evt_type, segments, chrom))\n", " for (_, s, e) in segments:\n", " all_coords.extend([s, e])\n", " \n", " if not parsed_events: return\n", " anchor = min(all_coords)\n", " primary_chrom = parsed_events[0][2]\n", " \n", " def shift(x): return x - anchor\n", " \n", " plt.figure(figsize=(18, 10))\n", " ax = plt.gca()\n", " \n", " y_offset = -1\n", " mxe_counter = 0\n", " legend_elements = {}\n", " parsed_events.sort(key=lambda x: 0 if x[0] == \"MXE\" else 1)\n", " \n", " for evt_type, segments, chrom in parsed_events:\n", " if evt_type == \"MXE\":\n", " color = COLORS[\"MXE\"][mxe_counter % len(COLORS[\"MXE\"])]\n", " label = f\"MXE-{mxe_counter + 1}\"\n", " mxe_counter += 1\n", " if label not in legend_elements:\n", " legend_elements[label] = mpatches.Patch(color=color, label=label)\n", " for (seg_type, s_raw, e_raw) in segments:\n", " ax.hlines(y_offset, shift(s_raw), shift(e_raw), color=color, linewidth=12)\n", " y_offset -= 1\n", " else:\n", " for (seg_type, s_raw, e_raw) in segments:\n", " color = COLORS.get(seg_type, \"#000000\")\n", " if seg_type not in legend_elements:\n", " legend_elements[seg_type] = mpatches.Patch(color=color, label=seg_type)\n", " ax.hlines(y_offset, shift(s_raw), shift(e_raw), color=color, linewidth=12)\n", " y_offset -= 1\n", "\n", " mxe_title = f\"MXE-1 to MXE-{mxe_counter}\" if mxe_counter > 0 else \"Shifted axis\"\n", " plt.title(f\"{gene} Alternative Splicing Architecture\\n(All Events, {mxe_title})\", fontsize=20, pad=20)\n", " plt.xlabel(f\"Shifted genomic coordinate (chr{primary_chrom}, {gene} locus)\", fontsize=14)\n", " plt.yticks([], [])\n", " plt.text(0, y_offset, f\"0 = chr{primary_chrom}:{anchor:,}\", fontsize=14)\n", " plt.axhline(0, color=\"black\", linewidth=1, linestyle=\"--\", alpha=0.5)\n", " plt.legend(handles=list(legend_elements.values()), title=\"Event legend\", loc=\"lower right\", fontsize=12, frameon=True)\n", " \n", " plt.tight_layout()\n", " \n", " # Save as PDF and PNG to avoid the Colab TIFF fetch error\n", " out_pdf = os.path.join(OUT_ARCH_DIR, f\"{gene}_Architecture_Map.pdf\")\n", " out_png = os.path.join(OUT_ARCH_DIR, f\"{gene}_Architecture_Map.png\")\n", " plt.savefig(out_pdf, dpi=300, format=\"pdf\", bbox_inches='tight')\n", " plt.savefig(out_png, dpi=300, format=\"png\", bbox_inches='tight')\n", " plt.close()\n", " print(f\"✅ Created Architecture Plot: {out_pdf}\")\n", "\n", "# ==========================================================\n", "# 4. ROBUST SASHIMI GRID GENERATOR\n", "# ==========================================================\n", "def create_sashimi_grid(gene, group_df):\n", " pdfs_to_plot = []\n", " \n", " for i, row in group_df.reset_index(drop=True).iterrows():\n", " evkey = str(row.get(\"rMATS_event_key\", \"\")).strip()\n", " evt_type = str(row.get(\"event_type\", \"Event\")).upper()\n", " coords_clean = str(row.get(\"Coordinate_Query_CLEAN\", \"Unknown Region\"))\n", " \n", " matches = []\n", " for root, dirs, sysfiles in os.walk(SASHIMI_DIR):\n", " if evkey in root:\n", " for f in sysfiles:\n", " if f.endswith('.pdf'):\n", " matches.append(os.path.join(root, f))\n", " \n", " if matches:\n", " pdf_path = matches[0]\n", " title = f\"{gene} {evt_type} {i+1}\\nRegion: {coords_clean}\"\n", " pdfs_to_plot.append((pdf_path, title))\n", "\n", " n_plots = len(pdfs_to_plot)\n", " if n_plots == 0:\n", " print(f\"⚠️ No Sashimi PDFs found for {gene} in {SASHIMI_DIR}\")\n", " return\n", "\n", " print(f\"Generating Sashimi Grid for {gene} ({n_plots} panels)...\")\n", "\n", " cols = 4 if n_plots > 6 else 3\n", " rows = math.ceil(n_plots / cols)\n", " \n", " fig, axes = plt.subplots(rows, cols, figsize=(8 * cols, 7 * rows))\n", " if n_plots > 1:\n", " axes = axes.flatten()\n", " else:\n", " axes = [axes]\n", " \n", " HIDE_HEADER_HEIGHT = 220\n", "\n", " for i, ax in enumerate(axes):\n", " if i < n_plots:\n", " pdf_path, title_text = pdfs_to_plot[i]\n", " try:\n", " images = convert_from_path(pdf_path, dpi=300)\n", " img = images[0]\n", " width, height = img.size\n", " \n", " ax.imshow(img)\n", " rect = patches.Rectangle((0, 0), width, HIDE_HEADER_HEIGHT,\n", " linewidth=0, edgecolor='white', facecolor='white')\n", " ax.add_patch(rect)\n", " \n", " ax.text(width/2, HIDE_HEADER_HEIGHT/2, title_text,\n", " ha='center', va='center', fontsize=18, fontweight='bold', color='black')\n", " \n", " label = string.ascii_uppercase[i] if i < 26 else str(i)\n", " ax.text(-0.05, 1.05, label, transform=ax.transAxes,\n", " fontsize=32, fontweight='bold', va='top', ha='right')\n", " ax.axis('off')\n", " except Exception as e:\n", " ax.axis('off')\n", " else:\n", " ax.axis('off')\n", "\n", " plt.tight_layout()\n", " out_file = os.path.join(OUT_GRID_DIR, f\"{gene}_Sashimi_Grid_Publication.pdf\")\n", " plt.savefig(out_file, dpi=300, bbox_inches='tight')\n", " plt.close()\n", " print(f\"✅ Created Sashimi Grid: {out_file}\")\n", "\n", "# ==========================================\n", "# 5. EXECUTION RUNNER\n", "# ==========================================\n", "if __name__ == \"__main__\":\n", " if not os.path.exists(INPUT_CSV):\n", " print(f\"ERROR: Cannot find {INPUT_CSV}\")\n", " else:\n", " df = pd.read_csv(INPUT_CSV)\n", " gene_col = next((c for c in df.columns if c.lower() in [\"gene\", \"genesymbol\", \"gene_name\"]), None)\n", " \n", " if gene_col:\n", " for gene_name in TARGET_GENES:\n", " group_data = df[df[gene_col] == gene_name]\n", " if not group_data.empty:\n", " # Generates 1 architecture plot and 1 sashimi grid per gene\n", " plot_architecture(gene_name, group_data)\n", " create_sashimi_grid(gene_name, group_data)\n", " else:\n", " print(f\"⚠️ No data found in CSV for {gene_name}\")\n", " \n", " print(\"\\n🎉 ALL TASKS COMPLETED. Zipping figures and initiating download...\")\n", " \n", " # Zip all publication plots\n", " subprocess.run(\"zip -r /content/Publication_Figures.zip /content/Publication_Architecture_Plots /content/Publication_Sashimi_Grids > /dev/null\", shell=True)\n", " \n", " # Trigger downloads directly to the browser\n", " files.download('/content/Publication_Figures.zip')\n", " \n", " else:\n", " print(\"ERROR: Invalid CSV format, 'gene' column missing.\")\n", "```\n", "\n" ], "metadata": { "id": "bSJwxThFHAUM" } } ] }