ORQAS fails at Ribomap step

[Alex-Nesta]

I have tried running both the run_ORQAS.sh script and running the pipeline step by step. I have tried using my ownrnaseq and riboseq data and the sample data. I keep getting the same error.

terminate called after throwing an instance of 'std::invalid_argument'
what(): stod
ORQAS/ribomap/scripts/run_ribomap.sh: line 329: 180610 Aborted

Do you have any experience with this error. I have been digging through this for a whole day now. I appreciate any insight!

I verified I am using salmon 0.7.2 for index and alignment as well as STAR 0.4.2j for the riboprof step.
If you suspect it is a software version incompatibility, would you consider making a docker or singularity container?

Thanks for your help!
Alex

(/projects/beck-lab/alex/miniconda/seqtk-env) [nestaa@sumner-log1 ORQAS]$ ORQAS/ribomap/scripts/run_ribomap.sh \
> --rnaseq_fq ${base}/RNA-Seq/GSM3611929_1.fastq \
> --riboseq_fq ${base}/Ribo-Seq/GSM3611926.fastq \
> --transcript_fa ${base}/Iso_BRCA/renamed_cds_sequences.fasta \
> --work_dir ${base}/stepwise/ \
> --offset ${base}/psite.txt \
> --softClipping false \
> --rnaUnstranded true \
> --sailfish_dir ${base}/stepwise/salmon_out
contaminant fasta not provided! Filter step skipped.
skipped filter read step.
aligning reads to transcriptome
aligning RNA_seq to transcriptome...
alignment skipped.
align_reads /pod/2/beck-lab/alex/Breast_Cancer/scripts/ORQAS/RNA-Seq/GSM3611929_1.fastq /pod/2/beck-lab/alex/Breast_Cancer/scripts/ORQAS/stepwise/alignment/GSM3611929_1_transcript_ /pod/2/beck-lab/alex/Breast_Cancer/scripts/ORQAS/stepwise/alignment/GSM3611929_1_transcript_Aligned.out.bam
aligning ribo_seq to transcriptome...
alignment skipped.
align_reads /pod/2/beck-lab/alex/Breast_Cancer/scripts/ORQAS/Ribo-Seq/GSM3611926.fastq /pod/2/beck-lab/alex/Breast_Cancer/scripts/ORQAS/stepwise/alignment/GSM3611926_transcript_ /pod/2/beck-lab/alex/Breast_Cancer/scripts/ORQAS/stepwise/alignment/GSM3611926_transcript_Aligned.out.bam
running riboprof...
cds_range file not provided, assume transcript fasta only contains cds sequences.
getting transcript info...
total number of transcripts in transcriptome: 94317
assigning ribo-seq reads...
constructing profile class...
terminate called after throwing an instance of 'std::invalid_argument'
  what():  stod
ORQAS/ribomap/scripts/run_ribomap.sh: line 329: 180610 Aborted                 /pod/2/beck-lab/alex/Breast_Cancer/scripts/ORQAS/ORQAS/ribomap/bin/riboprof ${options}
pipeline failed at ribosome profile generation!

[Alex-Nesta]

Just repeated with a fresh install of everything and using sample data. Same error. Sorry it's a lot to read through:

Collecting package metadata (current_repodata.json): done
Solving environment: failed with repodata from current_repodata.json, will retry with next repodata source.
Collecting package metadata (repodata.json): done
Solving environment: done


==> WARNING: A newer version of conda exists. <==
  current version: 4.10.0
  latest version: 4.10.1

Please update conda by running

    $ conda update -n base conda


conda activate orqas2-env

## Package Plan ##

  environment location: /projects/beck-lab/alex/miniconda/orqas2-env

  added / updated specs:
    - python=3.5


The following packages will be downloaded:

    package                    |            build
    ---------------------------|-----------------
    certifi-2018.8.24          |        py35_1001         139 KB  conda-forge
    pip-20.3.4                 |     pyhd8ed1ab_0         1.1 MB  conda-forge
    python-3.5.5               |       h5001a0f_2        20.2 MB  conda-forge
    setuptools-40.4.3          |           py35_0         572 KB  conda-forge
    sqlite-3.28.0              |       h8b20d00_0         1.9 MB  conda-forge
    ------------------------------------------------------------
                                           Total:        23.9 MB

