awk getline problem

yifangt · December 13, 2012, 11:41am

Hello,
I want to print out the DNA sequence entries (tens of thousand!) that are longer than certain value (i=200) from a file (FASTA file) as:

>S94D_ctg_8004 Average coverage: 402.95
ATAATGCCTGTGAATATGACATGTGTTCCTGTTTCTACATCAGACTACTATTCTTGCATA
TCCCGCATCAGTCACTTATGATAAGGCACCTTTAAGTCGCCACTGTGATATTGGTTGTGT
CTTGCATATGATTTCTAGCATGTACTGCTTCTAATTTGTTGAAACGCCATACAGTTTGAA
CTTGGCAGAGTGATATTGGTTGCATCTTTCATATGGTGTCTAGCTTGCAGTGCTTCTAAT
TT
>S94D_ctg_8022 Average coverage: 212.74
CGTTGTCTGTGTGACATTTAAGAGTGTTTTGTGCAGTGCAACAAAAGGTCACAATGGCAA
CGAGGTGAAGTAGATAACCCTGTTTACACAGAGGTGCAAAGGAAACAACTAGTGGTCAAT
AACCCTGTTTACCAGCACCATGTGCATCTACCGATGTACTGGGGTGATGCACAGTCTGTT
G
>S94D_ctg_8062 Average coverage: 710.87
ATGAATTTCAGACGAGTTTTGGATTTACTAGAATTTAAAAACCAGGCATCTCAATGTTTT
GCCGGCAATCAACGGTGCCCTGGTGTTTGAAATTCATTCCCATTTCTTGCATGGGACCTA
AGCATGCACCCAAGGACA
>S94D_ctg_9034 Average coverage: 37.84
AGATTATTTTGTCTGCCATGTATAATTTTGGTTGATGTTTAGCCTGTTGTGCTTAACATG
CTTCTCGACGTACCTACACAGGACAATTTGGGAACGACTGCTGTTTTCCATCGAGGTTAG
TTTCATCCCATGGCTTATATCTGCTCAATGTTCAGGATATCGGTAGCCGGTACCATATAG
GCCGGCGGCTGATAGGAGACTAATCGGTGAATCGGACT

I tried:

awk '! /^>S/ { next } { getline seq } (length(seq) >= i) { print $0 "\n" seq }' i=200 infile.fasta

Nothing was print out.
I noticed the DNA sequence line are 60bp wide, if my i is less than 60, only the header (the lines with ">" as separator) and the first sequence line were printed, but I am expecting all sequence lines under each header. If the i is bigger than 60, nothing was printed. For sure most of the entries in the infile.fasta are longer than 60bp. Did I miss anything?
Thanks a lot!

Corona688 · December 13, 2012, 11:46am

getline gets lines, it doesn't paste them together.

I don't think you need getline here in any case, awk doesn't need getline's help to read lines one at a time from the default source.

It'd be helpful to see what this program should be doing, not just what it isn't.

vgersh99 · December 13, 2012, 12:47pm

I *think* this is what you want.
awk -v len=200 -f yi.awk myFile
yi.awk:

function doPrint(r,h,   n)
{
     n=gsub(RS,RS,r)
     if (length(r)-n>len)
       print h RS r
}
/^>S/ {
   if (h) {
     doPrint(r,h)
     r=""
   }
   h=$0
   next
}
{r=(r)?r RS $0:$0}
END {
  doPrint(r,h)
}

yifangt · December 13, 2012, 2:15pm

The interesting part is the script works with another fasta file in which each sequence line/row has more 120bp wide.

>S94D_ctg_8004 Average coverage: 402.95
ATAATGCCTGTGAATATGACATGTGTTCCTGTTTCTACATCAGACTACTATTCTTGCATATCCCGCATCAGTCACTTATGATAAGGCACCTTTAAGTCGCCACTGTGATATTGGTTGTGT
CTTGCATATGATTTCTAGCATGTACTGCTTCTAATTTGTTGAAACGCCATACAGTTTGAACTTGGCAGAGTGATATTGGTTGCATCTTTCATATGGTGTCTAGCTTGCAGTGCTTCTAAT
TT
>S94D_ctg_8022 Average coverage: 212.74
CGTTGTCTGTGTGACATTTAAGAGTGTTTTGTGCAGTGCAACAAAAGGTCACAATGGCAACGAGGTGAAGTAGATAACCCTGTTTACACAGAGGTGCAAAGGAAACAACTAGTGGTCAAT
AACCCTGTTTACCAGCACCATGTGCATCTACCGATGTACTGGGGTGATGCACAGTCTGTTG
>S94D_ctg_8062 Average coverage: 710.87
ATGAATTTCAGACGAGTTTTGGATTTACTAGAATTTAAAAACCAGGCATCTCAATGTTTTGCCGGCAATCAACGGTGCCCTGGTGTTTGAAATTCATTCCCATTTCTTGCATGGGACCTA
AGCATGCACCCAAGGACA
>S94D_ctg_9034 Average coverage: 37.84
AGATTATTTTGTCTGCCATGTATAATTTTGGTTGATGTTTAGCCTGTTGTGCTTAACATGCTTCTCGACGTACCTACACAGGACAATTTGGGAACGACTGCTGTTTTCCATCGAGGTTAG
TTTCATCCCATGGCTTATATCTGCTCAATGTTCAGGATATCGGTAGCCGGTACCATATAGGCCGGCGGCTGATAGGAGACTAATCGGTGAATCGGACT

I thought the script may be related to the newline, but I did check and could not find any abnormal. Not sure the EOL is really the same of the two files. I will check it. Thank you both!

Corona688 · December 13, 2012, 2:18pm

