4

I need some help with a combination of awk & while loop. I have two simple files with columns (normal ones are very large), one representing simple intervals for an ID=10(of coding regions(exons),for chromosome 10 here):

#exons.bed
10  60005   60100   
10  61007   61130   
10  61200   61300   
10  61500   61650   
10  61680   61850   

and the other representing sequenced reads(=just intervals again but smaller) with an other value as last column, that I ll need later:

#reads.bed
10  60005   60010    34 
10  61010   61020    40
10  61030   61040    22
10  61065   61070    35 
10  61100   61105    41

So, I would like to search in a quick and efficient way and find which read intervals (of which line in the file) and how many, fall in one coding region:

exon 1(first interval of table 1) contains reads of line 1,2,3, etc. 
of   reads.file(2nd table)

so that I can get the value of 4th column of these lines later, for each exon.

I 've written a code,that probably needs some corrections on the while loop, since I cannot make it parse the reads lines one by one for each awk. Here it is:

while read chr a b cov; do  #for the 4-column file

#if <a..b> interval of read falls inside exon interval:
awk '($2<=$a && $b <= $3) {print NR}' exons.bed >> out_lines.bed

done < reads.bed

At the moment I can make the awk line running when I give manually a,b, but I want to make it run automatically for each a,b pair by file.

Any suggestion on changing syntax, or way of doing it, is highly appreciated!

FOLLOW UP

Finally I worked it out with this code:

    awk 'NR==FNR{
        a[NR]=$2; 
        b[NR]=$3;
        next; }
    {  #second file
    s[i]=0; m[i]=0;  k[i]=0;              # Add sum and mean calculation
    for (i in a){                                            
       if($2>=a[i] && $3<=b[i]){         # 2,3: cols of second file here
          k[i]+=1
          print k                      #Count nb of reads found in
          out[i]=out[i]" "FNR          # keep Nb of Line of read 
          rc[i]=rc[i]" "FNR"|"$4       #keep Line and cov value of $4th col
          s[i]= s[i]+$4                #sum over coverages for each exon
          m[i]= s[i]/k[i]             #Calculate mean (k will be the No or  
                                       #reads found on i-th exon)
     }}  
    }
    END{
       for (i in out){
          print "Exon", i,": Reads with their COV:",rc[i],\
          "Sum=",s[i],"Mean=",m[i] >> "MeanCalc.txt"

    }}' exons.bed  reads.bed

OUTPUT:

   Exon 2 : Reads with their COV:  2|40 3|22 4|35 5|41 Sum= 138  Mean= 34.5
   etc.
  • Please edit your question and show us the exact output you would like from your example input files. – terdon Mar 20 '15 at 11:41
  • 1
    (1) Some (!) of us here on Stack Exchange are not microbiologists; you might get better help if you stated more clearly what arithmetic/logical functionality you want, and relied less on your jargon. For instance, I've never heard of an "exon" and I don't know what you mean by "coding region". (2) I believe that your example is inconsistent. I believe that, for the files you show, exon 1 contains read 1, and exon 2 contains reads 2, 3, and 4. If I'm wrong, please explain. – Scott Mar 20 '15 at 12:00
  • Also, that last command — the one with the uniq in it — isn't a valid command; it seems to be missing a program name at the beginning, What is that supposed to be? – Scott Mar 20 '15 at 12:06
3

The first issue is that you can't use bash variables inside awk like that. $a within awk evaluates to field a but a is empty since it is not defined in awk, but in bash. One way around this is to use awk's -v option to define the variable

-v var=val
--assign var=val
   Assign the value val to the variable var,  before  execution  of
   the  program  begins.  Such variable values are available to the
   BEGIN rule of an AWK program.

So, you could do:

while read chr a b cov; do 
  awk -v a="$a" -v b="$b" '($2<=a && b <= $3) {print NR}' exons.bed > out$a$b 
done < reads.bed

You have another mistake there though. In order for a read to fall within an exon, the read's start position must be greater than the start position of the exon and its end position smaller than the end position of the exon. You are using $2<=a && b <= $3 which will select reads whose start is outside the exon's boundaries. What you want is $2>=a && $3<=b.

In any case, running this type of thing in a bash loop is very inefficient since it needs to read the input file once for every pair of a and b. Why not do the whole thing in awk?

awk 'NR==FNR{a[NR]=$2;b[NR]=$3; next} {
        for (i in a){
           if($2>=a[i] && $3<=b[i]){
            out[i]=out[i]" "FNR 
        }}}
        END{for (i in out){
                   print "Exon",i,"contains reads of line(s)"out[i],\
                   "of reads file" 
        }}' exons.bed reads.bed

The script above produces the following output if run on your example files:

Exon 1 contains reads of line(s) 1 of reads file
Exon 2 contains reads of line(s) 2 3 4 5 of reads file

