Hi Everyone, I have a potential bug question. The...
# rnasplice
m
Hi Everyone, I have a potential bug question. The execution is “successful.” But, not getting the results I want. The contrast I am using are these:
Copy code
contrast	treatment	control
HCD_FASTED	HCD	FASTED
HCD_MK2206_FASTED	HCD_MK2206	FASTED
HCD_SPHINX31_FASTED	HCD_SPHINX31	FASTED
HCD_MK2206_HCD	HCD_MK2206	HCD
HCD_SPHINX31_HCD	HCD_SPHINX31	HCD
For the rMATS, the folders are created (1st image attachment below). And, the files too, including “rmats_post” where the analysis takes place. However, the issue is that “HCD-FASTED” subdirectory has information, while the other 4 pairwise comparisons (like “HCD_MK2206-HCD” for example, below), their “.MATS.JC” and “.MATS.JCEC” text files are empty (2nd image attachment below). Has anyone come across this issue, that one out of [insert any number here] pairwise-comparisons contains valuable analysis and the rest empty, like for rMATS, the pipeline just picked one pairwise-comparison and leave the rest (for the “.MATS.JC” and “.MATS.JCEC” text files)?
Just in case clarity is an issue:
t
This hasn't happened to me, but I don't run this very often. If you have confirmed in the
.nextflow.log
that there really was no error, and if I had to guess, I would say there was a problem copying the results from the work directory to the results. (This would likely not be a pipeline-specific issue, but more a glitch somewhere between Nextflow and your system.) I would check the rMATS process work directories for the comparisons that are missing files to see if they are there (if not, have a look at the
.command.*
files there). If the files are indeed there, you could try rerunning the pipeline with
-resume
, in the hope that it finds the previous results and copies them properly this time. If all else fail, you can copy them manually.
m
When saying "This would likely not be a pipeline-specific issue, but more a glitch somewhere between Nextflow and your system," that "system" you're saying, is that referring to where I am running the pipeline? All of my bioinformatics run are done in my Temple Main Campus's HPC Cluster. Some stuff there are old (one person stated this in the email: "The Temple HPC resources are built on CentOS 7 which is very old and uses obsolete versions of the OpenSSL library (and other software). There is not much we can do about that. Since we support containers via apptainer / singularity, probably the easiest workaround is to download a suitable container image of a more recent Linux version and then start a shell from within that container and run your software."). Anyway, can that have something to do with it?
The other thing is that this is happening to all my rnasplice runs (5 in total). 1st time, with 3 contrasts, one contains details and the other 2 was empty, the “.MATS.JC” and “.MATS.JCEC” text files. When I increased to 5 contrast the 2nd time and the rest of the times, one out of 5 wasn't empty. It wasn't just one time for me. Can that still be that glitch?
t
... that "system" you're saying, is that referring to where I am running the pipeline?
Yes, I'm talking about the HPC in your case. I don't think it's necessarily an issue of age, I have had failures with writing files to disk completely (they were corrupted) when our cluster was overloaded (and that uses relatively new hardware and software). I'm not sure what the problem is, does the OS of the cluster not issue an appropriate error, or is Nextflow not detecting that error? Either way, something that goes wrong, and Nextflow does not learn about it, so it thinks that everything is fine (until it tries to use the corrupted file in a later process).
It wasn't just one time for me. Can that still be that glitch?
I honestly don't know. If you get the same (wrong) result repeatedly, that would point to a deterministic cause, and I would suspect wrong input or a bug in the code. But if could also be that every time you run this, at some point so many processes are running that load on the cluster increases to point where these I/O failures occur. And if the conditions are similar enough, that might occur always at roughly the same time, leading to the same problems. This does not sound super likely, but if you get only one result due to a bug, I would expect that there is an issue with the naming (i.e. everything is called the same) and newer results overwrite older ones. In your case, however, it is something different, since you are getting the directories (with correct, expected names), but they are empty ... So I'm not sure what would cause that, and it is tempting to blame some random glitch.
m
About the “bug in the code,” like what? So when I run the
rnasplice
pipeline, one of my five PBS/Shell scripts looks like this:
Copy code
#!/bin/bash

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

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

# 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

# Change directory to your pipeline directory
cd /home/tud03125/pipeline

# 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="$OUTPUT_DIR/work"

# Extract gene versions for MISO_SASHIMI
GTF_FILE="/home/tud03125/pipeline/Mus_musculus.GRCm39.112.gtf"
GENE_VERSIONS=$(awk -F '[\t ;"]+' '/gene_id "(ENSMUSG00000019836|ENSMUSG00000020887)"/ {for (i=1; i<=NF; i++) if ($i=="gene_id") {gid=$(i+1)} else if ($i=="gene_version") {print gid"."$(i+1)}}' $GTF_FILE | sort | uniq | paste -sd "," -)

# Run Nextflow with your local rnasplice pipeline and include the mouse liver data
echo "Starting Nextflow pipeline..."
./nextflow run nf-core/rnasplice -resume -profile singularity \
   --input /home/tud03125/pipeline/mouse_liver_2.csv \
   --fasta /home/tud03125/pipeline/Mus_musculus.GRCm39.dna_sm.primary_assembly.fa \
   --gtf /home/tud03125/pipeline/Mus_musculus.GRCm39.112.gtf \
   --save_reference \
   --outdir "$OUTPUT_DIR" \
   --star_index /home/tud03125/STAR_index_mouse \
   --salmon_index /home/tud03125/salmon_index_mouse_noversion \
   --rmats \
   --rmats_novel_splice_site \
   --suppa \
   --contrasts /home/tud03125/pipeline/contrast_mouse_liver_2.csv \
   --miso_genes "$GENE_VERSIONS"

# Debugging: Print completion message and end timestamp
echo "Processing complete."
echo "Job ended on $(date)"
Any potential bugs from there? From my experience, if there’s a flaw, like not right FASTA or GTF files, or indexes (STAR or salmon), I would get back an error execution report, pointing those things out. I did run into them. But, at least the pipeline stops, tells me those errors and I’d correct them. This one (above), after all the trials-and-errors, chiseling out any issues, runs successfully every time I run it. But, that just means I don’t get syntax errors or anything similar. So anyway, if I execute this Shell script, what kind of “bug” could I face? And, the “wrong input,” like wrong directory on where the FASTQ files are located?
t
I was just going through the things that could potentially cause an issue in general -- something in the code (I don't see anything obviously wrong in your code, I was more thinking about the pipeline), something in the data (which the code did not anticipate), or some (quasi-)random outside factor (an overloaded OS reacting unpredictably, a cosmic ray flipping a bit in RAM [on an HPC the RAM should have checksums, though], etc.) Issues with code and data should be reproducible. You can play around, testing possibilities by splitting the data up. E.g., what happens if you just run one comparison (one that gave no results)? Before you do that. though, find one of the work dirs from an "empty" process (either in the
.nextflow.log
, or maybe they are in the
execution_report_*.html
as well) and see what is there. Are there rMATS output files? Any interesting messages in
.command.*
files? Is there an
.exitcode
file with
0
in it?