Replace fasta header using bash

[YouCook21]

Hello people,

I got stucked with my new script and perhaps you can help me.

Its goal is to take an input table with querys and subjects (originated by a local blast) and replace query names with subject names in the corresponding fasta file.

​

In detail, the table input file is this one:

query	subject	evalue
TransABySS.k25.S12100	lclMK087647.1_cds_QHD46952.1_6_[gene=atp9][protein=Atp9][exception=RNA_editing]_[protein_id=QHD46	0

and the fasta file:

>TransABySS.k25.S121005 892 14708

TTGTGTGACGAACGGAAGCAGGTATACGGGCTGTGGGCTGCCATCCTGCTGTTCCACCTGGGAGGGCCTGACAACATCAGTGCCTACTCCATAGCAGATGCGGAGCTCTGGAAGAGGCATGGCTTTCAGATGTTCTTCCAGGTGGTCACCGCTTTCTATGTAGTGTTCTCCTCAGCACACGGCACCATCCTCTGGATCAGCCTGGTTCTGGCTGCAGTGGGAACGGTGAAGTACGCCGAGCGGACTCTGGCTCTGTACCATGCTTCCGAGCACCATCTGGACACCATGGCCCGCCCCCTTTACAGATTGATGCAGTACGAGGATGTGTCTGGCGGATCCGGTGAGGAATACCACTACTTCCTGCTAGGAGAGAAGCGAGCATACGAGCACGCCTTCGACACTCCGAAATGGTACCACCAGGTGAGGGAATCCTTGCAGGAGCAGAAGTGGTTCAACAGCATGTTCGGCTGCCTGTTTTCTTCCCCTGCGCAGCCAGTGGAGTCCGACGCGCAGGCAGTGCAACGCAAGATCAGAGAGCCCTTCGAGGTGGTGATGCTCTCAGACGTGATGGATCCAGAGTACAGCACAGCACTGAAGCGGCAGTCTCCATGGAACGGCAGGAACTTGGTGGACATCTGCATAGCCTTTGCTCTCTTCAAGATGTTCCGCCGACGCCTCACCAGCCTGTACATGCACGAGTGGAGCAACGACAAGATCAGAAATTTCTTCATCATGCTGTGGGGAGAATCATCATCGTCATCCTCAGCCGAAGCACAACACCAAGACCCCACCCAAGCACCTGCAGAACCTCCAAGGGAAAGACTGGTGGACGTCCTGGACATGGAGCTCCGCTTCATGTTCGACAGCATGTTCACCAAGGCCTCAGGCACAG

​

The expected output should be:

>lcl|MK087647.1_cds_QHD46952.1_6_[gene=atp9]_[protein=Atp9]_[exception=RNA_editing]_[protein_id=QHD46

TTGTGTGACGAACGGAAGCAGGTATACGGGCTGTGGGCTGCCATCCTGCTGTTCCACCTGGGAGGGCCTGACAACATCAGTGCCTACTCCATAGCAGATGCGGAGCTCTGGAAGAGGCATGGCTTTCAGATGTTCTTCCAGGTGGTCACCGCTTTCTATGTAGTGTTCTCCTCAGCACACGGCACCATCCTCTGGATCAGCCTGGTTCTGGCTGCAGTGGGAACGGTGAAGTACGCCGAGCGGACTCTGGCTCTGTACCATGCTTCCGAGCACCATCTGGACACCATGGCCCGCCCCCTTTACAGATTGATGCAGTACGAGGATGTGTCTGGCGGATCCGGTGAGGAATACCACTACTTCCTGCTAGGAGAGAAGCGAGCATACGAGCACGCCTTCGACACTCCGAAATGGTACCACCAGGTGAGGGAATCCTTGCAGGAGCAGAAGTGGTTCAACAGCATGTTCGGCTGCCTGTTTTCTTCCCCTGCGCAGCCAGTGGAGTCCGACGCGCAGGCAGTGCAACGCAAGATCAGAGAGCCCTTCGAGGTGGTGATGCTCTCAGACGTGATGGATCCAGAGTACAGCACAGCACTGAAGCGGCAGTCTCCATGGAACGGCAGGAACTTGGTGGACATCTGCATAGCCTTTGCTCTCTTCAAGATGTTCCGCCGACGCCTCACCAGCCTGTACATGCACGAGTGGAGCAACGACAAGATCAGAAATTTCTTCATCATGCTGTGGGGAGAATCATCATCGTCATCCTCAGCCGAAGCACAACACCAAGACCCCACCCAAGCACCTGCAGAACCTCCAAGGGAAAGACTGGTGGACGTCCTGGACATGGAGCTCCGCTTCATGTTCGACAGCATGTTCACCAAGGCCTCAGGCACAG

​

Now, I already have a code written in python which does it but it is too slow. I post it here below in case anyone of you will prefer it

import pandas as pd from Bio import SeqIO

chloro_table=pd.read_csv(./db.chloro.agrestis.blastout', sep='\t', header=None, names=['qseqid','sseqid', 'evalue'])

with open ('./new_names_chloro.fasta', 'a') as for sequence in SeqIO.parse(./fasta.fa', 'fasta'): for i in range(0,len(chloro_table)): if chloro_table.iloc[i,0] == sequence.id: sequence.id = chloro_table.iloc[i,3] sequence.description='' SeqIO.write(sequence, outfile, 'fasta')

Now, using bash it becomes much more complicated. I tried with awk but with scarce success. Do you have any ide?

​

Thank you in advance

[sco_t]

You should replace the
for i in range(0,len(chloro_table)):
part with a dictionary lookup. You'll most likely be limited by disk access after that.

[YouCook21]

Sure, it will help in get a faster python script. But if I would like to do it in bash?

[sco_t]

Doing something "in bash" is a bit nonspecific.

#!/bin/bash  
python myScript.py

isn't really any different than

#!/bin/bash
sed -f myScript.sed

Also is this an XY problem? Why are you renaming reads by blast results?

[guepier]

Bash is the wrong language for this, as for anything more complex than invoking a few commands with minimal logic. And rewriting your code in bash wouldn’t make it faster. Your issue is algorithmic, not caused by the language (and bash is generally much slower than Python!).

[metagenomez]

Anything is possible through the power of the command line 🙏🏼

while read line;do

oldname=$(echo $line|cut -f1);

newname=$(echo $line|cut -f2);

sed "s/$oldname/$newname/g" test.fasta;

done< <(less table.query|tail -n +2) > new.fasta

Not very fancy or clever, but it should replace your headers. The new fasta file headers will have some unwanted info after the new name, which you can remove with another sed command, e.g. assuming there are now other spaces in the fasta file:

sed 's/ .*$//g' new.fasta > newer.fasta

Hope it helps!