Question

Renaming fasta headers according to a matching name list

7

Entering edit mode

8.7 years ago

san.san ▴ 190

Hi all,

Can anyone advise me on how to rename the header in my fasta file? Say I have a seq.fa file with transcript sequences:

>TR1|c0_g1_i1
GTCGAGCATGGTCTTGGTCATCTTCCTTTCAAAGAA
>TR6|c0_g1_i1
GTGGAATATCGCCAGTGACCATCACTGATTAACCTG

I also have a file with contigs matching the transcripts - names.txt:

TR1|c0_g1_i1    scaf0432344_50037.734_wgs
TR6|c0_g1_i1    scaf0159424_10142.072_wgs

How to I add contig names to the fasta file headers so that the "scaf0..." identifier comes before the "TR..."?

Desired output:

>scaf0432344_50037.734_wgs|TR1|c0_g1_i1
GTCGAGCATGGTCTTGGTCATCTTCCTTTCAAAGAA
>scaf0159424_10142.072_wgs|TR6|c0_g1_i1
GTGGAATATCGCCAGTGACCATCACTGATTAACCTG

Cheers!

text-processing command-line sequence • 23k views

ADD COMMENT • link updated 2.7 years ago by Priya ▴ 20 • written 8.7 years ago by san.san ▴ 190

4

Entering edit mode

There are a lot of posts regarding this issue. I'd suggest you to have a look before asking because I think that this topic has been widely covered:

replace fasta headers with another name in a text file

ADD REPLY • link 8.7 years ago by iraun 6.2k

0

Entering edit mode

I have looked for a couple of hours now but didn't manage to find anything matching exactly what I need to do. I'm not trained so things that are not exactly what I need, unfortunately, were not a help to me. Since I don't need to replace the header but rather modify it, this particular post wasn't helpful to me.

ADD REPLY • link 8.7 years ago by san.san ▴ 190

1

Entering edit mode

Did you try anything? If yes post what you've tried which is not working?

ADD REPLY • link 8.7 years ago by venu 7.1k

0

Entering edit mode

No, sorry, I have no clue how to go about it.

ADD REPLY • link 8.7 years ago by san.san ▴ 190

0

Entering edit mode

Your fasta headers look like Trinity output. Did you run Trinity on linux? Or were you given the fasta?

ADD REPLY • link 8.3 years ago by st.ph.n ★ 2.7k

0

Entering edit mode

I am trying to use above mentioned perl script and modified it and trying to get a desirable output, but I am unable to get it. My input fasta file looks like-

>V091201769_1 H648IQM04I6TDD orig_bc=AGCACTGTAG
GACGAACGCTGGCGGCGTGCTTAACACATGCAAGTCGAACGAGAGGAGAAGGAAGCTTGCCT
>V091209099_2 H648IQM04JTWVX orig_bc=ACGAGTGCGT
GATAAACGCTGGCGGCGCACATAAGACATGCAAGTCGAACGAACTTAACCATTAGTTTAC
>V091201769_3 H648IQM04IRZLN orig_bc=AGCACTGTAG
GATGAACGCTGGCGGCGTGCTTAACACATGCAAGTCGAACGAAGCACTTTATTTGATTTCCTTCGGGACT

and the text file-

V091200036 TR01555
V091205377 TR01556
V091201769 TR01557
V091209099 TR01558

The expected outcome which i need is-

>TR01557_1 H648IQM04I6TDD orig_bc=AGCACTGTAG
GACGAACGCTGGCGGCGTGCTTAACACATGCAAGTCGAACGAGAGGAGAAGGAAGCTTGCCT
>TR01558_2 H648IQM04JTWVX orig_bc=ACGAGTGCGT
GATAAACGCTGGCGGCGCACATAAGACATGCAAGTCGAACGAACTTAACCATTAGTTTAC
>TR01557_3 H648IQM04IRZLN orig_bc=AGCACTGTAG
GATGAACGCTGGCGGCGTGCTTAACACATGCAAGTCGAACGAAGCACTTTATTTGATTTCCTTCGGGACT

I need to rename the header name before first underscore, and keep the rest of the header name as such.

Any suggestions?

ADD REPLY • link 7.9 years ago by shashank.gupta ▴ 10

6

Entering edit mode

Use SeqKit:

$ seqkit replace -p "^(.+?)_(\d+)" -r '{kv}_$2' -k names.txt seq.fa

ADD REPLY • link 7.9 years ago by shenwei356 8.7k

1

Entering edit mode

