1
$\begingroup$

I have a fasta file, like this:

>TRINITY_DN100000_c1_g1::TRINITY_DN100000_c1_g1_i3::g.3039::m.3039 TRINITY_DN100000_c1_g1::TRINITY_DN100000_c1_g1_i3::g.3039  ORF type:complete len:100 (-) TRINITY_DN100000_c1_g1_i3:1027-1326(-)
MVWIKFRGLHRVLTSTPLVKSGKTPSQTWAFLDISVELIVFLFLNVHKSPMPHFKIYSEA
FSEEWSLLWLQYSRHLIQKPKPWQIKIELLHLCCCNRLC*
>TRINITY_DN100000_c1_g6::TRINITY_DN100000_c1_g6_i2::g.84365::m.84365 TRINITY_DN100000_c1_g6::TRINITY_DN100000_c1_g6_i2::g.84365  ORF type:complete len:112 (-) TRINITY_DN100000_c1_g6_i2:379-714(-)
MEMMQEIIPFAREMLSARPSKGTMKVYLVGGTFAVLGIVSGMVEAACSLFPEQEESTLTK
LMEDCLTVTAQNQEPQTFIPEDDEQDAEMEAKAKDLPMFRQRRMSFRAHAS*

if I want to only keep the second header, like this (the amino acid sequence keep unchanged). which command could i USE?

>TRINITY_DN100000_c1_g1_i3
MVWIKFRGLHRVLTSTPLVKSGKTPSQTWAFLDISVELIVFLFLNVHKSPMPHFKIYSEA
FSEEWSLLWLQYSRHLIQKPKPWQIKIELLHLCCCNRLC*
>TRINITY_DN100000_c1_g6_i2
MEMMQEIIPFAREMLSARPSKGTMKVYLVGGTFAVLGIVSGMVEAACSLFPEQEESTLTK
LMEDCLTVTAQNQEPQTFIPEDDEQDAEMEAKAKDLPMFRQRRMSFRAHAS*

I know a command

> sed 's/::.*//' input > output 

can only keep the first header >TRINITY_DN100000_c1_g1 if I want to keep the second header with the isoform information >TRINITY_DN100000_c1_g1_i3, how should I correct this command?

$\endgroup$
  • $\begingroup$ You can also do this via eFetch or Biopython directly without regex-ing $\endgroup$ – Michael G. Oct 7 '19 at 9:37
2
$\begingroup$

Here's one way:

$ sed -E 's/>[^:]*::([^:]*)::.*/>\1/' file.fa 
>TRINITY_DN100000_c1_g1_i3
MVWIKFRGLHRVLTSTPLVKSGKTPSQTWAFLDISVELIVFLFLNVHKSPMPHFKIYSEA
FSEEWSLLWLQYSRHLIQKPKPWQIKIELLHLCCCNRLC*
>TRINITY_DN100000_c1_g6_i2
MEMMQEIIPFAREMLSARPSKGTMKVYLVGGTFAVLGIVSGMVEAACSLFPEQEESTLTK
LMEDCLTVTAQNQEPQTFIPEDDEQDAEMEAKAKDLPMFRQRRMSFRAHAS*

The idea is to match a >, then as many non: as possible ([^:]*), and then capture the string between the first :: and the second ::. Finally, you replace the entire line with the captured string.

Once you're satisfied this does what you need, you can edit the original file with -i:

sed -i -E 's/>[^:]*::([^:]*)::.*/>\1/' file.fa 

Alternatively, set awk's field separator to :: and print the second field only on lines that start with >:

$ awk -F'::' '/^>/ ? $0=">"$2 : 1' file.fa 
>TRINITY_DN100000_c1_g1_i3
MVWIKFRGLHRVLTSTPLVKSGKTPSQTWAFLDISVELIVFLFLNVHKSPMPHFKIYSEA
FSEEWSLLWLQYSRHLIQKPKPWQIKIELLHLCCCNRLC*
>TRINITY_DN100000_c1_g6_i2
MEMMQEIIPFAREMLSARPSKGTMKVYLVGGTFAVLGIVSGMVEAACSLFPEQEESTLTK
LMEDCLTVTAQNQEPQTFIPEDDEQDAEMEAKAKDLPMFRQRRMSFRAHAS*
$\endgroup$
1
$\begingroup$

A 2 step workaround is:

  1. Get rid of everything after the piece of info you want, that is anything starting with :: and goes on with g:

    sed "s/::g.*//" fasta > fasta2

>TRINITY_DN100000_c1_g1::TRINITY_DN100000_c1_g1_i3

  1. Get rid of the "first header" as you call it, after the first step, it is rather straightforward:

    sed "s/>TRINITY_.*::/>/" fasta2 > fasta3

>TRINITY_DN100000_c1_g1_i3

EDIT:

@terdon showed that multiple sed commands can be used when separated by ; in GNU sed (default on Linux):

sed "s/::g.*//; s/>TRINITY_.*::/>/" fasta > modified_fasta
$\endgroup$
  • $\begingroup$ it worked!! thank you! $\endgroup$ – lili Oct 6 '19 at 22:46
  • 1
    $\begingroup$ There is no need for multiple steps: sed 's/::g.*//; s/>TRINITY_.*::/>/'. Or, if you don't have GNU sed: sed -e 's/::g.*//' -e 's/>TRINITY_.*::/>/'. You could even pipe one sed into another, that would at least avoid the temp files: sed 's/::g.*//' | sed 's/>TRINITY_.*::/>/'. $\endgroup$ – terdon Oct 7 '19 at 9:30

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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