4

I have a sequence file:

$ cat file
CACCGTTGCCAAACAATG
TTAGAAGCCTGTCAGCCT
CATTGCTCTCAGACCCAC
GATGTACGTCACATTAGA
ACACGGAATCTGCTTTTT
CAGAATTCCCAAAGATGG

I want to calculate the longest stretch of C+T. I could only count total C+T, but I want the longest stretch.

$ cat file | awk '{ print $0, gsub(/[cCtT]/,"",$1)}'
CACCGTTGCCAAACAATG 9
TTAGAAGCCTGTCAGCCT 10
CATTGCTCTCAGACCCAC 12
GATGTACGTCACATTAGA 8
ACACGGAATCTGCTTTTT 11
CAGAATTCCCAAAGATGG 7

The Expected result would be to show the longest C+T stretch.

CACCGTTGCCAAACAATG 9 2
TTAGAAGCCTGTCAGCCT 10 3
CATTGCTCTCAGACCCAC 12 5
GATGTACGTCACATTAGA 8 2
ACACGGAATCTGCTTTTT 11 6
CAGAATTCCCAAAGATGG 7 5
5

FWIW here's a way to do it in perl, using max from List::Util

$ perl -MList::Util=max -lpe '$_ .= " " . max 0, map length, /[CT]+/gi' file
CACCGTTGCCAAACAATG 2
TTAGAAGCCTGTCAGCCT 3
CATTGCTCTCAGACCCAC 5
GATGTACGTCACATTAGA 2
ACACGGAATCTGCTTTTT 6
CAGAATTCCCAAAGATGG 5
4
$ awk '{ split($0, a, "[^CTct]+"); m=0
         for (i in a) {
             len=length(a[i])
             if (len > m) m=len
         }
         print $0, m  }' file
CACCGTTGCCAAACAATG 2
TTAGAAGCCTGTCAGCCT 3
CATTGCTCTCAGACCCAC 5
GATGTACGTCACATTAGA 2
ACACGGAATCTGCTTTTT 6
CAGAATTCCCAAAGATGG 5

This awk program splits each line on runs of anything that is not upper or lower case C or T. It then loops over the bits that the split results in, finding the longest one. It then prints the original line together with the found maximum length.

Since Roman looked at timings for various solutions, here's a quicker solution:

awk -F "[^CTct]+" '
    m = 0
    for (i = 1; i <= NF; ++i) {
        len = length($i)
        if (len > m) m = len
    }
    print m' file | paste file -

It's quicker since it only splits the lines once. The first code additionally tries to split the input lines on whitespaces.

Timing it using mawk shows 0.79s on 500000 lines. The first solution uses 1.69s on the same data, showing that it's probably the split operation that takes the most time.

  • Should be <= NF. Is that wall clock time or CPU time? That's significant as the two process variant would be impacted by the number of CPU cores. – Stéphane Chazelas Jun 29 '18 at 15:26
  • @StéphaneChazelas Wall clock. Dual processor system. Hyperthreading disabled. – Kusalananda Jun 29 '18 at 15:34
  • In the above example, can you please suggest how to print the starting and end point of the longest stretch within the given sequence – CN_229133 Oct 8 '18 at 22:24
3

With sed (assuming no more than 19 characters per line), just for the fun of it and using the greedy properties of RE matching:

sed '
  h;y/cCtT/xxxx/;x;H;s/./x/g;G
  s/^\(x*\).*\n.*\1.*\n/\1 /
  s/^x\{10\}/1/;s/$/:9876543210xxxxxxxxx/
  s/^\(1*\)\(x*\) \(.*\):.*\(.\).\{9\}\2$/\3 \1\4/'

A variation on @Kusalananda's solution:

awk -F '[^cCtT]+' '
  {
    max = 0
    for (i = 1; i <= NF; i++)
      if ((l = length($i)) > max)
        max = l
    print $0, max
  }'
3

Faster GNU awk solution:

awk -v FPAT='[ctCT]+' \
'{ 
     max_l = t_len = 0;
     for (i=1; i <= NF; i++) {
         len = length($i);
         if (len > max_l) max_l = len;
         t_len += len
     }
     print $0, t_len, max_l
 }' inputfile

The output:

CACCGTTGCCAAACAATG 9 2
TTAGAAGCCTGTCAGCCT 10 3
CATTGCTCTCAGACCCAC 12 5
GATGTACGTCACATTAGA 8 2
ACACGGAATCTGCTTTTT 11 6
CAGAATTCCCAAAGATGG 7 5

Time performance comparison (test inputfile has about 120000 lines):

$ time awk -v FPAT='[ctCT]+' '{ max_l = t_len = 0; for (i=1; i <= NF; i++) { len = length($i); if (len > max_l) max_l = len; t_len += len } print $0, t_len, max_l }' inputfile > /dev/null

real    0m1.018s
user    0m0.948s
sys 0m0.012s

$ time awk '{ split($0, a, "[^CTct]+"); m=0; for (i in a) { len=length(a[i]); if (len > m) m=len } print $0, m }' inputfile > /dev/null

real    0m1.802s
user    0m1.688s
sys 0m0.028s

$ time perl -MList::Util=max -lpe '$_ = "$_ " . max map { length $_ } /[CT]*/gi' inputfile > /dev/null

real    0m1.216s
user    0m1.160s
sys 0m0.016s

$ time sed 'h;y/cCtT/xxxx/;x;H;s/./x/g;G; s/^\(x*\).*\n.*\1.*\n/\1 /; s/^x\{10\}/1/;s/$/:9876543210xxxxxxxxx/; s/^\(1*\)\(x*\) \(.*\):.*\(.\).\{9\}\2$/\3 \1\4/' inputfile > /dev/null

real    1m4.165s
user    1m2.784s
sys 0m0.352s
  • 1
    It would be better to time a single execution on a large dataset rather than multiple executions on a small dataset. That would give a truer picture of which code is faster rather than how long it would take to start up the programs. – Kusalananda Jun 29 '18 at 13:02
  • 1
    @Kusalananda, thanks for the suggestion. I've made a test file with 120000 lines and rerun the "benchmark". My approach still remains the fastest one. – RomanPerekhrest Jun 29 '18 at 13:35
  • Yes, that look better. I've also modified my solution to split on [^CTct]+ (note the plus). That should make fewer bits to test against in my code. – Kusalananda Jun 29 '18 at 13:48
  • 1
    ok, that's better. See the extended "benchmark" results – RomanPerekhrest Jun 29 '18 at 13:56
  • I wonder how this compares in speed to a C solution? – qwr Jun 29 '18 at 15:50
0

Try also

awk '
        {T0 = $0
         while (match (T0, /[CTct]+/))  {if (RLENGTH > MX) MX = RLENGTH
                                         T0 = substr (T0, RSTART+RLENGTH)
                                        }
         print $0, MX
        }
' file

timing is a wee bit faster than the other awk proposals.

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.