2
$\begingroup$

I have a large fasta file with DNA sequences. I need to add a ; at the end of each header line (lines beginning with >) in the file.

For example, I have:

>Bacteria;Cyanobacteriales;Phormidiaceae;Lyngbya PCC-7419;Lyngbya sp. PCC 8106
GGAGAGUUUGAUCCUGGCUCAGGAUGAACGCUGGCGGUCUGCUUAACACAUGCAAGUCGAACGGACUCUUCGGAGUUAGU
>Bacteria;Cyanobacteriales;Nostocaceae;Nodularia PCC-9350;Nodularia spumigena CCY9414
GAGAGUUUGAUCCUGGCUCAGGAUGAACGCUGGCGGUAUGCUUAACACAUGCAAGUCGAACGGUCUCUUCGGAGAUAGUG

And, my desired output is:

>Bacteria;Cyanobacteriales;Phormidiaceae;Lyngbya PCC-7419;Lyngbya sp. PCC 8106;
GGAGAGUUUGAUCCUGGCUCAGGAUGAACGCUGGCGGUCUGCUUAACACAUGCAAGUCGAACGGACUCUUCGGAGUUAGU
>Bacteria;Cyanobacteriales;Nostocaceae;Nodularia PCC-9350;Nodularia spumigena CCY9414;
GAGAGUUUGAUCCUGGCUCAGGAUGAACGCUGGCGGUAUGCUUAACACAUGCAAGUCGAACGGUCUCUUCGGAGAUAGUG
$\endgroup$
1
  • $\begingroup$ Those sequences have U in them; they're RNA sequences, not DNA sequences $\endgroup$
    – gringer
    Apr 6 at 0:47

3 Answers 3

6
$\begingroup$

... using the legendary Perl pie

perl -p -i -e 's/^(>.*)/$1;/' mybacteria.fa

Input

>Bacteria;Cyanobacteriales;Phormidiaceae;Lyngbya PCC-7419;Lyngbya sp. PCC 8106
GGAGAGUUUGAUCCUGGCUCAGGAUGAACGCUGGCGGUCUGCUUAACACAUGCAAGUCGAACGGACUCUUCGGAGUUAGU
>Bacteria;Cyanobacteriales;Nostocaceae;Nodularia PCC-9350;Nodularia spumigena CCY9414
GAGAGUUUGAUCCUGGCUCAGGAUGAACGCUGGCGGUAUGCUUAACACAUGCAAGUCGAACGGUCUCUUCGGAGAUAGUG
A

Output

>Bacteria;Cyanobacteriales;Phormidiaceae;Lyngbya PCC-7419;Lyngbya sp. PCC 8106;
GGAGAGUUUGAUCCUGGCUCAGGAUGAACGCUGGCGGUCUGCUUAACACAUGCAAGUCGAACGGACUCUUCGGAGUUAGU
>Bacteria;Cyanobacteriales;Nostocaceae;Nodularia PCC-9350;Nodularia spumigena CCY9414;
GAGAGUUUGAUCCUGGCUCAGGAUGAACGCUGGCGGUAUGCUUAACACAUGCAAGUCGAACGGUCUCUUCGGAGAUAGUG

Details I copy pasted your first code into a file mybacteria.fa then ran the Perl pie script. Perl pie is very quick, there alternative regexs (regular expressions) that would be slightly more efficient, but Perl pie is fast regardless. In fact you don't need to use regex in this example and would be faster via a if $line == '>'; declaration (certainly true for Pythonic equivalents if line == '>') but it doesn't offer the ease and flexibility, which is combined with speed, that Perl pie does. The code will work in Linux (any flavour/flavor) or Mac OSX Darwin (i.e. Terminal app)

$\endgroup$
4
  • 1
    $\begingroup$ It will work on any flavour/flavor so long as Perl is installed. Note that not all distributions ship with Perl pre-installed. Also, the capture group is not necessary - all that's needed here is: perl -pe 's/$/;/ if /^>/' in.fa > out.fa. Note that using -i without an extension will modify the input file in-place without making a backup of the original. $\endgroup$
    – Steve
    Apr 6 at 2:50
  • $\begingroup$ In fairness I did point that out in the discussion below the code. The efficiencies are generally not noticeable. However, most regex manipulations take place within the string rather than at the end requested here, and capture can be one way to deal with different situations beyond this one. $\endgroup$
    – M__
    Apr 6 at 11:34
  • $\begingroup$ Fair point. Perl is included in most popular Linux versions including Ubuntu, Debian and Suse. Perl isn't preinstalled on RedHat or freeBSD after release 5 and nobody uses Kali in bioinformatics. Let me know if there's anything missed. $\endgroup$
    – M__
    Apr 6 at 12:15
  • 1
    $\begingroup$ Thank you very much -- this worked well and was simple to execute. $\endgroup$ Apr 7 at 0:46
