2

The problem:

I have multiple text files (.fas), which look like this:

file1.fas:

>species1
AICGICVIAGIAIYIAAICG
>species2
AICGIVVYICAGAYICAGCG

file2.fas:

>species1
AIG
>species2
GCI

I'm interested in counting the number of characters in the second line (it's the same across all the species in one file as they are aligned.

My current one liner:

for i in *.fas; do echo -n "$i," && grep -m 1 -A 1 '>' $i | tail -n 1 | wc -c; done;

This works to an extent but is counting one number higher than the actual number of characters, as it is counting the newline character. How do I fix this so it only counts the number of characters excluding newline characters?

Current output:

file1.fas,21
file2.fas,4

Desired output:

file1.fas,20
file2.fas,3
  • Are your protein sequences always on single lines? – Kusalananda Aug 26 at 9:56
  • Ideally, I can fix them to be if not, but if there is a way of counting until the next '>' that would be better. – Biomage Aug 26 at 9:58
  • And it's always the second sequence that you want the total length of? Or is it always the last one, or the one for species2 regardless of where in the file it occurs? – Kusalananda Aug 26 at 10:00
  • All the sequences in one file will be the same length, I was using grep -A 1 as that's the first sequence in the file. Species2 is here just to provide an example that there are multiple species in one file (but will always be the same length as species1) – Biomage Aug 26 at 10:02
  • 2
    Have you considered using samtools to create an index of the file and then simply read the length from the index file? Alternatively, infoseq from EMBOSS or one of the many other existing tools for working with Fasta files? – Kusalananda Aug 26 at 10:03
1
  • | wc -l prints the number of lines.
  • | wc -c prints the number of chars, including newline.
  • | wc -lc prints both (number of lines first).

So you can simply subtract them: (write it instead of |wc -c)

| wc -lc | awk '{print $2 - $1}'

If you only print your sequence on one line, you can substract 1 instead of the number of newlines.

Or you can use awk only, match the whole line and count its chars:

| awk '{match("[A-Z]*");print RLENGTH}'

RLENGTH is the length of the match (here, it's the whole line). Here, I assume you use only capital letters, else, use . instead of [A-Z].

On vim (text editor, so it could be incompatible with your script), visually select your area (one line or more) and run:

:'<,'>s/[A-Z]*//gn

You can also delete the newline characters with tr (also works for several lines) :

| tr -d '\n' | wc -c

Btw, there is probably lots of other ways to do it.

| improve this answer | |
  • wrt match("[A-Z]*") ITYM match($0,/[[:upper:]]+/) – Ed Morton Aug 26 at 17:35
1

Use find to find the .fas files and find command will push them onto awk's comnand line, as much as awk can take at a stretch. So the number of awk invocations is minimized.

$ find . -maxdepth 1 -type f -name '*.fas' -exec \
awk -v OFS=, 'FNR==2{print FILENAME, length();nextfile}' {} +

Output:

./file1.fas,20
./file2.fas,3
| improve this answer | |
  • 1
    Thanks, fixed now. find is necessary to prevent seeing too many files globbing error. – Rakesh Sharma Aug 26 at 18:00
  • Yeah, that's always true but there's no indication in the question that the OP has that situation and IMHO it doesn't make sense to use find... when writing code to open multiple files just in case it's one of those very rare situations where it ends up being needed. We can always just add it later if/when needed. – Ed Morton Aug 26 at 18:04
0

Use the following perl one-liner, which prints the file name and the length of the second line, comma-delimited. The newlines are stripped using -l flag:

for i in *.fas; do perl -lne 'if ( $. == 2) { print join ",", $ARGV, length $_; last LINE; }' $i ; done

Input:

cat > file1.fas <<EOF
>species1
AICGICVIAGIAIYIAAICG
>species2
AICGIVVYICAGAYICAGCG
EOF

cat > file2.fas <<EOF
>species1
AIG
>species2
GCI
EOF

Output:

file1.fas,20
file2.fas,3

The perl one-liner uses these command line flags:
-e : tells Perl to look for code in-line, instead of in a file.
-n : loop over the input one line at a time, assigning it to $_ by default.
-l : strip the input line separator ("\n" on *NIX by default) before executing the code in-line, and append it when printing.

$ARGV is the name of the current input file.
length $_ is the length of the current line, which has been selected as input line number 2 ($. == 2).
last LINE makes the code skip to the last input line after printing the line length, for speed.

SEE ALSO:

perlrun: command line switches

| improve this answer | |

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.