3

I have a dataset that I am working with made up of phylip files that I have been editing. Phylip format is a bioinformatics format that contains as a header the number of samples and the sequence length, followed by each sample and its sequence. for example:

5 10
sample_1 gaatatccga
sample_2 gaatatccga
sample_3 gaatatcgca
sample_4 caatatccga
sample_5 gaataagcga

My issue is that in trimming these datasets, the sample number in the header no longer is accurate (e.g. in above example might say five, but I've since trimmed to have only three samples). What I need to do is to replace that sample count with the new, accurate sample count but I cannot figure out how to do so without losing the sequence length number (e.g. the 10).

I have 550 files so simply doing this by hand is not an option. I can for-loop the wc but again I need to retain that sequence length information and somehow combine it with a new, accurate wc.

New contributor
erikusrex is a new contributor to this site. Take care in asking for clarification, commenting, and answering. Check out our Code of Conduct.
  • 3
    Please don't post pictures of text but post the actual text in a code block. – Jesse_b May 15 at 17:49
  • Also, it's unclear how the number should be changed. Always to 3? – choroba May 15 at 17:55
  • ok thank you, will do in the future. Not always to three as the sample number is different across files but most of the files have been edited and so the new number of samples is often fewer than what is currently stated in the header – erikusrex May 15 at 18:18
3

If I understand your requirement correctly you can use the following awk command:

awk -v samples="$(($(grep -c . input)-1))" 'NR == 1 { $1=samples }1' input

samples will be set to the number of lines in the input file minus one (since you aren't counting the header line).

awk will then change the first column of the first line to the new sample number and print everything.


$ cat input
5 10
sample_1 gaatatccga
sample_2 gaatatccga
sample_3 gaatatccga
$ awk -v samples="$(($(grep -c . input)-1))" 'NR == 1 { $1=samples }1' input
3 10
sample_1 gaatatccga
sample_2 gaatatccga
sample_3 gaatatccga

With GNU awk you can use the -i flag to modify the files in place but I would prefer to make a second set of modified files to ensure the correct changes have been made.

Something like:

for file in *.phy; do
    awk -v samples="$(($(grep -c . "$file")-1))" 'NR == 1 { $1=samples }1' "$file" > "${file}.new"
done
3

Another option would be to use ed (of course!):

for f in input*
do 
  printf '1s/[[:digit:]][[:digit:]]*/%d\nw\nq' $(( $(wc -l < "$f") - 1 )) | ed -s "$f"
done

This loops over the files (named, for example input-something) and sends a simple ed-script to ed:

  • on line 1, search and replace (s//) one or more digits at the beginning of the line with another number -- that replacement number being the result of computing the line length of the input minus one
  • after that, w write the file out and
  • then q quit ed
1

In Vim, run:

:execute '1s/^[0-9]\+/' . (line('$')-1) . '/'

(Thanks also to this answer for pointing me in the right direction.)

You can also do this in a loop, e.g. using :bufdo or just a shell for loop.

Your Answer

erikusrex is a new contributor. Be nice, and check out our Code of Conduct.

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.