The following NEW packages will be INSTALLED:

  _libgcc_mutex      conda-forge/linux-64::_libgcc_mutex-0.1-conda_forge
  _openmp_mutex      conda-forge/linux-64::_openmp_mutex-4.5-1_gnu
  bzip2              conda-forge/linux-64::bzip2-1.0.8-h7f98852_4
  ca-certificates    conda-forge/linux-64::ca-certificates-2020.12.5-ha878542_0
  certifi            conda-forge/linux-64::certifi-2018.8.24-py35_1001
  libffi             conda-forge/linux-64::libffi-3.2.1-he1b5a44_1007
  libgcc-ng          conda-forge/linux-64::libgcc-ng-9.3.0-h2828fa1_19
  libgomp            conda-forge/linux-64::libgomp-9.3.0-h2828fa1_19
  libstdcxx-ng       conda-forge/linux-64::libstdcxx-ng-9.3.0-h6de172a_19
  ncurses            conda-forge/linux-64::ncurses-6.2-h58526e2_4
  openssl            conda-forge/linux-64::openssl-1.0.2u-h516909a_0
  pip                conda-forge/noarch::pip-20.3.4-pyhd8ed1ab_0
  python             conda-forge/linux-64::python-3.5.5-h5001a0f_2
  readline           conda-forge/linux-64::readline-7.0-hf8c457e_1001
  setuptools         conda-forge/linux-64::setuptools-40.4.3-py35_0
  sqlite             conda-forge/linux-64::sqlite-3.28.0-h8b20d00_0
  tk                 conda-forge/linux-64::tk-8.6.10-h21135ba_1
  wheel              conda-forge/noarch::wheel-0.36.2-pyhd3deb0d_0
  xz                 conda-forge/linux-64::xz-5.2.5-h516909a_1
  zlib               conda-forge/linux-64::zlib-1.2.11-h516909a_1010


Proceed ([y]/n)? Invalid choice: conda activate orqas2-env
Proceed ([y]/n)? y


Downloading and Extracting Packages
python-3.5.5         | 20.2 MB   | ################################################################################################################################################################################################################################### | 100% 
certifi-2018.8.24    | 139 KB    | ################################################################################################################################################################################################################################### | 100% 
setuptools-40.4.3    | 572 KB    | ################################################################################################################################################################################################################################### | 100% 
sqlite-3.28.0        | 1.9 MB    | ################################################################################################################################################################################################################################### | 100% 
pip-20.3.4           | 1.1 MB    | ################################################################################################################################################################################################################################### | 100% 
Preparing transaction: done
Verifying transaction: done
Executing transaction: done
#
# To activate this environment, use
#
#     $ conda activate orqas2-env
#
# To deactivate an active environment, use
#
#     $ conda deactivate

(base) [nestaa@sumner030 ORQAS2]$ nano ~/.bash_profile 
(base) [nestaa@sumner030 ORQAS2]$ conda activate orqas2-env
(/projects/beck-lab/alex/miniconda/orqas2-env) [nestaa@sumner030 ORQAS2]$ ls
ORQAS
(/projects/beck-lab/alex/miniconda/orqas2-env) [nestaa@sumner030 ORQAS2]$ cd ORQAS/
(/projects/beck-lab/alex/miniconda/orqas2-env) [nestaa@sumner030 ORQAS]$ sh run_ORQAS.sh --rnaseq_fq sample_data/chr21_test_rna.fastq --riboseq_fq sample_data/chr21_test_ribo.fastq --cds_annotation sample_data/chr21_sample_CDS_annotation.fa --salmon_strand U --psites sample_data/sample_psites.txt --wd sample_data_test
Version Info: ### A newer version of Salmon is available. ####
###
The newest version, available at https://github.com/COMBINE-lab/salmon/releases 
contains new features, improvements, and minor bug fixes; please upgrade at your
earliest convenience.
###
index ["sample_data_test/sm_index/"] did not previously exist  . . . creating it
[2021-04-29 22:32:55.306] [jLog] [info] building index
RapMap Indexer

[Step 1 of 4] : counting k-mers
Elapsed time: 0.0516667s

