I found this interesting command:

grep -v '^>' test.fasta | tr -d '\n' | sed -e 's/\(.\)/\1\n/g' | sort | uniq -c | sort -rn

I have some grasp what it means (it counts letters from a text file), but my question is about this:

sed -e 's/\(.\)/\1\n/g'

I know that it is compose of three substitute commands. One is to substitute new lines (\n), one that matches any characters except newlines (\(.\)), but I am lost at /\1\?

  • 1
    It is fairly rough. The tr takes out all newlines, so the sed stdin is one long line (excluding the row headers, which start with >). Fasta files can be rather large (I see a topic in BioStars that quotes 54GB). sed is a line editor, so it reads the whole input into memory (as one line) before it processes it, at which point it doubles the size for the output (because every input char gets its own newline). A simple script (in Ruby or awk) would avoid six processes, five pipes, and two sorts (one of them with a line for every character). Oct 21 at 12:15
  • 1
    Your grep flags exclude lines starting with >, correct? So this grep code excludes the first (descriptive ) line in a .fasta file, and (I presume), the rest of the code is designed to count occurrences of individual nucleotides (or amino acids). Oct 21 at 22:55

The command

sed -e 's/\(.\)/\1\n/g'

is a single GNU sed substitution command that replaces every character with itself, followed by a newline character. The effect of this is to fold the input into a single column of single characters.

$ echo hello | sed -e 's/\(.\)/\1\n/g'

The \(.\) is a "capture group", capturing a single character. The \1 is a "back-reference" to the first capture group. Using \1 in the replacement text would insert whatever was captured by the first parentheses.

It could also be written without so many backslashes as

sed 's/./&\n/g'

where & simply means "whatever was matched by the expression".

The sed command requires GNU sed as standard sed can't insert newlines with \n like that.

To do that more efficiently with standard tools, use

fold -w 1

instead. This is more efficient as no regular expression matching is needed for each character in the input.

Using fold, your pipeline could be written

grep -v '^>' file | tr -d '\n' | fold -w 1 | sort | uniq -c | sort -rn

Alternatively, using awk to get rid of a few steps of that pipeline,

awk '!/^>/ { for (i = 1; i <= length; ++i) count[substr($0,i,1)]++ }
    END { for (ch in count) print count[ch], ch }' file |
sort -rn

The awk code counts the number of times each character has been seen. It does that by incrementing the value in the array count corresponding to each character in the input stream. At the end of input, a summary of the counts and characters counted are outputted.

  • For GNU awk 4.0 up you can eliminate the sort -rn by setting PROCINFO["sorted_in"]="@val_num_desc" Oct 22 at 4:26
  • 1
    @PauloSergioSchlogl I saw your edit but rejected it. The tr command is still needed if the input contains two (or more) newlines in a row.
    – they
    Oct 22 at 9:12

I hope this makes it clearer.

"I know that it is compose of three substitute commands"

It's just one substitute command (if you are referring to the sed command): s/<pattern to search>/<replacement>/, that will execute the following:

  • For every line search for <pattern> and substitute it with <replacement>.
  • The g flag means do it globally, because by default sed replaces only the first occurrence of <pattern>.

"but I am lost at /\1\"

You can capture a pattern by surroundig it with escaped parenthesis \(<pattern>\), or just parenthesis with the -E option, (<pattern>).

In the <replacement> section this captured pattern is referenced by backslash and a number, \<number>. The number refers to the position of the capture, since you can have several:

sed -E '/(<first capture>)(<second capture>)/\1\2/'

So the command sed -e 's/\(.\)/\1\n/g' means:

  • Capture every character \(.\) and replace it by itself and a new line \1\n.
  • With g, do it globally, don't stop at the first occurrence.

For example:

$ echo foo | sed -E 's/(.)/\1\n/g'

The -e options is not necessary here, unless you are concatenating several sed commands: sed -e '...' -e '...', etc.

You can find more information at Back-references and Subexpressions.

  • 1
    Functionally what the replacement is doing is putting each character on its own line, finding a single character and replacing it with that character plus \n. Oct 21 at 12:03
  • @JeffBreadner thanks, I'm still editing because I have power failures and I don't want to lose the changes. Oct 21 at 12:04

Using Raku (formerly known as Perl_6)

raku -e 'for lines.grep({ !/ ^ \> / }).join { .say for .comb.Bag.sort(*.values).reverse};'

Sample Input:

>sp|P01308|INS_HUMAN Insulin OS=Homo sapiens OX=9606 GN=INS PE=1 SV=1

Sample Output:

L => 20
G => 12
A => 10
E => 8
Q => 7
P => 6
C => 6
V => 6
R => 5
S => 5
Y => 4
F => 3
T => 3
N => 3
M => 2
D => 2
K => 2
I => 2
W => 2
H => 2

The code you present can be written in a number of languages (not only sed), any one of which may strike a chord with you. As an example, above your code has been re-written in Raku, a member of the Perl-family of languages.

Most of the Raku code should be fairly self-explanatory: lines are read in, and grep-ped for the ! absence of a ^ start-of-line > 'greater-than' angle, and join-ed. Joined lines are comb-ed (broken into individual characters), Bag-ged (each character present becomes a key and occurrences are counted/recorded as values), sort-ed in reverse to place highest # of occurrences first, then printed with say.



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.