EI-CoreBioinformatics/mikado

[main_samview] fail to read the header from "-".

NJeanray opened this issue · 4 comments

NJeanray commented

Hello,

I got some troubles running alignment with gsnap. Here is the error I get :

[Fri Aug  6 18:17:34 2021]
Job 6: Aligning RNAseq with gsnap (sample SNIPr_01_03_24_43_2009_9_16_VI_RP_1_S5_ME - run 0)
Reason: Missing output files: Daijin/2-alignments/output/gsnap-SNIPr_01_03_24_43_2009_9_16_VI_RP_1_S5_ME-0.bam; Input files updated by another job: /home/NJEANRAY/Daijin/Daijin/2-alignments/gsnap/gsnap_index.done

  gsnap --gunzip --bunzip2 --dir=Daijin/2-alignments/gsnap/index --db=Daijin  --novelsplicing=1 --localsplicedist=10000 --nthreads=40 --format=sam --npaths=20 Daijin/1-reads/SNIPr_01_03_24_43_2009_9_16_VI_RP_1_S5_ME.R1.fq.gz Daijin/1-reads/SNIPr_01_03_24_43_2009_9_16_VI_RP_1_S5_ME.R2.fq.gz 2> Daijin/2-alignments/logs/gsnap/gsnap-SNIPr_01_03_24_43_2009_9_16_VI_RP_1_S5_ME-0.log | samtools view -b -@ 40 - > Daijin/2-alignments/gsnap/SNIPr_01_03_24_43_2009_9_16_VI_RP_1_S5_ME-0/gsnap.bam && ln -sf ../gsnap/SNIPr_01_03_24_43_2009_9_16_VI_RP_1_S5_ME-0/gsnap.bam Daijin/2-alignments/output/gsnap-SNIPr_01_03_24_43_2009_9_16_VI_RP_1_S5_ME-0.bam && touch -h Daijin/2-alignments/output/gsnap-SNIPr_01_03_24_43_2009_9_16_VI_RP_1_S5_ME-0.bam
[main_samview] fail to read the header from "-".
[Fri Aug  6 18:17:34 2021]
Error in rule align_gsnap:
    jobid: 6
    output: Daijin/2-alignments/gsnap/SNIPr_01_03_24_43_2009_9_16_VI_RP_1_S5_ME-0/gsnap.bam, Daijin/2-alignments/output/gsnap-SNIPr_01_03_24_43_2009_9_16_VI_RP_1_S5_ME-0.bam
    log: Daijin/2-alignments/logs/gsnap/gsnap-SNIPr_01_03_24_43_2009_9_16_VI_RP_1_S5_ME-0.log (check log file(s) for error message)
    shell:
          gsnap --gunzip --bunzip2 --dir=Daijin/2-alignments/gsnap/index --db=Daijin  --novelsplicing=1 --localsplicedist=10000 --nthreads=40 --format=sam --npaths=20 Daijin/1-reads/SNIPr_01_03_24_43_2009_9_16_VI_RP_1_S5_ME.R1.fq.gz Daijin/1-reads/SNIPr_01_03_24_43_2009_9_16_VI_RP_1_S5_ME.R2.fq.gz 2> Daijin/2-alignments/logs/gsnap/gsnap-SNIPr_01_03_24_43_2009_9_16_VI_RP_1_S5_ME-0.log | samtools view -b -@ 40 - > Daijin/2-alignments/gsnap/SNIPr_01_03_24_43_2009_9_16_VI_RP_1_S5_ME-0/gsnap.bam && ln -sf ../gsnap/SNIPr_01_03_24_43_2009_9_16_VI_RP_1_S5_ME-0/gsnap.bam Daijin/2-alignments/output/gsnap-SNIPr_01_03_24_43_2009_9_16_VI_RP_1_S5_ME-0.bam && touch -h Daijin/2-alignments/output/gsnap-SNIPr_01_03_24_43_2009_9_16_VI_RP_1_S5_ME-0.bam
        (one of the commands exited with non-zero exit code; note that snakemake uses bash strict mode!)

Removing output files of failed job align_gsnap since they might be corrupted:
Daijin/2-alignments/gsnap/SNIPr_01_03_24_43_2009_9_16_VI_RP_1_S5_ME-0/gsnap.bam

Here are the versions of the tools used :

GSNAP version 2021-07-23 called with args: gsnap.avx2 --version
samtools 1.11
Using htslib 1.11

Could you please advice ?
Thanks in advance,

Best regards

ljyanesm commented

Dear @NJeanray,

This seems like an issue with the version of samtools. Probably a previous version will work for now.

lucventurini commented

Dear @NJeanray

Would you please be able to check the error log here:

Daijin/2-alignments/logs/gsnap/gsnap-SNIPr_01_03_24_43_2009_9_16_VI_RP_1_S5_ME-0.log

Another option would be to run the command:

gsnap --gunzip --bunzip2 --dir=Daijin/2-alignments/gsnap/index --db=Daijin --novelsplicing=1 --localsplicedist=10000 --nthreads=40 --format=sam --npaths=20 Daijin/1-reads/SNIPr_01_03_24_43_2009_9_16_VI_RP_1_S5_ME.R1.fq.gz Daijin/1-reads/SNIPr_01_03_24_43_2009_9_16_VI_RP_1_S5_ME.R2.fq.gz

as a standalone and check whether it actually functions? I cannot test locally at the moment (on a Mac machine and conda only has GSNAP 2021-05-27).

Things to look for: no errors and correct BAM header.

If there are errors we can take it from there.

Kind regards

NJeanray commented

Dear @lucventurini,

Please find below the content of

Daijin/2-alignments/logs/gsnap/gsnap-SNIPr_01_03_24_43_2009_9_16_VI_RP_1_S5_ME-0.log 

GSNAP version 2021-07-23 called with args: gsnap.avx2 --gunzip --bunzip2 --dir=Daijin/2-alignments/gsnap/index --db=Daijin --novelsplicing=1 --localsplicedist=10000 --nthreads=40 --format=sam --npaths=20 Daijin/1-reads/SNIPr_01_03_24_43_2009_9_16_VI_RP_1_S5_ME.R1.fq.gz Daijin/1-reads/SNIPr_01_03_24_43_2009_9_16_VI_RP_1_S5_ME.R2.fq.gz
Novel splicing (-N) turned on => assume reads are RNA-Seq
The novelend-splicedist 80000 is greater than the localsplicedist 10000.  Resetting novelend-splicedist to be 10000
Checking compiler assumptions for SSE2: 6B8B4567 327B23C6 xor=59F066A1
Checking compiler assumptions for SSE4.1: -103 -58 max=198 => compiler zero extends
Checking compiler assumptions for SSE4.2 options: 6B8B4567 __builtin_clz=1 __builtin_ctz=0 _mm_popcnt_u32=17 __builtin_popcount=17 
Finished checking compiler assumptions
This is a large genome of more than 2^32 (4 billion) bp.
You should run gsnapl instead.

Is this helping ?

Best regards

lucventurini commented

Dear @NJeanray

Yes this is helping ... the problem is that the genome is too big for a normal GMAP build and run. The logic in Daijin unfortunately does not support providing custom command line arguments to the building ... so there is currently no direct way to solve the issue.

We would need to amend the Daijin code to allow passing through building switches to GMAP for this to function.

Our oversight, we generally do not use Daijin with GSNAP or on such big genomes, so this bug slipped through.

CC @ljyanesm

Kind regards

Luca Venturini