https://nf-co.re logo
Join Slack
Powered by
Hey Everyone, One of my current, technical issues...
# rnasplice
m
Hey Everyone, One of my current, technical issues when running 
rnasplice
is this: 
Error executing process > 'NFCORE_RNASPLICE:RNASPLICE:VISUALISE_MISO:MISO_SASHIMI (1)'.
Specifically, this message:
Copy code
Workflow execution completed unsuccessfully!
The exit status of the task that caused the workflow execution to fail was: 1.
The full error message was:
Error executing process > 'NFCORE_RNASPLICE:RNASPLICE:VISUALISE_MISO:MISO_SASHIMI (1)'

Caused by:
 Process `NFCORE_RNASPLICE:RNASPLICE:VISUALISE_MISO:MISO_SASHIMI (1)` terminated with an error exit status (1)


Command executed:

 sashimi_plot --plot-event ENSG00000004961 index miso_settings.txt --output-dir sashimi
 
 cat <<-END_VERSIONS > versions.yml
 "NFCORE_RNASPLICE:RNASPLICE:VISUALISE_MISO:MISO_SASHIMI":
     python: $(python --version | sed "s/Python //g")
     misopy: $(python -c "import pkg_resources; print(pkg_resources.get_distribution('misopy').version)")
 END_VERSIONS

Command exit status:
 1

Command output:
 (empty)

Command error:
 /usr/local/lib/python2.7/site-packages/matplotlib/cbook/deprecation.py:107: MatplotlibDeprecationWarning: The mpl_toolkits.axes_grid module was deprecated in version 2.1. Use mpl_toolkits.axes_grid1 and mpl_toolkits.axisartist provies the same functionality instead.
   warnings.warn(message, mplDeprecation, stacklevel=1)
 Traceback (most recent call last):
   File "/usr/local/bin/sashimi_plot", line 11, in 
     sys.exit(main())
   File "/usr/local/lib/python2.7/site-packages/misopy/sashimi_plot/sashimi_plot.py", line 276, in main
     plot_label=plot_label)
   File "/usr/local/lib/python2.7/site-packages/misopy/sashimi_plot/sashimi_plot.py", line 142, in plot_event
     %(event_name, pickle_dir)
 Exception: Event ENSG00000004961 not found in pickled directory index. Are you sure this is the right directory for the event?

Work dir:
 /home/tud03125/pipeline/normal_human_liver_rnasplice/work/be/c27262d2dc58177e582cb86e3ea716

Tip: you can replicate the issue by changing to the process work dir and entering the command `bash .command.run`
I wonder if anyone had this before? My PBS job script I’ve used is this:
Copy code
#!/bin/bash
#PBS -N rnasplice_pipeline
#PBS -l nodes=1:ppn=16
#PBS -l mem=128gb
#PBS -l walltime=168:00:00 # Set to 168 hours
#PBS -o /home/tud03125/pipeline/normal_human_liver_rnasplice/$PBS_JOBID.out
#PBS -e /home/tud03125/pipeline/normal_human_liver_rnasplice/$PBS_JOBID.err
#PBS -M michael.levin@temple.edu
#PBS -m abe # 'a' for abort, 'b' for begin, 'e' for end

# Debugging: Print environment information
echo "Job started on $(date)"
echo "Running on $(hostname)"
echo "Current directory: $(pwd)"

# Output directory
OUTPUT_DIR="/home/tud03125/pipeline/normal_human_liver_rnasplice"

# Create the output directory if it doesn't exist
mkdir -p "$OUTPUT_DIR"

# Load necessary modules
module load gcc/11.2.0
module load singularity
module load python/3.10.2 # Load Python 3.10.2

# Create an alias for python to point to python3
alias python=python3

# Check that the alias works
python --version # Should output Python 3.10.2

# Set Nextflow version
export NXF_VER=24.04.4

# Install Nextflow if not already installed
if [ ! -f ./nextflow ]; then
   echo "Nextflow not found, downloading..."
   curl -s <https://get.nextflow.io> | bash
   chmod +x nextflow
fi

# Set Nextflow work directory to avoid quota issues
export NXF_WORK=/home/tud03125/pipeline/normal_human_liver_rnasplice/work

# Run Nextflow with your local rnasplice pipeline
echo "Starting Nextflow pipeline..."
./nextflow run /home/tud03125/pipeline/rnasplice -resume -profile singularity \
   --input /home/tud03125/pipeline/normal_human_liver_design_rnasplice_3.csv \
   --fasta /home/tud03125/pipeline/Homo_sapiens.GRCh38.dna_sm.primary_assembly.fa \
   --gtf /home/tud03125/pipeline/Homo_sapiens.GRCh38.108.gtf \
   --save_reference \
   --outdir /home/tud03125/pipeline/normal_human_liver_rnasplice \
   --star_index /home/tud03125/STAR_index \
   --salmon_index /home/tud03125/new_salmon_index \
   --rmats \
   --suppa \
   --ribo_database_manifest /home/tud03125/ribo_database_manifest.txt \
   --contrasts /home/tud03125/pipeline/contrast_NORMAL_2.csv 