Replaced 105 non-ATCG nucleotides
Clipped poly-A tails from 0 transcripts
Building rank-select dictionary and saving to disk done
Elapsed time: 0.000251107s
Writing sequence data to file . . . done
Elapsed time: 0.000917063s
[info] Building 32-bit suffix array (length of generalized text is 1105278)
Building suffix array . . . success
saving to disk . . . done
Elapsed time: 0.00291992s
done
Elapsed time: 0.0864539s
processed 1000000 positions
khash had 341627 keys
saving hash to disk . . . done
Elapsed time: 0.0260355s
[2021-04-29 22:32:55.875] [jLog] [info] done building index
run_ORQAS.sh: line 73: error_msg: command not found
Version Info: ### A newer version of Salmon is available. ####
###
The newest version, available at https://github.com/COMBINE-lab/salmon/releases 
contains new features, improvements, and minor bug fixes; please upgrade at your
earliest convenience.
###
### salmon (mapping-based) v0.7.2
### [ program ] => salmon 
### [ command ] => quant 
### [ unmatedReads ] => { sample_data/chr21_test_rna.fastq }
### [ index ] => { sample_data_test/sm_index/ }
### [ libType ] => { U }
### [ output ] => { sample_data_test/sm_quant/ }
Logs will be written to sample_data_test/sm_quant/logs
[2021-04-29 22:32:56.198] [jointLog] [info] parsing read library format
[2021-04-29 22:32:56.198] [jointLog] [info] There is 1 library.
[2021-04-29 22:32:56.342] [jointLog] [info] Loading Quasi index
[2021-04-29 22:32:56.342] [jointLog] [info] Loading 32-bit quasi index
[2021-04-29 22:32:56.377] [jointLog] [info] done
[2021-04-29 22:32:56.377] [jointLog] [info] Index contained 834 targets
[2021-04-29 22:32:56.343] [stderrLog] [info] Loading Suffix Array 
[2021-04-29 22:32:56.345] [stderrLog] [info] Loading Transcript Info 
[2021-04-29 22:32:56.346] [stderrLog] [info] Loading Rank-Select Bit Array
[2021-04-29 22:32:56.346] [stderrLog] [info] There were 834 set bits in the bit array
[2021-04-29 22:32:56.347] [stderrLog] [info] Computing transcript lengths
[2021-04-29 22:32:56.347] [stderrLog] [info] Waiting to finish loading hash
[2021-04-29 22:32:56.377] [stderrLog] [info] Successfully loaded position hash
[2021-04-29 22:32:56.377] [stderrLog] [info] Done loading index








[2021-04-29 22:33:09.985] [jointLog] [info] Computed 679 rich equivalence classes for further processing
[2021-04-29 22:33:09.985] [jointLog] [info] Counted 229935 total reads in the equivalence classes 






[2021-04-29 22:33:10.052] [jointLog] [warning] 5.7467% of fragments were shorter than the k used to build the index (31).
If this fraction is too large, consider re-building the index with a smaller k.
The minimum read size found was 25.


[2021-04-29 22:33:10.053] [jointLog] [warning] Only 229935 fragments were mapped, but the number of burn-in fragments was set to 5000000.
The effective lengths have been computed using the observed mappings.

[2021-04-29 22:33:10.053] [jointLog] [info] Mapping rate = 92.8775%

[2021-04-29 22:33:10.053] [jointLog] [info] finished quantifyLibrary()
[2021-04-29 22:33:10.056] [jointLog] [info] Starting optimizer
[2021-04-29 22:33:10.058] [jointLog] [info] Marked 0 weighted equivalence classes as degenerate
[2021-04-29 22:33:10.058] [jointLog] [info] iteration = 0 | max rel diff. = 8.58943
[2021-04-29 22:33:10.071] [jointLog] [info] iteration = 100 | max rel diff. = 0.112244
[2021-04-29 22:33:10.100] [jointLog] [info] iteration = 200 | max rel diff. = 0.0547118
[2021-04-29 22:33:10.144] [jointLog] [info] iteration = 300 | max rel diff. = 0.0383001
[2021-04-29 22:33:10.190] [jointLog] [info] iteration = 400 | max rel diff. = 0.0273019
[2021-04-29 22:33:10.230] [jointLog] [info] iteration = 500 | max rel diff. = 0.019537
[2021-04-29 22:33:10.276] [jointLog] [info] iteration = 600 | max rel diff. = 0.0119184
[2021-04-29 22:33:10.318] [jointLog] [info] iteration = 700 | max rel diff. = 0.0116674
[2021-04-29 22:33:10.362] [jointLog] [info] iteration = 800 | max rel diff. = 0.0117435
[2021-04-29 22:33:10.410] [jointLog] [info] iteration = 900 | max rel diff. = 0.0117991
[2021-04-29 22:33:10.455] [jointLog] [info] iteration = 1000 | max rel diff. = 0.0118409
[2021-04-29 22:33:10.500] [jointLog] [info] iteration = 1100 | max rel diff. = 0.0118732
[2021-04-29 22:33:10.538] [jointLog] [info] iteration = 1200 | max rel diff. = 0.0118985
[2021-04-29 22:33:10.573] [jointLog] [info] iteration = 1286 | max rel diff. = 0.0083639
[2021-04-29 22:33:10.573] [jointLog] [info] Finished optimizer
[2021-04-29 22:33:10.573] [jointLog] [info] writing output 