Here's the same thing in a less condensed form for clarity

#!/usr/bin/awk -f

## While we're reading the 1st file, exons.bed
NR==FNR{
    ## Save the start position in array a and the end 
    ## in array b. The keys of the arrays are the line numbers.
    a[NR]=$2;
    b[NR]=$3; 
    ## Move to the next line, without continuing
    ## the script.
    next;
}
 ## Once we move on to the 2nd file, reads.bed
 {
     ## For each set of start and end positions
     for (i in a){
         ## If the current line's 2nd field is greater than
         ## this start position and smaller than this end position,
         ## add this line number (FNR is the current file's line number)
         ## to the list of reads for the current value of i. 
         if($2>=a[i] && $3<=b[i]){
             out[i]=out[i]" "FNR 
         }
     }
 }
 ## After both files have been processed
 END{
     ## For each exon in the out array
     for (i in out){
         ## Print the exon name and the redas it contains
         print "Exon",i,"contains reads of line(s)"out[i],
             "of reads file" 
        }
  • Thanks terdon for your suggestion ! I used it and adapted it a while on what I need further,so you can check next,your code with some changes-for using 4th col info that I need in the end, and a further question,for what I didn't manage to do. I would like actually to calculate for each exon with its reads, the mean value of this 'c' value that corresponds to the reads of one exon, on $4th column or read.bed. – MariaK Mar 22 '15 at 22:56
  • @MariaK please post follow-up questions separately. Comments are ephemeral and don't serve for this type of thing. – terdon Mar 23 '15 at 12:56
2

I know it's not quite what you're after, but personally - I don't get along with awk and so would suggest having a go in perl.

Something like this:

#!/usr/bin/perl

#REALLY GOOD IDEA at the start of any perl code
use strict;
use warnings;

#open some files for input
open( my $exons, "<", 'exons.bed' ) or die $!;

#record where our exons start and finish. 
my %start_of;
my %end_of;

#read line by line our exons file. 
#extract the 3 fields and save 'start' and 'end' in a hash table. 
while (<$exons>) {
    my ( $something, $start, $end ) = split;

    my $exon_id = $.;    #line number;
    $start_of{$exon_id} = $start;
    $end_of{$exon_id}   = $end;
}
close ( $exons );

my %exons;
#run through 'reads' line by line, extracting the files. 

open( my $reads, "<", 'reads.bed' ) or die $!;
while (<$reads>) {
    my ( $thing, $read_start, $read_end, $value ) = split;

    #cycle through each exon. 
    foreach my $exon_id ( keys %start_of ) {

        #check if _this_ 'read' is within the start and end ranges. 
        if (    $read_start >= $start_of{$exon_id}
            and $read_end <= $end_of{$exon_id} )
        {
            #store the line number in our hash %exons. 
            push( @{ $exons{$exon_id} }, $. );
        }
    }
}
close ( $reads ); 

#cycle through %exons - in 'id' order. 
foreach my $exon_id ( sort keys %exons ) {
    #print any matches. 
    print "exon ",$exon_id, " (", $start_of{$exon_id}, " - ", $end_of{$exon_id},
        ") contains reads of line:", join( ",", @{ $exons{$exon_id} } ), "\n";
}

Which given your sample data gives:

exon 1 (60005 - 60100) contains reads of line:1
exon 2 (61007 - 61130) contains reads of line:2,3,4,5

You should be able to extend this to do some more complicated range checking/validation trivially too!

  • Just a nitpick: neither use strict nor use warnings are in any way mandatory. Including them is a good idea and promotes writing clean code but is not necessary. – terdon Mar 20 '15 at 15:31
  • Duly noted. I've amended it. The code will work without, but anyone you ask about a problem will probably point that out as the very first thing that's "wrong" – Sobrique Mar 20 '15 at 15:56
  • No, there's nothing "wrong". There are perfectly legitimate reasons to remove the use warnings and strict is just to promote clean code. Nothing wrong with not using them, it's just good practice. A Really Good Idea® as you said. They might point out that you haven't closed your open file descriptors though :) – terdon Mar 20 '15 at 15:58
  • They don't warn over that. 's a fair point, but one of the advantages of lexically scoped filehandles in this particular style is they close automatically when they leave scope. But in the interests of good example, I'll add that too :) – Sobrique Mar 20 '15 at 16:12
  • 1
    Oh yes, absolutely. It was just the "mandatory" that I was taking issue at. And I'd upvoted your answer even before you edited. By the way, if you scroll down to the very end of Perl's documentation for open() , you'll find this: "Always close the handle yourself and inspect the return value". I think they're thinking about fringe cases and I admit to leaving mine open in my simple scripts but. . . All I know is that Perl gurus shout at me when I don't close mine. – terdon Mar 20 '15 at 16:35

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.