0

I checked the different threads but couldn't find a solution to my problem. I have a file containing string of letters A, T, C, G (representing DNA sequences) as follows:

$ cat species_gene

>genus_1_species_1
TTGGATATTGAGTTACATACATATGATTTGGGAATGGAGAACCGAGATAAGACTGAAGATCAAGTCACAATTGATTGTGC
TAATGCTGTCAAGAAATACAATGTCGGCATCAAATGTGCAACTATTACTCCTGATGAAAATCGAGTTGAGGAGTTCAAGC

>tribe_1_species_1
CTGGATATAGAGTTGCACACATACGATCTGGGAATGGAAAACCGCGACAAGACTGATGACCAAGTCACAATCGATTGCGC
CAACGCGATTAAAAAATACAACGTTGGAATTAAGTGCGCGACCATCACTCCCGACGAGAAGAGAGTCGAGGAATTCAAGC

>species_2
CTTGATATTGAACTACACACCTATGATTTGGGAATGGAATACCGCGATAAAACTGATGACCAAGTCACAATTGACTGCGC
KAATGCTATTAAGAAATACAATGTTGGTATTAAGTGCGCTACCATAACTCCAGATGAAAAAAGAGTTGAAGAATTCAAAC

I would like to count the number of letters in the strings following the header ">genus_1_species_1" until the next header. In this case, the length of the following strings, equal to 160:

TTGGATATTGAGTTACATACATATGATTTGGGAATGGAGAACCGAGATAAGACTGAAGATCAAGTCACAATTGATTGTGC TAATGCTGTCAAGAAATACAATGTCGGCATCAAATGTGCAACTATTACTCCTGATGAAAATCGAGTTGAGGAGTTCAAGC

I would like to use awk and use ">" as Record Separator and line-breaks as Field Separator , and subsequently count the number of fields (i. e. lines), starting from the second line of the record. I used the following command:

$ awk 'BEGIN {RS=">"} NR==2 {print length}' species_gene

I obtained 180 (the number of characters of the first record). The output I'd like to get is 160 (the length of the first record, starting from the second line). I tried the following one to discard the header:

$ awk 'BEGIN {RS=">"; FS="\n"} NR==2 {if (NF!=1) {length}}' species_gene

There is something I don't get here and I would be glad if someone could help me.

  • I edited your question in an attempt to make it clearer, but, since I don’t understand what you’re saying, I probably didn’t get it right. Please check it and fix it. Also, please clarify your question. What is a “sequence”? (Remember, you’re talking to computer people, not fellow geneticists.) … (Cont’d) – Scott Sep 13 '16 at 3:30
  • (Cont’d) … The fact that you have six lines with the same number of characters would make it impossible to figure out exactly what you’re asking for, even if you said what value you’re expecting for your example data, which you don’t.  (Do you want an output of 160?)  For example, are you saying that you want the number of characters on the lines between the first > line and the second one? … … … … … … … … … P.S. Using > as your shell prompt when you have >s in your data is really confusing. – Scott Sep 13 '16 at 3:37
1

IMHO it would be mre natural to use the blank line as record separator. In this 'paragraph mode', awk includes newlines as field separators, so in your case each line of the record becomes a field. You can then remove the first field by assigning the empty string to it; this has a nice side effect of re-evaluating the whole record with the given output field separator: by setting this to the empty string as well, the remaining fields are concatenated. You can then simply print the length of the record.

awk 'BEGIN {RS="";OFS=""} NR==1 {$1=""; print length}' species_gene

Testing with your given input

$ awk 'BEGIN {RS="";OFS=""} NR==1 {$1=""; print length}' species_gene
160
  • I don't want to use the blank line as record separator because I sometimes have files where there are no blank lines between the records. I used your command and with some editing I found this one to do the job: awk 'BEGIN {RS=">"; OFS=""} NR==2 {$1=""; print length}' species_gene – jibbah Sep 15 '16 at 5:38
1

length is not only a variable in Awk, but a function:

$ awk 'BEGIN {RS=">"} NR==2 {print length($2) + length($3)}' species_gene
160

That's it: add the lengths of the two fields, rather than using the length of the whole record.

Shorter:

$ awk 'BEGIN {RS=">"} NR==2 {print length($2 $3)}' species_gene
160

Equally long:

$ awk 'BEGIN {RS=">"} NR==2 && $0 = length($2 $3)' species_gene
160
0

A possible minor improvement over your improvement to steeldriver’s answer:

awk -vRS=">" -vOFS="" '$1=="genus_1_species_1" {$1=""; print length;}' species_gene

The -vRS=">" is equivalent to BEGIN {RS=">"}, and likewise with the -vOFS="".  My contribution is that this version lets you specify the header of the record whose length, so you don’t have to count records.

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.