CAN YOU PLZ EXPLAIN THIS COMMAND, giving explanation with parts of it's functions so that i can use it according to my task , bcz i need to just add unique set of strings that i created right after the accession number with an #. for better understanding i eed my files like this-

V091201769_1#TR01555 H648IQM04I6TDD orig_bc=AGCACTGTAG GACGAACGCTGGCGGCGTGCTTAACACATGCAAGTCGAACGAGAGGAGAAGGAAGCTTGCCT V091209099_2#TR01556 H648IQM04JTWVX orig_bc=ACGAGTGCGT GATAAACGCTGGCGGCGCACATAAGACATGCAAGTCGAACGAACTTAACCATTAGTTTAC V091201769_3#TR01557 H648IQM04IRZLN orig_bc=AGCACTGTAG GATGAACGCTGGCGGCGTGCTTAACACATGCAAGTCGAACGAAGCACTTTATTTGATTTCCTTCGGGACT

ADD REPLY • link 2.7 years ago by Priya ▴ 20

0

Entering edit mode

That works. Thanks!

ADD REPLY • link 7.9 years ago by shashank.gupta ▴ 10

0

Entering edit mode

seqkit is a life saver

ADD REPLY • link 4.5 years ago by 1769mkc ★ 1.2k

0

Entering edit mode

Dear Shenwei, As I'm a biologist who is totally new to any programming task would you please tell me what would be the Seqkit's command if I want to do the following.

My input fasta is like

>AFA46815.1 chromosomal replication initiator protein DnaA [Acetobacterium woodii DSM 1030]
MDTSFKKIWDAALEIIKPDISPTSFNTWFLKIKPINYVNNTYYFLSENEFEKGILESRHIPLITNALAEVTGKTGQVKIV
LKEEDAMIAQNADSPNVSITNVEEKISNYQNNSMNPKYTFDSFVIGENNRFAHAACVAVSEAPSERYNPLFIYGGVGLGK
THLMQAIGHYILSYAPHKKVVYVSCEKFTNDFIDAIQNKTNISFRNKYRSADILLIDDIQFLAKKEGTQEEFFHTFNTLH
QENKQIVISSDRPPKEIPTLEERLRSRFESGLITDIYAPNFETRIAIIRKKAESFSDEIPNEVLSFIADSIHSNIRELEG
ALTTVFAYSKLHNAPINLDSAKNALKDIFRKKDDIVITSEYIKEVTAKYFNITVEDMNSKKRTRSISLPRQVAMYITREI
TDLSLPKIGDEFGGRDHSTVIHACQKITEEMTINTDFKNLILRIEREINGK
>AFA46816.1 DNA polymerase III subunit beta [Acetobacterium woodii DSM 1030]
MKFNCSKESLMAALTIAQKAVSSKTTLTILEGILFSAADNKLLLRSTDLEIGVEVTIPAEVEKEGEIVLTASIIGELIRK
MSGSDIFFESNENHQIKIECLLSNFTLKGLSSEEFPAFPEIIDDHIFSIDAGVLKELIKGTLFSVATNENIPVLTGVKIE
LHADDINLIALDGYRLALKSGKIKNNLTEDLSVIIPSKSLSEINKVLSNYSGDVTVKFSKNQIFFEMDNTQFTSRLLEGD
FINYKQIIPVEKATQVKINRRLLLESSERAALLAREGKNNLIKMDFNQDQLILTSNADIGDVFEVIPIENSGESLKIAFN
SKFLIEALRVIEGDELMIDMTTSVGPGVLLPIEGNNFIYLILPVRIAEEN
>AFA46817.1 RNA-binding protein containing S4 domain [Acetobacterium woodii DSM 1030]
MQIIEINTEFIKIDQLLKYAGIVGNGSDVKFMILDGLIRVNGELCTQRNKKIRDGDLVEIEDYEPLQVKGITTTCLSKR

and the text file is:

222528058
222528059
222528060

The expected output is:

