1
$\begingroup$

I try to pass a .fasta file with headers where tabs are used as separators as input file to NCBI BLAST+ segmasker:
>A0A059VAR9 IPR022879 V-type ATP synthase regulatory subunit B/beta

Unfortunately, the segmasker output's headers are shortened to the first tab:
>A0A059VAR9

Is there a way to get whole headers in the output file (without replacing all the tabs in the input file with white spaces)?

$\endgroup$
3
  • $\begingroup$ Could you supply a small data set please? You can fill is with dumbie amino acids ... I'll look at NCBI'ssegmasker local version because thats what I've got installed. $\endgroup$
    – M__
    Nov 26, 2022 at 12:17
  • $\begingroup$ I think I also use local version. Here's a sample data: link $\endgroup$
    – maciejwww
    Nov 26, 2022 at 12:26
  • $\begingroup$ Thanks for the file, answered below. Sorry about the delay. This question isn't simplistic. $\endgroup$
    – M__
    Nov 28, 2022 at 9:50

1 Answer 1

1
$\begingroup$

I used a local version of segmasker *. One way to solve this is,

perl -p -i -e 's/ //g' bioinfo.fa
perl -p -e 's/\t/^/g' bioinfo.fa | segmasker | perl -p -e 's/^/\t/g'  > bioinfo2.sm

Edit: I changed , for ^ because the OPs code already has , and will insert an extra tab (\t). Same change in the line below

Here the input file you presented in the link is called bioinfo.fa. If you need the pipe explaining let me know, but it's fairly straight forward.

If you want to preserve space characters, i.e. ' ' or "\s" - because in the code above I delete them - then (not tested):

perl -p -e 's/\t/^/g' bioinfo.fa | perl -p -e 's/ /_/g' | segmasker | perl -p -e 's/_/ /g' | perl -p -e 's/^/\t/g'  > bioinfo2.sm

The issue the OP encounter is not unique to tabs (\t) but also occurs with spaces (\s), at least with makebasedb.

input

>A0A000 IPR004839   Aminotransferase, class I/classII   PF00155 41  381
GVSAIDGMNTLVRGQKLPIFSSSGQPHSELAMQIARQASVPEEEEGGDDEEGSEFAVIFGAMGITQEEANEFMEDFERTGALERSVVFMNLADDPAVERTVTPRMVLTTAEYLAFEKDYHVLVILTDMTNYCEALREIGAAREEVPGRRGYPG
>A0A000 IPR010961   Tetrapyrrole biosynthesis, 5-aminolevulinic acid synthase   TIGR01821   12  391
...

output

>A0A000 IPR004839       Aminotransferase        class_I/classII        PF00155 41      381
41 - 53
>A0A000 IPR010961       Tetrapyrrole_biosynthesis       _5aminolevulinic_acid_synthase  TIGR01821       12      391
41 - 53
>A0A000 IPR015421       Pyridoxal phosphate dependent transferase        major domainG3DSA:3.40.640.10  48      288
106 - 126
...

Please don't ask about solving this via an upstream makeblastdb because that gets complicated. The rate limiting step, i.e. why it is ignoring tabs is because of makeblastdb. Thats a consequence of the way its written ... its a developer issue (see blurb)


* Local segmasker installation ...

  1. Install conda via Anaconda or miniconda
  2. conda create -n blast3 python=3.11 (creating an environment ... yep 3.11 is out (wow!, if you don't understand why WOW! ask in chat)
  3. conda activate blast3 (this is the conda environment, you can call it anything)
  4. conda install -c bioconda blast

perl | segmasker | perl will now work. The Python installation is just good practice.


blurb about NCBI's blast The issue is makeblastdb, what happens is that long headers are a complete pain in algorithms because they use loads of RAM and you've got to carry them through reiterative loops (convergence etc ...). The way to solve it is via recoding (thats how I've done it). Essentially, you have a mini database forcing header recoding which activates at the start of the algorithm and reassembles the original header at the end of the algorithm. The user never sees this step, but it's really important for efficiency. Thus within the reiteration the header isn't consuming RAM, because it's an internal code, beyond base 10 to squash the overhead even further. The thing is header can become a lot bigger than the sequence data and thats not cool when you've loads and loads of data, e.g. a blast database. NCBI just haven't bothered, like if you can't see it why do it attitude. What they've said is that 50 characters is the limitation of the header for makeblastdb... what NCBI is looking for is a unique id at the start of the header it then truncates everything beyond that, i.e. tabs, spaces ... anything it can grab. NCBI have carried this through into segmasker it's not needed, but makeblastdb and segmasker would normally be used together - thats why. In the OPs case this results in duplicate headers, so the file format falls over. This isn't the OPs fault.

$\endgroup$
3
  • 1
    $\begingroup$ Thank you for the exhaustive answer! $\endgroup$
    – maciejwww
    Nov 28, 2022 at 20:46
  • $\begingroup$ I agree, Python 3.11 has many great features. What's your biggest WOW there? $\endgroup$
    – maciejwww
    Nov 28, 2022 at 20:48
  • 1
    $\begingroup$ @maciejwww The criticism of Python is its slow. 3.11 is the major upgrade to focus on speed - its the fight back and may there be many more. Answer, you are welcome, note was a slight bug, I switched , for ^. $\endgroup$
    – M__
    Nov 29, 2022 at 1:06

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.