This is a test version of Biostars. For the public version, visit https://www.biostars.org.
How to limit fasta header to 40 characters?

I have FASTA headers with long annotation names, but the program it will be run through for proteomics has a limit of roughly 40 characters or else it will crash.

The file starts off like this:

>TRINITY_DN0_c0_g1_i1.p1 - RecName: Full=E3 ubiquitin-protein ligase CIP8; AltName: Full=COP1-interacting protein 8; AltName: Full=RING-type E3 ubiquitin transferase CIP8
SEQUENCE

>TRINITY_DN10003_c0_g1_i12.p1 - RecName: Full=Polycomb group protein FIE1; AltName: Full=Protein  FERTILIZATION-INDEPENDENT ENDOSPERM 2; Short=OsFIE2; AltName: Full=WD40 repeat-containing protein 153; Short=OsWD40-153
SEQUENCE

Ideally I want the FASTA to look like this:

>DN0_c0_g1_i1.p1 - E3 ubiquitin-protein 
SEQUENCE

>DN10003_c0_g1_i12.p1 - Polycomb group
SEQUENCE

I used sed and seqkit to cut out the repetitive parts

sed 's/>.*Y_/>/' proteome.fasta
seqkit replace -p " RecName: Full=" -r ' ' proteome.fasta > proteome2.fasta

The fasta looks like this now:

>DN0_c0_g1_i1.p1 - E3 ubiquitin-protein ligase CIP8; AltName: Full=COP1-interacting protein 8; AltName: Full=RING-type E3 ubiquitin transferase CIP8
SEQUENCE

>DN10003_c0_g1_i12.p1 - Polycomb group protein FIE1; AltName: Full=Protein  FERTILIZATION-INDEPENDENT ENDOSPERM 2; Short=OsFIE2; AltName: Full=WD40 repeat-containing protein 153; Short=OsWD40-153
SEQUENCE

What can I do to limit the header length? Can I do it with seqkit?

unix seqkit fasta
• 4,885 views • link
updated by jena • jena • written by fishingline37 • fishingline37

How about this (starting with your last example file)

$ awk '(">" ~ $0); {print substr($0,0,40)}' fasta_file
>DN0_c0_g1_i1.p1 - E3 ubiquitin-protein 
SEQUENCE
SEQUENCE2


>DN10003_c0_g1_i12.p1 - Polycomb group p
SEQUENCE
• 0 views • link
written by GenoMax GenoMax

OP also wants TRINITY_ removed so a sed might be required before the awk. However, given that OP has already figured out the sed, maybe use proteome.fasta instead of fasta_file in your code so OP knows not to replace their entire code with your awk.

• 0 views • link
written by Ram Ram

I used the last example that OP showed above. fasta_file is just a place holder for file name.

• 0 views • link
written by GenoMax GenoMax

Only keep the sequence identifiers:

seqkit seq -i proteome.fasta -o proteome2.fasta

Where,

-i, --only-id                   print IDs instead of full headers

Or

$ seqkit replace -p ' - RecName: Full=(.+?);.+' -r ' $1' proteome.fasta  
>TRINITY_DN0_c0_g1_i1.p1 E3 ubiquitin-protein ligase CIP8
SEQUENCE
>TRINITY_DN10003_c0_g1_i12.p1 Polycomb group protein FIE1
SEQUENCE

But the number of the characters may exceed the limit.

• 0 views • link
written by shenwei356 shenwei356

This is a neat use of seqkit, I will definitely keep this in mind for future projects.

• 0 views • link
written by fishingline37 • fishingline37

Perhaps you should also check that after trimming the names you don't get duplicate IDs. These two commands should give the same output (not checked):

grep '>' before_trimming.fasta | sort | uniq | wc -l
grep '>' after_trimming.fasta | sort | uniq | wc -l
• 0 views • link
written by dariober dariober

3 answers

You can do sth like this-

cat proteome.fasta
>TRINITY_DN0_c0_g1_i1.p1 - RecName: Full=E3 ubiquitin-protein ligase CIP8; AltName: Full=COP1-interacting protein 8; AltName: Full=RING-type E3 ubiquitin transferase CIP8
SEQUENCE

>TRINITY_DN10003_c0_g1_i12.p1 - RecName: Full=Polycomb group protein FIE1; AltName: Full=Protein  FERTILIZATION-INDEPENDENT ENDOSPERM 2; Short=OsFIE2; AltName: Full=WD40 repeat-containing protein 153; Short=OsWD40-153
SEQUENCE

cut -d ' ' -f1-5 proteome.fasta |sed 's/RecName: Full=//g' |sed 's/TRINITY_//g'
>DN0_c0_g1_i1.p1 - E3 ubiquitin-protein
SEQUENCE

>DN10003_c0_g1_i12.p1 - Polycomb group
SEQUENCE
• 0 views • link
written by bk11 • bk11

This won't work if the non-blank-space separated words are too long. Given that the 40 char limit is the primary concern, the cut should be something like -c 1-40 and it should come after the sed operations.

• 0 views • link
written by Ram Ram

Please check this too if it works.

sed 's/TRINITY_//g' proteome.fasta |sed 's/RecName: Full=//g' |awk '/^>/{print substr($0, 1, 39)} !/^>/{print}'
• 0 views • link
written by bk11 • bk11

This is the most complete and concise answer. The two seds instead of using seqkit is great.

The awk for the character limit also works like a charm. If you make this the main post, I will accept it as the answer. Thanks for the nice solution!

• 0 views • link
written by fishingline37 • fishingline37

Ok. I am making it as the main post. Please go ahead and accept it as answer. Thanks!

• 0 views • link
written by bk11 • bk11

Please check this, it will work-

sed 's/TRINITY_//g' proteome.fasta |sed 's/RecName: Full=//g' |awk '/^>/{print substr($0, 1, 39)} !/^>/{print}'
• 0 views • link
written by bk11 • bk11

Looks like a job for bioawk. For example, if you just need the IDs (like in the seqkit answer above), you can do this:

bioawk -c fastx '{print ">"$name; print $seq}' test.fa
>DN0_c0_g1_i1.p1
SEQUENCE
>DN10003_c0_g1_i12.p1
SEQUENCE

This works because bioawk -c fastx parses FASTA/FASTQ into a table of 3 or 4 fields for you (FASTA is missing $qual of course). From bioawk -c help you can see the column names:

...
fastx:
    1:name 2:seq 3:qual 4:comment

You can also use $comment in substr or gsub functions to extract other info and add it after your ID ($name). Comments are separated from names by the first whitespace in the header.

• 0 views • link
written by jena • jena

Log in to answer this question.

More posts like this
  • Blast command line
    written by doinelpierrot • doinelpierrot
    Hello all, * I would like to run a blastx on a fasta file in order to perform a quick functionnal annotation on swissprot. * …