6
$\begingroup$

One way, using and a regex address to select only the header lines and apply a substitution to append a semicolon to the end of each line:

sed '/^>/s/$/;/' in.fa > out.fa

Or, to edit the file in place:

sed -i '/^>/s/$/;/' in.fa

Or one way using and a ternary construct to print the line followed by a semicolon or an empty string, depending on whether or not the line is a header line:

awk '{ print $0 (/^>/ ? ";" : "") }' in.fa > out.fa

Alternatively, you may prefer:

awk '/^>/ { print $0 ";"; next } { print $0 }' in.fa > out.fa

Which is the same as:

awk '/^>/ { print $0 ";"; next }1' in.fa > out.fa

If you have installed, all you need is:

perl -pe 's/$/;/ if /^>/' in.fa > out.fa
$\endgroup$
3
$\begingroup$

I recommend to use Biopython as a clearer way.

from Bio import SeqIO

records = []
for record in SeqIO.parse("file.fas", "fasta"):
    record.description += ";"
    records.append(record)
SeqIO.write(records, "file.fas", "fasta")

Before:

>OK087550.1 Procyon lotor isolate Plo1 cytochrome b gene, partial cds; mitochondrial
TCAAACATCTCAGCATGATGAAATTTCGGCTCCCTCCTCGGAATTTGTTTGCTTCTACAGATCGCAACAG
GTTTATTCTTAGCCATGCACTATACACCAGATACAACCACAGCTTTCTCATCAGTGACCCACATTTGCCG
AGATGTAAATTATGGCTGAATTATCCGATATATACACGCTAACGGAGCTTCTATATTCTTTATATGCCTA
TTCTTACACGTAGGACGAGGCTTATACTATGGCTCCTATACATTCTCTGAAACATGAAATATTGGTGTTA
TTTTACTACTCACAGTAATAGCAACCGCTTTCATAGGTTACGTCCTCCCATGAGGACAAATATCATTCTG
AGGAGCAACAGTAATTACAAATCTCCTGTCAGCTATCCCCTATATCGGATCTAACCTTGTAGAATGAATT
TGAGGAGGGTTTTCAGTAGACAAAGCCACCCTAACACGATTCTTCGCATTCCACTTCATTCTACCATTTA
TTATCACAACGCTAGCAATAATTCACCTGCTATTCCTACACGAAACAGGATCCAATAATCCTGCTGGAAT
TACATCAGAATCTGACAAAATTCCATTTCACCCATACTACACCATTAAAGACATTCTGGGAATCCTATTC
CTTATTTTTATACTTATAGGTTTAGTGCTATTTACGCCAGACCTACTAGGAGACCCAGATAACTACACAC
CCGCCAACCCCTTAAGCACCCCACCTCACATTAAACCTGAATGATATTTTCTATTCGCCTACGCAATTCT
ACGTTCCATTCCCAACAAACTAGGAGGAGTCCTAGCCCTAGTCCTCTCCATCTTAATCCTAATCATCATC
CCACTCCTACATACCTCAAAACAACGAAGCATAATATTTCGGCCACTTAGCCAATGTTCATTCTGATTCC
TAGTAGCAGACCTCCTCGTCCTAACATGAATTGGAGGTCAACCAGTAGAATATCCCTTCATCATCATCGG
CCAACTAGCCTCCATCTTCTACTTTATAATCCTCCTGGTCCT

