1
$\begingroup$

I'm using efetch to dig out a middle "envelope" protein within a viral genome specifically within fasta format. The code below works but may not be robust for an entire database. Thus I want to replace the two "perl" lines with formal code within the Entrez package. I tried using efilter -query but couldn't get it to work. Genome positions are one approach but are not absolute from one virus to another. Rewriting this may be the best option.

ncbi2=ACF49259
efetch -db protein -format gpc -id "$ncbi2" |\
xtract -pattern INSDSeq -element INSDQualifier_value |\
perl -p -e 's/^(.+?)\s+taxon.*envelope protein.*\s+\b([A-Z]{300}[A-Z]+).*nonstructural protein 1.*/>$1\n$2/'|\
perl -p -e "s/[\t \;\/]+/_/g"

Output

>Dengue_virus_1_VR1254_1_US_Hawaii_1944...
MRCVGIGNRDFVEGLSGATWVDVVL.....

I will provide a more detailed explanation shortly. Any ideas welcome.

$\endgroup$
  • 1
    $\begingroup$ In the meantime, you could streamline the perl a little to | perl -pe 's/^(.+?)\s+taxon.*envelope protein.*\s+\b([A-Z]{300}[A-Z]+).*nonstructural protein 1.*/>$1\n$2/; s/[\t ;\/]+/_/g. There's no need to process the file twice and you also don't need to escape the ; inside a character class. $\endgroup$ – terdon May 14 at 12:33
  • $\begingroup$ Hi @terdon, yes and thanks now you've mentioned it, it's obvious. $\endgroup$ – Michael G. May 14 at 13:16
  • 1
    $\begingroup$ Not exactly what you were looking for but may be you can use something like this? efetch -db protein -format gpc -id "$ncbi2" | xtract -pattern INSDSeq_feature-table -group INSDFeature -if INSDFeature_key -equals 'mat_peptide' -tab '\n' -element INSDFeature_location,INSDQualifier_value | grep 'envelope' $\endgroup$ – vkkodali May 15 at 0:56
  • $\begingroup$ Perfect it works a treat @vkkodali. Huge thanks. A small amount of regex will finish the fasta format. Maybe you would like to provide a formal answer so I can sign this question as complete? $\endgroup$ – Michael G. May 15 at 11:48
2
$\begingroup$

Not exactly what you were looking for but may be you can use something like this?

ncbi2=ACF49259
efetch -db protein -format gpc -id "$ncbi2" \
    | xtract -pattern INSDSeq_feature-table \
        -group INSDFeature \
        -if INSDFeature_key -equals 'mat_peptide' \
        -tab '\n' \
        -element INSDFeature_location,INSDQualifier_value \
    | grep 'envelope'

For explorations of the INSDSet XMLs, xtract has -insd option that you can use. In this particular case, you can do something like this:

efetch -db nuccore -id "$ncbi2" -format gpc \
    | xtract -insd mat_peptide product peptide

For more information, checkout EDirect documentation and some example queries

$\endgroup$
  • $\begingroup$ Thank you its perfect. Stacked together in a small subroutine and parsed with a bit of regex its fantastic. "esearch -query" is used to obtain the accession numbers. $\endgroup$ – Michael G. May 16 at 12:22

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.