4
$\begingroup$

I have a FASTA header that looks like this:

>7c8250ef-c89f-4d42-9d48-12c8fe245fb2 runid=606f271fc97598006ba5a922136a2c304cef75a5 sampleid=Pool12-1 read=19008 ch=301 start_time=2021-07-03T08:48:18Z barcode=barcode01

And I am able to change it to the desired output here:

>7c8250ef-c89f-4d42-9d48-12c8fe245fb2_001

Using this command:

grep '^>' 001_old.fasta | cut -d ' ' -f 1,8 | sed '/^>/s/$/_001/' > 001_new.fasta

However it completely deletes the DNA sequence for that read and only outputs the names of the headers to the new fasta file:

>7c8250ef-c89f-4d42-9d48-12c8fe245fb2_001
>815440d5-e033-446c-b7dd-5fad4105184a_001
>d36a0fea-101d-4bd9-bee4-adec7fe4d43c_001
>38ca9eb8-bbd2-4d0a-ab06-79d134ad46f1_001
...

I think it has something to do with the cut command I'm using. I tried just doing 1 but that gave the same output so I figured if I included the 8th position of text after a space it would work but again no dice. If someone could help me to include the DNA sequence under the FASTA header in this new file that would be great.

$\endgroup$
1
  • $\begingroup$ Do yourself a favor and use Bipython for parsing and modifying FASTA files. $\endgroup$
    – Supertech
    Jun 29 at 1:03

2 Answers 2

5
$\begingroup$

The grep, cut, sed pipeline could also be replaced with the following command:

awk '/^>/ { print $1 "_001"; next }1' file.fa

The general syntax of an AWK command is condition { action } where condition specifies a condition or regular expression, and action specifies a block of code or set of commands to be run when the condition is met. Both condition and action are optional. Here, we use a regular expression to match lines beginning with a ">" symbol, and print the first whitespace-separated field followed by "_001". We then use next to skip the remaining actions and move to the next line, which avoids any further processing of the current line. The 1 on the end is actually another condition that always evaluates to true. Since there's no action associated with this condition, AWK performs the default action which is to print the entire line. This effectively prints all lines that do not begin with a ">" symbol.

$\endgroup$
1
  • 1
    $\begingroup$ thank you this worked perfectly! and thank you for the explanation too $\endgroup$
    – rimo
    Jun 28 at 18:54
3
$\begingroup$

The problem is in the grep command you're using, which only sends the fasta file headers to the output. Here's my guess at what the commands in your pipeline do:

  1. grep '^>' 001_old.fasta - extract lines starting with '>' from 001_old.fasta

  2. cut -d ' ' -f 1,8 - assume file contains space-delimited records; extract only the first and eighth field from each line

  3. sed '/^>/s/$/_001/' - in lines that begin with a '>', replace the end of the line with '_001'

  4. > 001_new.fasta - write result out to 001_new.fasta

Of all of these, I think only the sed and output parts are important, but it needs to be slightly modified to work as intended:

sed '/^>/s/ .*$/_001/' 001_old.fasta > 001_new.fasta

Assuming I've got the syntax right, this should search for lines beginning with '>', and replace everything from the first space up to the end of the line with '_001'.

I'm more familiar with perl, so here's a perl one-liner that should do the same thing:

perl -pe '/^>/ && s/ .*$/_001/' 001_old.fasta > 001_new.fasta
$\endgroup$
1
  • 1
    $\begingroup$ thank you so much! this makes a lot of sense to me as well. $\endgroup$
    – rimo
    Jun 28 at 18:54

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge that you have read and understand our privacy policy and code of conduct.

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