# Debugging: Print completion message and end timestamp
echo "Processing complete."
echo "Job ended on $(date)"
The 
miso_sashimi.nf
script from 
/home/tud03125/pipeline/rnasplice/modules/local/
directory is this:
Copy code
process MISO_SASHIMI {
   label 'process_single'

   conda "conda-forge::python=2.7 bioconda::misopy=0.5.4"
   container "${ workflow.containerEngine == 'singularity' && !task.ext.singularity_pull_docker_container ?
       '<https://depot.galaxyproject.org/singularity/misopy:0.5.4--py27h9801fc8_5>' :
       'biocontainers/misopy:0.5.4--py27h9801fc8_5' }"

   input:
   path index_path
   tuple path(miso_settings), val(miso_gene)
   path ("bam_files/*")  // Need bam and bai in working directory
   path ("miso_data/*")  // Need Miso files in working directory

   output:
   path "sashimi/*"          , emit: sashimi
   path "versions.yml"       , emit: versions

   when:
   task.ext.when == null || task.ext.when

   script:
   """
   sashimi_plot --plot-event $miso_gene $index_path $miso_settings --output-dir sashimi

   cat <<-END_VERSIONS > versions.yml
   "${task.process}":
       python: \$(python --version | sed "s/Python //g")
       misopy: \$(python -c "import pkg_resources; print(pkg_resources.get_distribution('misopy').version)")
   END_VERSIONS
   """

}
I’ve noticed that for 
misopy
, they require outdated Python version: 
2.7
. I definitely have samples like 
ENSG00000004961
, in fact 62,703 of them. It’s just, for some reason, they couldn’t identify it. Could that potentially create issues or is there a way to bypass 
Python 2.7
version while using 
misopy
? If not, then can it be skipped (ChatGPT suggested that if I use 
rMAT
or 
SUPPA
, which’s in my parameters, then I’d probably don’t need 
misopy
)?
t
Yes, this a common issue. The problem is that the default genes (it's looking for the exact version of the annotation) to be used for the Sashimi plot are not in the annotation you happen to be using. If you are interested in specific genes, pick them from your GTF and use them with 
--miso_genes
, otherwise use 
--sashimi_plot false
.
m
@Thomas Danhorn Thanks for pointing it out! And, I thought there’s something fishy about it. So, I'm planning to run this, instead:
Copy code
#!/bin/bash
#PBS -N rnasplice_pipeline
#PBS -l nodes=1:ppn=16
#PBS -l mem=128gb
#PBS -l walltime=168:00:00 # Set to 168 hours
#PBS -o /home/tud03125/pipeline/normal_human_liver_rnasplice/$PBS_JOBID.out
#PBS -e /home/tud03125/pipeline/normal_human_liver_rnasplice/$PBS_JOBID.err
#PBS -M <mailto:michael.levin@temple.edu|michael.levin@temple.edu>
#PBS -m abe # 'a' for abort, 'b' for begin, 'e' for end

# Debugging: Print environment information
echo "Job started on $(date)"
echo "Running on $(hostname)"
echo "Current directory: $(pwd)"

# Output directory
OUTPUT_DIR="/home/tud03125/pipeline/normal_human_liver_rnasplice"

# Create the output directory if it doesn't exist
mkdir -p "$OUTPUT_DIR"

# Load necessary modules
module load gcc/11.2.0
module load singularity
module load python/3.10.2 # Load Python 3.10.2

# Create an alias for python to point to python3
alias python=python3

# Check that the alias works
python --version # Should output Python 3.10.2

# Set Nextflow version
export NXF_VER=24.04.4

# Install Nextflow if not already installed
if [ ! -f ./nextflow ]; then
   echo "Nextflow not found, downloading..."
   curl -s <https://get.nextflow.io> | bash
   chmod +x nextflow
fi

# Set Nextflow work directory to avoid quota issues
export NXF_WORK=/home/tud03125/pipeline/normal_human_liver_rnasplice/work

# Run Nextflow with your local rnasplice pipeline
echo "Starting Nextflow pipeline..."
./nextflow run /home/tud03125/pipeline/rnasplice -resume -profile singularity \
   --input /home/tud03125/pipeline/normal_human_liver_design_rnasplice_3.csv \
   --fasta /home/tud03125/pipeline/Homo_sapiens.GRCh38.dna_sm.primary_assembly.fa \
   --gtf /home/tud03125/pipeline/Homo_sapiens.GRCh38.108.gtf \
   --save_reference \
   --outdir /home/tud03125/pipeline/normal_human_liver_rnasplice \
   --star_index /home/tud03125/STAR_index \
   --salmon_index /home/tud03125/new_salmon_index \
   --rmats \
   --suppa \
   --ribo_database_manifest /home/tud03125/ribo_database_manifest.txt \
   --contrasts /home/tud03125/pipeline/contrast_NORMAL_2.csv \
   --sashimi_plot \
   --miso_genes_file /home/tud03125/gene_list.txt

# Debugging: Print completion message and end timestamp
echo "Processing complete."
echo "Job ended on $(date)"
It’s just that instead of 
--miso_genes
, I use 
--miso_genes_file
, which’s part of its list of available parameters: https://nf-co.re/rnasplice/1.0.4/parameters/#miso. That 
gene_list.txt
I've created contains all 62,703 of them. ChatGPT says that it shouldn’t matter if I use 
--miso_genes
or 
--miso_genes_file
. But, what do you think?
t
It makes sense to use 
--miso_genes_file
(I hope that overrides 
--miso_genes
and you don't have to "unset" that one), but using that with >60k genes seems a bit over the top. The point of these plots is that you view them one-by-one, maybe put them in a publication to highlight something. You are not realistically looking at 60k genes, and most of these are not even expressed. If you really care about some pathway, you could put some key genes from that, etc.; certainly apply some filter (e.g. based on expression or splicing events). My approach would be to skip the plots in the first pass, look at the results from all the tools to figure out what is going on (and what you care about), and then maybe generate some plots to highlight key findings.
m
That's interesting. I'm actually more focused on intron retention analysis. I was looking into 
rMAT
and 
SUPPA
for that task, and found out that 
nf-core/rnasplice
uses 
rMAT
and 
SUPPA
and some others. So, I thought "this's useful. Could help with our intron retention analysis. So, why not? Let's just try it out to see what kind of info I'll get." So, the part about "most of these are not even expressed," that's actually the plan because my next experiment (which I've wrote recently in #C0493FQQZ3P) will involve analyzing Fasting vs. 3 different inhibitors (HCD and combinations). The goal is to analyze dsRNA, or to find out why some mRNA strands are not fully spliced out and capture those ones. Since we're looking for introns, specifically, we're looking more splicing events and less on expressions. Right now, I'm just doing a test run to see if I can get this pipeline working in our Temple HPC server. Getting the technicality work first before delving into the real, biological study of dsRNA or not-fully-spliced mRNA.
👍 1
t
For testing you could pick a couple of genes, just to see what you would get for the Sashimi plots.
👍 1
m
@Thomas Danhorn I’ve just re-ran the PBS script with 
--miso_gene_file
parameter we were discussing (using 
gene_list.txt
as that input). Still getting 
MISO_SASHIMI
errors. While I haven’t yet received the execution report (HTML), inside 
execution_trace_2024-08-27_00-26-17.txt
file, I see these:
Copy code
300       46/fc5354          191572             NFCORE_RNASPLICE:RNASPLICE:VISUALISE_MISO:MISO_SASHIMI (2)        FAILED 1           2024-08-27 21:32:14.553            6.5s      6.5s      -             -            -            -            -
301       10/a7710b        191733             NFCORE_RNASPLICE:RNASPLICE:VISUALISE_MISO:MISO_SASHIMI (3)        FAILED 1           2024-08-27 21:32:14.760            6.4s      6.4s      -             -            -            -            -
299       0e/9b3910        191671             NFCORE_RNASPLICE:RNASPLICE:VISUALISE_MISO:MISO_SASHIMI (1)        FAILED 1           2024-08-27 21:32:14.704            6.6s      6.5s      -             -            -            -            -
303       cb/40472c         191774             NFCORE_RNASPLICE:RNASPLICE:VISUALISE_MISO:MISO_SASHIMI (5)        FAILED 1           2024-08-27 21:32:14.788            6.7s      6.6s      -             -            -            -            -
302       0a/8e7ddb        191654             NFCORE_RNASPLICE:RNASPLICE:VISUALISE_MISO:MISO_SASHIMI (4)        FAILED 1           2024-08-27 21:32:14.680            6.9s      6.8s      -             -            -            -            -
304       09/b82862        191807             NFCORE_RNASPLICE:RNASPLICE:VISUALISE_MISO:MISO_SASHIMI (6)        FAILED 1           2024-08-27 21:32:14.811            7.2s      7.1s      -             -            -            -            -
306       47/35dedc         196325             NFCORE_RNASPLICE:RNASPLICE:VISUALISE_MISO:MISO_SASHIMI (8)        FAILED 1           2024-08-27 21:32:21.098            5.3s      5.3s      -             -            -            -            -
305       b7/b2cb33        196369             NFCORE_RNASPLICE:RNASPLICE:VISUALISE_MISO:MISO_SASHIMI (7)        FAILED 1           2024-08-27 21:32:21.174            5.5s      5.4s      -             -            -            -            -
Before was 3 of them instead of 8. But still, same problem even with including 
gene_list.txt
as gene identifiers. Just in case, this is an example of what 
gene_list.txt
looks like:
Copy code
ENSG00000000003
ENSG00000000005
ENSG00000000419
ENSG00000000457
ENSG00000000460
ENSG00000000938
ENSG00000000971
ENSG00000001036
ENSG00000001084
ENSG00000001167
ENSG00000001460
ENSG00000001461
ENSG00000001497
ENSG00000001561
ENSG00000001617
ENSG00000001626
ENSG00000001629
ENSG00000001630
ENSG00000001631
ENSG00000002016
ENSG00000002079
ENSG00000002330
ENSG00000002549
ENSG00000002586
ENSG00000002587
ENSG00000002726
ENSG00000002745
ENSG00000002746
ENSG00000002822
ENSG00000002834
ENSG00000002919
ENSG00000002933
ENSG00000003056
ENSG00000003096
ENSG00000003137
ENSG00000003147
ENSG00000003249
ENSG00000003393
ENSG00000003400
ENSG00000003402
ENSG00000003436
ENSG00000003509
ENSG00000003756
…
…
t
Just to clarify -- you ran with (at least) dozens of genes, and eight of them failed? Does that mean the others succeeded, or were they not attempted? Also, is the error message the same about the gene not being found in the "pickled" index? If so, are these the exact gene IDs from your GTF, and was the MISO index generated from that GTF (because that would be very strange)?
m
@Thomas Danhorn Just to clarify -- you ran with (at least) dozens of genes, and eight of them failed? • The whole 
gene_list.txt
was the whole Ensemble list, about ~62,000 of them. Does that mean the others succeeded, or were they not attempted? • I don't know why first time, without 
gene_list.txt
, shows 3 of them and then when adding 
gene_list.txt
shows 8 of them. But either one, after it tried 3 or 8 of them, it just stopped. So, probably the rest was "not attempted" after the first few failed. But, another interesting thing, I got 
root
email about execution terminated around 3:13 AM, and those 
FAILED
ones were in 
21:32:14.553
. The last line in 
trace
report was 
281	96/c58bfe	205849	NFCORE_RNASPLICE:RNASPLICE:DEXSEQ_DEU:DEXSEQ_EXON (1)	FAILED	-	2024-08-27 11:12:43.738	16h	16h	-	-	-	-	-
, which I can't tell what time was that and if that will be another 
FAILED
report I would get. Also, is the error message the same about the gene not being found in the "pickled" index? • Yes, the very same one! If so, are these the exact gene IDs from your GTF, and was the MISO index generated from that GTF (because that would be very strange)? • The whole 
gene_list.txt
was gathered from GTF. I even just ran the test check today about it: “1. Extract Gene IDs from GTF File,” “2. Compare Gene IDs with `gene_list.txt`:,” and then “3. Inspect Results:”
Copy code
(base) [tud03125@compute ~]$ awk -F 'gene_id "|";' '/gene_id/ {print $2}' /home/tud03125/pipeline/Homo_sapiens.GRCh38.108.gtf | sort | uniq > gtf_gene_ids.txt 
(base) [tud03125@compute ~]$ comm -23 <(sort gene_list.txt) <(sort gtf_gene_ids.txt) > missing_genes.txt 
(base) [tud03125@compute ~]$ cat missing_genes.txt (base) [tud03125@compute ~]$
Since you’ve found this one also very strange, should I also post this one in GitHub’s bug report?
@Thomas Danhorn when investigating this, while my 
gene_list.txt
got all the gene IDs from my GTF file, when checking gene IDs from miso index (I've made 
miso_gene_ids.txt
for that, to compare to my GTF's gene IDs), I found that the miso index has 0 gene IDs in it. And, when running tests on any issues with GTF, I've used these commands, copy-and-paste the outputs to ChatGPT for analysis:
Copy code
head /home/tud03125/pipeline/Homo_sapiens.GRCh38.108.gtf
grep 'gene_id' /home/tud03125/pipeline/Homo_sapiens.GRCh38.108.gtf | head -n 20
awk '$9 !~ /gene_id/ {print}' /home/tud03125/pipeline/Homo_sapiens.GRCh38.108.gtf | head -n 20
awk -F '\t' 'NF != 9' /home/tud03125/pipeline/Homo_sapiens.GRCh38.108.gtf | head -n 20
gffread /home/tud03125/pipeline/Homo_sapiens.GRCh38.108.gtf -o- | head -n 20
When ChatGPT checks it, it says “*Your GTF file appears to be correctly formatted, and there don't seem to be any obvious issues.* However, since your 
miso_gene_ids.txt
file is empty, it might indicate that the MISO_INDEX step didn't generate the expected gene entries, leading to the missing genes issue.” Can that be compatibility or configuration script issue, or something else?
t
How did you make that 
miso_gene_ids.txt
file? If the 
MISO_INDEX
step didn't use any genes, then I would have expected it to fail, killing the pipeline then and there ...
m
From this command: 
awk -F 'gene_id "|";' '/gene_id/ {print $2}' /home/tud03125/pipeline/normal_human_liver_rnasplice/misopy/index/*/genes.gff | sort | uniq > miso_gene_ids.txt
. I had to update 
/home/tud03125/pipeline/normal_human_liver_rnasplice/misopy/index/*/genes.gff
one (technically 
/home/tud03125/pipeline/normal_human_liver_rnasplice/misopy/index/genes.gff
) in order to properly find it. But, it's got nothing to do with the pipeline. It just extracts all the 
gene_id
from the 
/home/tud03125/pipeline/normal_human_liver_rnasplice/misopy/index/genes.gff
folder and then make 
miso_gene_ids.tx
file out of it. While 
gene_list.txt
extracts 
gene_id
from GTF file, 
miso_gene_ids.tx
one is from 
/home/tud03125/pipeline/normal_human_liver_rnasplice/misopy/index/genes.gff
that the pipeline automatically creates. Anyway, I'm not sure what I should do here. I could 1) check the configuration script and see what's wrong (though that would involve modyfing), 2) make "Bug Report" in GitHub, or 2) 
--sashimi_plot false
just to skip it and move on to the next one?
t
I doubt that exact same 
awk
command would work on a GFF file, the format is different from GTF. Before you file a bug report, let's find the actual problem. Maybe it does not warrant a report, and even if it does, it will help the devs fix it.
Here is an a way to extract gene IDs from a GFF and GTF:
Copy code
# GFF3:
grep -oPe '(?<=gene_id=)[^;]+' "$mygff3" | sort -u > gff3_geneIDs.txt
# GTF (a.k.a. GFF2.2):
grep -oPe '(?<=gene_id ")[^"]+' "$mygtf" | sort -u > gtf_geneIDs.txt
# mgff3 and mygtf should contain the respective file names
m
@Thomas Danhorn Those are good points! In fact, I’ve just learned the 
wc -l
trick. It's actually very useful, esp. with lots of lines in the command prompt. So instead of visually looking in Command Line, this command will count the total lines for me. So, using your methods, I did these:
Copy code
grep -oPe '(?<=gene_id=)[^;]+' "/home/tud03125/pipeline/normal_human_liver_rnasplice/misopy/Homo_sapiens.GRCh38.108_genes.gff3" | sort -u > gff3_geneIDs.txt
#or
grep -oPe '(?<=gene_id=)[^;]+' "/home/tud03125/pipeline/normal_human_liver_rnasplice/misopy/index/genes.gff" | sort -u > gff3_index_geneIDs.txt
grep -oPe '(?<=gene_id ")[^"]+' "/home/tud03125/pipeline/Homo_sapiens.GRCh38.108.gtf" | sort -u > gtf_geneIDs.txt
comm -23 <(sort gtf_geneIDs.txt) <(sort gff3_geneIDs.txt) > missing_geneIDs.txt
#or
comm -23 <(sort gtf_geneIDs.txt) <(sort gff3_index_geneIDs.txt) > missing_index_geneIDs.txt
wc -l missing_geneIDs.txt  # Either one of these two choices
And, this is the output:
Copy code
(gffread_env) [tud03125@compute ~]$ grep -oPe '(?<=gene_id=)[^;]+' "/home/tud03125/pipeline/normal_human_liver_rnasplice/misopy/Homo_sapiens.GRCh38.108_genes.gff3" | sort -u > gff3_geneIDs.txt
(gffread_env) [tud03125@compute ~]$ grep -oPe '(?<=gene_id ")[^"]+' "/home/tud03125/pipeline/Homo_sapiens.GRCh38.108.gtf" | sort -u > gtf_geneIDs.txt
(gffread_env) [tud03125@compute ~]$ comm -23 <(sort gtf_geneIDs.txt) <(sort gff3_geneIDs.txt) > missing_geneIDs.txt
(gffread_env) [tud03125@compute ~]$ wc -l missing_geneIDs.txt
62703 missing_geneIDs.txt
(gffread_env) [tud03125@compute ~]$ grep -oPe '(?<=gene_id=)[^;]+' "/home/tud03125/pipeline/normal_human_liver_rnasplice/misopy/index/genes.gff" | sort -u > gff3_index_geneIDs.txt
(gffread_env) [tud03125@compute ~]$ comm -23 <(sort gtf_geneIDs.txt) <(sort gff3_index_geneIDs.txt) > missing_index_geneIDs.txt
(gffread_env) [tud03125@compute ~]$ wc -l missing_geneIDs.txt
62703 missing_geneIDs.txt
(gffread_env) [tud03125@compute ~]$
So, when comparing GTF to GFF, 62,703 missing gene IDs, which’s actually what’s the total in my 
gene_list.txt
. Huh!? 😕 I also did 
awk
commands too such as 
awk -F'\t' '$3 == "gene" {print $9}' /home/tud03125/pipeline/Homo_sapiens.GRCh38.108.gtf | awk -F'; ' '{print $1}' | sed 's/gene_id "//; s/"//g' > gtf_gene_ids.txt
and 
awk -F'\t' '$3 == "gene" {print $9}' /home/tud03125/pipeline/normal_human_liver_rnasplice/misopy/index/genes.gff | awk -F';' '{for(i=1;i<=NF;i++){if($i ~ /^ID=/) print $i}}' | sed 's/ID=//g' > pipeline_index_gff_gene_ids.txt
or 
awk -F'\t' '$3 == "gene" {print $9}' /home/tud03125/pipeline/normal_human_liver_rnasplice/misopy/Homo_sapiens.GRCh38.108_genes.gff3 | awk -F';' '{for(i=1;i<=NF;i++){if($i ~ /^ID=/) print $i}}' | sed 's/ID=//g' > pipeline_gff3_gene_ids.txt
, then followed by 
comm
and 
wc -l
commands. Same “62,703 missing gene IDs.” So, I was trying to investigate why that’s happening. I was testing 
gffread
command for this. When initially doing 
gffread /home/tud03125/pipeline/Homo_sapiens.GRCh38.108.gtf -T -o Homo_sapiens.GRCh38.108_genes_full.gff3
, 62,703 missing gene IDs still persist. It’s not until I used this that it worked: 
gffread /home/tud03125/pipeline/Homo_sapiens.GRCh38.108.gtf -O -F -o Homo_sapiens.GRCh38.108_genes_full_with_all_features.gff3
. When comparing, using awk command (though that’s before your grep suggestions):
Copy code
# Step 1: Extract Gene IDs from the GTF file
awk -F'\t' '$3 == "gene" {print $9}' /home/tud03125/pipeline/Homo_sapiens.GRCh38.108.gtf | awk -F'; ' '{print $1}' | sed 's/gene_id "//; s/"//g' > gtf_gene_ids.txt

# Step 2: Extract Gene IDs from the GFF3 file
awk -F'\t' '$3 == "gene" {print $9}' Homo_sapiens.GRCh38.108_genes_full_with_all_features.gff3 | awk -F';' '{for(i=1;i<=NF;i++){if($i ~ /^ID=/) print $i}}' | sed 's/ID=//g' > gff3_gene_ids.txt

# Step 3: Compare the lists and find missing gene IDs
comm -23 <(sort gtf_gene_ids.txt) <(sort gff3_gene_ids.txt) > missing_gene_ids.txt

# Step 4: Count the missing gene IDs
wc -l missing_gene_ids.txt
The result shows this:
Copy code
(gffread_env) [tud03125@compute ~]$ awk -F'\t' '$3 == "gene" {print $9}' /home/tud03125/pipeline/Homo_sapiens.GRCh38.108.gtf | awk -F'; ' '{print $1}' | sed 's/gene_id "//; s/"//g' > gtf_gene_ids.txt
(gffread_env) [tud03125@compute ~]$ awk -F'\t' '$3 == "gene" {print $9}' Homo_sapiens.GRCh38.108_genes_full_with_all_features.gff3 | awk -F';' '{for(i=1;i<=NF;i++){if($i ~ /^ID=/) print $i}}' | sed 's/ID=//g' > gff3_gene_ids.txt
(gffread_env) [tud03125@compute ~]$ comm -23 <(sort gtf_gene_ids.txt) <(sort gff3_gene_ids.txt) > missing_gene_ids.txt (gffread_env)
[tud03125@compute ~]$ wc -l missing_gene_ids.txt 
0 missing_gene_ids.txt
(gffread_env) [tud03125@compute ~]$
@Thomas Danhorn If that's not "bug," could that be some sort of flaw inside modules configuration script, somewhere?
@Thomas Danhorn Using just one gene, e.g. 
ENSG00000000003
, ChatGPT found this critical info: Key Differences 1. Feature Types: ◦ The GFF3 file only includes gene and mRNA annotations, while the GTF file contains a more granular breakdown with exon, CDS, start_codon, stop_codon, and UTR annotations. 2. Missing Annotations: ◦ The absence of these specific features in the GFF3 file indicates that the conversion process from GTF to GFF3 might be filtering out or not including these detailed features.
t
GFF3 is very liberal with its naming conventions. I used 
gene_id
for the gene IDs, because that's what Ensembl uses, but MISO uses its own format. I no longer have the MISO index for my run, so I didn't know, but I found a different one, where genes are stored as 
ensgene_id=gene:ENSG000...
. You can use 
less
to look at text files (like GTFs, GFFs, etc., use 
zless
if they are gzipped) to see what they look like (and which conventions they use or features they contain. You can easily modify the 
grep
or 
awk
commands accordingly. Instead of relying on Chat GPT, I would suggest reading actual specifications and documentations. GFF3s can contain all the same features like GTFs, but if you think about it, we are looking at splicing, so the only relevant ones are exons. The GFF3 generated by MISO is not a straight "conversion" of the GTF, it was built for a very particular purpose, and it's not containing a gene annotation (hence no exons, let alone CDRa, etc.), but the intervals for the exons/introns used in particular splicing events (e.g. the SE index has all the annotated exon skipping events). The gene IDs here are just used as annotation. If you forget about the 
--miso_gene_file
for a moment (keeping things we change at a minimum) and grep one (or a few) of your genes (from the MISO GFF3) and use them with 
--miso_genes
, does that work? Before doing that, make sure the gene actually has enough exons (minimum 3 for skipping, 2 for RI) to have the potential for events. If it still complains about that gene not being in the pickled index, can you search the MISO index for that file, to make sure it's somewhere there, not just in the GFF3? The index has a directory for each event class (SE, RI, ...), and in that a subdirectory for each contig. You can search like this: 
grep -lFe ENG00... /path/to/SE.*_index/*/*.pickle
(adjust the path and use a full gene ID). This should give you at least one name of a file containing the gene ID.
m
@Thomas Danhorn. So, based on what you’ve suggested regarding 
grep -lFe ENG00... /path/to/SE.*_index/*/*.pickle
, I made this python script to search for two IDs: 
ENSG00000005302
and 
ENSG00000004961
based on the error messages in the HTML report I’m getting:
Copy code
import sys
sys.path.append('/path/to/misopy')

import pickle
import os

# Directory containing the pickle files
pickle_dir = "/home/tud03125/pipeline/normal_human_liver_rnasplice/misopy/index"

# List of gene IDs you want to search for
gene_ids_to_search = ["ENSG00000005302", "ENSG00000004961"]

# Function to search for gene IDs in a pickle file
def search_pickle_for_gene_ids(pickle_file_path, gene_ids):
   with open(pickle_file_path, 'rb') as f:
       data = pickle.load(f)
       for gene_id in gene_ids:
           if gene_id in str(data): # Convert data to string to search
               print("Found {0} in {1}".format(gene_id, pickle_file_path))

# Walk through the directory and search in each pickle file
for root, dirs, files in os.walk(pickle_dir):
   for file in files:
       if file.endswith(".pickle"):
           pickle_file_path = os.path.join(root, file)
           search_pickle_for_gene_ids(pickle_file_path, gene_ids_to_search)
After that run, for some time, this’s the output:
Copy code
((miso_env) [tud03125@compute ~]$ python search_pickle_files.py
Found ENSG00000005302 in /home/tud03125/pipeline/normal_human_liver_rnasplice/misopy/index/chrX/ENSG00000005302.19.pickle
Found ENSG00000004961 in /home/tud03125/pipeline/normal_human_liver_rnasplice/misopy/index/chrX/ENSG00000004961.15.pickle
(miso_env) [tud03125@compute ~]$
Can that mean anything to you?
@Thomas Danhorn regarding 
--miso_gene_file
and 
--miso_genes
, the first run's without 
--miso_gene_file
and the 2nd run's with it. Both lead me to that same error. I haven't tested 
--miso_genes
one, though.
Also, when looking into MISO, it apparently is only compatible with Python 2.7 version. That and that there's "deprecated modules in 
Matplotlib
," the 
MatplotlibDeprecationWarning
and Matplotlib's version is 2.2.5, could that play role here?
One more message for tonight. So, I was doing this run, to see if 
ENSG00000005302
exists in the `shelve`:
Copy code
(miso_env) [tud03125@compute ~]$ python -c "
> import shelve
> shelve_file = '/home/tud03125/pipeline/normal_human_liver_rnasplice/work/45/bed4975c16bbed45c4df4bfd0e3614/index/genes_to_filenames.shelve'
> db = shelve.open(shelve_file)
> try:
>    if 'ENSG00000005302' in db:
>        print('Pickle filename for ENSG00000005302:', db['ENSG00000005302'])
>    else:
>        print('ENSG00000005302 does not exist in the shelve.')
> finally:
>    db.close()
> "
ENSG00000005302 does not exist in the shelve.
(miso_env) [tud03125@compute ~]$
So, 
genes_to_filenames.shelve
doesn’t have it despite 
ENSG00000005302
existing in my pickle directories (BTW the message above was involving 
/home/tud03125/pipeline/normal_human_liver_rnasplice/misopy/
directory, while 
work
directory has more of them). Not only that, but I also found this in one of the 
.command.sh
lists:
Copy code
(base) [tud03125@compute ~]$ cat /home/tud03125/pipeline/normal_human_liver_rnasplice/work/46/fc53548fb4b289d78a609eb2388218/.command.sh 
#!/bin/bash -euo pipefail 
sashimi_plot --plot-event ENSG00000005302 index miso_settings.txt --output-dir sashimi 
cat <<-END_VERSIONS > versions.yml "NFCORE_RNASPLICE:RNASPLICE:VISUALISE_MISO:MISO_SASHIMI": python: $(python --version | sed "s/Python //g") misopy: $(python -c "import pkg_resources; 
print(pkg_resources.get_distribution('misopy').version)") 
END_VERSIONS
Currently, I don’t know which script 
rnasplice
used to create 
genes_to_filenames.shelve
, yet. And also, I've test run 
sashimi_plot --plot-event ENSG00000005302 index miso_settings.txt --output-dir sashimi
command in my prompt, using their 
sashimi_plot.py
script, and I'm still getting this Traceback error: 
Exception: Event ENSG00000004961 not found in pickled directory index. Are you sure this is the right directory for the event?
. So, given this situation, do you think it’s worth making a bug report about it in GitHub (and then proceed with 
--sashimi_plot false
) or leave it alone?
t
• Your Python script found the genes in the index, that was all we need to know (
grep
would have given the same result faster, though). So I don't understand why you are getting the error -- is that index the one being used by the pipeline?? • The first error you got is because 
--miso_genes
is by default set something that does not work if you use your own GTF (I consider this a bug, but that has been submitted a number of times already and there is no consensus). So I want to try overriding this without using yet another parameter (of which I don't know if it normally works). • MISO is old and unmaintained, hence the need for Python 2.7. Similarly with the deprecation warning, but it is just a warning (not an error), and it means that some syntax or feature will no longer be available in a future version; it does not affect this version (it's a hint for the programmer to update the code in time before it breaks). • The 
sashimi_plot --plot-event ENSG00000005302
is expected, that's how the plot is being made (if it works), and I think that's the process throwing the error. • It makes sense that you can reproduce the error outside the pipeline. I don't understand why though -- was the index not made properly? We have to find the root of the problem, otherwise a report is not terribly useful. • It is interesting that the gene is not in the 
genes_to_filenames.shelve
, but I honestly don't know if (and how) it is supposed to. Is there any gene that works (if the pipeline has not produces any plot already, you could run the Sashimi plot command outside the pipeline with various gene IDs and see if any of them works)? If so, can you find that ID in the shelve file? • In case there is a problem with the index, it might be worthwhile to recreate it (there is some process that does that, or you can probably find it in the MISO documentation) and run the command with the new index to test it. Sashimi plots work in general, we have to find out what is wrong in your specific case.
m
• You’ve mentioned “index” a lot. Just to check, like “is that index the one being used by the pipeline” or “I don't understand why though -- was the index not made properly? We have to find the root of the problem…” Just to check, you’re talking about 
miso_index
? Would that be the “root of the problem” you’re referring to? • About the “The 
sashimi_plot --plot-event ENSG00000005302
is expected, that's how the plot is being made (if it works), and I think that's the process throwing the error,” I agree. • About “if the pipeline has not produces any plot already, you could run the Sashimi plot command outside the pipeline with various gene IDs and see if any of them works,” no Sashimi plots made yet, ever since that error pops up. But, I would like to know if there’s a way to bypass these errors, like some if/else statement, like “if XXX gene ID is not found, continue or ‘move forward” or something like that. But for now, I’ll try to identify which gene IDs are found in pickle directory. • Some useful info, on one of the 
work
indexes:
Copy code
(base) [tud03125@compute ~]$ cd /home/tud03125/pipeline/normal_human_liver_rnasplice/work/45/bed4975c16bbed45c4df4bfd0e3614/index/
(base) [tud03125@compute index]$ ls -l
total 27104
drwxr-xr-x 2 tud03125 eng    5661 Aug 27 02:09 chr1
drwxr-xr-x 2 tud03125 eng    2430 Aug 27 02:10 chr10
drwxr-xr-x 2 tud03125 eng    3472 Aug 27 02:10 chr11
drwxr-xr-x 2 tud03125 eng    3146 Aug 27 02:10 chr12
drwxr-xr-x 2 tud03125 eng    1460 Aug 27 02:10 chr13
drwxr-xr-x 2 tud03125 eng    2342 Aug 27 02:10 chr14
drwxr-xr-x 2 tud03125 eng    2323 Aug 27 02:10 chr15
drwxr-xr-x 2 tud03125 eng    2675 Aug 27 02:10 chr16
drwxr-xr-x 2 tud03125 eng    3166 Aug 27 02:10 chr17
drwxr-xr-x 2 tud03125 eng    1267 Aug 27 02:10 chr18
drwxr-xr-x 2 tud03125 eng    3078 Aug 27 02:10 chr19
drwxr-xr-x 2 tud03125 eng    4345 Aug 27 02:09 chr2
drwxr-xr-x 2 tud03125 eng    1504 Aug 27 02:09 chr20
drwxr-xr-x 2 tud03125 eng     899 Aug 27 02:09 chr21
drwxr-xr-x 2 tud03125 eng    1447 Aug 27 02:09 chr22
drwxr-xr-x 2 tud03125 eng    3305 Aug 27 02:09 chr3
drwxr-xr-x 2 tud03125 eng    2734 Aug 27 02:09 chr4
drwxr-xr-x 2 tud03125 eng    3076 Aug 27 02:09 chr5
drwxr-xr-x 2 tud03125 eng    3183 Aug 27 02:09 chr6
drwxr-xr-x 2 tud03125 eng    3149 Aug 27 02:09 chr7
drwxr-xr-x 2 tud03125 eng    2542 Aug 27 02:10 chr8
drwxr-xr-x 2 tud03125 eng    2418 Aug 27 02:10 chr9
drwxr-xr-x 2 tud03125 eng       3 Aug 27 02:10 chrGL000009.2
drwxr-xr-x 2 tud03125 eng       4 Aug 27 02:09 chrGL000194.1
drwxr-xr-x 2 tud03125 eng       4 Aug 27 02:09 chrGL000195.1
drwxr-xr-x 2 tud03125 eng       3 Aug 27 02:09 chrGL000205.2
drwxr-xr-x 2 tud03125 eng       4 Aug 27 02:09 chrGL000213.1
drwxr-xr-x 2 tud03125 eng       3 Aug 27 02:10 chrGL000216.2
drwxr-xr-x 2 tud03125 eng       3 Aug 27 02:09 chrGL000218.1
drwxr-xr-x 2 tud03125 eng       3 Aug 27 02:09 chrGL000219.1
drwxr-xr-x 2 tud03125 eng       6 Aug 27 02:10 chrGL000220.1
drwxr-xr-x 2 tud03125 eng       3 Aug 27 02:10 chrGL000225.1
drwxr-xr-x 2 tud03125 eng       4 Aug 27 02:10 chrKI270442.1
drwxr-xr-x 2 tud03125 eng       3 Aug 27 02:10 chrKI270711.1
drwxr-xr-x 2 tud03125 eng       6 Aug 27 02:09 chrKI270713.1
drwxr-xr-x 2 tud03125 eng       6 Aug 27 02:10 chrKI270721.1
drwxr-xr-x 2 tud03125 eng       4 Aug 27 02:09 chrKI270726.1
drwxr-xr-x 2 tud03125 eng       9 Aug 27 02:09 chrKI270727.1
drwxr-xr-x 2 tud03125 eng       8 Aug 27 02:10 chrKI270728.1
drwxr-xr-x 2 tud03125 eng       4 Aug 27 02:09 chrKI270731.1
drwxr-xr-x 2 tud03125 eng       6 Aug 27 02:09 chrKI270733.1
drwxr-xr-x 2 tud03125 eng       6 Aug 27 02:10 chrKI270734.1
drwxr-xr-x 2 tud03125 eng       3 Aug 27 02:10 chrKI270744.1
drwxr-xr-x 2 tud03125 eng       3 Aug 27 02:10 chrKI270750.1
drwxr-xr-x 2 tud03125 eng      39 Aug 27 02:10 chrMT
drwxr-xr-x 2 tud03125 eng    2483 Aug 27 02:10 chrX
drwxr-xr-x 2 tud03125 eng     555 Aug 27 02:10 chrY
-rw-r--r-- 1 tud03125 eng       0 Aug 27 02:11 compressed_ids_to_genes.shelve.bak
-rw-r--r-- 1 tud03125 eng       0 Aug 27 02:11 compressed_ids_to_genes.shelve.dat
-rw-r--r-- 1 tud03125 eng       0 Aug 27 02:11 compressed_ids_to_genes.shelve.dir
-rw-r--r-- 1 tud03125 eng 7718907 Aug 27 02:11 genes.gff
-rw-r--r-- 1 tud03125 eng 2250012 Aug 27 02:12 genes_to_filenames.shelve.bak
-rw-r--r-- 1 tud03125 eng 32103470 Aug 27 02:12 genes_to_filenames.shelve.dat
-rw-r--r-- 1 tud03125 eng 2250012 Aug 27 02:12 genes_to_filenames.shelve.dir
(base) [tud03125@compute index]$ wc -l genes_to_filenames.shelve.dir
62703 genes_to_filenames.shelve.dir
(base) [tud03125@compute index]$ wc -l genes_to_filenames.shelve.bak
62703 genes_to_filenames.shelve.bak
(base) [tud03125@compute index]$ wc -l genes_to_filenames.shelve.dat
125406 genes_to_filenames.shelve.dat
(base) [tud03125@compute index]$ wc -l genes.gff
62703 genes.gff
(base) [tud03125@compute index]$
• And, this info too:
Copy code
(base) [tud03125@compute ~]$ wc -l Homo_sapiens.GRCh38.108.gtf
3409311 Homo_sapiens.GRCh38.108.gtf
(base) [tud03125@compute ~]$
t
• Yes, I'm talking about the 
miso_index
-- the error message complains that the gene(s) in question cannot be found there, even though they are in the GFF3 in the iindex, as well as the 
*.pickle
files. ◦ The "root of the problem" is the reason why these genes can't be found in the MISO index -- that's what we want to know. • You can "bypass errors" by setting 
errorStrategy = 'ignore'
for the process (since this is a terminal process in the sense that nothing depends on it, this should have no implications for the rest of the pipeline). Doing that would be interesting, as all genes in your list/parameter would have a chance to get plotted, so we can see if any of them succeed. Though given what we have seen so far, I doubt it ... • The 
*shelve
files are binary, so 
wc -l
(which basically counts newline characters) is not a meaningful operation on them. (It is for text files, like GTF/GFF3 etc.)
2 Views
Open in Slack
PreviousNext
Powered by