CODE CELL 1-METADATA¶
The genomes were downloaded along with their phenotype data with the GEO accession number. Though initially 22 control samples of normal esophgeal cell genome and 22 test samples of genomes of esophgeal adenocarcinoma cell lines were taken, due to poor quality issues and increased percentage of multimapped reads finally only 15 each were selected . The GEO accession numbers and the metadata of the 15 selected in the control and test arm has been mentioned in the table below | Run | Assay Type | AvgSpotLen | Bases | BioProject | BioSample | Bytes | cell_line | cell_type | Center Name | Consent | DATASTORE filetype | DATASTORE provider | DATASTORE region | Experiment | GEO_Accession (exp) | Instrument | LibraryLayout | LibrarySelection | LibrarySource | Organism | Platform | ReleaseDate | create_date | version | Sample Name | source_name | SRA Study | |:------------|:-------------|-------------:|------------:|:-------------|:-------------|-----------:|:-----------------------------|:---------------------------|:--------------|:----------|:---------------------|:---------------------|:-------------------------------------|:-------------|:----------------------|:----------------------|:----------------|:-------------------|:----------------|:-------------|:-----------|:---------------------|:---------------------|----------:|:--------------|:------------------------------|:------------| | SRR12076665 | RNA-Seq | 300 | 6813593700 | PRJNA641420 | SAMN15352566 | 2019498862 | OACP4 C (RRID: CVCL_1843) | Oesophageal Adenocarcinoma | GEO | public | fastq,run.zq,sra | gs,ncbi,s3 | gs.us-east1,ncbi.public,s3.us-east-1 | SRX8604046 | GSM4634039 | Illumina NovaSeq 6000 | PAIRED | cDNA | TRANSCRIPTOMIC | Homo sapiens | ILLUMINA | 2022-06-02T00:00:00Z | 2020-06-23T17:13:00Z | 1 | GSM4634039 | Gastric cardia adenocarcinoma | SRP268506 | | SRR12076666 | RNA-Seq | 300 | 21084481200 | PRJNA641420 | SAMN15352566 | 6552777357 | OACP4 C (RRID: CVCL_1843) | Oesophageal Adenocarcinoma | GEO | public | fastq,run.zq,sra | gs,ncbi,s3 | gs.us-east1,ncbi.public,s3.us-east-1 | SRX8604046 | GSM4634039 | Illumina NovaSeq 6000 | PAIRED | cDNA | TRANSCRIPTOMIC | Homo sapiens | ILLUMINA | 2022-06-02T00:00:00Z | 2020-06-23T17:31:00Z | 1 | GSM4634039 | Gastric cardia adenocarcinoma | SRP268506 | | SRR12076667 | RNA-Seq | 300 | 23223490800 | PRJNA641420 | SAMN15352566 | 7028831354 | OACP4 C (RRID: CVCL_1843) | Oesophageal Adenocarcinoma | GEO | public | fastq,run.zq,sra | gs,ncbi,s3 | gs.us-east1,ncbi.public,s3.us-east-1 | SRX8604046 | GSM4634039 | Illumina NovaSeq 6000 | PAIRED | cDNA | TRANSCRIPTOMIC | Homo sapiens | ILLUMINA | 2022-06-02T00:00:00Z | 2020-06-23T17:35:00Z | 1 | GSM4634039 | Gastric cardia adenocarcinoma | SRP268506 | | SRR12076668 | RNA-Seq | 300 | 1991410500 | PRJNA641420 | SAMN15352565 | 705923129 | OE19 (RRID: CVCL_1622) | Oesophageal Adenocarcinoma | GEO | public | fastq,run.zq,sra | gs,ncbi,s3 | gs.us-east1,ncbi.public,s3.us-east-1 | SRX8604047 | GSM4634040 | Illumina NovaSeq 6000 | PAIRED | cDNA | TRANSCRIPTOMIC | Homo sapiens | ILLUMINA | 2022-06-02T00:00:00Z | 2020-06-23T17:04:00Z | 1 | GSM4634040 | Gastric cardia adenocarcinoma | SRP268506 | | SRR12076669 | RNA-Seq | 300 | 7935932700 | PRJNA641420 | SAMN15352565 | 2521697106 | OE19 (RRID: CVCL_1622) | Oesophageal Adenocarcinoma | GEO | public | fastq,run.zq,sra | gs,ncbi,s3 | gs.us-east1,ncbi.public,s3.us-east-1 | SRX8604047 | GSM4634040 | Illumina NovaSeq 6000 | PAIRED | cDNA | TRANSCRIPTOMIC | Homo sapiens | ILLUMINA | 2022-06-02T00:00:00Z | 2020-06-23T17:24:00Z | 1 | GSM4634040 | Gastric cardia adenocarcinoma | SRP268506 | | SRR12076670 | RNA-Seq | 300 | 13922438700 | PRJNA641420 | SAMN15352565 | 4308927619 | OE19 (RRID: CVCL_1622) | Oesophageal Adenocarcinoma | GEO | public | fastq,run.zq,sra | gs,ncbi,s3 | gs.us-east1,ncbi.public,s3.us-east-1 | SRX8604047 | GSM4634040 | Illumina NovaSeq 6000 | PAIRED | cDNA | TRANSCRIPTOMIC | Homo sapiens | ILLUMINA | 2022-06-02T00:00:00Z | 2020-06-23T17:23:00Z | 1 | GSM4634040 | Gastric cardia adenocarcinoma | SRP268506 | | SRR12076671 | RNA-Seq | 300 | 14874900900 | PRJNA641420 | SAMN15352564 | 4434421815 | FLO-1 (RRID: CVCL_2045) | Oesophageal Adenocarcinoma | GEO | public | fastq,run.zq,sra | gs,ncbi,s3 | gs.us-east1,ncbi.public,s3.us-east-1 | SRX8604048 | GSM4634041 | Illumina NovaSeq 6000 | PAIRED | cDNA | TRANSCRIPTOMIC | Homo sapiens | ILLUMINA | 2022-06-02T00:00:00Z | 2020-06-23T17:28:00Z | 1 | GSM4634041 | Gastric cardia adenocarcinoma | SRP268506 | | SRR12076673 | RNA-Seq | 300 | 12507345900 | PRJNA641420 | SAMN15352564 | 3856103738 | FLO-1 (RRID: CVCL_2045) | Oesophageal Adenocarcinoma | GEO | public | fastq,run.zq,sra | gs,ncbi,s3 | gs.us-east1,ncbi.public,s3.us-east-1 | SRX8604048 | GSM4634041 | Illumina NovaSeq 6000 | PAIRED | cDNA | TRANSCRIPTOMIC | Homo sapiens | ILLUMINA | 2022-06-02T00:00:00Z | 2020-06-23T17:19:00Z | 1 | GSM4634041 | Gastric cardia adenocarcinoma | SRP268506 | | SRR12076674 | RNA-Seq | 300 | 5925931200 | PRJNA641420 | SAMN15352563 | 1768626380 | SK-GT-4 (RRID: CVCL_2195) | Oesophageal Adenocarcinoma | GEO | public | fastq,run.zq,sra | gs,ncbi,s3 | gs.us-east1,ncbi.public,s3.us-east-1 | SRX8604049 | GSM4634042 | Illumina NovaSeq 6000 | PAIRED | cDNA | TRANSCRIPTOMIC | Homo sapiens | ILLUMINA | 2022-06-02T00:00:00Z | 2020-06-23T17:15:00Z | 1 | GSM4634042 | Gastric cardia adenocarcinoma | SRP268506 | | SRR12076675 | RNA-Seq | 300 | 11800185900 | PRJNA641420 | SAMN15352563 | 3711645365 | SK-GT-4 (RRID: CVCL_2195) | Oesophageal Adenocarcinoma | GEO | public | fastq,run.zq,sra | gs,ncbi,s3 | gs.us-east1,ncbi.public,s3.us-east-1 | SRX8604049 | GSM4634042 | Illumina NovaSeq 6000 | PAIRED | cDNA | TRANSCRIPTOMIC | Homo sapiens | ILLUMINA | 2022-06-02T00:00:00Z | 2020-06-23T17:24:00Z | 1 | GSM4634042 | Gastric cardia adenocarcinoma | SRP268506 | | SRR12076677 | RNA-Seq | 300 | 1248840600 | PRJNA641420 | SAMN15352562 | 388361482 | ESO-26 (RRID: CVCL_2035) | Oesophageal Adenocarcinoma | GEO | public | fastq,run.zq,sra | gs,ncbi,s3 | gs.us-east1,ncbi.public,s3.us-east-1 | SRX8604050 | GSM4634043 | Illumina NovaSeq 6000 | PAIRED | cDNA | TRANSCRIPTOMIC | Homo sapiens | ILLUMINA | 2022-06-02T00:00:00Z | 2020-06-23T17:03:00Z | 1 | GSM4634043 | Gastric cardia adenocarcinoma | SRP268506 | | SRR12076681 | RNA-Seq | 300 | 5930880300 | PRJNA641420 | SAMN15352561 | 1894524755 | OE33 (RRID: CVCL_0471) | Oesophageal Adenocarcinoma | GEO | public | fastq,run.zq,sra | gs,ncbi,s3 | gs.us-east1,ncbi.public,s3.us-east-1 | SRX8604051 | GSM4634044 | Illumina NovaSeq 6000 | PAIRED | cDNA | TRANSCRIPTOMIC | Homo sapiens | ILLUMINA | 2022-06-02T00:00:00Z | 2020-06-23T17:12:00Z | 1 | GSM4634044 | Gastric cardia adenocarcinoma | SRP268506 | | SRR12076685 | RNA-Seq | 300 | 16033146000 | PRJNA641420 | SAMN15352560 | 4961667986 | ESO-51 (RRID: CVCL_2036) | Oesophageal Adenocarcinoma | GEO | public | fastq,run.zq,sra | gs,ncbi,s3 | gs.us-east1,ncbi.public,s3.us-east-1 | SRX8604052 | GSM4634045 | Illumina NovaSeq 6000 | PAIRED | cDNA | TRANSCRIPTOMIC | Homo sapiens | ILLUMINA | 2022-06-02T00:00:00Z | 2020-06-23T17:33:00Z | 1 | GSM4634045 | Gastric cardia adenocarcinoma | SRP268506 | | SRR12076687 | RNA-Seq | 300 | 12086426100 | PRJNA641420 | SAMN15352559 | 3873887481 | OAC,5.1 C (RRID: CVCL_1842) | Oesophageal Adenocarcinoma | GEO | public | fastq,run.zq,sra | gs,ncbi,s3 | gs.us-east1,ncbi.public,s3.us-east-1 | SRX8604053 | GSM4634046 | Illumina NovaSeq 6000 | PAIRED | cDNA | TRANSCRIPTOMIC | Homo sapiens | ILLUMINA | 2022-06-02T00:00:00Z | 2020-06-23T17:19:00Z | 1 | GSM4634046 | Gastric cardia adenocarcinoma | SRP268506 | | SRR12076691 | RNA-Seq | 300 | 10687040400 | PRJNA641420 | SAMN15352558 | 3322884050 | JH-EsoAd1 (RRID: CVCL_8098) | Oesophageal Adenocarcinoma | GEO | public | fastq,run.zq,sra | gs,ncbi,s3 | gs.us-east1,ncbi.public,s3.us-east-1 | SRX8604054 | GSM4634047 | Illumina NovaSeq 6000 | PAIRED | cDNA | TRANSCRIPTOMIC | Homo sapiens | ILLUMINA | 2022-06-02T00:00:00Z | 2020-06-23T17:21:00Z | 1 | GSM4634047 | Gastric cardia adenocarcinoma | SRP268506 | THe table showing the list of controls selected are as follows | Run | Assay Type | AvgSpotLen | Bases | BioProject | BioSample | Bytes | Center Name | Consent | DATASTORE filetype | DATASTORE provider | DATASTORE region | Experiment | GEO_Accession (exp) | Instrument | LibraryLayout | LibrarySelection | LibrarySource | Organism | Platform | ReleaseDate | create_date | version | Sample Name | sample_type | source_name | SRA Study | tissue_type | |:------------|:-------------|-------------:|-----------:|:-------------|:-------------|-----------:|:--------------|:----------|:---------------------|:---------------------|:-------------------------------------|:-------------|:----------------------|:--------------------|:----------------|:-------------------|:----------------|:-------------|:-----------|:---------------------|:---------------------|----------:|:--------------|:--------------|:----------------------------|:------------|:--------------| | SRR8931991 | RNA-Seq | 202 | 9156592330 | PRJNA533799 | SAMN11467208 | 5212344969 | GEO | public | fastq,run.zq,sra | gs,s3,ncbi | gs.us-east1,s3.us-east-1,ncbi.public | SRX5712732 | GSM3731538 | Illumina HiSeq 2000 | PAIRED | cDNA | TRANSCRIPTOMIC | Homo sapiens | ILLUMINA | 2019-10-29T00:00:00Z | 2019-04-19T14:40:00Z | 1 | GSM3731538 | normal | ESCC adjacent normal tissue | SRP193095 | ESCC | | SRR8931993 | RNA-Seq | 202 | 7944461636 | PRJNA533799 | SAMN11467206 | 4620123535 | GEO | public | run.zq,fastq,sra | s3,gs,ncbi | ncbi.public,s3.us-east-1,gs.us-east1 | SRX5712734 | GSM3731540 | Illumina HiSeq 2000 | PAIRED | cDNA | TRANSCRIPTOMIC | Homo sapiens | ILLUMINA | 2019-10-29T00:00:00Z | 2019-04-19T14:39:00Z | 1 | GSM3731540 | normal | ESCC adjacent normal tissue | SRP193095 | ESCC | | SRR8931995 | RNA-Seq | 202 | 8374573368 | PRJNA533799 | SAMN11467204 | 4281649348 | GEO | public | fastq,sra,run.zq | s3,gs,ncbi | gs.us-east1,s3.us-east-1,ncbi.public | SRX5712736 | GSM3731542 | Illumina HiSeq 2000 | PAIRED | cDNA | TRANSCRIPTOMIC | Homo sapiens | ILLUMINA | 2019-10-29T00:00:00Z | 2019-04-19T14:44:00Z | 1 | GSM3731542 | normal | ESCC adjacent normal tissue | SRP193095 | ESCC | | SRR8931997 | RNA-Seq | 202 | 7780259876 | PRJNA533799 | SAMN11467202 | 5207013296 | GEO | public | run.zq,fastq,sra | gs,ncbi,s3 | gs.us-east1,ncbi.public,s3.us-east-1 | SRX5712738 | GSM3731544 | Illumina HiSeq 2000 | PAIRED | cDNA | TRANSCRIPTOMIC | Homo sapiens | ILLUMINA | 2019-10-29T00:00:00Z | 2019-04-19T14:38:00Z | 1 | GSM3731544 | normal | ESCC adjacent normal tissue | SRP193095 | ESCC | | SRR8931999 | RNA-Seq | 202 | 7547482146 | PRJNA533799 | SAMN11467200 | 5076859591 | GEO | public | fastq,run.zq,sra | gs,s3,ncbi | s3.us-east-1,ncbi.public,gs.us-east1 | SRX5712740 | GSM3731546 | Illumina HiSeq 2000 | PAIRED | cDNA | TRANSCRIPTOMIC | Homo sapiens | ILLUMINA | 2019-10-29T00:00:00Z | 2019-04-19T14:38:00Z | 1 | GSM3731546 | normal | ESCC adjacent normal tissue | SRP193095 | ESCC | | SRR8932001 | RNA-Seq | 202 | 7898539158 | PRJNA533799 | SAMN11467198 | 4468647504 | GEO | public | fastq,sra,run.zq | ncbi,s3,gs | s3.us-east-1,gs.us-east1,ncbi.public | SRX5712742 | GSM3731548 | Illumina HiSeq 2000 | PAIRED | cDNA | TRANSCRIPTOMIC | Homo sapiens | ILLUMINA | 2019-10-29T00:00:00Z | 2019-04-19T14:41:00Z | 1 | GSM3731548 | normal | ESCC adjacent normal tissue | SRP193095 | ESCC | | SRR8932011 | RNA-Seq | 202 | 8929573216 | PRJNA533799 | SAMN11467222 | 5038274084 | GEO | public | run.zq,fastq,sra | gs,ncbi,s3 | ncbi.public,s3.us-east-1,gs.us-east1 | SRX5712752 | GSM3731558 | Illumina HiSeq 2000 | PAIRED | cDNA | TRANSCRIPTOMIC | Homo sapiens | ILLUMINA | 2019-10-29T00:00:00Z | 2019-04-19T14:42:00Z | 1 | GSM3731558 | normal | ESCC adjacent normal tissue | SRP193095 | ESCC | | SRR10173245 | RNA-Seq | 202 | 8245201458 | PRJNA533799 | SAMN12828211 | 4142319993 | GEO | public | sra,run.zq,fastq | gs,s3,ncbi | ncbi.public,gs.us-east1,s3.us-east-1 | SRX6894900 | GSM4094341 | Illumina HiSeq 2000 | PAIRED | cDNA | TRANSCRIPTOMIC | Homo sapiens | ILLUMINA | 2019-10-29T00:00:00Z | 2019-09-24T12:16:00Z | 1 | GSM4094341 | normal | ESCC adjacent normal tissue | SRP193095 | ESCC | | SRR10173247 | RNA-Seq | 202 | 6886854074 | PRJNA533799 | SAMN12828241 | 3502400021 | GEO | public | sra,fastq,run.zq | ncbi,s3,gs | s3.us-east-1,ncbi.public,gs.us-east1 | SRX6894902 | GSM4094343 | Illumina HiSeq 2000 | PAIRED | cDNA | TRANSCRIPTOMIC | Homo sapiens | ILLUMINA | 2019-10-29T00:00:00Z | 2019-09-24T11:57:00Z | 1 | GSM4094343 | normal | ESCC adjacent normal tissue | SRP193095 | ESCC | | SRR10173249 | RNA-Seq | 202 | 6984377452 | PRJNA533799 | SAMN12828236 | 3563133111 | GEO | public | sra,fastq,run.zq | gs,s3,ncbi | ncbi.public,s3.us-east-1,gs.us-east1 | SRX6894904 | GSM4094345 | Illumina HiSeq 2000 | PAIRED | cDNA | TRANSCRIPTOMIC | Homo sapiens | ILLUMINA | 2019-10-29T00:00:00Z | 2019-09-24T12:03:00Z | 1 | GSM4094345 | normal | ESCC adjacent normal tissue | SRP193095 | ESCC | | SRR10173251 | RNA-Seq | 202 | 6626240138 | PRJNA533799 | SAMN12828234 | 3401621901 | GEO | public | fastq,run.zq,sra | ncbi,gs,s3 | ncbi.public,s3.us-east-1,gs.us-east1 | SRX6894906 | GSM4094347 | Illumina HiSeq 2000 | PAIRED | cDNA | TRANSCRIPTOMIC | Homo sapiens | ILLUMINA | 2019-10-29T00:00:00Z | 2019-09-24T11:56:00Z | 1 | GSM4094347 | normal | ESCC adjacent normal tissue | SRP193095 | ESCC | | SRR10173255 | RNA-Seq | 202 | 6740870088 | PRJNA533799 | SAMN12828226 | 3772474902 | GEO | public | sra,run.zq,fastq | gs,s3,ncbi | gs.us-east1,ncbi.public,s3.us-east-1 | SRX6894910 | GSM4094351 | Illumina HiSeq 2000 | PAIRED | cDNA | TRANSCRIPTOMIC | Homo sapiens | ILLUMINA | 2019-10-29T00:00:00Z | 2019-09-24T11:56:00Z | 1 | GSM4094351 | normal | ESCC adjacent normal tissue | SRP193095 | ESCC | | SRR10173257 | RNA-Seq | 202 | 6542993918 | PRJNA533799 | SAMN12828224 | 3623421871 | GEO | public | fastq,sra,run.zq | gs,ncbi,s3 | ncbi.public,s3.us-east-1,gs.us-east1 | SRX6894912 | GSM4094353 | Illumina HiSeq 2000 | PAIRED | cDNA | TRANSCRIPTOMIC | Homo sapiens | ILLUMINA | 2019-10-29T00:00:00Z | 2019-09-24T11:56:00Z | 1 | GSM4094353 | normal | ESCC adjacent normal tissue | SRP193095 | ESCC | | SRR10173259 | RNA-Seq | 202 | 7030409818 | PRJNA533799 | SAMN12828222 | 3787094622 | GEO | public | fastq,run.zq,sra | gs,ncbi,s3 | s3.us-east-1,ncbi.public,gs.us-east1 | SRX6894914 | GSM4094355 | Illumina HiSeq 2000 | PAIRED | cDNA | TRANSCRIPTOMIC | Homo sapiens | ILLUMINA | 2019-10-29T00:00:00Z | 2019-09-24T12:05:00Z | 1 | GSM4094355 | normal | ESCC adjacent normal tissue | SRP193095 | ESCC |
CODE CELL 2-METADATA¶
Public RNA-seq data for esophageal adenocarcinoma cell lines were obtained from GEO accession GSE153096, submitted by Shorthouse et al. at the MRC Cancer Unit, University of Cambridge. GEO accession: GSE153096. RNA-seq of Oesophageal Adenocarcinoma Cell Lines. Available at: https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE153096
At the time of writing, no formal publication was linked to this dataset
Cell lines included OE19, OE33, ESO26, FLO‑1, SK-GT‑4, OACP4C, and OACM5.1C, which are commonly accepted models of human EAC. Although some metadata fields (e.g., source_name) referenced “gastric cardia adenocarcinoma,” Cellosaurus (RRID: CVCL_1843) and the literature confirm that these lines originate from esophageal tumors, particularly those arising near the gastroesophageal junction.
Public rna seq data for normal esophageal mucosa was obtained from GEO accession GSE130078, submitted by Nam J et al at the Department of Life Science Seongdong-Gu Hangdang-dong Seoul.
CODE CELL 3-COLAB SET UP¶
The analysis was done in xterm in google colab with proplus subscription so that there is enough computing power and memory. The code to load xterm that was used is as follows
# This is formatted as code
!pip install colab-xterm
%load_ext colabxterm
%xterm
google drive with 2tb memory subscription was mounted on colab with the following code in order to store the results of each analysis
# This is formatted as code
from google.colab import drive
drive.mount('/content/drive')
The pipeline for the analysis is mentioned below
Step 1¶
Downloading fastq with prefetch, fasterqdump from sra toolkit, doing a fastqc to check for quality and then run fastp to clean the files and save it as cleaned_fastq. sratoolkit.3.2.1 was downloaded in drive and the fastq files were saved in the folder esodata with eso subfolder for tumor samples and healthy subfolder for controls. fastqc and fastp were downloaded in colab to do the qualtiy analysis and cleaning prior to further analysis
# This is formatted as code
cd /content/drive/MyDrive/sratools/sratoolkit.3.2.1-ubuntu64
export PATH=$PATH:$(pwd)/bin
chmod +x $(pwd)/bin/*
cd /content/drive/MyDrive/esodata/esodata1/eso
#!/bin/bash
# Replace this list with your own SRR IDs
SRR_LIST=("SRR16927859" "SRR16927860" "SRR16927861" "SRR16927856" "SRR16927857" "SRR16927858" "SRR16927853" "SRR16927854" "SRR16927855" "SRR16927850"
) # <-- add your SRRs here
for SRR in "${SRR_LIST[@]}"
do
echo "Processing $SRR..."
# Create directory for this SRR and go into it
mkdir -p "$SRR"
cd "$SRR"
# Download SRA file using prefetch (stored in default SRA path)
prefetch "$SRR"
# Dump FASTQ to current directory
fasterq-dump --split-files --threads 4 --outdir . "$SRR"
# Run FastQC and output to subdir
mkdir -p qc_reports
fastqc *.fastq -o qc_reports
echo "FastQC summary for $SRR:"
grep -H '>>' qc_reports/*_fastqc/summary.txt
# Run fastp on 1 or 2 files depending on single/paired
R1="${SRR}_1.fastq"
R2="${SRR}_2.fastq"
if [ -f "$R2" ]; then
# Paired-end
fastp -i "$R1" -I "$R2" -o "${SRR}_1.cleaned.fastq" -O "${SRR}_2.cleaned.fastq" -h "${SRR}_fastp.html" -j "${SRR}_fastp.json"
else
# Single-end
fastp -i "$R1" -o "${SRR}.cleaned.fastq" -h "${SRR}_fastp.html" -j "${SRR}_fastp.json"
fi
echo "Finished processing $SRR."
cd ..
done
CODE CELL 4 STAR ALIGNMENT¶
STAR version 2.7.11b was donwloaded in drive.
The human genome annotation file (GENCODE release 40)and the primary genome assembly of Homo sapiens (GRCh38.p13) was downloaded from the GENCODE FTP server using the following command:
# This is formatted as code
!wget ftp://ftp.ebi.ac.uk/pub/databases/gencode/Gencode_human/release_40/GRCh38.primary_assembly.genome.fa.gz
!wget ftp://ftp.ebi.ac.uk/pub/databases/gencode/Gencode_human/release_40/gencode.v40.annotation.gtf.gz
This was used for STAR genome index generation using the following command after gunzipping the annotation and fasta files
# This is formatted as code
STAR --runThreadN 6 \
--runMode genomeGenerate \
--genomeDir /content/drive/MyDrive/esodata/esodata1/STAR/index \
--genomeFastaFiles /content/drive/MyDrive/esodata/esodata1/STAR/GRCh38.primary_assembly.genome.fa \
--sjdbGTFfile /content/drive/MyDrive/esodata/esodata1/STAR/gencode.v40.annotation.gtf \
--sjdbOverhang 50
Then we aligned the control and test samples with the following command files
# This is formatted as code
#!/bin/bash -l
cd /content/drive/MyDrive/M1/data/STAR
export PATH="/content/drive/MyDrive/M1/data/STAR/source:$PATH"
chmod +x /content/drive/MyDrive/M1/data/STAR/source/*
cd /content
BASEDIR=/content/drive/MyDrive/esodata/esodata1
DATADIR=${BASEDIR}/healthy
WORKDIR=${BASEDIR}/STAR
STAR2IDX=${BASEDIR}/STAR/index
for i in SRR8931987 SRR8931989 SRR8931991 SRR8931993 SRR8931995 SRR8931997 SRR8931999 SRR8932001 SRR8932003 SRR8932005 SRR8932007 SRR8932009 SRR8932011 SRR10173247
do
time \
STAR \
--limitBAMsortRAM 1200000000 \
--runThreadN 6 \
--genomeDir $STAR2IDX/ \
--readFilesIn \
"${DATADIR}/${i}/${i}_1.cleaned.fastq" \
"${DATADIR}/${i}/${i}_2.cleaned.fastq" \
--outSAMtype BAM SortedByCoordinate \
--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 "${WORKDIR}/${i}_" \
&> "/content.star.log"
done
wait
exit;
the command file for the star alignment of the test samples is as follows
# This is formatted as code
#!/bin/bash -l
set -euo pipefail
cd /content/drive/MyDrive/M1/data/STAR
export PATH="/content/drive/MyDrive/M1/data/STAR/source:$PATH"
chmod +x /content/drive/MyDrive/M1/data/STAR/source/*
cd /content
BASEDIR=/content/drive/MyDrive/esodata/esodata1
DATADIR=${BASEDIR}/eso
WORKDIR=${BASEDIR}/STAR
STAR2IDX=${BASEDIR}/STAR/index
# Use DISK for STAR temp to avoid /dev/shm truncation
DISK_TMP_BASE=/content/star_tmp
mkdir -p "$DISK_TMP_BASE"
for i in SRR12076665 SRR12076666 SRR12076667 SRR12076668 SRR12076669 SRR12076670 SRR12076671 SRR12076673 SRR12076674 SRR12076675 SRR12076681 SRR12076685 SRR12076687 SRR12076691 ;
do
# Use ~60% of currently available RAM for sort (safe on Colab)
LM=$(awk '/MemAvailable/ {printf "%.0f",$2*1024*0.60}' /proc/meminfo)
TMPDIR="${DISK_TMP_BASE}/${i}_STARtmp"
# STAR requires per-sample tmp dir to NOT exist before start
[ -d "$TMPDIR" ] && rm -rf "$TMPDIR"
echo "[$(date)] STAR ${i} --limitBAMsortRAM=${LM} --outTmpDir=${TMPDIR}"
# Clean temp even if interrupted
trap 'rm -rf "$TMPDIR"' EXIT
time STAR \
--limitBAMsortRAM "$LM" \
--outBAMsortingBinsN 300 \
--outBAMcompression 10 \
--runThreadN 6 \
--genomeDir "$STAR2IDX/" \
--readFilesIn \
"${DATADIR}/${i}/${i}_1.cleaned.fastq" \
"${DATADIR}/${i}/${i}_2.cleaned.fastq" \
--outSAMtype BAM SortedByCoordinate \
--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 "${WORKDIR}/${i}_" \
--outTmpDir "$TMPDIR" \
&>> "${WORKDIR}/${i}_STAR.log"
# Index the sorted BAM
samtools index -@ 8 "${WORKDIR}/${i}_Aligned.sortedByCoord.out.bam"
# Cleanup temp
rm -rf "$TMPDIR"
trap - EXIT
done
wait
echo "All STAR runs complete."
CODE CELL 5 - CUFFDIFF CODE¶
to run cuffdiff, mounting the bamfiles from gdrive onto colab was not successful as it was too slow , so a cloud bucket was created in google cloud and the gdrive genetic data were transferred to the cloudbucket with the following code
# This is formatted as code
from google.colab import auth
auth.authenticate_user()
!gsutil -m cp -r "/content/drive/MyDrive/esodata" "gs://esodatacloud/"
Also, i had to move to a customisable vm with a 300 gb disk space to hold the bam files and therefore changed to colab enterprise with a e2-standard-32 runtime. For the cuffdiff analysis, 12 controls and 13 test samples were selected excluding those which were found to have errors after running STAR. Thereafter i ran cuffdiff with the following commands breaking it into individual chrmomosomes because of the sheer bulk of the analysis and later clubbing it into one and exporting the results back into drive using rsync. The code is shown below
# This is formatted as code
#!/usr/bin/env bash
set -euo pipefail
export LC_ALL=C LANG=C
umask 0022
# =========[ 0) CONFIGURE YOUR RUN HERE ]=======================================
# Project / location (only needed if you use gcloud/gcsfuse steps)
PROJECT_ID="YOUR_PROJECT_ID"
REGION="us-central1"
ZONE="us-central1-a"
# Canonical bucket & prefixes
GCS_BUCKET="esodatacloud"
GCS_ROOT="gs://${GCS_BUCKET}/esodata/esodata1"
# Input BAMs on GCS (two groups: control & test). Keep them as GCS prefixes:
CTRL_PREFIX="${GCS_ROOT}/bam/control" # e.g. gs://.../bam/control/*.bam
TEST_PREFIX="${GCS_ROOT}/bam/test" # e.g. gs://.../bam/test/*.bam
# Reference files (GTF & genome FASTA) hosted on GCS
GTF_GCS="${GCS_ROOT}/ref/annotation.gtf"
FA_GCS="${GCS_ROOT}/ref/genome.fa" # optional; cuffdiff doesn’t need fasta
# Labels & library type for cuffdiff
LABELS="c,t"
LIBTYPE="fr-unstranded" # change if needed (e.g., fr-firststrand / fr-secondstrand)
# Per-chromosome scheduling (HDD-friendly defaults)
THREADS_PER_CHR="${THREADS_PER_CHR:-3}" # threads per cuffdiff process
CAP_ACTIVE="${CAP_ACTIVE:-4}" # max concurrent active chromosomes
EXCLUDE_CHR="chr14" # example from your earlier runs
# Run name (timestamped) and local working root
RUN_TAG="$(date -u +run_%Y%m%d_%H%M%S)"
WORK_ROOT="/content/cuffdiff_out"
RUN_DIR="${WORK_ROOT}/${RUN_TAG}"
LOG_DIR="${RUN_DIR}/logs"
CHUNKS="${RUN_DIR}/chunks"
RERUN="${RUN_DIR}/chunks_rerun"
REF_DIR="${RUN_DIR}/ref"
GTF_CHUNKS="${REF_DIR}/gtf_chunks"
# Where to publish results on GCS
GCS_RESULTS_BASE="${GCS_ROOT}/cuffdiffresultsfast"
GCS_RUN_DIR="${GCS_RESULTS_BASE}/${RUN_TAG}" # per-chr outputs
GCS_NEOFINAL="${GCS_RESULTS_BASE}/neofinal/${RUN_TAG}" # stitched classic outputs
# Optional: Google Drive publish via rclone (set to 1 to enable)
PUBLISH_TO_DRIVE="${PUBLISH_TO_DRIVE:-0}"
DRIVE_DST_BASE="mydrive:esodata/esodata1/cuffdiffresultsfast"
# Service account JSON for GCS (if you use rclone for Drive publishing)
SA_JSON="/content/esodata-3691b52fd16f.json" # change if different
# ==============================================================================
echo "[cfg] RUN_TAG=${RUN_TAG}"
echo "[cfg] THREADS_PER_CHR=${THREADS_PER_CHR} CAP_ACTIVE=${CAP_ACTIVE}"
echo "[cfg] GCS_ROOT=${GCS_ROOT}"
mkdir -p "${RUN_DIR}" "${LOG_DIR}" "${CHUNKS}" "${RERUN}" "${REF_DIR}" "${GTF_CHUNKS}"
# =========[ 1) (OPTIONAL) CREATE VM & ATTACH DISK (documentation only) ]=======
cat <<'DOC'
# If you want to reproduce the VM creation via gcloud (run on your workstation or Cloud Shell):
gcloud config set project YOUR_PROJECT_ID
gcloud compute instances create cd-e2-32 \
--zone=us-central1-a \
--machine-type=e2-standard-32 \
--create-disk=size=300GB,type=pd-standard,auto-delete=yes \
--scopes=cloud-platform \
--image-family=debian-12 --image-project=debian-cloud
# SSH and then continue with the steps below on that VM (or run in Colab).
DOC
# =========[ 2) STAGE INPUTS: SYMLINK BAMs from GCS to local short paths ]======
# We’ll use gsutil to list GCS BAMs and symlink them under /dev/shm/cd/{ctrl,test}/
# (fast short paths help some tools and keep logs concise)
echo "[stage] Listing BAMs from GCS..."
mkdir -p /dev/shm/cd/ctrl /dev/shm/cd/test
rm -f /dev/shm/cd/ctrl/* /dev/shm/cd/test/* || true
# Requires gcloud SDK/gsutil to be authenticated (user or service account).
# If needed, uncomment one of the auth snippets:
# gcloud auth login --no-launch-browser && gcloud config set project "${PROJECT_ID}"
# gcloud auth activate-service-account --key-file="${SA_JSON}" && gcloud config set project "${PROJECT_ID}"
mapfile -t CTRL_BAMS < <(gsutil ls "${CTRL_PREFIX}/*.bam")
mapfile -t TEST_BAMS < <(gsutil ls "${TEST_PREFIX}/*.bam")
echo "[stage] Control BAMs: ${#CTRL_BAMS[@]} Test BAMs: ${#TEST_BAMS[@]}"
((${#CTRL_BAMS[@]})) || { echo "[ERR] No control BAMs found at ${CTRL_PREFIX}"; exit 1; }
((${#TEST_BAMS[@]})) || { echo "[ERR] No test BAMs found at ${TEST_PREFIX}"; exit 1; }
echo "[stage] Downloading BAMs to /content/bams (this can take time)..."
mkdir -p /content/bams/ctrl /content/bams/test
for uri in "${CTRL_BAMS[@]}"; do gsutil -m cp -n "$uri" /content/bams/ctrl/; done
for uri in "${TEST_BAMS[@]}"; do gsutil -m cp -n "$uri" /content/bams/test/; done
# Symlink to short paths in shared memory namespace
i=1; for f in /content/bams/ctrl/*.bam; do ln -sf "$f" "/dev/shm/cd/ctrl/$(printf 'c%02d.bam' "$i")"; i=$((i+1)); done
i=1; for f in /content/bams/test/*.bam; do ln -sf "$f" "/dev/shm/cd/test/$(printf 't%02d.bam' "$i")"; i=$((i+1)); done
CTRL_CSV="$(ls -1v /dev/shm/cd/ctrl/*.bam | paste -sd, -)"
TEST_CSV="$(ls -1v /dev/shm/cd/test/*.bam | paste -sd, -)"
echo "[stage] CTRL_CSV=${CTRL_CSV}"
echo "[stage] TEST_CSV=${TEST_CSV}"
# =========[ 3) PREP REFERENCE: split GTF by chromosome ]=======================
echo "[ref] fetching GTF..."
gsutil -m cp -n "${GTF_GCS}" "${REF_DIR}/annotation.gtf"
echo "[ref] splitting GTF to ${GTF_CHUNKS}"
awk '
BEGIN{OFS="\t"}
/^#/ {next}
{chr=$1; gsub(/\r$/,"",chr); out=sprintf("%s/%s.gtf","'"${GTF_CHUNKS}"'", chr); print > out}
' "${REF_DIR}/annotation.gtf"
# Normalize filenames like "chr 14" -> "chr14"
if [ -f "${GTF_CHUNKS}/chr 14.gtf" ] && [ ! -f "${GTF_CHUNKS}/chr14.gtf" ]; then
mv "${GTF_CHUNKS}/chr 14.gtf" "${GTF_CHUNKS}/chr14.gtf"
fi
# Build chromosome universe (exclude any you want)
mapfile -t ALL_CHRS < <(find "${GTF_CHUNKS}" -maxdepth 1 -type f -name 'chr*.gtf' -printf '%f\n' \
| sed 's/\.gtf$//' | sort -V | grep -v "^${EXCLUDE_CHR}$" || true)
((${#ALL_CHRS[@]})) || { echo "[ERR] No chr*.gtf in ${GTF_CHUNKS}"; exit 1; }
echo "[ref] chromosomes: ${#ALL_CHRS[@]} (excluded: ${EXCLUDE_CHR:-none})"
# =========[ 4) RUN: per-chromosome cuffdiff, HDD-friendly scheduler ]==========
cat >"${RUN_DIR}/launch_pending_strict.sh" <<'SH'
#!/usr/bin/env bash
set -euo pipefail
export LC_ALL=C LANG=C
RUN_DIR="${1:?Usage: launch_pending_strict.sh RUN_DIR}"
THREADS_PER_CHR="${THREADS_PER_CHR:-3}"
CAP_ACTIVE="${CAP_ACTIVE:-4}"
LABELS="${LABELS:-c,t}"
LIBTYPE="${LIBTYPE:-fr-unstranded}"
CHUNKS="$RUN_DIR/chunks"
RERUN="$RUN_DIR/chunks_rerun"
GTF_DIR="$RUN_DIR/ref/gtf_chunks"
LOG_DIR="$RUN_DIR/logs"
LOCK_DIR="$RUN_DIR/locks"
mkdir -p "$RERUN" "$LOG_DIR" "$LOCK_DIR"
CTRL_CSV="$(ls -1v /dev/shm/cd/ctrl/*.bam 2>/dev/null | paste -sd, - || true)"
TEST_CSV="$(ls -1v /dev/shm/cd/test/*.bam 2>/dev/null | paste -sd, - || true)"
[ -n "$CTRL_CSV" ] && [ -n "$TEST_CSV" ] || { echo "[ERR] Missing /dev/shm/cd/{ctrl,test}"; exit 1; }
mapfile -t CHRS < <(find "$GTF_DIR" -maxdepth 1 -type f -name 'chr*.gtf' -printf '%f\n' | sed 's/\.gtf$//' | sort -V)
has_all(){ d="$1"; [ -s "$d/gene_exp.diff" ] && [ -s "$d/isoform_exp.diff" ] && [ -s "$d/tss_group_exp.diff" ] && [ -s "$d/cds_exp.diff" ]; }
dir_nonempty(){ [ -d "$1" ] && find "$1" -mindepth 1 -maxdepth 1 -quit >/dev/null 2>&1; }
# never (re)start if we have outputs/log/lock
should_skip(){ c="$1";
[ -f "$LOCK_DIR/$c.lock" ] && return 0
has_all "$CHUNKS/$c" && return 0
has_all "$RERUN/$c" && return 0
dir_nonempty "$CHUNKS/$c" && return 0
dir_nonempty "$RERUN/$c" && return 0
[ -s "$LOG_DIR/cuffdiff_${c}.log" ] && return 0
pgrep -af "/cufflinks.*/cuffdiff" | grep -q "gtf_chunks/${c}\.gtf" && return 0
return 1
}
# current active/paused counts
active_n(){ ps -eo stat,cmd | awk "/cuffdiff .*gtf_chunks\\/chr/ && \$1 !~ /T/ {c++} END{print c+0}"; }
for c in "${CHRS[@]}"; do
# Skip excluded
if [ -f "$GTF_DIR/${c}.skip" ]; then echo "[skip] $c (marked skip)"; continue; fi
if should_skip "$c"; then echo "[skip] $c"; continue; fi
# throttle active jobs
while [ "$(active_n)" -ge "${CAP_ACTIVE}" ]; do
sleep 5
done
echo $$ > "$LOCK_DIR/$c.lock"
GTF="$GTF_DIR/$c.gtf"
OUT="$RERUN/$c"
LOG="$LOG_DIR/cuffdiff_${c}.log"
mkdir -p "$OUT"
echo "[start] $c (-p ${THREADS_PER_CHR})"
( export OMP_NUM_THREADS="${THREADS_PER_CHR}" MALLOC_ARENA_MAX=2 TMPDIR="/dev/shm/cd/tmp_${c}";
stdbuf -oL -eL cuffdiff --no-update-check \
-o "$OUT" -p "${THREADS_PER_CHR}" --library-type "$LIBTYPE" \
-u -L "$LABELS" "$GTF" "$CTRL_CSV" "$TEST_CSV" >> "$LOG" 2>&1 || true
# mark done if 4 tables exist
if has_all "$OUT"; then touch "$LOCK_DIR/$c.done"; fi
rm -f "$LOCK_DIR/$c.lock" ) &
done
wait
echo "[done] launch_pending_strict: all submitted jobs have finished."
SH
chmod +x "${RUN_DIR}/launch_pending_strict.sh"
echo "[run] launching per-chromosome cuffdiff..."
THREADS_PER_CHR="${THREADS_PER_CHR}" CAP_ACTIVE="${CAP_ACTIVE}" \
"${RUN_DIR}/launch_pending_strict.sh" "${RUN_DIR}"
# =========[ 5) STITCH: build classic cuffdiff outputs (neofinal) ]==============
NEOFINAL="/content/cuffdiff_out/neofinal/${RUN_TAG}"
mkdir -p "${NEOFINAL}"
echo "[merge] building neofinal (gene/isoform/tss_group/cds diffs + tracking/read-group/counts + run.info)"
declare -a BASES=(
gene_exp.diff isoform_exp.diff tss_group_exp.diff cds_exp.diff
genes.fpkm_tracking isoforms.fpkm_tracking tss_groups.fpkm_tracking cds.fpkm_tracking
genes.count_tracking isoforms.count_tracking tss_groups.count_tracking cds.count_tracking
genes.read_group_tracking isoforms.read_group_tracking tss_groups.read_group_tracking cds.read_group_tracking
run.info
)
# Choose source per chr (your earlier rule: use chunks_rerun by default, but you could special-case chr14)
src_dir_for_chr(){
local c="$1"
if [ -d "${RERUN}/${c}" ] && find "${RERUN}/${c}" -mindepth 1 -maxdepth 1 -quit >/dev/null 2>&1; then
echo "${RERUN}/${c}"
elif [ -d "${CHUNKS}/${c}" ] && find "${CHUNKS}/${c}" -mindepth 1 -maxdepth 1 -quit >/dev/null 2>&1; then
echo "${CHUNKS}/${c}"
else
echo ""
fi
}
# chromosome list as used for the run
mapfile -t CHRS_NOW < <(find "${GTF_CHUNKS}" -maxdepth 1 -type f -name 'chr*.gtf' -printf '%f\n' | sed 's/\.gtf$//' | sort -V | grep -v "^${EXCLUDE_CHR}$" || true)
stitch_one(){
local base="$1" out="${NEOFINAL}/${base}"
: >"$out"
local wrote=0 used=0
for c in "${CHRS_NOW[@]}"; do
sd="$(src_dir_for_chr "$c")"; [ -n "$sd" ] || continue
src="${sd}/${base}"; [ -s "$src" ] || continue
if [ $wrote -eq 0 ]; then
cat "$src" >> "$out"; wrote=1
else
if [[ "$base" == "run.info" ]]; then
{ echo ""; echo "### ${c} ---"; cat "$src"; } >> "$out"
else
tail -n +2 "$src" >> "$out"
fi
fi
used=$((used+1))
done
if [ $used -eq 0 ]; then
rm -f "$out"; echo "[skip] $base (no sources)"
else
# de-dup identical rows for tabular (keeps first header)
if [[ "$base" != "run.info" ]]; then
awk 'NR==1{print; next} !seen[$0]++' "$out" > "${out}.tmp" && mv "${out}.tmp" "$out"
fi
echo "[made] $base (from $used chromosomes)"
fi
}
for b in "${BASES[@]}"; do stitch_one "$b"; done
# Manifest of per-chr sources used
{
echo -e "chr\tsource_dir"
for c in "${CHRS_NOW[@]}"; do sd="$(src_dir_for_chr "$c")"; [ -n "$sd" ] && echo -e "$c\t$sd"; done
} > "${NEOFINAL}/_manifest_sources.tsv"
# =========[ 6) PUBLISH: rsync per-chr outputs + neofinal back to GCS ]=========
echo "[publish] syncing per-chr outputs to ${GCS_RUN_DIR}"
gsutil -m rsync -r "${RUN_DIR}" "${GCS_RUN_DIR}"
echo "[publish] syncing neofinal to ${GCS_NEOFINAL}"
gsutil -m rsync -r "${NEOFINAL}" "${GCS_NEOFINAL}"
# =========[ 7) (OPTIONAL) PUBLISH to Google Drive via rclone ]=================
if [ "${PUBLISH_TO_DRIVE}" = "1" ]; then
echo "[drive] publishing to Google Drive via rclone..."
# lightweight rclone install into $HOME/bin if needed
if ! command -v rclone >/dev/null 2>&1; then
work="/tmp/rclone_dl"; rm -rf "$work"; mkdir -p "$work"
curl -fsSL https://downloads.rclone.org/rclone-current-linux-amd64.zip -o "$work/rclone.zip"
unzip -q "$work/rclone.zip" -d "$work"
pkg="$(find "$work" -maxdepth 1 -type d -name 'rclone-*-linux-amd64' | head -n1)"
mkdir -p "$HOME/bin"; cp "$pkg/rclone" "$HOME/bin/rclone"; chmod +x "$HOME/bin/rclone"
export PATH="$HOME/bin:$PATH"
fi
# GCS remote using service account; Drive remote must be pre-auth’d once (see notes)
rclone listremotes | grep -qx 'mygcs:' || rclone config create mygcs gcs service_account_file="${SA_JSON}" --non-interactive
rclone listremotes | grep -qx 'mydrive:' || rclone config create mydrive drive scope=drive file_scope=false --non-interactive || true
# If mydrive not yet authorized, run on your laptop: rclone authorize "drive"
# then set token here: rclone config update mydrive token '<JSON>' --non-interactive
rclone copy -P "${NEOFINAL}" "${DRIVE_DST_BASE}/neofinal/${RUN_TAG}" --create-empty-src-dirs
rclone copy -P "${RUN_DIR}" "${DRIVE_DST_BASE}/${RUN_TAG}" --create-empty-src-dirs
fi
# =========[ 8) SUMMARY ]=======================================================
echo
echo "======================= SUMMARY ======================="
echo "[run] ${RUN_DIR}"
echo "[per-chr] ${GCS_RUN_DIR}"
echo "[neofinal] ${NEOFINAL}"
echo "[gcs out] ${GCS_NEOFINAL}"
if [ "${PUBLISH_TO_DRIVE}" = "1" ]; then
echo "[drive] ${DRIVE_DST_BASE}/neofinal/${RUN_TAG}"
fi
echo "======================================================="
##CODE CELL -6- BIAS CORRECTION AND RESIDUALISATION ================================¶
Raw → Residualized DE (combined)¶
================================¶
Required inputs in /content/¶
- gene_exp.diff¶
- isoform_exp.diff¶
- genes.read_group_tracking (per-replicate FPKM with condition)¶
¶
Outputs (all in /content/):¶
- TopUpregulated_Raw_Cuffdiff.csv¶
- PerGene_Isoform_Metrics.csv¶
- Residualized_DE_splicingBurden.csv¶
- TopUpregulated_Residualized_Splicing.csv¶
- Residualization_Combined_Table.csv¶
- TopUpregulated_Residualized_SampleLevel.csv¶
- Top10_Raw_vs_Residualized_BothModels.csv¶
- TopUpregulated_BiasCorrected_Combined.csv¶
- Residualization_geneLevel_OLS_Summary.txt¶
bold text !pip install statsmodels import pandas as pd, numpy as np import statsmodels.api as sm from pathlib import Path
---------- Paths ----------¶
GENE_DIFF = "/content/gene_exp.diff" ISO_DIFF = "/content/isoform_exp.diff" READ_GROUP = "/content/genes.read_group_tracking"
pd.set_option("display.max_columns", 200) pd.set_option("display.width", 180)
def safe_log2fc(v2, v1, eps=1e-9): v2 = pd.to_numeric(v2, errors="coerce") v1 = pd.to_numeric(v1, errors="coerce") return np.log2((v2 + eps) / (v1 + eps))
============================================================¶
1) RAW CUFFDIFF: top upregulated genes (gene_exp.diff)¶
============================================================¶
gene = pd.read_csv(GENE_DIFF, sep="\t") gene["significant"] = gene["significant"].astype(str).str.lower()
gene_sig = gene[gene["significant"].eq("yes")].copy() gene_sig["log2FC_raw"] = safe_log2fc(gene_sig["value_2"], gene_sig["value_1"]) gene_sig = ( gene_sig .replace([np.inf, -np.inf], np.nan) .dropna(subset=["log2FC_raw"]))
top_raw = ( gene_sig .sort_values("log2FC_raw", ascending=False) [["gene","gene_id","log2FC_raw","p_value","q_value","value_1","value_2"]] .head(50).reset_index(drop=True)) top_raw.insert(0, "Rank_Raw", np.arange(1, len(top_raw)+1)) top_raw.to_csv("/content/TopUpregulated_Raw_Cuffdiff.csv", index=False)
print("\n=== Top 10 Upregulated (RAW Cuffdiff) ===") print(top_raw.head(10).to_string(index=False))
============================================================¶
2) ISOFORM SUMMARY → per-gene splicing burden metrics¶
============================================================¶
iso = pd.read_csv(ISO_DIFF, sep="\t") iso["significant"] = iso["significant"].astype(str).str.lower().eq("yes")
log2(fold_change)¶
if "log2(fold_change)" in iso.columns: iso["log2fc"] = pd.to_numeric(iso["log2(fold_change)"], errors="coerce") else: iso["log2fc"] = safe_log2fc(iso.get("value_2"), iso.get("value_1"))
abundance filter (mean FPKM >= 1)¶
if {"value_1","value_2"}.issubset(iso.columns): iso["mean_fpkm"] = ( pd.to_numeric(iso["value_1"], errors="coerce").fillna(0) + pd.to_numeric(iso["value_2"], errors="coerce").fillna(0) ) / 2.0 iso = iso[iso["mean_fpkm"] >= 1.0].copy() else: iso["mean_fpkm"] = np.nan # no filter if missing
grp = iso.groupby("gene", dropna=False) per_gene = pd.DataFrame({ "n_isoforms": grp["test_id"].count(), "n_sig_iso": grp["significant"].sum(), "n_sig_up": grp.apply( lambda d: (pd.to_numeric(d["log2fc"], errors="coerce") > 0).sum(), include_groups=False ), "n_sig_down": grp.apply( lambda d: (pd.to_numeric(d["log2fc"], errors="coerce") < 0).sum(), include_groups=False ), }).reset_index().rename(columns={"gene":"Gene"})
Splicing-burden index (up - down)¶
per_gene["SBI"] = per_gene["n_sig_up"] - per_gene["n_sig_down"] per_gene.to_csv("/content/PerGene_Isoform_Metrics.csv", index=False)
============================================================¶
3) GENE-LEVEL RESIDUALIZATION vs splicing burden¶
log2FC_raw ~ n_sig_iso + SBI (OLS)¶
============================================================¶
de = gene_sig.rename(columns={"gene":"Gene"})[ ["Gene","gene_id","log2FC_raw","q_value","p_value"] ] merged_genelevel = de.merge(per_gene, on="Gene", how="left")
for c in ["n_isoforms","n_sig_iso","n_sig_up","n_sig_down","SBI"]: merged_genelevel[c] = merged_genelevel[c].fillna(0).astype(int)
X = sm.add_constant( merged_genelevel[["n_sig_iso","SBI"]].astype(float).fillna(0.0) ) y = merged_genelevel["log2FC_raw"].astype(float)
model_gene = sm.OLS(y, X, missing="drop").fit(cov_type="HC3") merged_genelevel["log2FC_resid_splicing"] = model_gene.resid
Save splicing-residualized table & top list¶
merged_genelevel.to_csv("/content/Residualized_DE_splicingBurden.csv", index=False)
top_resid_splice = ( merged_genelevel .sort_values("log2FC_resid_splicing", ascending=False) [["Gene","gene_id","log2FC_resid_splicing","log2FC_raw","n_sig_iso","SBI"]] .head(50).reset_index(drop=True) ) top_resid_splice.insert( 0, "Rank_Residual_Splicing", np.arange(1, len(top_resid_splice)+1) ) top_resid_splice.to_csv( "/content/TopUpregulated_Residualized_Splicing.csv", index=False )
print("\n=== Top 10 Upregulated (Residualized vs Splicing Burden) ===") print(top_resid_splice.head(10).to_string(index=False))
Save OLS summary (for supplement)¶
Path("/content/Residualization_geneLevel_OLS_Summary.txt").write_text( str(model_gene.summary()) )
============================================================¶
4) SAMPLE-LEVEL RESIDUALIZATION using genes.read_group_tracking¶
log2FPKM ~ condition + replicate (per gene)¶
============================================================¶
rg = pd.read_csv(READ_GROUP, sep="\t", low_memory=False)
Normalise condition labels¶
rg["condition"] = rg["condition"].astype(str).str.strip() cond_map = { "c": "normal", "control": "normal", "t": "tumor", "test": "tumor", } rg = rg[rg["condition"].isin(cond_map.keys())].copy() rg["condition_norm"] = rg["condition"].map(cond_map)
Build per-row table¶
if "FPKM" not in rg.columns: raise ValueError("Expected column 'FPKM' not found in genes.read_group_tracking")
M = rg[["tracking_id","condition_norm","replicate","FPKM"]].copy() M = M.rename(columns={ "tracking_id": "gene_id", "condition_norm": "condition" })
log2(FPKM + 1e-6)¶
M["log2FPKM"] = np.log2(pd.to_numeric(M["FPKM"], errors="coerce") + 1e-6)
Add Gene symbol from gene_exp.diff¶
gene_symbol_map = gene[["gene_id","gene"]].drop_duplicates() M = M.merge(gene_symbol_map, on="gene_id", how="left") M = M.rename(columns={"gene":"Gene"})
Design matrix: condition (tumor vs normal) + replicate dummies¶
M["cond_Tumor"] = (M["condition"] == "tumor").astype(int) M["replicate"] = M["replicate"].astype(str) rep_dummies = pd.get_dummies(M["replicate"], prefix="rep", drop_first=True) M = pd.concat([M, rep_dummies], axis=1)
design_cols = ["cond_Tumor"] + [c for c in M.columns if c.startswith("rep_")]
coef_rows = [] for g, df_g in M.groupby("Gene", dropna=False): # need both tumor and normal for this gene if df_g["cond_Tumor"].nunique() < 2: continue Z = df_g[design_cols].copy().astype(float).fillna(0.0) Z = sm.add_constant(Z) yy = df_g["log2FPKM"].astype(float) # Ensure enough degrees of freedom if len(df_g) <= Z.shape[1]: continue try: fit = sm.OLS(yy, Z, missing="drop").fit(cov_type="HC3") coef_rows.append({ "Gene": g, "log2FC_resid_sample": float(fit.params.get("cond_Tumor", np.nan)) }) except Exception: # silently skip genes that fail to fit continue
sample_resid = pd.DataFrame(coef_rows)
============================================================¶
5) MERGE ALL + BUILD COMBINED BIAS-CORRECTED SCORE¶
============================================================¶
final = merged_genelevel.merge(sample_resid, on="Gene", how="left")
def zseries(s): s = pd.to_numeric(s, errors="coerce") m = s.mean(skipna=True) sd = s.std(skipna=True, ddof=0) if not np.isfinite(sd) or sd == 0: return pd.Series(np.nan, index=s.index) return (s - m) / sd
final["z_splicing"] = zseries(final["log2FC_resid_splicing"]) final["z_sample"] = zseries(final["log2FC_resid_sample"])
def combined_row(row): vals = [] if pd.notna(row["z_splicing"]): vals.append(row["z_splicing"]) if pd.notna(row["z_sample"]): vals.append(row["z_sample"]) if not vals: return np.nan return float(np.mean(vals))
final["CombinedScore"] = final.apply(combined_row, axis=1)
Save master table¶
final.to_csv("/content/Residualization_Combined_Table.csv", index=False)
Top by sample-residual¶
top_resid_sample = ( final.dropna(subset=["log2FC_resid_sample"]) .sort_values("log2FC_resid_sample", ascending=False) [["Gene","gene_id","log2FC_resid_sample"]] .head(50).reset_index(drop=True) ) top_resid_sample.insert( 0, "Rank_Residual_Sample", np.arange(1, len(top_resid_sample)+1) ) top_resid_sample.to_csv( "/content/TopUpregulated_Residualized_SampleLevel.csv", index=False )
Top by combined bias-corrected score¶
top_combined = ( final.dropna(subset=["CombinedScore"]) .sort_values("CombinedScore", ascending=False) [["Gene","gene_id","CombinedScore", "log2FC_resid_sample","log2FC_resid_splicing","log2FC_raw"]] .head(50).reset_index(drop=True) ) top_combined.insert(0, "Rank_Combined", np.arange(1, len(top_combined)+1)) top_combined.to_csv( "/content/TopUpregulated_BiasCorrected_Combined.csv", index=False )
Side-by-side comparison (top-10)¶
t10_raw = (top_raw.head(10)[["Rank_Raw","gene","log2FC_raw"]] .rename(columns={"gene":"Gene","log2FC_raw":"Raw_log2FC"})) t10_spl = top_resid_splice.head(10)[ ["Rank_Residual_Splicing","Gene","log2FC_resid_splicing"] ] t10_smp = top_resid_sample.head(10)[ ["Rank_Residual_Sample","Gene","log2FC_resid_sample"] ] cmp = (t10_raw .merge(t10_spl, on="Gene", how="outer") .merge(t10_smp, on="Gene", how="outer")) cmp.to_csv("/content/Top10_Raw_vs_Residualized_BothModels.csv", index=False)
print("\n=== Top 10 Upregulated (Residualized, Sample-level) ===") print(top_resid_sample.head(10).to_string(index=False))
print("\n=== Top 10 Upregulated (Combined bias-corrected score) ===") print(top_combined.head(10).to_string(index=False))
print("\nFiles written:") for p in [ "/content/TopUpregulated_Raw_Cuffdiff.csv", "/content/PerGene_Isoform_Metrics.csv", "/content/Residualized_DE_splicingBurden.csv", "/content/TopUpregulated_Residualized_Splicing.csv", "/content/Residualization_Combined_Table.csv", "/content/TopUpregulated_Residualized_SampleLevel.csv", "/content/Top10_Raw_vs_Residualized_BothModels.csv", "/content/TopUpregulated_BiasCorrected_Combined.csv", "/content/Residualization_geneLevel_OLS_Summary.txt", ]:
print(" -", p)¶
CODE CELL 7-RMATS CODE Uploaded the bam files and gtf annotation files in /content in colab enterprise with a customised runtime e232 and 200 gb disk space and ran rmats with the following code
At first , rMATS-Turbo (v4.1.2) was installed within a dedicated Conda environment with all required dependencies pinned, enabling consistent detection and quantification of alternative splicing events across all RNA-seq samples.
# This is formatted as code
# ============================================================
# 1) Install Miniconda (if not already installed)
# ============================================================
wget --quiet https://repo.anaconda.com/miniconda/Miniconda3-latest-Linux-x86_64.sh \
-O miniconda.sh
bash miniconda.sh -b -p /content/miniconda
eval "$(/content/miniconda/bin/conda shell.bash hook)"
# ============================================================
# 2) Create conda environment for rMATS
# ============================================================
conda create -y -n rmats_env python=3.8
conda activate rmats_env
# ============================================================
# 3) Install dependencies exactly as earlier
# ============================================================
conda install -y -c bioconda samtools
conda install -y -c bioconda star
conda install -y -c bioconda rmats=4.1.2
# 4.1.2 is what your folder structure / previous logs matched
# ============================================================
# 4) Verify installation
# ============================================================
rmats.py --version
samtools --version
star --version
echo "rMATS conda environment ready!"
Next we ran rmats with the following code
# This is formatted as code
rm -rf /content/rmats_tmp_full /content/rmats_output_full
rmats.py \
--b1 /content/b1.txt \
--b2 /content/b2.txt \
--gtf /content/gencode.v40.annotation.gtf \
--od /content/rmats_output_full \
--tmp /content/rmats_tmp_full \
--readLength 126 \
--variable-read-length \
--libType fr-unstranded \
--novelSS \
--allow-clipping \
--nthread 32
CODE CELL-8 RUNNING MNTJULIP In order to furhter analyse intron level driver changes between the 2 sets of bam files, i downloaded MntJULiP latest github releease 1.5.2.zip and uploaded it in colab enterprise notebook and ran MntJULiP with the following code
# This is formatted as code
%%bash
# ================== Reproducible MntJULiP 1.5.2 Run (Colab) ==================
set -euo pipefail
cd /content
# --- Inputs (edit paths if needed) ---
ZIP="/content/MntJULiP-1.5.2.zip"
SRC="/content/MntJULiP-1.5.2"
BAM_LIST="/content/bam_list.txt"
GTF="/content/gencode.v40.annotation.gtf"
OUT="/content/mntjulip_out"
THREADS=6 # safer than 32 on Colab
# --- Fresh environment ---
rm -rf "$OUT" /content/mntjulip_env
python3.8 -V >/dev/null 2>&1 || (apt-get -qq update && apt-get -qq install -y python3.8 python3.8-venv python3.8-dev >/dev/null)
apt-get -qq install -y samtools build-essential >/dev/null
python3.8 -m venv /content/mntjulip_env
source /content/mntjulip_env/bin/activate
python -V
pip -q install --upgrade pip
# --- Pin compatible Python deps (fixes np.int issue) ---
# NumPy 1.23.5 still has np.int; we also pin pandas & cython & pystan versions used earlier.
pip -q install "numpy==1.23.5" "pandas==1.3.5" "cython==0.29.36" "pystan==2.19.1.1" \
"statsmodels==0.12.2" "scikit-learn" "pysam" "tqdm" "joblib" "psutil" "cloudpickle"
echo "[ENV] Versions:"
python - <<'PY'
import sys, numpy, pandas, pystan, statsmodels
print("Python:", sys.version.split()[0])
print("NumPy:", numpy.__version__)
print("Pandas:", pandas.__version__)
print("PyStan:", pystan.__version__)
print("statsmodels:", statsmodels.__version__)
PY
# --- Install MntJULiP from your uploaded ZIP (clean restore) ---
rm -rf "$SRC"
unzip -q "$ZIP" -d /content
cd "$SRC"
python setup.py install >/dev/null
cd /content
# --- Patch deprecated NumPy aliases in models.py (future-proof) ---
python - <<'PY'
from pathlib import Path
import re
p = Path("/content/MntJULiP-1.5.2/models.py")
s = p.read_text()
s = re.sub(r"dtype\s*=\s*np\.int\b", "dtype=int", s)
s = re.sub(r"dtype\s*=\s*np\.float\b", "dtype=float", s)
s = re.sub(r"dtype\s*=\s*np\.bool\b", "dtype=bool", s)
p.write_text(s)
print("[PATCH] models.py: replaced np.int/np.float/np.bool with builtins.")
PY
# --- Thread caps (stability on Colab) ---
export OMP_NUM_THREADS=2
export OPENBLAS_NUM_THREADS=2
export MKL_NUM_THREADS=2
export NUMEXPR_NUM_THREADS=2
export VECLIB_MAXIMUM_THREADS=2
export BLIS_NUM_THREADS=2
export STAN_NUM_THREADS=1
# --- Run MntJULiP (reuses splice files automatically if rerun with --save-tmp true) ---
mkdir -p "$OUT"
set -x
python "$SRC/run.py" \
--bam-list "$BAM_LIST" \
--anno-file "$GTF" \
--out-dir "$OUT" \
--num-threads "$THREADS"
set +x
# --- Record exact environment for your methods section (reproducibility) ---
pip freeze | sed 's/@ file:.*//g' > "$OUT/requirements_exact.txt"
python - <<'PY' > "$OUT/run_metadata.txt"
import platform, sys, numpy, pandas, pystan, statsmodels, json, os
meta = {
"python": sys.version.split()[0],
"platform": platform.platform(),
"numpy": numpy.__version__,
"pandas": pandas.__version__,
"pystan": getattr(pystan, "__version__", "n/a"),
"statsmodels": statsmodels.__version__,
"threads": {
"num_threads_flag": int(os.environ.get("THREADS", "NA")) if os.environ.get("THREADS") else "NA",
"OMP_NUM_THREADS": os.environ.get("OMP_NUM_THREADS",""),
"OPENBLAS_NUM_THREADS": os.environ.get("OPENBLAS_NUM_THREADS",""),
"MKL_NUM_THREADS": os.environ.get("MKL_NUM_THREADS",""),
}
}
print(json.dumps(meta, indent=2))
PY
# --- Zip outputs and auto-download (works in Colab) ---
STAMP=$(date +%Y%m%d_%H%M%S)
ZIP_OUT="/content/mntjulip_out_${STAMP}.zip"
(cd /content && zip -qr "$ZIP_OUT" "$(basename "$OUT")")
echo "[ZIP] Created: $ZIP_OUT"
python - <<PY
try:
from google.colab import files
files.download("$ZIP_OUT")
print("[DL] Initiated download of $ZIP_OUT")
except Exception as e:
print("[DL] Skipped auto-download (not Colab or permission). File at: $ZIP_OUT")
PY
echo "[DONE] MntJULiP run complete. Outputs in: $OUT"
# ==============================================================================
The versions of the software that were used to run mntjulip in colab enterprise
# This is formatted as code
# Core environment
python==3.8.x
pip==24.0
setuptools==70.0.0
wheel==0.43.0
# MntJULiP 1.5.2 dependencies
numpy==1.23.5
pandas==1.3.5
cython==0.29.36
pystan==2.19.1.1
statsmodels==0.12.2
scikit-learn==1.3.0
pysam==0.22.0
tqdm==4.66.4
joblib==1.3.2
psutil==5.9.8
cloudpickle==3.0.0
# System tools
samtools==1.17 (via apt)
build-essential==12.9
runtime metadata used
# This is formatted as code
{
"python": "3.8.20",
"platform": "Linux-5.15.120+-x86_64-with-glibc2.31",
"numpy": "1.23.5",
"pandas": "1.3.5",
"pystan": "2.19.1.1",
"statsmodels": "0.12.2",
"threads": {
"num_threads_flag": 6,
"OMP_NUM_THREADS": "2",
"OPENBLAS_NUM_THREADS": "2",
"MKL_NUM_THREADS": "2"
}
}
CODE CELL 9-SPLICING FREQUENCY RANKING, INTEGRATION, CROSS COHORT VALIDATION¶
- Overview
All splicing analyses were performed using a unified, fully reproducible Python pipeline integrating rMATS (exon-level alternative splicing) and MntJULiP (intron usage and differential intron regulation), with structural and cross-cohort validation from GENCODE v40 and Snaptron (GTEx + TCGA). This pipeline was executed as a single code block in Google Colab, requiring no manual intervention or external dependencies.
The workflow ensures deterministic, reproducible identification of splicing events and the correspondence between exon-defined and intron-defined signatures across genes.
- Computational Environment
Analyses were performed on Google Colab Enterprise:
Runtime: e2-standard-32 (32 vCPUs, ~128 GB RAM, ~300 GB SSD disk)
Operating System: Ubuntu 22.04 LTS
Python: 3.10+
Core Packages:
pandas 2.1+
numpy 1.26+
requests 2.31+
Standard libraries: gzip, zipfile, io, json
All analyses are containerized within this environment for reproducibility.
- Input Data and Annotation Sources RNA-seq Alignment
28 tumor and matched normal RNA-seq datasets were aligned with:
STAR v2.7.x
Reference: GRCh38
Annotation: GENCODE v40 GTF
Splicing Quantification Tools
rMATS-Turbo v4.1.2
Event types: SE, RI, MXE, A5SS, A3SS
Both JC and JCEC modes were used
Significant events defined by:
FDR ≤ 0.05
|ΔPSI| ≥ 0.20
MntJULiP v1.5.2
DSA (diff_introns.txt): base intron significance
DSR (diff_spliced_introns.txt): intron clusters with regulatory grouping
Significant events defined by q ≤ 0.05
Reference Gene Models
GENCODE v40 transcripts were used for exon and intron validation.
External Cohort Validation
Using Snaptron:
GTEx tissues:
Esophagus – Mucosa
Esophagus – Gastroesophageal Junction (GEJ)
TCGA Pan-Cancer:
Used as TCGA-level junction support for ESCA-like introns.
- Unified MntJULiP Integration (DSR→DSA Fallback)
MntJULiP provides two sources:
DSA (diff_introns.txt)
Contains all introns tested
Includes q-values
Does not include regulatory group annotation
DSR (diff_spliced_introns.txt)
Contains intron clusters with:
group_id
group_loc (cluster anchor)
group_structure ("o"/"i")
group_q
Integration Strategy
The pipeline constructs a unified intron table as:
If intron exists in DSR: Use DSR annotations (group_id, group_structure, group_loc, group_q) MntJULiP_source = "DSR" Else: Use DSA annotations MntJULiP_source = "DSA"
This ensures maximal coverage (from DSA) while preserving regulatory structure (from DSR).
Filtering applied:
q ≤ 0.05
|Δ intron usage| ≥ 0.20 (if available)
Optional: restrict to non-coding / pseudogene loci using GENCODE biotypes
- Anchor-Based Event–Intron Matching Logic (Generalized Algorithm)
(GAS5 name removed; the algorithm is universal.)
Pairing exon-defined rMATS events with MntJULiP introns follows a deterministic anchor-based framework.
Step 1 — rMATS Junction Segment Reconstruction
For each rMATS event type (SE, RI, MXE, A5SS, A3SS), all splice junction segments are reconstructed from the rMATS coordinate fields.
Step 2 — Event Center Position
The event center (evctr) is computed as the mean of junction midpoints.
Step 3 — Event-Specific Anchor Points
Anchor positions represent the biological decision point of a splice event and include:
RI → upstream exon start
MXE → start positions of both alternative exons
A5SS → short-isoform donor site
Other event types inherit anchor positions from their rMATS-defined segment boundaries when relevant.
Step 4 — Identify Overlapping MntJULiP Introns
All MntJULiP introns overlapping any rMATS junction segment by ≥1 bp are retained as candidates.
Step 5 — Hierarchical Selection of Best-Matching Intron
Candidate introns are ranked through strict, ordered criteria:
Anchor co-localization (≤1 bp between group_loc and anchor)
Group structure preference: group_structure = "o"
Lower q-value (higher significance)
Higher overlap length (bp)
Minimal distance between intron midpoint and event center (evctr)
Minimal genomic start coordinate (final deterministic tie-breaker)
This framework ensures system-wide reproducibility and removes subjectivity from event–intron pairing.
- GENCODE v40 Structural Support
Each selected intron is evaluated against GENCODE transcript models:
Exact intron match (boundaries match a GENCODE intron)
Retained intron within exon (intron falls entirely inside a GENCODE exon)
Unsupported / novel intron
Transcript IDs supporting each match are recorded, enabling structural classification of intron usage.
- External Cohort Junction Support (GTEx and TCGA)
For every intron, Snaptron is queried to obtain:
Number of samples with non-zero junction read counts
Mean junction read count among positive samples
Presence/absence in:
GTEx Esophagus – Mucosa
GTEx Esophagus – GEJ
TCGA Pan-Cancer (ESCA-like)
This allows classification into:
normal-esophagus-supported introns
tumor-specific introns
widely-used housekeeping introns
- Gene-Level Splicing Burden and Regulation Mode
All rMATS and MntJULiP events are unified per gene.
Splicing burden (per gene)
Total unique rMATS + MntJULiP events per gene
Sorted to produce SplicingFrequencyRanking_byGene.csv
Regulation Mode
A gene is labeled as:
both (exon + intron regulated)
exon-level-only
intron-only
neither
Saved as SplicingRegulationModes_TopMostSpliced.csv
Gene selection
The pipeline automatically selects:
Top N most spliced genes
Forced inclusion list (e.g., specific lncRNAs)
- Per-Gene Supplementary Output Tables
For each selected gene, two outputs are produced:
Minimal, clean table suitable for direct supplementation
Contains rMATS event, matched intron, q-values, overlap, match type, source (DSR/DSA)
Full structural and cohort annotation
Includes GENCODE status and GTEx/TCGA support
A combined table for all selected genes is saved as:
TopGenes_Combined_With_Annotation_and_CohortPresence.csv
- Reproducibility
The entire workflow executes from fixed ZIP archives (rMATS, MntJULiP, GENCODE).
All operations are deterministic and use strict ordering rules.
No randomness, external state, or user interaction affects the output.
Re-running the Colab code block yields identical byte-for-byte CSV outputs, satisfying journal reproducibility standards
# This is formatted as code
# ================================================================
# Unified rMATS ↔ MntJULiP pipeline (DSR → DSA fallback)
# + GENCODE intron/exon annotation
# + Snaptron GTEx (EM, GEJ) + TCGA (PanCan ESCA-like) junction usage
#
# Outputs:
# /content/SplicingFrequencyRanking_byGene.csv
# /content/SplicingRegulationModes_TopMostSpliced.csv
# /content/<GENE>_PublicationView.csv
# /content/<GENE>_With_Annotation_and_CohortPresence.csv
# /content/TopGenes_Combined_With_Annotation_and_CohortPresence.csv
#
# Key behaviour:
# - MntJULiP introns come from BOTH:
# * diff_spliced_introns.txt (DSR) → carries group_id / group_loc / group_q
# * diff_introns.txt (DSA) → full intron set, with q
# - If a given intron has DSR info, it is labeled MntJULiP_source = "DSR"
# and group_* columns are filled.
# - If only DSA has that intron, it's MntJULiP_source = "DSA"
# and group_* columns are NaN.
# - GAS5 anchor logic is preserved.
# - GENCODE intron/exon support and Snaptron cohorts are added per event.
# ================================================================
import os, re, io, gzip, zipfile, math, json, shutil, time
from pathlib import Path
import pandas as pd
import numpy as np
import requests
pd.set_option("display.max_columns", 200)
pd.set_option("display.width", 260)
# ---------------- USER PATHS ----------------
RMATS_ZIP = "/content/rmats_output_full.zip"
MNT_ZIP = "/content/mntjulip_out.zip"
GTF_ZIP = "/content/gencode.v40.annotation.gtf.zip"
# ---------------- OUTPUTS -------------------
OUT_FREQ_RANK = "/content/SplicingFrequencyRanking_byGene.csv"
OUT_REGMODES = "/content/SplicingRegulationModes_TopMostSpliced.csv"
OUT_COMBINED_ANNOT = "/content/TopGenes_Combined_With_Annotation_and_CohortPresence.csv"
# ---------------- PARAMS --------------------
ALPHA_FDR = 0.05
MIN_ABS_DPSI = 0.20
MNT_Q = 0.05 # synonym of ALPHA_FDR for MntJULiP
NONCODING_ONLY = True # keep True for lncRNAs/pseudogenes focus
TOP_PICK_N = 3
FORCE_INCLUDE = ["GAS5", "PVT1", "MUC20-OT1", "MALAT1"]
COLOC_TOL = 1
GENCODE_SLOP = 0 # exact intron match by default
GTEX_TISSUES = {
"EM": "Esophagus - Mucosa",
"GEJ": "Esophagus - Gastroesophageal Junction"
}
SNAP_GTX = "https://snaptron.cs.jhu.edu/gtex/snaptron"
SNAP_TCGA = "https://snaptron.cs.jhu.edu/tcga/snaptron"
# --------------- UTILITIES ------------------
def _find_first(cols, *cands):
low = {c.lower(): c for c in cols}
for c in cands:
if c.lower() in low:
return low[c.lower()]
return None
def _gi(v):
try:
if v is None or (isinstance(v, float) and np.isnan(v)):
return None
return int(float(str(v).strip()))
except Exception:
return None
def _span_overlap(a, b):
s = max(a[0], b[0]); e = min(a[1], b[1])
return max(0, e - s)
def _normalize_chr(x):
return str(x).replace("chr", "")
EVENT_RE = re.compile(r"(?:^|/|\.)(SE|RI|A3SS|A5SS|MXE|TXI)\.MATS\.(?:JC|JCEC)\.txt$", re.IGNORECASE)
def _largest_entries_by_key(zf, key_fn):
bykey = {}
for info in zf.infolist():
key = key_fn(info)
if key is None:
continue
if (key not in bykey) or (info.file_size > bykey[key].file_size):
bykey[key] = info
return bykey
# ----------- GTF open & noncoding filter set -----------
def _open_gtf_any_from_zip(zip_path):
with zipfile.ZipFile(zip_path, 'r') as z:
gtf_name = next((n for n in z.namelist() if n.endswith(".gtf")), None)
if gtf_name is None:
raise FileNotFoundError("No .gtf inside GTF_ZIP")
return io.TextIOWrapper(z.open(gtf_name), encoding="utf-8", errors="ignore")
def load_noncoding_symbol_set(gtf_zip):
if not NONCODING_ONLY:
return None
nc_types = {
"lncRNA","lincRNA","macro_lncRNA","bidirectional_promoter_lncRNA",
"antisense","sense_intronic","sense_overlapping","3prime_overlapping_ncRNA",
"processed_transcript","non_coding","ncRNA",
"miRNA","snRNA","snoRNA","rRNA","tRNA","scaRNA","srpRNA","vaultRNA","Y_RNA","Mt_rRNA","Mt_tRNA",
"pseudogene","transcribed_pseudogene","processed_pseudogene","unprocessed_pseudogene",
"rRNA_pseudogene","tRNA_pseudogene","snRNA_pseudogene","snoRNA_pseudogene",
"unitary_pseudogene","polymorphic_pseudogene","artifact","TEC"
}
nc = set()
with _open_gtf_any_from_zip(gtf_zip) as fh:
for line in fh:
if not line or line.startswith("#"):
continue
p = line.rstrip("\n").split("\t")
if len(p) < 9 or p[2] != "gene":
continue
attrs = p[8]
def grab(k):
m = re.search(rf'{k}\s+"([^"]+)"', attrs)
return m.group(1) if m else None
gname = grab("gene_name")
biotype = grab("gene_type") or grab("gene_biotype")
if gname and biotype and biotype in nc_types:
nc.add(gname)
return nc
NONCODING_SET = load_noncoding_symbol_set(GTF_ZIP) if NONCODING_ONLY else None
# ----------- rMATS loaders (all + significant) -----------
def _read_tsv(zf, info):
with zf.open(info) as f:
return pd.read_csv(f, sep="\t")
def _event_key_from_row(evt, row, colnames):
if evt == "SE":
cols = ["exonStart_0base","exonEnd","upstreamES","upstreamEE","downstreamES","downstreamEE"]
elif evt == "RI":
cols = ["exonStart_0base","exonEnd","upstreamEE","downstreamES"]
elif evt in ("A3SS","A5SS"):
cols = ["longExonStart_0base","longExonEnd","flankingES","flankingEE","shortES","shortEE"]
elif evt == "MXE":
cols = ["1stExonStart_0base","1stExonEnd","2ndExonStart_0base","2ndExonEnd",
"upstreamES","upstreamEE","downstreamES","downstreamEE"]
else:
cols = []
cc = _find_first(colnames, "chr","chrom","chromosome") or "chr"
ss = _find_first(colnames, "strand") or "strand"
parts = [str(row.get(cc, "")).replace("chr", "").strip(),
str(row.get(ss, "")).strip()]
for c in cols:
v = row.get(c, None)
parts.append("" if v is None else ("" if _gi(v) is None else str(_gi(v))))
parts.append(evt)
return "_".join(parts)
def load_all_rmats_events(zip_path):
with zipfile.ZipFile(zip_path, 'r') as z:
def k(info):
m = EVENT_RE.search(info.filename)
if not m:
return None
evt = m.group(1).upper()
variant = "JCEC" if info.filename.endswith(".MATS.JCEC.txt") else "JC"
return (evt, variant)
chosen = _largest_entries_by_key(z, k)
rows = []
for (evt, variant), info in sorted(chosen.items()):
df = _read_tsv(z, info)
df.columns = [c.strip() for c in df.columns]
gcol = _find_first(df.columns, "geneSymbol","GeneSymbol","gene_name","gene","Gene","GeneID","gene_id")
if gcol is None:
continue
if NONCODING_SET is not None:
df = df[df[gcol].astype(str).isin(NONCODING_SET)]
if df.empty:
continue
colnames = list(df.columns)
df["__event_key"] = df.apply(lambda r: _event_key_from_row(evt, r, colnames), axis=1)
df["gene"] = df[gcol].astype(str)
df["__evt_type"] = evt
df["__variant"] = variant
keep = ["gene","__event_key","__evt_type","__variant",
"chr","Chromosome","chrom","strand",
"exonStart_0base","exonEnd","upstreamES","upstreamEE","downstreamES","downstreamEE",
"longExonStart_0base","longExonEnd","flankingES","flankingEE","shortES","shortEE",
"1stExonStart_0base","1stExonEnd","2ndExonStart_0base","2ndExonEnd"]
rows.append(df[[c for c in keep if c in df.columns]].drop_duplicates())
if not rows:
return pd.DataFrame(columns=["gene","__event_key","__evt_type","__variant"])
return pd.concat(rows, ignore_index=True).drop_duplicates(["__event_key"])
def load_sig_rmats_events(zip_path, alpha=0.05, min_abs_dpsi=0.20):
with zipfile.ZipFile(zip_path, 'r') as z:
def k(info):
m = EVENT_RE.search(info.filename)
if not m:
return None
evt = m.group(1).upper()
variant = "JCEC" if info.filename.endswith(".MATS.JCEC.txt") else "JC"
return (evt, variant)
chosen = _largest_entries_by_key(z, k)
parts = []
for (evt, variant), info in sorted(chosen.items()):
df = _read_tsv(z, info)
df.columns = [c.strip() for c in df.columns]
gcol = _find_first(df.columns, "geneSymbol","GeneSymbol","gene_name","gene","Gene","GeneID","gene_id")
fdr_col = _find_first(df.columns, "FDR","fdr","qvalue","q_value","adjp","adj_p")
dpsi_col= _find_first(df.columns, "IncLevelDifference","incLevelDifference","deltaPSI","dPSI","delta_psi")
if gcol is None or fdr_col is None or dpsi_col is None:
continue
if NONCODING_SET is not None:
df = df[df[gcol].astype(str).isin(NONCODING_SET)]
if df.empty:
continue
df[fdr_col] = pd.to_numeric(df[fdr_col], errors="coerce")
df[dpsi_col]= pd.to_numeric(df[dpsi_col], errors="coerce")
df = df[(df[fdr_col] <= alpha) & (df[dpsi_col].abs() >= min_abs_dpsi)]
if df.empty:
continue
colnames = list(df.columns)
df["__event_key"] = df.apply(lambda r: _event_key_from_row(evt, r, colnames), axis=1)
df["gene"] = df[gcol].astype(str)
df["__evt_type"] = evt
df["__variant"] = variant
parts.append(df[["gene","__event_key","__evt_type","__variant",fdr_col,dpsi_col]].rename(
columns={fdr_col:"fdr", dpsi_col:"delta_psi"}))
if not parts:
return pd.DataFrame()
sig = pd.concat(parts, ignore_index=True)
sig["__var_rank"] = sig["__variant"].map({"JCEC":0,"JC":1}).fillna(2)
sig = (sig.sort_values(["__event_key","__var_rank","fdr"])
.drop_duplicates("__event_key")
.drop(columns="__var_rank"))
return sig
# ------------- MntJULiP: detect group columns -------------
def detect_group_columns(df):
"""
Heuristic detector for MntJULiP group_id / group_structure / group_loc / group_q.
"""
def pick(regex, numeric=None, prefer_vals=None):
cand = [c for c in df.columns if re.search(regex, c.lower())]
if prefer_vals is not None:
scored = []
for c in cand:
vals = df[c].dropna().astype(str).head(200)
hit = any(v in {"o","i"} for v in vals)
scored.append((2 if hit else 1, c))
cand = [c for _, c in sorted(scored, reverse=True)]
if numeric is True:
cand = [c for c in cand if pd.api.types.is_numeric_dtype(df[c])]
if numeric is False:
cand = [c for c in cand if not pd.api.types.is_numeric_dtype(df[c])]
return cand[0] if cand else None
gid = pick(r"(group.*id|^gid$|cluster)", numeric=False)
gstruct = pick(r"(group.*structure|structure|orient|^group$)", prefer_vals={"o","i"})
gloc = pick(r"(group.*loc|^loc$|locus|anchor|junction.*pos)", numeric=True)
gq = pick(r"(group.*q|posterior|post_prob|prob|score)$", numeric=True)
if gid is None:
for c in df.columns:
s = df[c].dropna().astype(str).head(200)
if s.str.match(r"g\d+").any():
gid = c
break
if gstruct is None:
for c in df.columns:
s = df[c].dropna().astype(str).head(200)
if s.isin(["o","i"]).any():
gstruct = c
break
if gloc is None:
for c in df.columns:
if pd.api.types.is_integer_dtype(df[c]):
gloc = c
break
if gq is None:
for c in df.columns:
if pd.api.types.is_float_dtype(df[c]):
gq = c
break
return gid, gstruct, gloc, gq
def _read_mnt(zf, info):
with zf.open(info) as f:
if info.filename.endswith(".gz"):
return pd.read_csv(gzip.open(io.BytesIO(f.read()), "rt"), sep="\t")
return pd.read_csv(f, sep="\t")
def load_mnt_sig_introns_union(zip_path, alpha=0.05, min_abs_dpsi=0.20):
"""
UNION of DSR + DSA:
- Start from diff_introns.txt (DSA) → full intron set, with q.
- Annotate with group_id / group_loc / group_q from diff_spliced_introns.txt (DSR) when available.
- MntJULiP_source = "DSR" if group_id present, else "DSA".
- Apply q-value and (if available) |delta_psi| filter.
- Apply noncoding gene filter via NONCODING_SET.
"""
with zipfile.ZipFile(zip_path, 'r') as z:
p_dsr = None
p_dsa = None
for info in z.infolist():
name = info.filename
if name.endswith("diff_spliced_introns.txt"):
p_dsr = info
elif name.endswith("diff_introns.txt"):
p_dsa = info
if p_dsa is None:
return pd.DataFrame()
# --- DSA (diff_introns) ---
dall = _read_mnt(z, p_dsa)
dall.columns = [c.strip() for c in dall.columns]
chr_a = _find_first(dall.columns, "chrom","chr","chromosome")
st_a = _find_first(dall.columns, "start","loc_start","s")
en_a = _find_first(dall.columns, "end","loc_end","e")
sd_a = _find_first(dall.columns, "strand")
q_a = _find_first(dall.columns, "q_value","qvalue","q","FDR","fdr")
g_a = _find_first(dall.columns, "gene_name","gene","Gene")
if not all([chr_a, st_a, en_a, sd_a, q_a]):
return pd.DataFrame()
dall_proc = dall[[chr_a, st_a, en_a, sd_a, q_a] + ([g_a] if g_a else [])].copy()
dall_proc = dall_proc.rename(columns={
chr_a:"chrom", st_a:"start", en_a:"end", sd_a:"strand", q_a:"MntJULiP_q", (g_a or "gene"):"gene_raw"
})
dall_proc["chrom"] = dall_proc["chrom"].astype(str).str.replace("chr","")
dall_proc["start"] = pd.to_numeric(dall_proc["start"], errors="coerce")
dall_proc["end"] = pd.to_numeric(dall_proc["end"], errors="coerce")
dall_proc["strand"] = dall_proc["strand"].astype(str)
if g_a:
dall_proc["gene_dsa"] = dall_proc["gene_raw"].astype(str)
else:
dall_proc["gene_dsa"] = np.nan
dpsi_dsa_col = _find_first(dall.columns, "delta_psi","dpsi","deltaPSI","IncLevelDifference")
if dpsi_dsa_col:
dall_proc["delta_psi_dsa"] = pd.to_numeric(dall[dpsi_dsa_col], errors="coerce")
else:
dall_proc["delta_psi_dsa"] = np.nan
# --- DSR (diff_spliced_introns) ---
intr_proc = None
if p_dsr is not None:
intr = _read_mnt(z, p_dsr)
intr.columns = [c.strip() for c in intr.columns]
chr_i = _find_first(intr.columns, "chrom","chr","chromosome")
st_i = _find_first(intr.columns, "start","loc_start","s")
en_i = _find_first(intr.columns, "end","loc_end","e")
sd_i = _find_first(intr.columns, "strand")
g_i = _find_first(intr.columns, "gene_name","gene","Gene")
if all([chr_i, st_i, en_i, sd_i, g_i]):
intr_proc = intr[[chr_i, st_i, en_i, sd_i, g_i]].copy()
intr_proc = intr_proc.rename(columns={
chr_i:"chrom", st_i:"start", en_i:"end", sd_i:"strand", g_i:"gene_dsr"
})
intr_proc["chrom"] = intr_proc["chrom"].astype(str).str.replace("chr","")
intr_proc["start"] = pd.to_numeric(intr_proc["start"], errors="coerce")
intr_proc["end"] = pd.to_numeric(intr_proc["end"], errors="coerce")
intr_proc["strand"] = intr_proc["strand"].astype(str)
# detect group columns on full intr
gid_col, gstruct_col, gloc_col, gq_col = detect_group_columns(intr)
if gid_col is not None:
intr_proc["group_id"] = intr[gid_col].astype(str)
else:
intr_proc["group_id"] = np.nan
if gstruct_col is not None:
intr_proc["group_structure"] = intr[gstruct_col].astype(str)
else:
intr_proc["group_structure"] = np.nan
if gloc_col is not None:
intr_proc["group_loc"] = pd.to_numeric(intr[gloc_col], errors="coerce")
else:
intr_proc["group_loc"] = np.nan
if gq_col is not None:
intr_proc["group_q"] = pd.to_numeric(intr[gq_col], errors="coerce")
else:
intr_proc["group_q"] = np.nan
dpsi_dsr_col = _find_first(intr.columns, "delta_psi","dpsi","deltaPSI","IncLevelDifference")
if dpsi_dsr_col:
intr_proc["delta_psi_dsr"] = pd.to_numeric(intr[dpsi_dsr_col], errors="coerce")
else:
intr_proc["delta_psi_dsr"] = np.nan
# --- merge DSA base with DSR annotations ---
if intr_proc is not None:
m = dall_proc.merge(
intr_proc,
on=["chrom","start","end","strand"],
how="left",
suffixes=("","_dsr")
)
else:
m = dall_proc.copy()
m["gene_dsr"] = np.nan
m["group_id"] = np.nan
m["group_structure"] = np.nan
m["group_loc"] = np.nan
m["group_q"] = np.nan
m["delta_psi_dsr"] = np.nan
# unify gene
m["gene"] = m["gene_dsr"].fillna(m["gene_dsa"])
# unify delta_psi
m["delta_psi"] = m["delta_psi_dsr"]
m.loc[m["delta_psi"].isna(), "delta_psi"] = m.loc[m["delta_psi"].isna(), "delta_psi_dsa"]
# MntJULiP_source: DSR if group_id present, else DSA
m["MntJULiP_source"] = np.where(m["group_id"].notna(), "DSR", "DSA")
# q-value filter
m["MntJULiP_q"] = pd.to_numeric(m["MntJULiP_q"], errors="coerce")
m = m[m["MntJULiP_q"] <= alpha]
# optional delta_psi filter
if "delta_psi" in m.columns:
m = m[m["delta_psi"].abs() >= min_abs_dpsi]
# noncoding filter
if NONCODING_SET is not None:
m = m[m["gene"].astype(str).isin(NONCODING_SET)]
m = m.dropna(subset=["start","end"]).copy()
m["chrom"] = m["chrom"].astype(str)
m["strand"] = m["strand"].astype(str)
m["start"] = m["start"].astype(int)
m["end"] = m["end"].astype(int)
m["src"] = m["MntJULiP_source"]
m["__key"] = (
m["chrom"].astype(str) + ":" +
m["start"].astype(int).astype(str) + "-" +
m["end"].astype(int).astype(str) + ":" +
m["strand"].astype(str)
)
return m.drop_duplicates(["gene","__key"])
# ---------- Rank genes by frequency (events per gene) ----------
rmats_all_df = load_all_rmats_events(RMATS_ZIP)
mnt_all_df = load_mnt_sig_introns_union(MNT_ZIP, ALPHA_FDR, MIN_ABS_DPSI)
print(f"[INFO] rMATS events: {len(rmats_all_df)}, MntJULiP introns (DSR+DSA): {len(mnt_all_df)}")
master = pd.concat([
rmats_all_df.assign(source="rMATS", __event_key=rmats_all_df["__event_key"]),
mnt_all_df.assign(source="MntJULiP", __event_key=mnt_all_df["__key"])
], ignore_index=True)
master = master.dropna(subset=["gene"]).drop_duplicates(["gene","__event_key","source"], keep="first")
freq = (master.groupby("gene")
.agg(unique_events=("__event_key","nunique"),
total_rows=("__event_key","size"))
.reset_index()
.sort_values(["unique_events","total_rows","gene"], ascending=[False,False,True])
.reset_index(drop=True))
freq["splicing_rank"] = np.arange(1, len(freq)+1)
freq.to_csv(OUT_FREQ_RANK, index=False)
# classify regulation mode
rmats_sig = load_sig_rmats_events(RMATS_ZIP, ALPHA_FDR, MIN_ABS_DPSI)
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"])
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"])
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"))
reg["rmats_sig_unique"] = reg["rmats_sig_unique"].fillna(0).astype(int)
reg["mntjulip_any_sig_unique"] = reg["mntjulip_any_sig_unique"].fillna(0).astype(int)
def _mode(row):
r = row["rmats_sig_unique"] > 0
m = row["mntjulip_any_sig_unique"] > 0
if r and m:
return "both"
if r and not m:
return "exon_level_only"
if (not r) and m:
return "intron_only"
return "neither"
reg["regulation_mode"] = reg.apply(_mode, axis=1)
reg = reg.sort_values(["unique_events","total_rows","gene"], ascending=[False,False,True]).reset_index(drop=True)
reg.to_csv(OUT_REGMODES, index=False)
# ---------- pick top N + force include list ----------
top_auto = reg.head(TOP_PICK_N)["gene"].tolist()
selected_genes = list(dict.fromkeys(top_auto + [g for g in FORCE_INCLUDE if g in reg["gene"].values]))
print("[INFO] Selected genes:", selected_genes)
# ---------- rMATS per-event junction segments ----------
def rmats_junction_segments(evt, row):
segs = []
gi = _gi
if evt == "SE":
uEE = gi(row.get("upstreamEE")); exS = gi(row.get("exonStart_0base"))
exE = gi(row.get("exonEnd")); dES = gi(row.get("downstreamES"))
if uEE is not None and exS is not None: segs.append((min(uEE, exS), max(uEE, exS)))
if exE is not None and dES is not None: segs.append((min(exE, dES), max(exE, dES)))
elif evt == "RI":
uEE = gi(row.get("upstreamEE")); dES = gi(row.get("downstreamES"))
if uEE is not None and dES is not None: segs.append((min(uEE, dES), max(uEE, dES)))
elif evt in ("A5SS","A3SS"):
lS = gi(row.get("longExonStart_0base")); lE = gi(row.get("longExonEnd"))
sS = gi(row.get("shortES")); sE = gi(row.get("shortEE"))
fS = gi(row.get("flankingES")); fE = gi(row.get("flankingEE"))
if fE is not None and lS is not None: segs.append((min(fE,lS), max(fE,lS)))
if lE is not None and fS is not None: segs.append((min(lE,fS), max(lE,fS)))
if fE is not None and sS is not None: segs.append((min(fE,sS), max(fE,sS)))
if sE is not None and fS is not None: segs.append((min(sE,fS), max(sE,fS)))
elif evt == "MXE":
uEE = gi(row.get("upstreamEE")); e1S = gi(row.get("1stExonStart_0base")); e1E = gi(row.get("1stExonEnd"))
e2S = gi(row.get("2ndExonStart_0base")); e2E = gi(row.get("2ndExonEnd")); dES = gi(row.get("downstreamES"))
if uEE is not None and e1S is not None: segs.append((min(uEE,e1S), max(uEE,e1S)))
if e1E is not None and dES is not None: segs.append((min(e1E,dES), max(e1E,dES)))
if uEE is not None and e2S is not None: segs.append((min(uEE,e2S), max(uEE,e2S)))
if e2E is not None and dES is not None: segs.append((min(e2E,dES), max(e2E,dES)))
return sorted({(s,e) for (s,e) in segs if s is not None and e is not None and e > s})
# ---------- GENCODE intron/exon catalog ----------
def build_gencode_catalog_for_gene(gtf_zip, gene_name):
exons = []
with _open_gtf_any_from_zip(gtf_zip) as f:
for line in f:
if not line or line.startswith("#"):
continue
parts = line.rstrip("\n").split("\t")
if len(parts) < 9:
continue
chrom, source, feature, start, end, score, strand, frame, attrs = parts
if feature != "exon":
continue
d = {m.group(1): m.group(2) for m in re.finditer(r'(\S+)\s+"(.*?)";', attrs)}
gname = d.get("gene_name") or d.get("gene_id","")
if gname != gene_name and gene_name not in gname:
continue
tx = d.get("transcript_id","")
exons.append({"chr": ("chr"+chrom if not str(chrom).startswith("chr") else chrom),
"start": int(start), "end": int(end), "strand": strand, "transcript_id": tx})
exons_df = pd.DataFrame(exons)
introns = []
for tx, sub in exons_df.groupby("transcript_id"):
sub = sub.sort_values("start")
for (_, e1), (_, e2) in zip(sub.iloc[:-1].iterrows(), sub.iloc[1:].iterrows()):
if e2["start"] > e1["end"]:
introns.append({"chr": e1["chr"], "start": e1["end"], "end": e2["start"], "strand": e1["strand"], "transcript_id": tx})
introns_df = pd.DataFrame(introns)
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 {}
exon_index = {(c,s): sub[["start","end","transcript_id"]].to_records(index=False)
for (c,s), sub in exons_df.groupby(["chr","strand"])}
return exons_df, introns_df, exon_index, intron_index
def gencode_intron_match(ch, st, a, b, intron_index, slop=0):
key = (ch, st)
recs = intron_index.get(key, [])
hits = []
for (s,e,tx) in recs:
if (s - slop) <= a <= (s + slop) and (e - slop) <= b <= (e + slop):
hits.append(tx)
return hits
# ---------- Snaptron (GTEx tissue-filtered, TCGA all) ----------
def snaptron_query_counts(base_url, chrom, start, end, tissue=None, timeout=25):
region = f"{chrom}:{int(start)}-{int(end)}"
params = {"regions": region, "dtype": "s"}
if tissue:
params["sqs"] = f'tissue:"{tissue}"'
try:
r = requests.get(base_url, params=params, timeout=timeout)
r.raise_for_status()
text = r.text.strip()
if not text:
return {}
per_sample = {}
for line in text.splitlines():
parts = line.rstrip("\n").split("\t")
if len(parts) < 2:
continue
sid, val = parts[0], parts[1]
try:
c = float(val); per_sample[sid] = c
except Exception:
continue
return per_sample
except Exception:
return {}
def summarize_counts(d):
if not d:
return 0, 0.0
vals = [v for v in d.values() if v > 0]
if not vals:
return 0, 0.0
return len(vals), float(sum(vals)/len(vals))
# ---------- Anchor helpers (GAS5-style) ----------
def anchor_positions(evt, row):
"""
Consistent with your GAS5 script:
- RI : upstream exon START (upES)
- MXE : 1stExonStart_0base, 2ndExonStart_0base
- A5SS: shortES
"""
if evt == "RI":
upES = _gi(row.get("upstreamES"))
return [upES] if upES is not None else []
if evt == "MXE":
s1 = _gi(row.get("1stExonStart_0base"))
s2 = _gi(row.get("2ndExonStart_0base"))
return [x for x in [s1, s2] if x is not None]
if evt == "A5SS":
s = _gi(row.get("shortES"))
return [s] if s is not None else []
return []
def coloc_to_anchors(group_loc, anchors, tol=1):
if pd.isna(group_loc) or not anchors:
return False
try:
gl = int(group_loc)
except Exception:
return False
return any(abs(gl - a) <= tol for a in anchors)
def min_anchor_dist(group_loc, anchors):
if pd.isna(group_loc) or not anchors:
return np.inf
try:
gl = int(group_loc)
except Exception:
return np.inf
return min(abs(gl - a) for a in anchors)
def max_overlap_with_segments(segs, start, end):
best = 0
for (s, e) in segs:
best = max(best, _span_overlap((s, e), (start, end)))
return best
# ---------- Per-gene core (generalized GAS5 logic) ----------
def per_gene_pipeline(gene, rmats_all_df, rmats_sig, mnt_all_df):
r_any = rmats_all_df[rmats_all_df["gene"] == gene].copy()
if r_any.empty:
return None
chr_col = next((c for c in ["chr","chrom","Chromosome"] if c in r_any.columns), None)
sd_col = "strand" if "strand" in r_any.columns else None
if chr_col is None or sd_col is None:
return None
r_any["chrom"] = r_any[chr_col].astype(str).map(_normalize_chr)
r_any["strand"] = r_any[sd_col].astype(str)
coords = []
for _, row in r_any.iterrows():
for s,e in rmats_junction_segments(row["__evt_type"], row):
coords += [s,e]
if not coords:
return None
if rmats_sig.empty:
return None
coord_cols = [c for c in ["chr","Chromosome","chrom","strand",
"exonStart_0base","exonEnd","upstreamES","upstreamEE","downstreamES","downstreamEE",
"longExonStart_0base","longExonEnd","flankingES","flankingEE","shortES","shortEE",
"1stExonStart_0base","1stExonEnd","2ndExonStart_0base","2ndExonEnd"]
if c in rmats_all_df.columns]
r_coords = rmats_all_df[["__event_key","gene","__evt_type","__variant",*coord_cols]].drop_duplicates("__event_key")
r_sig = rmats_sig.merge(r_coords, on=["__event_key","gene","__evt_type"], how="left")
r_sig = r_sig[(r_sig["gene"] == gene)].copy()
if r_sig.empty:
return None
m = mnt_all_df[mnt_all_df["gene"] == gene].copy()
if m.empty:
return None
m["chrom"] = m["chrom"].astype(str)
m["strand"] = m["strand"].astype(str)
r_sig["chrom"] = r_sig[coord_cols[0]].astype(str).map(_normalize_chr) if coord_cols else ""
r_sig["strand"] = r_sig["strand"].astype(str)
pub_rows = []
for _, row in r_sig.iterrows():
evt = row["__evt_type"]
segs = rmats_junction_segments(evt, row)
if not segs:
continue
evctr = float(np.mean([(s+e)/2.0 for (s,e) in segs]))
anchors = anchor_positions(evt, row)
msub = m[(m["chrom"] == row["chrom"]) & (m["strand"] == row["strand"])].copy()
if msub.empty:
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, "MntJULiP_intron_end": np.nan,
"MntJULiP_q": np.nan, "Overlap_bp": 0, "Overlap_type": "None",
"MntJULiP_group_id": np.nan, "MntJULiP_group_structure": np.nan,
"MntJULiP_group_loc": np.nan, "MntJULiP_group_q": np.nan,
"MntJULiP_source": "None",
"note": "No overlapping significant intron"
})
continue
msub["Overlap_bp"] = msub.apply(
lambda r: max_overlap_with_segments(segs, int(r["start"]), int(r["end"])),
axis=1
)
msub = msub[msub["Overlap_bp"] > 0].copy()
if msub.empty:
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, "MntJULiP_intron_end": np.nan,
"MntJULiP_q": np.nan, "Overlap_bp": 0, "Overlap_type": "None",
"MntJULiP_group_id": np.nan, "MntJULiP_group_structure": np.nan,
"MntJULiP_group_loc": np.nan, "MntJULiP_group_q": np.nan,
"MntJULiP_source": "None",
"note": "No overlapping significant intron"
})
continue
# Anchor-based tie-breaking, GAS5-style
if anchors and "group_loc" in msub.columns:
msub["anchor_dist"] = msub["group_loc"].apply(lambda gl: min_anchor_dist(gl, anchors))
msub["CoLoc_ANCHOR"] = msub["group_loc"].apply(lambda gl: coloc_to_anchors(gl, anchors, tol=COLOC_TOL))
else:
msub["anchor_dist"] = np.inf
msub["CoLoc_ANCHOR"] = False
msub["pref_o"] = (msub.get("group_structure","").astype(str) == "o").astype(int)
msub["q_sort"] = pd.to_numeric(msub.get("MntJULiP_q", 1.0), errors="coerce").fillna(1.0)
msub["ctr_diff"] = (msub[["start","end"]].mean(axis=1) - evctr).abs()
if msub["CoLoc_ANCHOR"].any():
c_sel = msub[msub["CoLoc_ANCHOR"]].copy()
else:
c_sel = msub
c_sel = c_sel.sort_values(
["anchor_dist","pref_o","q_sort","Overlap_bp","ctr_diff","start","end"],
ascending=[True, False, True, False, True, True, True]
)
best = c_sel.iloc[0]
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": int(best["start"]),
"MntJULiP_intron_end": int(best["end"]),
"MntJULiP_q": float(best["MntJULiP_q"]) if pd.notna(best["MntJULiP_q"]) else np.nan,
"Overlap_bp": int(best["Overlap_bp"]),
"Overlap_type": "Strong" if int(best["Overlap_bp"]) >= 20 else "Partial",
"MntJULiP_group_id": best.get("group_id", np.nan),
"MntJULiP_group_structure": best.get("group_structure", np.nan),
"MntJULiP_group_loc": best.get("group_loc", np.nan),
"MntJULiP_group_q": best.get("group_q", np.nan),
"MntJULiP_source": best.get("MntJULiP_source", "DSA/DSR"),
"note": "Matched MntJULiP significant intron"
})
pub = pd.DataFrame(pub_rows).drop_duplicates(["gene","rMATS_event_key"])
if pub.empty:
return None
# -------- GENCODE structural support --------
exons_df, introns_df, exon_index, intron_index = build_gencode_catalog_for_gene(GTF_ZIP, gene)
def gencode_check(row):
if pd.isna(row["MntJULiP_intron_start"]) or pd.isna(row["MntJULiP_intron_end"]):
return "No", "no_match", ""
ch = row["chr"]
st = row["strand"]
a = int(row["MntJULiP_intron_start"])
b = int(row["MntJULiP_intron_end"])
hit = gencode_intron_match(ch, st, a, b, intron_index, slop=GENCODE_SLOP)
if hit:
return "Yes", ("exact_intron±%d" % GENCODE_SLOP), ";".join(sorted(set(hit)))
exrecs = exon_index.get((ch, st), [])
exon_hits = []
for (s,e,tx) in exrecs:
if s <= a and e >= b:
exon_hits.append(tx)
if exon_hits:
return "Yes", "retained_intron_exon", ";".join(sorted(set(exon_hits)))
return "No", "no_match", ""
g_ann = pub.apply(gencode_check, axis=1, result_type="expand")
pub[["annotated_in_gencode","gencode_match_type","gencode_transcripts"]] = g_ann
# -------- Snaptron presence (GTEx EM/GEJ, TCGA) --------
def cohort_presence(row):
ch = row["chr"]; a = row["MntJULiP_intron_start"]; b = row["MntJULiP_intron_end"]
if pd.isna(a) or pd.isna(b):
return pd.Series({
"GTEx_EM_present":"No","GTEx_EM_n":0,"GTEx_EM_mean":0.0,
"GTEx_GEJ_present":"No","GTEx_GEJ_n":0,"GTEx_GEJ_mean":0.0,
"TCGA_present":"No","TCGA_n":0,"TCGA_mean":0.0
})
em = snaptron_query_counts(SNAP_GTX, ch, a, b, tissue=GTEX_TISSUES["EM"])
em_n, em_mean = summarize_counts(em)
gej = snaptron_query_counts(SNAP_GTX, ch, a, b, tissue=GTEX_TISSUES["GEJ"])
gej_n, gej_mean = summarize_counts(gej)
tc = snaptron_query_counts(SNAP_TCGA, ch, a, b, tissue=None)
tc_n, tc_mean = summarize_counts(tc)
return pd.Series({
"GTEx_EM_present":"Yes" if em_n>0 else "No","GTEx_EM_n":int(em_n),"GTEx_EM_mean":float(em_mean),
"GTEx_GEJ_present":"Yes" if gej_n>0 else "No","GTEx_GEJ_n":int(gej_n),"GTEx_GEJ_mean":float(gej_mean),
"TCGA_present":"Yes" if tc_n>0 else "No","TCGA_n":int(tc_n),"TCGA_mean":float(tc_mean)
})
cohort_df = pub.apply(cohort_presence, axis=1)
out = pd.concat([pub, cohort_df], axis=1)
# write per gene
out_pub = f"/content/{gene}_PublicationView.csv"
out_annot = f"/content/{gene}_With_Annotation_and_CohortPresence.csv"
out.to_csv(out_annot, index=False)
pub.to_csv(out_pub, index=False)
return out
# -------- run per-gene on the chosen genes --------
combined = []
for g in selected_genes:
print(f"[RUN] {g}")
res = per_gene_pipeline(g, rmats_all_df, rmats_sig, mnt_all_df)
if res is not None and not res.empty:
combined.append(res)
if combined:
comb = pd.concat(combined, ignore_index=True)
comb.to_csv(OUT_COMBINED_ANNOT, index=False)
print(f"[OK] Combined table → {OUT_COMBINED_ANNOT}")
else:
print("[WARN] No per-gene outputs were generated. Check inputs/filters.")
CODE CELL -10 CREATING RMATS SASHIMI PLOT¶
the event files to create rmats sashimi plots of the significant events were generated using the following code
# New Section
# This is formatted as code
import os
from pathlib import Path
import pandas as pd
import numpy as np
# ---------- USER PATHS ----------
TABLE3 = "/Table3_TopGenes_With_Annotation_CohortPresence.csv"
RMATS_DIR = Path("/mnt/gs/rmats/rmats_output_full") # directory with *.MATS.*.txt
OUT_DIR = Path("/content/Table3_rMATS_Events_A5SSFIX")
OUT_DIR.mkdir(parents=True, exist_ok=True)
print(f"Using Table3: {TABLE3}")
print(f"Using rMATS: {RMATS_DIR}")
print(f"Output: {OUT_DIR}")
# ---------- HELPERS ----------
def _find_first(cols, *cands):
low = {c.lower(): c for c in cols}
for c in cands:
if c.lower() in low:
return low[c.lower()]
return None
def _gi(v):
"""Safe int converter used when building keys."""
try:
if v is None or (isinstance(v, float) and np.isnan(v)):
return None
s = str(v).strip()
if s == "" or s.lower() == "nan":
return None
return int(float(s))
except Exception:
return None
def build_event_key(evt, row, cols):
"""
Rebuild the SAME rMATS_event_key used in the main pipeline that generated Table 3.
IMPORTANT:
- For A3SS/A5SS we MUST use:
[longExonStart_0base, longExonEnd,
flankingES, flankingEE,
shortES, shortEE]
because that is what was used in _event_key_from_row earlier.
- MXE/SE/RI orders are unchanged.
"""
chr_col = _find_first(cols, "chr", "chrom", "chromosome")
strand_col = _find_first(cols, "strand")
chrom = str(row[chr_col]).replace("chr", "").strip()
strand = str(row[strand_col]).strip()
def gi(colname):
if colname not in cols:
return ""
g = _gi(row[colname])
return "" if g is None else str(g)
evt = evt.upper()
parts = [chrom, strand]
if evt == "SE":
for c in [
"exonStart_0base", "exonEnd",
"upstreamES", "upstreamEE",
"downstreamES", "downstreamEE",
]:
parts.append(gi(c))
elif evt == "RI":
# RI: exonStart_0base/exonEnd + upstream/downstream anchors
for c in [
"exonStart_0base", "exonEnd",
"upstreamES", "upstreamEE",
"downstreamES", "downstreamEE",
]:
parts.append(gi(c))
elif evt in ("A3SS", "A5SS"):
# 🔴 CRITICAL FIX: match the ORIGINAL pipeline:
# [longExonStart_0base, longExonEnd,
# flankingES, flankingEE,
# shortES, shortEE]
for c in [
"longExonStart_0base", "longExonEnd",
"flankingES", "flankingEE",
"shortES", "shortEE",
]:
parts.append(gi(c))
elif evt == "MXE":
for c in [
"1stExonStart_0base", "1stExonEnd",
"2ndExonStart_0base", "2ndExonEnd",
"upstreamES", "upstreamEE",
"downstreamES", "downstreamEE",
]:
parts.append(gi(c))
parts.append(evt)
return "_".join(parts)
# coordinate columns that MUST NOT be empty for each event type
REQ_COORDS = {
"SE": [
"exonStart_0base", "exonEnd",
"upstreamES", "upstreamEE",
"downstreamES", "downstreamEE",
],
"RI": [
"exonStart_0base", "exonEnd",
"upstreamES", "upstreamEE",
"downstreamES", "downstreamEE",
],
"A5SS": [
"longExonStart_0base", "longExonEnd",
"flankingES", "flankingEE",
"shortES", "shortEE",
],
"A3SS": [
"longExonStart_0base", "longExonEnd",
"flankingES", "flankingEE",
"shortES", "shortEE",
],
"MXE": [
"1stExonStart_0base", "1stExonEnd",
"2ndExonStart_0base", "2ndExonEnd",
"upstreamES", "upstreamEE",
"downstreamES", "downstreamEE",
],
}
# ---------- LOAD TABLE 3 ----------
t3 = pd.read_csv(TABLE3, dtype=str)
if "event_type" not in t3.columns or "rMATS_event_key" not in t3.columns:
raise SystemExit("Table 3 must contain 'event_type' and 'rMATS_event_key' columns.")
t3_ok = t3[
t3["event_type"].notna() &
t3["rMATS_event_key"].notna()
].copy()
print(f"Table3 usable rows: {len(t3_ok)}")
# ---------- LOAD rMATS MATS FILES & BUILD INDEX ----------
rmats_index = {} # event_key -> (evt, variant, src_fname, row)
for evt in ["SE", "RI", "A3SS", "A5SS", "MXE"]:
for variant in ["JC", "JCEC"]:
fname = f"{evt}.MATS.{variant}.txt"
path = RMATS_DIR / fname
if not path.exists():
continue
print(f"Loaded: {fname}")
df = pd.read_csv(path, sep="\t", dtype=str)
df_cols = list(df.columns)
chr_col = _find_first(df_cols, "chr", "chrom", "chromosome")
strand_col = _find_first(df_cols, "strand")
if chr_col is None or strand_col is None:
print(f" ⚠ Skipping {fname}: no chr/strand columns detected")
continue
# Build event keys exactly as in Table 3 pipeline
keys = []
for _, r in df.iterrows():
keys.append(build_event_key(evt, r, df_cols))
df["__event_key"] = keys
# Store best representative per key (prefer JCEC over JC)
for _, r in df.iterrows():
key = r["__event_key"]
if not key:
continue
prev = rmats_index.get(key)
if prev is None:
rmats_index[key] = (evt, variant, fname, r)
else:
_, prev_variant, _, _ = prev
if prev_variant == "JCEC":
continue
if variant == "JCEC":
rmats_index[key] = (evt, variant, fname, r)
print(f"Indexed events: {len(rmats_index)}")
# ---------- WRITE EVENT FILES ----------
written = 0
missing = 0
for _, row in t3_ok.iterrows():
gene = str(row.get("gene", row.get("geneSymbol", "GENE"))).strip()
evt_tbl = str(row["event_type"]).upper().strip()
evkey = str(row["rMATS_event_key"]).strip()
hit = rmats_index.get(evkey)
if hit is None:
# Special message for RI (we accept this as MntJULiP-only)
if evt_tbl == "RI":
print(f"❗ RI skipped (MntJULiP-only, rMATS has no matching RI): {evkey}")
else:
print(f"❌ Missing rMATS: {evkey}")
missing += 1
continue
evt0, variant, src_fname, mrow = hit
df_evt = evt0 # event type from rMATS (should agree with evt_tbl)
cols = list(mrow.index)
# Ensure required coordinate columns are not empty, or EventCoor will break
req = REQ_COORDS.get(df_evt, [])
mrow = mrow.copy()
for c in req:
if c in cols:
val = str(mrow.get(c, "")).strip()
if val == "" or val.lower() == "nan":
# fallback to -1 so int() succeeds inside EventCoor
mrow[c] = "-1"
# filename uses the Table3 rMATS_event_key for traceability
fname = f"{gene}_{evkey}.event.txt"
out_path = OUT_DIR / fname
with open(out_path, "w") as f:
f.write("\t".join(cols) + "\n")
f.write("\t".join(str(mrow.get(c, "")) for c in cols) + "\n")
print(f"[OK] {fname} (from {src_fname}, {variant})")
written += 1
print("\nDONE.")
print(f"Written: {written}")
print(f"Missing: {missing}")
print(f"Output: {OUT_DIR}")
Once the event files were generated they were fed in to rmats sashimi to produce the plots using the following code . It was necessary to use python 3.10 here for compatibility for running rmats sashimiplot. The rMATS sashimi plots shown in this study were generated using the official rmats2sashimiplot package (Xinglab), converted to full Python 3.10 compatibility to match the environment of Colab Enterprise pipeline and ensure reproducible visualization of splice‐junction usage for each significant event.
# This is formatted as code
%%bash
set -e
# --- STEP 1: INSTALL PYTHON 3.10 ---
echo "=== 1. INSTALLING PYTHON 3.10 ==="
sudo apt-get update -y > /dev/null 2>&1
sudo apt-get install -y python3.10 python3.10-dev python3.10-distutils > /dev/null 2>&1
curl -sS https://bootstrap.pypa.io/get-pip.py | python3.10 > /dev/null 2>&1
# --- STEP 2: INSTALL DEPENDENCIES ---
echo "=== 2. INSTALLING LIBRARIES ==="
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
# --- STEP 3: DOWNLOAD CODE ---
echo "=== 3. DOWNLOADING CODE ==="
cd /content
rm -rf rmats2sashimiplot_fixed
git clone https://github.com/Xinglab/rmats2sashimiplot.git rmats2sashimiplot_fixed > /dev/null 2>&1
cd rmats2sashimiplot_fixed
# --- STEP 4: CONVERT CODE TO PYTHON 3 ---
echo "=== 4. CONVERTING CODE TO PYTHON 3 ==="
cat << 'EOF' > convert.py
import sys
from lib2to3.main import main
sys.argv = ['2to3', '-w', '-n', '--no-diffs', 'src']
sys.exit(main("lib2to3.fixes"))
EOF
python3.10 convert.py
echo "✅ Code successfully converted."
# --- STEP 5: RUN THE PLOT ---
echo "=== 5. RUNNING SASHIMI PLOT ==="
B1="/fixed_b1.txt"
B2="/fixed_b2.txt"
OUT="/content/Finalplotgas5_mxe8"
EVENT_FILE="/content/Table3_rMATS_Events_A5SSFIX/GAS5_1_-_173866527_173866581_173866760_173866796_173866176_173866206_173866990_173867043_MXE.event.txt"
rm -rf "$OUT"
mkdir -p "$OUT"
export PYTHONPATH=$PYTHONPATH:/content/rmats2sashimiplot_fixed/src
python3.10 ./src/rmats2sashimiplot/rmats2sashimiplot.py \
-o "$OUT" \
--l1 Control \
--l2 Tumor \
--event-type MXE \
-e "$EVENT_FILE" \
--b1 "$B1" \
--b2 "$B2" \
--keep-event-chr-prefix \
--min-counts 0
echo ""
echo "=== RESULTS ==="
find "$OUT" -name "*.pdf"
CODE CELL-11- Creating final image files¶
code to create image -1
# This is formatted as code
import matplotlib.pyplot as plt
import matplotlib.patches as mpatches
# -----------------------------------------------------------
# FIXED ANCHOR FOR SHIFT
# -----------------------------------------------------------
ANCHOR = 173864075
def shift(x): return x - ANCHOR
# -----------------------------------------------------------
# MXE structure (each MXE has two parts → must stay same row)
# -----------------------------------------------------------
MXEs = {
1: [(173864256,173864506), (173865228,173865282)],
2: [(173864479,173864704), (173865228,173865282)],
3: [(173865509,173865547), (173864256,173864304)],
4: [(173865856,173865894), (173866760,173866796)],
5: [(173866176,173866206), (173865856,173865894)],
6: [(173866176,173866567), (173866760,173866796)],
7: [(173866527,173866567), (173865856,173866206)],
8: [(173866527,173866581), (173866176,173866206)]
}
# -----------------------------------------------------------
# Distinct MXE colours
# -----------------------------------------------------------
COLOR_MXE = [
"#7b61ff", "#ff6b6b", "#4ecdc4", "#ffa600",
"#1e90ff", "#ff1493", "#32cd32", "#8b4513"
]
COLOR_A5SS_LONG = "#3b75ff"
COLOR_A5SS_SHORT = "#38e5d2"
COLOR_RI = "#ff4a4a"
# -----------------------------------------------------------
# Assign one horizontal lane per MXE
# MXE-1 at y=1, MXE-2 at y=2, ... MXE-8 at y=8
# -----------------------------------------------------------
mxe_y = {m: m for m in MXEs} # simple and clean
plt.figure(figsize=(18,10))
ax = plt.gca()
# -----------------------------------------------------------
# Draw all MXEs cleanly, one lane each
# -----------------------------------------------------------
for mxe_id, parts in MXEs.items():
y = mxe_y[mxe_id]
color = COLOR_MXE[mxe_id-1]
for (s_raw, e_raw) in parts:
s, e = shift(s_raw), shift(e_raw)
ax.hlines(y, s, e, color=color, linewidth=12)
# -----------------------------------------------------------
# Draw A5SS + RI lanes below MXEs
# -----------------------------------------------------------
ax.hlines(-1, shift(173864256), shift(173865547), color=COLOR_A5SS_LONG, linewidth=12)
ax.hlines(-2, shift(173865282), shift(173865509), color=COLOR_A5SS_SHORT, linewidth=12)
ax.hlines(-3, shift(173865282), shift(173865470), color=COLOR_RI, linewidth=12)
# -----------------------------------------------------------
# Labels & formatting
# -----------------------------------------------------------
plt.title("GAS5 alternative splicing architecture\n"
"(All MXEs + A5SS + RI, shifted axis, MXE-1 to MXE-8)",
fontsize=20, pad=20)
plt.xlabel("Shifted genomic coordinate (chr1, GAS5 locus)")
plt.yticks([], [])
plt.text(0, -4, f"0 = chr1:{ANCHOR:,}", fontsize=14)
plt.axhline(0, color="black", linewidth=1)
# -----------------------------------------------------------
# Legend bottom right
# -----------------------------------------------------------
legend_items = []
for i,c in enumerate(COLOR_MXE,1):
legend_items.append(mpatches.Patch(color=c, label=f"MXE-{i}"))
legend_items.append(mpatches.Patch(color=COLOR_A5SS_LONG, label="A5SS-long"))
legend_items.append(mpatches.Patch(color=COLOR_A5SS_SHORT, label="A5SS-short"))
legend_items.append(mpatches.Patch(color=COLOR_RI, label="RI (MntJULiP-only)"))
plt.legend(handles=legend_items,
title="Event legend",
loc="lower right",
fontsize=12,
frameon=True)
plt.tight_layout()
plt.savefig("/content/GAS5_MXE_simple_clean.png", dpi=300)
plt.show()
print("Saved figure: /content/GAS5_MXE_simple_clean.png")
Code to create image -2
# This is formatted as code
# 1. INSTALL TOOLS
!apt-get install -y poppler-utils > /dev/null
!pip install pdf2image > /dev/null
import matplotlib.pyplot as plt
import matplotlib.patches as patches
from pdf2image import convert_from_path
import os
import string
# ==========================================
# CONFIGURATION
# ==========================================
# 1. FILES (Exact requested order)
# Row 1: pdf1, pdf2, pdfa5ss
# Row 2: pdf3, pdf4, pdf5
# Row 3: pdf6, pdf7, pdf8
pdf_files = [
'pdf1.pdf', # MXE 1
'pdf2.pdf', # MXE 2
'pdfa5ss.pdf', # A5SS
'pdf3.pdf', # MXE 3
'pdf4.pdf', # MXE 4
'pdf5.pdf', # MXE 5
'pdf6.pdf', # MXE 6
'pdf7.pdf', # MXE 7
'pdf8.pdf' # MXE 8
]
# 2. DETAILED CUSTOM TITLES (Highlighting the Differences)
custom_titles = [
# 1. MXE 1
"GAS5 MXE 1\nVar Exons: 173,864,256 | 173,865,228",
# 2. MXE 2
"GAS5 MXE 2\nVar Exons: 173,864,479 | 173,865,228",
# 3. A5SS
"GAS5 A5SS\nAlt Splice: 173,865,470 vs 173,865,509",
# 4. MXE 3
"GAS5 MXE 3\nVar Exons: 173,865,509 | 173,866,176",
# 5. MXE 4
"GAS5 MXE 4\nVar Exons: 173,865,856 | 173,866,760",
# 6. MXE 5
"GAS5 MXE 5\nVar Exon 1: 173,866,176 (30bp)",
# 7. MXE 6 (Note the longer exon)
"GAS5 MXE 6\nVar Exon 1: 173,866,176 (391bp)",
# 8. MXE 7
"GAS5 MXE 7\nVar Exon 1: 173,866,527 (40bp)",
# 9. MXE 8
"GAS5 MXE 8\nVar Exon 1: 173,866,527 (54bp)"
]
# 3. SETTINGS
HIDE_HEADER_HEIGHT = 220
# ==========================================
# GENERATION SCRIPT
# ==========================================
fig, axes = plt.subplots(3, 3, figsize=(24, 21))
axes = axes.flatten()
print("Generating detailed publication figure...")
for i, ax in enumerate(axes):
if i < len(pdf_files):
base_name = pdf_files[i].strip().replace("'", "")
# Smart Find
possible_paths = [base_name, f"/{base_name}", f"/content/{base_name}"]
found_path = None
for p in possible_paths:
p = p.replace("//", "/")
if os.path.exists(p):
found_path = p
break
if found_path:
try:
images = convert_from_path(found_path, dpi=300)
img = images[0]
width, height = img.size
ax.imshow(img)
# White-Out Old Title
rect = patches.Rectangle((0, 0), width, HIDE_HEADER_HEIGHT,
linewidth=0, edgecolor='white', facecolor='white')
ax.add_patch(rect)
# Write New Detailed Title
title_text = custom_titles[i] if i < len(custom_titles) else base_name
ax.text(width/2, HIDE_HEADER_HEIGHT/2, title_text,
ha='center', va='center', fontsize=18, fontweight='bold', color='black')
# Panel Label
label = string.ascii_uppercase[i]
ax.text(-0.05, 1.05, label, transform=ax.transAxes,
fontsize=32, fontweight='bold', va='top', ha='right')
ax.axis('off')
print(f"✅ Panel {label}: {title_text.splitlines()[0]}")
except Exception as e:
print(f"❌ Error {found_path}: {e}")
ax.axis('off')
else:
print(f"⚠️ MISSING: {base_name}")
ax.axis('off')
else:
ax.axis('off')
plt.tight_layout()
output_file = "GAS5_Sashimi_Publication_Detailed.pdf"
plt.savefig(output_file, dpi=300, bbox_inches='tight')
print(f"\nSUCCESS! Download below:")
from google.colab import files
files.download(output_file)