Question: extract same all similar sequences in FASTA based on the header
0
gravatar for akhilvbioinfo
3.1 years ago by
akhilvbioinfo140
India, chennai
akhilvbioinfo140 wrote:

Hello I downloaded more than 10k cds from NCBI and Headers look like this

lcl|DS180867.1_cds_EDM56551.1_2 [protein=hemolysin] [protein_id=EDM56551.1] [location=complement(<1..1012)]

lcl|DS180866.1_cds_EDM56342.1_3 [protein=phosphogluconate dehydratase] [protein_id=EDM56342.1] [location=complement(279..>893)]

lcl|DS180865.1_cds_EDM56120.1_4 [protein=3-hydroxyisobutyrate dehydrogenase] [frame=3] [protein_id=EDM56120.1] [location=complement(162..>559)]

lcl|DS180863.1_cds_EDM56465.1_5 [protein=hemolysin] [protein_id=EDM56465.1] [location=218..>977]

lcl|DS180862.1_cds_EDM56350.1_6 [protein=phosphogluconate dehydratase] [protein_id=EDM56350.1] [location=complement(<1..857)]

i want to extract same all similar sequences (based on the header) into a new fasta file based on protein name (example hemolysin)

please suggest any tool or programme for this

alignment sequence fasta • 1.4k views
ADD COMMENTlink modified 3.1 years ago by shenwei3565.0k • written 3.1 years ago by akhilvbioinfo140

Have you tried anything?

ADD REPLYlink written 3.1 years ago by venu6.3k

Yah i tried with AWK but its not working properly

ADD REPLYlink written 3.1 years ago by akhilvbioinfo140
2

can you post your awk code, so that people here can fix the problem in it. I feel that this would be a good idea instead of directly getting a working solution.

ADD REPLYlink written 3.1 years ago by venu6.3k

awk -F "|" '/^>/ {F = $3".fasta"} {print > F}' input.fasta>out.fasta

ADD REPLYlink written 3.1 years ago by akhilvbioinfo140

How do you plan to define similarity based on the header? By the [protein=...] tag?

ADD REPLYlink written 3.1 years ago by cschu1811.9k
1
gravatar for shenwei356
3.1 years ago by
shenwei3565.0k
China
shenwei3565.0k wrote:

Hi @akhilvbioinfo, remember SeqKit in your last post?

In which, we sort FASTA sequences by protein names? Here we can also split FASTA records by protein names again.

Same strategy: use flag --id-regexp to specify the protein name as the sequence identifier.

$ seqkit split --by-id --id-regexp "\[protein=(.+?)\]" seqs.fa
[INFO] split by ID. idRegexp: \[protein=(.+?)\]
[INFO] read sequences ...
[INFO] read 5 sequences
[INFO] write 2 sequences to file: seqs.fa.split/seqs.id_hemolysin.fa
[INFO] write 2 sequences to file: seqs.fa.split/seqs.id_phosphogluconate dehydratase.fa
[INFO] write 1 sequences to file: seqs.fa.split/seqs.id_3-hydroxyisobutyrate dehydrogenase.fa

Very simple, right?

See other 18 sub commands provided by SeqKit, you may need it again :)

ADD COMMENTlink modified 3.1 years ago • written 3.1 years ago by shenwei3565.0k

its very simple sir but i got an error like this

[ERRO] open seqs.fna.split/seqs.id_ketopantoate reductase PanE/ApbA family protein.fna: no such file or directory

ADD REPLYlink modified 3.1 years ago • written 3.1 years ago by akhilvbioinfo140

I changed headers using this comment ( Seq tool kit wont accept "/" characters)

sed -i -e 's/\//-/g'seqs.fna

Now its working perfectly

Thank you

ADD REPLYlink written 3.1 years ago by akhilvbioinfo140
1

it's a bug and already fixed since v0.3.8, please download the latest version.

ADD REPLYlink written 3.1 years ago by shenwei3565.0k
Please log in to add an answer.

Help
Access

Use of this site constitutes acceptance of our User Agreement and Privacy Policy.
Powered by Biostar version 2.3.0
Traffic: 1065 users visited in the last hour