>gi|222528058|ref|AFA46815.1| chromosomal replication initiator protein DnaA [Acetobacterium woodii DSM 1030]
MDTSFKKIWDAALEIIKPDISPTSFNTWFLKIKPINYVNNTYYFLSENEFEKGILESRHIPLITNALAEVTGKTGQVKIV
LKEEDAMIAQNADSPNVSITNVEEKISNYQNNSMNPKYTFDSFVIGENNRFAHAACVAVSEAPSERYNPLFIYGGVGLGK
THLMQAIGHYILSYAPHKKVVYVSCEKFTNDFIDAIQNKTNISFRNKYRSADILLIDDIQFLAKKEGTQEEFFHTFNTLH
QENKQIVISSDRPPKEIPTLEERLRSRFESGLITDIYAPNFETRIAIIRKKAESFSDEIPNEVLSFIADSIHSNIRELEG
ALTTVFAYSKLHNAPINLDSAKNALKDIFRKKDDIVITSEYIKEVTAKYFNITVEDMNSKKRTRSISLPRQVAMYITREI
TDLSLPKIGDEFGGRDHSTVIHACQKITEEMTINTDFKNLILRIEREINGK
>gi|222528059|ref|AFA46816.1| chromosomal replication initiator protein DnaA [Acetobacterium woodii DSM 1030]
MKFNCSKESLMAALTIAQKAVSSKTTLTILEGILFSAADNKLLLRSTDLEIGVEVTIPAEVEKEGEIVLTASIIGELIRK
MSGSDIFFESNENHQIKIECLLSNFTLKGLSSEEFPAFPEIIDDHIFSIDAGVLKELIKGTLFSVATNENIPVLTGVKIE
LHADDINLIALDGYRLALKSGKIKNNLTEDLSVIIPSKSLSEINKVLSNYSGDVTVKFSKNQIFFEMDNTQFTSRLLEGD
FINYKQIIPVEKATQVKINRRLLLESSERAALLAREGKNNLIKMDFNQDQLILTSNADIGDVFEVIPIENSGESLKIAFN
SKFLIEALRVIEGDELMIDMTTSVGPGVLLPIEGNNFIYLILPVRIAEEN
>gi|222528060|ref|AFA46817.1| chromosomal replication initiator protein DnaA [Acetobacterium woodii DSM 1030]
MQIIEINTEFIKIDQLLKYAGIVGNGSDVKFMILDGLIRVNGELCTQRNKKIRDGDLVEIEDYEPLQVKGITTTCLSKR

I will greatly appreciate if you could help

ADD REPLY • link 6.6 years ago by na.poorin • 0

1

Entering edit mode

Firstly preparing the mapping file of accession and GI:

If you have Unix/Linux, it's simple:

paste <(seqkit seq -n -i seqs.fasta) gi.txt
AFA46815.1      222528058
AFA46816.1      222528059
AFA46817.1      222528060

If not, you may need help of csvtk which has windows version:

seqkit seq -n -i seqs.fasta | csvtk sample -H -t -p 1 -n > n-acc.tsv

csvtk sample -H -t -p 1 -n gi.txt > n-gi.tsv

csvtk join -H -t n-acc.tsv n-gi.tsv | csvtk -H -t cut -f 2,3 > acc2gi.tsv

Then you can replace AFA46815.1 with gi|222528058|ref|AFA46815.1|:

seqkit replace -k acc2gi.tsv -p "^(.+?) (.+)$" -r "gi|{kv}|ref|\$1| \$2" seqs.fasta > seqs.renamed.fasta

seqkit seq -n seqs.renamed.fasta 
gi|222528058|ref|AFA46815.1| chromosomal replication initiator protein DnaA [Acetobacterium woodii DSM 1030]
gi|222528059|ref|AFA46816.1| DNA polymerase III subunit beta [Acetobacterium woodii DSM 1030]
gi|222528060|ref|AFA46817.1| RNA-binding protein containing S4 domain [Acetobacterium woodii DSM 1030]

ADD REPLY • link 6.6 years ago by shenwei356 8.7k

0

Entering edit mode

Hi Shenwei

seqkit replace -k acc2gi.tsv -p "^(.+?) (.+)$" -r "gi|{kv}|ref|\$1| \$2" seqs.fasta > seqs.renamed.fasta

Can you modify this command for me?

This is my fasta header

>orf05188   draft|P. aeruginosa 908|Pseudomonas|orf05188  Scaffold_1  5479915 5478419  len=1497

My output header should look like

>geneA

Thanks alot

ADD REPLY • link updated 3.7 years ago by GenoMax 147k • written 3.7 years ago by Optimist ▴ 190

0

Entering edit mode

What is the relationship between current fasta header and expected one?

ADD REPLY • link 3.7 years ago by GenoMax 147k

0

Entering edit mode

I want to replace orf05188 with geneA

ADD REPLY • link 3.7 years ago by Optimist ▴ 190

score 5 · Answer 1 · 2016-03-14

5

Entering edit mode

8.7 years ago

iraun 6.2k

Using the solution proposed in replace fasta headers with another name in a text file by @Kenosis just with a little modification, you'll get the desired output:

  use strict;
  use warnings;

    my @arr;

    while (<>) {
        chomp;
        my @a = split /\t/;
        push @arr, $a[1] . "|" . $a[0] if length;
        last if eof; }

    while (<>) {
        print /^>/ ? ">" . shift(@arr) . "\n" : $_; }