>OK087549.1 Procyon lotor isolate Plo17 cytochrome b gene, partial cds; mitochondrial
TCAAACATCTCAGCATGATGAAATTTCGGCTCCCTCCTCGGAATTTGTTTGCTTCTACAGATCGCAACAG
GTTTATTCTTAGCCATGCACTATACACCAGATACAACCACAGCTTTCTCATCAGTGACCCACATTTGCCG
AGATGTAAATTATGGCTGAATTATCCGATATATACACGCTAACGGAGCTTCTATATTCTTTATATGCCTA
TTCTTACACGTAGGACGAGGCTTATACTATGGCTCCTATACATTCTCTGAAACATGAAATATTGGTGTTA
TTTTACTACTCACAGTAATAGCAACCGCTTTCATAGGTTACGTCCTCCCATGAGGACAAATATCATTCTG
AGGAGCAACAGTAATTACAAATCTCCTGTCAGCTATCCCCTATATCGGATCTAACCTTGTAGAATGAATT
TGAGGAGGGTTTTCAGTAGACAAAGCCACCCTAACACGATTCTTCGCATTCCACTTCATTCTACCATTTA
TTATCACAACGCTAGCAATAATTCACCTGCTATTCCTACACGAAACAGGATCCAATAATCCTGCTGGAAT
TACATCAGAATCTGACAAAATTCCATTTCACCCATACTACACCATTAAAGACATTCTGGGAATCCTATTC
CTTATTTTTATACTTATAGGTTTAGTGCTATTTACGCCAGACCTACTAGGAGACCCAGATAACTACACAC
CCGCCAACCCCTTAAGCACCCCACCTCACATTAAACCTGAATGATATTTTCTATTCGCCTACGCAATTCT
ACGTTCCATTCCCAACAAACTAGGAGGAGTCCTAGCCCTAGTCCTCTCCATCTTAATCCTAATCATCATC
CCACTCCTACATACCTCAAAACAACGAAGCATAATATTTCGGCCACTTAGCCAATGTTCATTCTGATTCC
TAGTAGCAGACCTCCTCGTCCTAACATGAATTGGAGGTCAACCAGTAGAATATCCCTTCATCATCATCGG
CCAACTAGCCTCCATCTTCTACTTTATAATCCTCCTGGTCCT

>OK087548.1 Procyon lotor isolate Plo11 cytochrome b gene, partial cds; mitochondrial
TCAAACATCTCAGCATGATGAAATTTCGGCTCCCTCCTTGGAATTTGTTTGCTTCTACAGATCGCAACAG
GTTTATTCTTAGCCATGCACTATACACCAGATACAACCACAGCTTTCTCATCAGTGACCCACATTTGCCG
AGATGTAAATTATGGCTGAATTATCCGATATATACACGCTAACGGAGCTTCTATATTCTTTATATGCCTA
TTCTTACACGTAGGACGAGGCTTATACTATGGCTCCTATACATTCTCTGAAACATGAAATATTGGTGTTA
TTTTACTACTCACAGTAATAGCAACCGCTTTCATAGGTTACGTCCTCCCATGAGGACAAATATCATTCTG
AGGAGCAACAGTAATTACAAATCTCCTGTCAGCTATCCCCTATATCGGATCTAACCTTGTAGAATGAATT
TGAGGAGGGTTTTCAGTAGACAAAGCCACCCTAACACGATTCTTCGCATTCCACTTCATTCTACCATTTA

After:

