This is a test version of Biostars. For the public version, visit https://www.biostars.org.
Any Script To Parse Fasta Headers?

I have a 47GB file to parse. The sequences are in following format:

>TSCS_00041 gene0EA_12345_rframe2_ORF
MLAATHYYKFAIRRLFPLLKDTICASYSISIKHHENFMALSNMPKIWEDVEVDGNNMQWTRFQTTPVMPVYFIAAGVFNLSFITNWNTKLLYRKDILPYMTFAYNVAKNIAWFLSHIRKTKITNHI
>TSCS_00044 gene0EA_12341_rframe2_ORF
MTICASYSISIKHHENFMAIKHHENFMALSNMPKIWEDV

I simply want to format this file like:

>TSCS_00041
MLAATHYYKFAIRRLFPLLKDTICASYSISIKHHENFMALSNMPKIWEDVEVDGNNMQWTRFQTTPVMPVYFIAAGVFNLSFITNWNTKLLYRKDILPYMTFAYNVAKNIAWFLSHIRKTKITNHI
>TSCS_00044 
MTICASYSISIKHHENFMAIKHHENFMALSNMPKIWEDV

Could anyone share the script

fasta next-gen perl
• 12,068 views • link
updated by Ram Ram • written by SK • SK

can this be done with cut only? the OP seems to shorten the fasta header not the other lines

• 0 views • link
written by Istvan Albert Istvan Albert

cut -d" " -f 1 will work as long as no spaces in sequence.

• 122 views • link
written by brentp brentp

and that's how homework is done ;)

• 0 views • link
written by Jorge Amigo Jorge Amigo

4 answers

Assuming that you do not have spaces in your sequences you can try that:

sed -e '/^>/ s/ .*//' mybigfile

Awk should work too (even if your sequences contain spaces):

awk '{print /^>/ ? $1 : $0}' mybigfile

and just for fun, a pure bash version:

while read l ; do echo "${l%% *}" ; done < mybigfile

but a simple cut command would be faster:

cut -d " " -f 1 mybigfile
• 211 views • link
written by Frédéric Mahé Frédéric Mahé

nice, can you explain the bash version?

• 0 views • link
written by brentp brentp

Sure, "${l%% *}" is a parameter expansion. It says remove the longest suffix string containing a space from the variable l. So, for all lines in mybigfile, remove anything right of a space (if any) and write the result on the standard output.

• 187 views • link
written by Frédéric Mahé Frédéric Mahé

cheers .

• 0 views • link
written by brentp brentp

Try the following:

use strict;
use warnings;

while (<>) {
    s/>\S+\K.+//;
    print;
}

Results:

>TSCS_00041
MLAATHYYKFAIRRLFPLLKDTICASYSISIKHHENFMALSNMPKIWEDVEVDGNNMQWTRFQTTPVMPVYFIAAGVFNLSFITNWNTKLLYRKDILPYMTFAYNVAKNIAWFLSHIRKTKITNHI
>TSCS_00044
MTICASYSISIKHHENFMAIKHHENFMALSNMPKIWEDV

Usage: perl script.pl inFile >outFile

As the fasta inFile is read line-by-line, the substituting regex will keep the > and the non-whitespace characters up to the first whitespace. The >outFile notation directs the printing to the file outFile.

Hope this helps!

• 179 views • link
written by Kenosis Kenosis

the equivalent perl-one-liner:

  perl -pe 's/>\S+\K.+//' < file.fasta > new_file.fasta
• 161 views • link
written by JC JC

Actually I am not sure why you always use -plane. -p and -n contradict each other, while -a enables autosplit, which causes overhead and is not necessary in this case. -l is useless here, too. The proper one should be perl -pe 's/^>(\S+).*/>$1/' or something similar. Sed also works.

• 158 views • link
written by lh3 lh3

You (as usual) are right, I removed the unnecessary flags. Thanks.

• 0 views • link
written by JC JC

Put the following code in a file called parse.py

import sys

for line in open(sys.argv[1]):
    if line.startswith(">"):
        line = line.split()
        print line[0]
    else:
        print line.strip()

Then assuming your file is named "myfile.fa", you type:

python parse.py myfile.fa
• 18 views • link
written by KCC KCC

if you remove "if line.startswith(">"):" from the code, the script will still work.

• 0 views • link
written by Geparada Geparada

Fair enough. A little more slick than I was going for, I guess. I guess I would argue the program as written is more robust, explicit and easier to tweak if you need something slightly different done later.

• 0 views • link
written by KCC KCC

If you keep the "if line.startswith('>'), you can remove the first strip() function. You do not need to remove the newline character (or any trailing spaces) since you only output the first item.

• 0 views • link
written by Frédéric Mahé Frédéric Mahé

Thanks. I modified it.

• 0 views • link
written by KCC KCC

This is good!

Can also use a substitution regex here w/a capture:

import sys
import re

for line in open(sys.argv[1]):
    print re.sub(r'(>\S+)\s+.+', r'\1', line) or line
• 0 views • link
written by Kenosis Kenosis

Nice. Much more suave than what I did. :)

• 0 views • link
written by KCC KCC

This is good!

Can also use a substutition w/capture here:

import sys, re

for line in open(sys.argv[1]):
    print re.sub(r'(>\S+)\s+.+', r'\1', line) or line,

Sorry, George, about deleting my earlier posting. I didn't know what I was doing...

• 239 views • link
written by Kenosis Kenosis

The Biopieces www.biopieces.org) solution:

read_fasta -i in.fna | split_vals -k SEQ_NAME -d ' ' | rename_keys -k SEQ_NAME_0,SEQ_NAME | write_fasta -o out.fna -x
• 0 views • link
written by Martin A Hansen Martin A Hansen

Log in to answer this question.

More posts like this