I have multifasta sequence and i am trying to rename the header and print the first id only
>gene_1|id=abc|protein=capA|ABD5672.4
ATGCGTCGACGTCGTACGGGTTTT
CGTACGGGTTATGCGTCGACGTCG
TACGGGTTTT
desired output:
>gene_1
ATGCGTCGACGTCGTACGGGTTTT
CGTACGGGTTATGCGTCGACGTC
GTACGGGTTTT
I tried following using sed but didnt work it seems sed -r 's/^([^|]+\|){0}//;s/\|//'any improvement on this? or any other alternative ?
•
written
by
####
8 answers
there is no 4th Id you have only 2 here
1 = gene_1
and
2 = id=abc;protein=capA
also your sequence is multiline format meanwhile the command above assumes that is is one line in this case you can linearize it first
then
cat gene.fasta | awk '{ split($1,a,"|");print a[4] "\n" $2}' > result.fa
but this will be too much work the above solution is much better
Linearizing the input is not necessary. You can simply awk for lines that are sequence headers — lines which start with > — and those lines which are not headers (sequences). You can then modify the header lines by whatever logic you want and print out the modified header. In the case of non-header lines — sequences — you just print out the line as-is.
Looking for a rigorous solution?
As @genomax2 said, "Fasta format is not a strictly defined standard". So you need to check FASTA format very carefully.
Here's a fake but may existing FASTA file. Note that the last record gene_6 is appending to
sequence of gene_5, this often happens by cating too FASTA files where the
first one is not ending with new-line-character \n.
Similarly, concatenating files could result in the space in
front of >gene_5.
$ cat seqs.fa
>gene_1|id=abc;protein=capA
ATGCGTCGACGTCGTACGGGTTTT
CGTACGGGTTATGCGTCGACGT
>gene_2 id=abc;protein=capA
ATGCGTCGACGTCGTACGGGTTTT
CGTACGGGTTATGCGTCGAC
>gene_3
ATGCGTCGACGTCGTACGGGTTTT
CGTACGGGTTATGCGTCG
>gene_4|
ATGCGTCGACGTCGTACGGGTTTT
CGTACGGGTTATGCGT
>gene_5|
ATGCGTCGACGTCGTACGGGTTTT
CGTACGGGTTATGC>gene_6|id=123;protein=fake
ACTG
Let's do this step by step.
Counting sequences with different methods
$ grep -c '>' seqs.fa 6 $ grep -c '^>' seqs.fa 4 $ seqtk seq seqs.fa | grep -c '^>' 4 $ seqkit stat seqs.fa file format type num_seqs sum_len min_len avg_len max_len seqs.fa FASTA DNA 4 251 42 62.8 119Seems some bad formats are there. Let's check the sequences:
$ seqkit seq -v seqs.fa [ERRO] error when parsing seq: gene_4| (seq: invalid DNAredundant letter: >)Retrieving
gene_4, there they are!!!$ seqkit grep -p "gene_4|" seqs.fa >gene_4| ATGCGTCGACGTCGTACGGGTTTTCGTACGGGTTATGCGT >gene_5| ATGCGTCGAC GTCGTACGGGTTTTCGTACGGGTTATGC>gene_6|id=123;protein=fakeACTGChecking
gene_5andgene_6in original file$ grep "gene_5|" -A 1 -B 1 seqs.fa CGTACGGGTTATGCGT >gene_5| ATGCGTCGACGTCGTACGGGTTTT $ grep "gene_6|" -A 1 -B 1 seqs.fa ATGCGTCGACGTCGTACGGGTTTT CGTACGGGTTATGC>gene_6|id=123;protein=fake ACTGFix it:
$ sed -i 's/ >gene_5/>gene_5/; s/>gene_6/\n>gene_6/' seqs.faValidate again
$ seqkit seq -v seqs.fa | seqkit stat file format type num_seqs sum_len min_len avg_len max_len - FASTA DNA 6 214 4 35.7 46Format is OK now, you may want to remove gaps in sequences.
$ seqkit seq -g seqs.fa| seqkit stat file format type num_seqs sum_len min_len avg_len max_len - FASTA DNA 6 214 4 35.7 46Well, no gaps.
Finaly, we can handle the sequence IDs.
$ seqkit seq -n -i seqs.fa gene_1|id=abc;protein=capA gene_2 gene_3 gene_4| gene_5| gene_6|id=123;protein=fakeOh no! some headers do not have
|. The Methods our friends kindlly provide may not work.Let's use a general way provided by seqkit, i.e., using regular expression to retrieve sequence id.
$ seqkit seq -n -i seqs.fa --id-regexp '^([^| ]+)[\| ]?' gene_1 gene_2 gene_3 gene_4 gene_5 gene_6^([^| ]+)[\| ]?captures leading non-space and non-|characters.$ seqkit seq -i seqs.fa --id-regexp '^([^| ]+)[\| ]?' >gene_1 ATGCGTCGACGTCGTACGGGTTTTCGTACGGGTTATGCGTCGACGT >gene_2 ATGCGTCGACGTCGTACGGGTTTTCGTACGGGTTATGCGTCGAC >gene_3 ATGCGTCGACGTCGTACGGGTTTTCGTACGGGTTATGCGTCG >gene_4 ATGCGTCGACGTCGTACGGGTTTTCGTACGGGTTATGCGT >gene_5 ATGCGTCGACGTCGTACGGGTTTTCGTACGGGTTATGC >gene_6 ACTG
Learn more about the SeqKit.
An introductory post on Biostars: SeqKit: a cross-platform and ultrafast toolkit for FASTA/Q file manipulation in Golang
$ pip install pyfaidx
$ faidx -e "lambda x: x.split('|')[0]" genes.fa
>gene_1
ATGCGTCGACGTCGTACGGGTTTT
CGTACGGGTTATGCGTCGACGTC
GTACGGGTTTT
...
You can just use awk directly. No need to use cat:
$ awk '{ if ($0 ~ /^>/) { n = split($0, a, "|"); printf("%s\n", a[1]); } else { print $0; } }' in.fa > out.fa
Really no need to use slow parsers, either, until you do something more complex. This is a very basic text-editing exercise that all informaticists should be able to do.
You don't need cat with the sed solutions either, if you want to use inplace editing...my answer could be re-written sed -i 's/|.*//' filename but I figured simply printing the output avoids the OP making mistakes and irreversibly overwriting their files till they are sure what it'll do ;)
Using sed works too, definitely. I tend to shy away from sed because of differences between BSD and GNU versions (but that's solved on OS X with Homebrew, I guess).
Do you think so? :P
I actually find a cat - sed - redirect, while not computationally efficient maybe because of an extra system call, to be a bit more elegant because you don't have to create extra messy intermediate files and you can test/build the command as you go.
Personal preference though!
A million solutions :p
True. I'm usually building pipelines as I go, and I start off with cat file. Then I hit up to go to the previous command, and I append a sed to it. It is easier building that way.
But once I build it the first time, I optimize it for later re-use. At this point, all the unnecessary cats are thrown out.
Same, and piping the end of the pipeline to a head to look how it's doing (if the file is big).
Even simpler, just 2 calls one after the other, but your way would totally work too - I wrapped it in a slightly more complicated if so that I could pass a number of lines to print if I wanted more than default 10
ends(){
if [ "$2" == "" ]; then
head $1
tail $1
elif [ "$2" != "" ] ; then
head -"${2}" $1
tail -"${2}" $1
else
echo "Parameter not recognised"
fi
}
** It's not written as a oneliner in bashrc, I've just copied it though, so there are some missing ";" as a oneliner FYI.
Please just use a proper FASTA parser. awk/grep/sed cannot parse FASTA reliably. Why do people even bother writing FASTA/FASTQ parsers if no one is going to use them? Why do people bother making custom bioinformatic file formats if people are going to treat them like a weird CSV? Convert all your data to CSV if you want CSV. But if you want to stick with FASTA - for whatever reason - use a FASTA parsing API like BioPython/Perl/Ruby/etc.
Everything else is loaded with assumptions about the rest of your data which we cannot possible know.
For complex input data, or situations where you don't know your own files that well I very much agree. However, for working with files on the commandline that don't require you to have installed Biopython (for example, can be a right pain in the ass to install), I think knowing the necessary oneliners with sed/cut/grep/awk and so on is an absolute requirement...
Is using pip that painful? :-p
My comment above was meant to be in jest. I have modified it to more accurately reflect that.
By using a standard file parser one is able to account for even edge cases and get the desired output. You are not able to post your entire file (due to space constraints) and the command line solutions people suggest consider only those snippets. Fasta format is not a strictly defined standard so if your file does not follow the same exact structure throughout then you could end up with silent corruption of data that you may not discover until some other downstream tool throws an error.
Please be careful on how you comment ...
I'd word that differently - it comes across as too aggressive/confrontational.
And it's not about logic, it's about using a robust and scalable solution. I'd use cut or awk or sed only after ensuring all known edge cases, and even then I'm risking unknown/unpredictable edge cases. An actively maintained tool would bypass that.
That's OK Ram - i'm not high value i'm just chatty :) I didn't get my points on Biostars like you, Genomax, Devon or Pierre by being helpful - I just made lots of posts :P More over, I think changing your opinion is like losing a friend. Denial, Anger, Bargaining, Depression, Acceptance. Anger is an unhelpful but very important part of that. If people aren't getting angry, they're not confronting their own assumptions. So all i'm saying is that while being rude/insulting is never OK, being grumpy is fine by me. And nothing will make me less chatty... unfortunately ;) haha
John, it's about credibility within a group. I've found that no matter how great you are in a group, when you move out, you have to start rebuilding that reputation.
I believe in that philosophy, and that's how I look at users - not by points, but by how credible they are and how much effort goes into their input.
With OP, for example, it is hard to imagine someone not googling "edit FASTA header" or "parse FASTA". OP's expertise is in NGS analysis, which needs a lot of experience with UNIX, and OP has difficulties with a simple sed here. That against your post history is a stark contrast in credibility, and that's why I made the comment on earning the right to criticize. Critiquing in great, but criticizing needs the person to have some weight attached to their words.
This is my personal opinion, and I might be wrong, so please treat it with all relevant disclaimers.
I'm really uncomfortable with attacking users. Can we please either answer the question or move on? Thanks.
I feel like cleaning up the off topic comments on the answer of John, but they're mixed with good/helpful bits.
If this keeps up, I might start doing some cleaning, myself.
I cleaned up the unpleasant off-topic discussion which originated from this reply.
A summary of the useful bits of replies on this is below. I apologise for any context that has been lost by this operation.
jrj.healey
For complex input data, or situations where you don't know your own files that well I very much agree. However, for working with files on the commandline that don't require you to have installed Biopython (for example, can be a right pain in the ass to install), I think knowing the necessary oneliners with sed/cut/grep/awk and so on is an absolute requirement...
the point he's getting at is that you shouldn't necessarily treat a fasta as a file that can be coerced in to something which is simple a delimited text file, when parsers already exist, which is what a couple of the suggestions in this thread do.
genomax2
By using a standard file parser one is able to account for even edge cases and get the desired output. You are not able to post your entire file (due to space constraints) and the command line solutions people suggest consider only those snippets. Fasta format is not a strictly defined standard so if your file does not follow the same exact structure throughout then you could end up with silent corruption of data that you may not discover until some other downstream tool throws an error.
Ram
I'd use cut or awk or sed only after ensuring all known edge cases, and even then I'm risking unknown/unpredictable edge cases. AN actively maintained tool would bypass that.
Why one would use cut or awk with CSV files without running a few checks beforehand) when they are insensitive to quoted commas is beyond me :) Also, encoding, quoted new lines, vertical tabs, etc. Takeaway: Not even CSV file are as simple as they seem.
John
That's OK Ram - i'm not high value i'm just chatty
Thank you for keeping Biostars a nice place for all users!
Log in to answer this question.
-
Lifting over from Hg37 to Hg38 coordinates doesn't work
written by Jamie Watson •
-
Problem renaming headers in multi fasta file using certain scripts
written by lukhanyomakhabane •
-
How to edit fasta headers to keep only ID and organism?
written by A_heath •
-
Rename fasta headers using a variable already present in the header?
written by FreshBio •
-
Reformat sam header issue
written by pg_canada •
-
Best way to processing huge fasta files using python.
written by elisheva
-
Trim FASTA header in multifasta file
written by j.m.carsten •
-
I want to make reference table so that I can annotate the antibiotic resistance gene hits with anti…
written by ghanbari.msc •
-
Adapt perl script
written by weissert •
-
[Fastacmd] How To Retrieve Sequence From Blast Db
written by Leszek
Second ID is not gene_1. Do you want
abcinstead?No in this case I want first id ..i.e. gene_1
Why'd you say "second id" in the question then?
Hi ####
It's perfectly possible that the discussion in this thread didn't make you happy, but deleting the thread isn't the right way to cope with that. I've reversed that now. There are many interesting ideas and solutions, so there is no reason to wipe this all away "forever".
Cheers, Wouter
Hello ####!
We believe that this post does not fit the main topic of this site.
Please close this post as it is solved ..It's too much now
For this reason we have closed your question. This allows us to keep the site focused on the topics that the community can help with.
If you disagree please tell us why in a reply below, we'll be happy to talk about it.
Cheers!