A program which doesn't do what you want really doesn't help explain what you do want, and you still haven't showed output demonstrating what you actually want. So we're still a bit in the dark, here.

yifangt · December 13, 2012, 2:34pm

Hi, Corona:
For example, infile:

>S94D_ctg_8004 Average coverage: 402.95
ATAATGCCTGTGAATATGACATGTGTTCCTGTTTCTACATCAGACTACTATTCTTGCATA
TCCCGCATCAGTCACTTATGATAAGGCACCTTTAAGTCGCCACTGTGATATTGGTTGTGT
CTTGCATATGATTTCTAGCATGTACTGCTTCTAATTTGTTGAAACGCCATACAGTTTGAA
CTTGGCAGAGTGATATTGGTTGCATCTTTCATATGGTGTCTAGCTTGCAGTGCTTCTAAT
TT
>S94D_ctg_8022 Average coverage: 212.74
CGTTGTCTGTGTGACATTTAAGAGTGTTTTGTGCAGTGCAACAAAAGGTCACAATGGCAA
CGAGGTGAAGTAGATAACCCTGTTTACACAGAGGTGCAAAGGAAACAACTAGTGGTCAAT
AACCCTGTTTACCAGCACCATGTGCATCTACCGATGTACTGGGGTGATGCACAGTCTGTT
G
>S94D_ctg_8062 Average coverage: 710.87
ATGAATTTCAGACGAGTTTTGGATTTACTAGAATTTAAAAACCAGGCATCTCAATGTTTT
GCCGGCAATCAACGGTGCCCTGGTGTTTGAAATTCATTCCCATTTCTTGCATGGGACCTA
AGCATGCACCCAAGGACA
>S94D_ctg_9034 Average coverage: 37.84
AGATTATTTTGTCTGCCATGTATAATTTTGGTTGATGTTTAGCCTGTTGTGCTTAACATG
CTTCTCGACGTACCTACACAGGACAATTTGGGAACGACTGCTGTTTTCCATCGAGGTTAG
TTTCATCCCATGGCTTATATCTGCTCAATGTTCAGGATATCGGTAGCCGGTACCATATAG
GCCGGCGGCTGATAGGAGACTAATCGGTGAATCGGACT

I want all the entries with sequence less than 200bp are filtered to get outfile as:

>S94D_ctg_8004 Average coverage: 402.95
ATAATGCCTGTGAATATGACATGTGTTCCTGTTTCTACATCAGACTACTATTCTTGCATA
TCCCGCATCAGTCACTTATGATAAGGCACCTTTAAGTCGCCACTGTGATATTGGTTGTGT
CTTGCATATGATTTCTAGCATGTACTGCTTCTAATTTGTTGAAACGCCATACAGTTTGAA
CTTGGCAGAGTGATATTGGTTGCATCTTTCATATGGTGTCTAGCTTGCAGTGCTTCTAAT
TT
>S94D_ctg_9034 Average coverage: 37.84
AGATTATTTTGTCTGCCATGTATAATTTTGGTTGATGTTTAGCCTGTTGTGCTTAACATG
CTTCTCGACGTACCTACACAGGACAATTTGGGAACGACTGCTGTTTTCCATCGAGGTTAG
TTTCATCCCATGGCTTATATCTGCTCAATGTTCAGGATATCGGTAGCCGGTACCATATAG
GCCGGCGGCTGATAGGAGACTAATCGGTGAATCGGACT

There are some bioperl script for this job, but I like awk that can do it in a flash. Thanks!

Corona688 · December 13, 2012, 2:43pm

Do you have GNU awk? You can change the meaning of 'record' with RS, split on > instead of \n, and count the length more directly:

awk -v RS="\n>" -F"\n" '{ V=$0; gsub(/\n/, "", V); if((length(V)-length($1)) > 200) print ">"$0 }'

vgersh99 · December 13, 2012, 2:44pm

doesn't my awk script produces the correct output based on your sample above?
Is there still an issue?

rdrtx1 · December 13, 2012, 3:10pm

another awk suggestion:

awk -v i=200 '
/>/ {if (length(s)>i)printf h s; c=0; s=""; h=$0"\n" ; next;}
{ c+=length ; s=s $0 "\n";}
END {if (length(s)>i)printf h s}
' infile.fasta

yifangt · December 13, 2012, 3:16pm

Thanks! All worked like charm!
Can I ask what was wrong with my own script? I thought it should work. Actually it did sometime, yet sometime did not!
By the way, I believe the scripts will be used by a lot of "bioinformaticians" for this routine job with DNA sequences!
Thank you all again!

Corona688 · December 13, 2012, 3:24pm

I couldn't tell what your original script was even trying to do.

getline profoundly doesn't do what what you thought. It reads lines, but doesn't squeeze them together.

Your control logic was also kind of odd. That 'next' guarantees no lines containing > would ever be printed unless you read them manually with getline. And when you didn't see a >, 'getline' ran every single line -- so you were reading double-speed, then throwing half your lines away. Doing a manual getline bypasses your /^>S/ completely, which is the only way any of those lines could have been printed at all.

I can only guess that some of your input had very long lines, and that work only meant "work". Or perhaps you weren't reading the output you thought you were.

yifangt · December 13, 2012, 3:36pm

Thanks Corona!
That's the part I was afraid of, which can lead to a totally wrong result! I am trying to catch more details about awk and hope use it correctly.
Thank you again!

Corona688 · December 13, 2012, 3:38pm

Understand that awk has a built-in while loop. This code:

awk '/asdf/'

amounts to:

while( (getline) > 0)
{
        if($0 ~ /asdf/) print
}

Putting a getline after that 'if' would mean it wouldn't get checked by the if-statement.