contaminant fasta not provided! Filter step skipped.
skipped filter read step.
aligning reads to transcriptome
building transcriptome index...
Apr 29 22:33:10 ..... Started STAR run
Apr 29 22:33:10 ... Starting to generate Genome files
Apr 29 22:33:11 ... starting to sort  Suffix Array. This may take a long time...
Apr 29 22:33:11 ... sorting Suffix Array chunks and saving them to disk...
Apr 29 22:33:12 ... loading chunks from disk, packing SA...
Apr 29 22:33:12 ... writing Suffix Array to disk ...
Apr 29 22:33:13 ... Finished generating suffix array
Apr 29 22:33:13 ... starting to generate Suffix Array index...
Apr 29 22:33:13 ... writing SAindex to disk
Apr 29 22:33:13 ..... Finished successfully
aligning RNA_seq to transcriptome...
Apr 29 22:33:13 ..... Started STAR run
Apr 29 22:33:13 ..... Started mapping
Apr 29 22:33:29 ..... Finished successfully
align_reads sample_data/chr21_test_rna.fastq sample_data_test/alignment/chr21_test_rna_transcript_ sample_data_test/alignment/chr21_test_rna_transcript_Aligned.out.bam
aligning ribo_seq to transcriptome...
Apr 29 22:33:29 ..... Started STAR run
Apr 29 22:33:29 ..... Started mapping
Apr 29 22:33:49 ..... Finished successfully
align_reads sample_data/chr21_test_ribo.fastq sample_data_test/alignment/chr21_test_ribo_transcript_ sample_data_test/alignment/chr21_test_ribo_transcript_Aligned.out.bam
running riboprof...
cds_range file not provided, assume transcript fasta only contains cds sequences.
getting transcript info...
total number of transcripts in transcriptome: 834
assigning ribo-seq reads...
constructing profile class...
terminate called after throwing an instance of 'std::invalid_argument'
  what():  stod
./ribomap/scripts/run_ribomap.sh: line 329: 28354 Aborted                 riboprof ${options}
pipeline failed at ribosome profile generation!
2021-04-29 22:33:49,946 - INFO: Accessing dir: sample_data_test/
2021-04-29 22:33:49,949 - INFO: File sample_data_test/sm_quant/quant.sf closed.
2021-04-29 22:33:49,950 - INFO: Error on line 227
2021-04-29 22:33:49,950 - INFO: not enough values to unpack (expected 1, got 0)
2021-04-29 22:33:49,950 - INFO: Error on line 60
2021-04-29 22:33:49,951 - INFO: 1
Transcript to cds assignment not provided... calculating it.
2021-04-29 22:33:50,017 - INFO: File sample_data/chr21_sample_CDS_annotation.fa closed.
2021-04-29 22:33:50,020 - INFO: Process ended succesfully.
2021-04-29 22:33:50,077 - INFO: Reading file sample_data_test/sample_data_test.Abundance.txt
2021-04-29 22:33:50,077 - INFO: Error on line 109
2021-04-29 22:33:50,077 - INFO: [Errno 2] No such file or directory: 'sample_data_test/sample_data_test.Abundance.txt'
2021-04-29 22:33:50,077 - INFO: Error on line 63
2021-04-29 22:33:50,077 - INFO: 1
DONE! Outputs are in sample_data_test
To obtain periodicity and uniformity values run validateORF with a single sample or a pool of samples from the same condition
(/projects/beck-lab/alex/miniconda/orqas2-env) [nestaa@sumner030 ORQAS]$ 

