Allow tabs in fasta header when creating decoys for salmon index

[paoloAngelino]

Description of the bug

A fasta header can contain comments together with the name of the contig. Example:
>HLA-DRB1*16:02:01 HLA00878

The corresponding line in the decoy.txt file to be passed to salmon index would be
HLA-DRB1*16:02:01 HLA00878

the problem is that the comment is interpreted by salmon as an extra decoy, while creating the index, and it stops with an error. In my case:

[2022-09-05 14:58:36.028] [puff::index::jointLog] [critical] The decoy file contained the names of 3892 decoy sequences, but 3367 were matched by sequences in the reference file provided. To prevent unintentional errors downstream, please ensure that the decoy file exactly matches with the fasta file that is being indexed.
[2022-09-05 14:58:36.424] [puff::index::jointLog] [error] The fixFasta phase failed with exit code 1

An additional cleaning step in rnaseq/modules/nf-core/modules/salmon/index/main.nf would fix the issue. What I propose is to replace line 31:
sed -i.bak -e 's/>//g' decoys.txt
with

mv decoys.txt decoys.txt.bak
awk '{print $1}' decoys.txt.bak > decoys.txt 
sed -i -e 's/>//g' > decoys.txt

Command used and terminal output

No response

Relevant files

No response

System information

No response

[drpatelh]

Hi @paoloAngelino ! Thanks for reporting.

I am unable to reproduce this if I run:

$ echo ">HLA-DRB1*16:02:01      HLA00878" | grep '^>' | cut -d ' ' -f 1
>HLA-DRB1*16:02:01

Which is essentially running this bit of code in the module that splits on spaces and only takes the first value. As you can see above this prints the contig name without the comment delimited by spaces.

[drpatelh]

Closing this for now but feel free to re-open if the issue persists or if you are able to identify another reason for this failure.

[paoloAngelino]

@drpatelh Thanks for your answer.
I've actually found that in the reference that I'm using, the comments in the contig header are separated by a TAB.

grep '^>' hg38.fa  | cat -T
>HLA-DRB1*15:03:01:01^IHLA00870 11567 bp
>HLA-DRB1*15:03:01:02^IHLA03454 11569 bp
>HLA-DRB1*16:02:01^IHLA00878 11005 bp

This is the reason why, in my case, cut -d ' ' -f 1 doesn't split the header.
I don't know if having TAB as separator is a special case, or if it could happen for other fasta references. But, maybe including also TAB as a separator could be more general? This would work for me:
grep '^>' hg38.fa | cut -d ' ' -f 1 | cut -d $'\t' -f 1

[drpatelh]

Cool. Will re-open and fix in the next release. In the meantime, you will have to clean the fasta before passing to the pipeline.