Call the script in the following way: perl script.pl textFile fastaFile [>outFile]

I'd strongly suggest you to try to reach the solution yourself, otherwise you'll never learn. And it is recommended to show up your unsuccessful attempts to get the solution, otherwise people can think that you are asking just to avoid the effort of doing it yourself. Even if your attempts are not the best ones (all of us have been through that), I'd suggest you to explain to the community what you've tried.

ADD COMMENT • link 8.7 years ago by iraun 6.2k

1

Entering edit mode

Thank you so much. Unfortunately learning in my own time isn't an option since there's pressure for me to produce results asap. Which really isn't helping anyone in the long run, but I don't have a say in it. But I get your point and next time will explain that I'm not just trying to bum off someone.

ADD REPLY • link 8.7 years ago by san.san ▴ 190

0

Entering edit mode

I am trying to use above mentioned perl script and modified it and trying to get a desirable output, but I am unable to get it. My input fasta file looks like-

V091201769_1 H648IQM04I6TDD orig_bc=AGCACTGTAG GACGAACGCTGGCGGCGTGCTTAACACATGCAAGTCGAACGAGAGGAGAAGGAAGCTTGCCT V091209099_2 H648IQM04JTWVX orig_bc=ACGAGTGCGT GATAAACGCTGGCGGCGCACATAAGACATGCAAGTCGAACGAACTTAACCATTAGTTTAC V091201769_3 H648IQM04IRZLN orig_bc=AGCACTGTAG GATGAACGCTGGCGGCGTGCTTAACACATGCAAGTCGAACGAAGCACTTTATTTGATTTCCTTCGGGACT

and the text file-

V091200036 TR01555 V091205377 TR01556 V091201769 TR01557 V091209099 TR01558

The expected outcome which i need is-

TR01557_1 H648IQM04I6TDD orig_bc=AGCACTGTAG GACGAACGCTGGCGGCGTGCTTAACACATGCAAGTCGAACGAGAGGAGAAGGAAGCTTGCCT TR01558_2 H648IQM04JTWVX orig_bc=ACGAGTGCGT GATAAACGCTGGCGGCGCACATAAGACATGCAAGTCGAACGAACTTAACCATTAGTTTAC TR01557_3 H648IQM04IRZLN orig_bc=AGCACTGTAG GATGAACGCTGGCGGCGTGCTTAACACATGCAAGTCGAACGAAGCACTTTATTTGATTTCCTTCGGGACT

I need to rename the header name before first underscore, and keep the rest of the header name as such.

Any suggestions?

ADD REPLY • link 7.9 years ago by shashank.gupta ▴ 10

score 5 · Answer 2 · 2016-03-14

5

Entering edit mode

8.7 years ago

san.san ▴ 190

This is another solution that was proposed to me, if anyone finds it helpful:

awk 'FNR==NR{
  a[">"$1]=$2;next
}
$1 in a{
  sub(/>/,">"a[$1]"|",$1)
}1' names.txt seq.fa

ADD COMMENT • link 8.7 years ago by san.san ▴ 190

0

Entering edit mode

Good one-liner. Thanks

ADD REPLY • link 4.4 years ago by JMMM ▴ 10

0

Entering edit mode

It worked perfeclty fine .Thanks alot

ADD REPLY • link 3.3 years ago by serfraz.saad • 0

score 4 · Answer 3 · 2016-08-18

SeqKit supports this case now by subcommand replace, usage, download .

seqs:

$ cat seq.fa 
>TR1|c0_g1_i1
GTCGAGCATGGTCTTGGTCATCTTCCTTTCAAAGAA
>TR6|c0_g1_i1
GTGGAATATCGCCAGTGACCATCACTGATTAACCTG

names.txt should be tab-delimited:

$ cat names.txt 
TR1|c0_g1_i1    scaf0432344_50037.734_wgs
TR6|c0_g1_i1    scaf0159424_10142.072_wgs

replace

$ seqkit replace -p "(.+)" -r '{kv}|$1' -k names.txt seq.fa 
[INFO] read key-value file: names.txt
[INFO] 2 pairs of key-value loaded
>scaf0432344_50037.734_wgs|TR1|c0_g1_i1
GTCGAGCATGGTCTTGGTCATCTTCCTTTCAAAGAA
>scaf0159424_10142.072_wgs|TR6|c0_g1_i1
GTGGAATATCGCCAGTGACCATCACTGATTAACCTG