[Alex-Nesta]

UPDATE: Tried with a fresh ubuntu VM instead of using my university cluster. Same error.

Getting closer... Definitely an issue with the riboprof step of ribomap:

I believe I am failing on line 147 of ribomap in read_mapper.cpp

https://github.com/Kingsford-Group/ribomap/blob/e96093d97c8fd2bf2f6abe453b9fc7abaad00a9c/src/read_mapper.cpp

bool readmapper_pipeline(const transcript_info& tinfo, const char* mRNA_bam, const char* ribo_bam, int lmin, int lmax, bool useSecondary, bool useRC, const string& offset, const char* tabd_fname, const string& tabd_type, double tabd_cutoff, const string& log_prefix)
{
  cout<<"assigning ribo-seq reads..."<<endl;
  cout<<"constructing profile class..."<<endl;
  ribo_profile rprofile(tinfo, tabd_fname, tabd_type, tabd_cutoff);
  cout<<"number of transcripts in profile class: "<<rprofile.number_of_transcripts()<<endl;
  cout<<"loading reads from bam..."<<endl;

Specifically, I think one of the function arguments for rprofile is being provided as a string, when a double is expected, leading to the following error:

terminate called after throwing an instance of 'std::invalid_argument'
  what():  stod

I think that would be one of these arguments when running riboprof, although they are all numeric....:

--min_fplen 25 
--max_fplen 33
--offset 11
--tabd_cutoff 0.01  

I think the next step is to cout<< all of the arguments and see what is getting hung up.

(/projects/beck-lab/alex/miniconda/orqas2-env) [nestaa@sumner030 ORQAS]$ riboprof \
> --fasta /pod/2/beck-lab/alex/Breast_Cancer/scripts/ORQAS2/ORQAS/sample_data/chr21_sample_CDS_annotation.fa \
> --cds_range /pod/2/beck-lab/alex/Breast_Cancer/scripts/ORQAS2/ORQAS/sample_data_test/sample_data_test.ENSTtoCDS.txt \
> --mrnabam /pod/2/beck-lab/alex/Breast_Cancer/scripts/ORQAS2/ORQAS/sample_data_test/alignment/chr21_test_rna_transcript_Aligned.out.bam \
> --ribobam /pod/2/beck-lab/alex/Breast_Cancer/scripts/ORQAS2/ORQAS/sample_data_test/alignment/chr21_test_ribo_transcript_Aligned.out.bam \
> --min_fplen 25 \
> --max_fplen 33 \
> --offset 11 \
> --sf /pod/2/beck-lab/alex/Breast_Cancer/scripts/ORQAS2/ORQAS/sample_data_test/sm_quant/quant.sf \
> --tabd_cutoff 0.01  \
> --out /pod/2/beck-lab/alex/Breast_Cancer/scripts/ORQAS2/ORQAS/sample_data_test/outputs/ribomapout
getting transcript info...
total number of transcripts in transcriptome: 834
assigning ribo-seq reads...
constructing profile class...
terminate called after throwing an instance of 'std::invalid_argument'
  what():  stod
Aborted
(/projects/beck-lab/alex/miniconda/orqas2-env) [nestaa@sumner030 ORQAS]$ 

[Alex-Nesta]

I wrote a hack to get it working for the most part.

The issue seems to be that the line headers of the quant.sf are read into riboprof and cause it to crash. deleting the header of quant.sf allows riboprof to complete successfully.

I replaced step 2 or run_ORQAS.sh starting at line 99 with the following code.
Essentially, the run_ribomap.sh is run once, fails, then the first line of quant.sf is deleted, and then run_ribomap.sh is run again.

while this is not elegant, I ran into other issues when I tried to delete the first line of quant.sf before running run_ribomap.sh once.

I hope this information helps you implement a proper solution, but maybe helps others in the meantime.

####################
# 2: RUN RIBOMAP   #
####################

if [[ ${salmon_strand} == *"U"* ]]; then
     rnaseq_strand="true"
else
     rnaseq_strand="false"
fi

	
if [[ -e  "${wd}/outputs/" ]]; then
      echo "Ribomap output already exists. Skipping..."
else 
      ${orqas_dir}/ribomap/scripts/run_ribomap.sh \
      --rnaseq_fq ${rnaseq_fq} \
      --riboseq_fq ${riboseq_fq} \
      --transcript_fa ${cds_annotation} \
      --work_dir "${wd}/" \
      --offset ${psites} \
      --softClipping false \
      --rnaUnstranded ${rnaseq_strand}  \
      --sailfish_dir ${salmon_output}    
      
      sed '1d' ${salmon_output}/quant.sf > ${salmon_output}/quant.sftmpfile; mv ${salmon_output}/quant.sftmpfile ${salmon_output}/quant.sf
      
      ${orqas_dir}/ribomap/scripts/run_ribomap.sh \
      --rnaseq_fq ${rnaseq_fq} \
      --riboseq_fq ${riboseq_fq} \
      --transcript_fa ${cds_annotation} \
      --work_dir "${wd}/" \
      --offset ${psites} \
      --softClipping false \
      --rnaUnstranded ${rnaseq_strand}  \
      --sailfish_dir ${salmon_output}    
fi

[noahpieta]

Hi Alex-Nesta,

I'm also trying to use ORQAS and what you have reported is really useful for me. Really appreciate!

May I ask if you or any other has ever meet another issue which reported like this:
run_ribomap.sh: line 331: 45925 Segmentation fault (core dumped) /riboprof ${options}

Please let me know if you know how to deal with it. Thank you!

I wrote a hack to get it working for the most part.

The issue seems to be that the line headers of the quant.sf are read into riboprof and cause it to crash. deleting the header of quant.sf allows riboprof to complete successfully.

I replaced step 2 or run_ORQAS.sh starting at line 99 with the following code. Essentially, the run_ribomap.sh is run once, fails, then the first line of quant.sf is deleted, and then run_ribomap.sh is run again.

while this is not elegant, I ran into other issues when I tried to delete the first line of quant.sf before running run_ribomap.sh once.

I hope this information helps you implement a proper solution, but maybe helps others in the meantime.

####################
# 2: RUN RIBOMAP   #
####################

if [[ ${salmon_strand} == *"U"* ]]; then
     rnaseq_strand="true"
else
     rnaseq_strand="false"
fi

	
if [[ -e  "${wd}/outputs/" ]]; then
      echo "Ribomap output already exists. Skipping..."
else 
      ${orqas_dir}/ribomap/scripts/run_ribomap.sh \
      --rnaseq_fq ${rnaseq_fq} \
      --riboseq_fq ${riboseq_fq} \
      --transcript_fa ${cds_annotation} \
      --work_dir "${wd}/" \
      --offset ${psites} \
      --softClipping false \
      --rnaUnstranded ${rnaseq_strand}  \
      --sailfish_dir ${salmon_output}    
      
      sed '1d' ${salmon_output}/quant.sf > ${salmon_output}/quant.sftmpfile; mv ${salmon_output}/quant.sftmpfile ${salmon_output}/quant.sf
      
      ${orqas_dir}/ribomap/scripts/run_ribomap.sh \
      --rnaseq_fq ${rnaseq_fq} \
      --riboseq_fq ${riboseq_fq} \
      --transcript_fa ${cds_annotation} \
      --work_dir "${wd}/" \
      --offset ${psites} \
      --softClipping false \
      --rnaUnstranded ${rnaseq_strand}  \
      --sailfish_dir ${salmon_output}    
fi

[EduEyras]

I cc Marina in case she can indicate what could be the cause of this ribomap issue E.

…

On Wed, 30 Mar 2022 at 21:39, Noah Yu Zhang ***@***.***> wrote: Hi Alex-Nesta, I'm also trying to use ORQAS and what you have reported is really useful for me. Really appreciate! May I ask if you or any other has ever meet another issue which reported like this: run_ribomap.sh: line 331: 45925 Segmentation fault (core dumped) /riboprof ${options} Please let me know if you know how to deal with it. Thank you! — Reply to this email directly, view it on GitHub <#1 (comment)>, or unsubscribe <https://github.com/notifications/unsubscribe-auth/ADCZKB34Q57ZHZSRXB4MQ7TVCQVPZANCNFSM433DVDIA> . You are receiving this because you are subscribed to this thread.Message ID: ***@***.***>

-- Prof. E Eyras EMBL Australia Group Leader The John Curtin School of Medical Research - Australian National University https://github.com/comprna http://scholar.google.com/citations?user=LiojlGoAAAAJ

[Alex-Nesta]

I believe the issue is the output file column names. Not sure if it was a software version mismatch? Hopefully that narrows things down. A docker/singularity container might help eliminate such issues in the future. On Oct 11, 2022, at 4:09 AM, Eduardo Eyras ***@***.***> wrote:  I cc Marina in case she can indicate what could be the cause of this ribomap issue E.

…

On Wed, 30 Mar 2022 at 21:39, Noah Yu Zhang ***@***.***> wrote: Hi Alex-Nesta, I'm also trying to use ORQAS and what you have reported is really useful for me. Really appreciate! May I ask if you or any other has ever meet another issue which reported like this: run_ribomap.sh: line 331: 45925 Segmentation fault (core dumped) /riboprof ${options} Please let me know if you know how to deal with it. Thank you! — Reply to this email directly, view it on GitHub <#1 (comment)>, or unsubscribe <https://github.com/notifications/unsubscribe-auth/ADCZKB34Q57ZHZSRXB4MQ7TVCQVPZANCNFSM433DVDIA> . You are receiving this because you are subscribed to this thread.Message ID: ***@***.***>

-- Prof. E Eyras EMBL Australia Group Leader The John Curtin School of Medical Research - Australian National University https://github.com/comprna http://scholar.google.com/citations?user=LiojlGoAAAAJ — Reply to this email directly, view it on GitHub<#1 (comment)>, or unsubscribe<https://github.com/notifications/unsubscribe-auth/AJ6DWRDTRXYCYBFFRA26H5LWCUODVANCNFSM433DVDIA>. You are receiving this because you authored the thread.Message ID: ***@***.***> --- The information in this email, including attachments, may be confidential and is intended solely for the addressee(s). If you believe you received this email by mistake, please notify the sender by return email as soon as possible.

[Kratos12138]

I wrote a hack to get it working for the most part.

The issue seems to be that the line headers of the quant.sf are read into riboprof and cause it to crash. deleting the header of quant.sf allows riboprof to complete successfully.

I replaced step 2 or run_ORQAS.sh starting at line 99 with the following code. Essentially, the run_ribomap.sh is run once, fails, then the first line of quant.sf is deleted, and then run_ribomap.sh is run again.

while this is not elegant, I ran into other issues when I tried to delete the first line of quant.sf before running run_ribomap.sh once.

I hope this information helps you implement a proper solution, but maybe helps others in the meantime.

####################
# 2: RUN RIBOMAP   #
####################

if [[ ${salmon_strand} == *"U"* ]]; then
     rnaseq_strand="true"
else
     rnaseq_strand="false"
fi

	
if [[ -e  "${wd}/outputs/" ]]; then
      echo "Ribomap output already exists. Skipping..."
else 
      ${orqas_dir}/ribomap/scripts/run_ribomap.sh \
      --rnaseq_fq ${rnaseq_fq} \
      --riboseq_fq ${riboseq_fq} \
      --transcript_fa ${cds_annotation} \
      --work_dir "${wd}/" \
      --offset ${psites} \
      --softClipping false \
      --rnaUnstranded ${rnaseq_strand}  \
      --sailfish_dir ${salmon_output}    
      
      sed '1d' ${salmon_output}/quant.sf > ${salmon_output}/quant.sftmpfile; mv ${salmon_output}/quant.sftmpfile ${salmon_output}/quant.sf
      
      ${orqas_dir}/ribomap/scripts/run_ribomap.sh \
      --rnaseq_fq ${rnaseq_fq} \
      --riboseq_fq ${riboseq_fq} \
      --transcript_fa ${cds_annotation} \
      --work_dir "${wd}/" \
      --offset ${psites} \
      --softClipping false \
      --rnaUnstranded ${rnaseq_strand}  \
      --sailfish_dir ${salmon_output}    
fi

Hi, Alex, lots of thanks, this helps me a lot!
Now I'm facing the difficulty of make a p-sites file by myself, can I get some help or take a look of your p-sites file? Thank you!