>OK087550.1 Procyon lotor isolate Plo1 cytochrome b gene, partial cds; mitochondrial;
TCAAACATCTCAGCATGATGAAATTTCGGCTCCCTCCTCGGAATTTGTTTGCTTCTACAG
ATCGCAACAGGTTTATTCTTAGCCATGCACTATACACCAGATACAACCACAGCTTTCTCA
TCAGTGACCCACATTTGCCGAGATGTAAATTATGGCTGAATTATCCGATATATACACGCT
AACGGAGCTTCTATATTCTTTATATGCCTATTCTTACACGTAGGACGAGGCTTATACTAT
GGCTCCTATACATTCTCTGAAACATGAAATATTGGTGTTATTTTACTACTCACAGTAATA
GCAACCGCTTTCATAGGTTACGTCCTCCCATGAGGACAAATATCATTCTGAGGAGCAACA
GTAATTACAAATCTCCTGTCAGCTATCCCCTATATCGGATCTAACCTTGTAGAATGAATT
TGAGGAGGGTTTTCAGTAGACAAAGCCACCCTAACACGATTCTTCGCATTCCACTTCATT
CTACCATTTATTATCACAACGCTAGCAATAATTCACCTGCTATTCCTACACGAAACAGGA
TCCAATAATCCTGCTGGAATTACATCAGAATCTGACAAAATTCCATTTCACCCATACTAC
ACCATTAAAGACATTCTGGGAATCCTATTCCTTATTTTTATACTTATAGGTTTAGTGCTA
TTTACGCCAGACCTACTAGGAGACCCAGATAACTACACACCCGCCAACCCCTTAAGCACC
CCACCTCACATTAAACCTGAATGATATTTTCTATTCGCCTACGCAATTCTACGTTCCATT
CCCAACAAACTAGGAGGAGTCCTAGCCCTAGTCCTCTCCATCTTAATCCTAATCATCATC
CCACTCCTACATACCTCAAAACAACGAAGCATAATATTTCGGCCACTTAGCCAATGTTCA
TTCTGATTCCTAGTAGCAGACCTCCTCGTCCTAACATGAATTGGAGGTCAACCAGTAGAA
TATCCCTTCATCATCATCGGCCAACTAGCCTCCATCTTCTACTTTATAATCCTCCTGGTC
CT
>OK087549.1 Procyon lotor isolate Plo17 cytochrome b gene, partial cds; mitochondrial;
TCAAACATCTCAGCATGATGAAATTTCGGCTCCCTCCTCGGAATTTGTTTGCTTCTACAG
ATCGCAACAGGTTTATTCTTAGCCATGCACTATACACCAGATACAACCACAGCTTTCTCA
TCAGTGACCCACATTTGCCGAGATGTAAATTATGGCTGAATTATCCGATATATACACGCT
AACGGAGCTTCTATATTCTTTATATGCCTATTCTTACACGTAGGACGAGGCTTATACTAT
GGCTCCTATACATTCTCTGAAACATGAAATATTGGTGTTATTTTACTACTCACAGTAATA
GCAACCGCTTTCATAGGTTACGTCCTCCCATGAGGACAAATATCATTCTGAGGAGCAACA
GTAATTACAAATCTCCTGTCAGCTATCCCCTATATCGGATCTAACCTTGTAGAATGAATT
TGAGGAGGGTTTTCAGTAGACAAAGCCACCCTAACACGATTCTTCGCATTCCACTTCATT
CTACCATTTATTATCACAACGCTAGCAATAATTCACCTGCTATTCCTACACGAAACAGGA
TCCAATAATCCTGCTGGAATTACATCAGAATCTGACAAAATTCCATTTCACCCATACTAC
ACCATTAAAGACATTCTGGGAATCCTATTCCTTATTTTTATACTTATAGGTTTAGTGCTA
TTTACGCCAGACCTACTAGGAGACCCAGATAACTACACACCCGCCAACCCCTTAAGCACC
CCACCTCACATTAAACCTGAATGATATTTTCTATTCGCCTACGCAATTCTACGTTCCATT
CCCAACAAACTAGGAGGAGTCCTAGCCCTAGTCCTCTCCATCTTAATCCTAATCATCATC
CCACTCCTACATACCTCAAAACAACGAAGCATAATATTTCGGCCACTTAGCCAATGTTCA
TTCTGATTCCTAGTAGCAGACCTCCTCGTCCTAACATGAATTGGAGGTCAACCAGTAGAA
TATCCCTTCATCATCATCGGCCAACTAGCCTCCATCTTCTACTTTATAATCCTCCTGGTC
CT
>OK087548.1 Procyon lotor isolate Plo11 cytochrome b gene, partial cds; mitochondrial;
TCAAACATCTCAGCATGATGAAATTTCGGCTCCCTCCTTGGAATTTGTTTGCTTCTACAG
ATCGCAACAGGTTTATTCTTAGCCATGCACTATACACCAGATACAACCACAGCTTTCTCA
TCAGTGACCCACATTTGCCGAGATGTAAATTATGGCTGAATTATCCGATATATACACGCT
AACGGAGCTTCTATATTCTTTATATGCCTATTCTTACACGTAGGACGAGGCTTATACTAT
GGCTCCTATACATTCTCTGAAACATGAAATATTGGTGTTATTTTACTACTCACAGTAATA
GCAACCGCTTTCATAGGTTACGTCCTCCCATGAGGACAAATATCATTCTGAGGAGCAACA
$\endgroup$
10
  • 1
    $\begingroup$ Clearer than what? Why is using a complex package requiring several lines of code "clearer" than using standard *nix tools and one liners? $\endgroup$
    – terdon
    Apr 6 at 8:35
  • 1
    $\begingroup$ @terdon For the same reasons fasta, nexus, &c. are explicit. It is easier to read at a glance. $\endgroup$
    – Vovin
    Apr 6 at 9:29
  • 2
    $\begingroup$ @terdon As someone who has taught programming to juniors I would bet a significant amount of money that a larger percentage of 'people in IT' will correctly "decode" what this does compared to any regex based solution. The strength of regexes is their flexibility, not their cleanness or readability. Appending a character at the end of a line really isn't playing to the strengths of regexes. $\endgroup$ Apr 6 at 12:30
  • 2
    $\begingroup$ Steve's point deserves emphasis. I love Biopython but its implementation of the id versus description attributes on SeqRecord objects never fails to trip me up. $\endgroup$
    – Jesse
    Apr 6 at 14:41
  • 2
    $\begingroup$ This script loads all sequences into memory, which wouldn't work for huge input. $\endgroup$
    – user172818
    Apr 6 at 23:21

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.