Aninhado 'awk' em um loop 'while', analisa dois arquivos linha por linha e compara valores de colunas

Question

Preciso de ajuda com uma combinação de awk & while loop. Eu tenho dois arquivos simples com colunas (os normais são muito grandes), um representando intervalos simples para um ID = 10 (de regiões de codificação (exons), para o cromossomo 10 aqui):

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

e o outro representando leituras sequenciadas (= apenas intervalos novamente, mas menores) com um outro valor como última coluna, que eu precisarei mais tarde:

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

Então, eu gostaria de pesquisar de forma rápida e eficiente e encontrar quais intervalos de leitura (de qual linha no arquivo) e quantos, cair em uma região de codificação:

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

para que eu possa obter o valor da quarta coluna dessas linhas mais tarde, para cada exon.

Eu escrevi um código, que provavelmente precisa de algumas correções no loop while, já que não posso fazer com que ele analise as linhas de leitura uma por uma para cada awk. Aqui está:

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

No momento eu posso fazer a linha awk rodando quando eu der manualmente a, b, mas eu quero fazer com que ela rode automaticamente para cada par a, b por arquivo.

Qualquer sugestão de mudar a sintaxe ou a maneira de fazê-lo é muito apreciada!

FOLLOW UP

Por fim, desenvolvi esse código:

    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.

    

Answer 1

A primeira questão é que você não pode usar variáveis bash dentro de awk dessa forma. $a dentro de awk é avaliado como campo a , mas a está vazio, pois não está definido em awk , mas em bash . Uma maneira de contornar isso é usar a opção awk de -v para definir a variável

-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.

Então, você poderia fazer:

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

Você tem outro erro aí. Para que uma leitura caia dentro de um exon, a posição inicial da leitura deve ser maior que a posição inicial do exon e sua posição final menor que a posição final do exon. Você está usando $2<=a && b <= $3 , que selecionará leituras cujo início está fora dos limites do exon. O que você quer é $2>=a && $3<=b .

Em qualquer caso, executar esse tipo de coisa em um loop bash é muito ineficiente, já que ele precisa ler o arquivo de entrada uma vez para cada par de a e b . Por que não fazer a coisa toda em 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

O script acima produz a seguinte saída, se executada nos seus arquivos de exemplo:

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

Aqui está a mesma coisa em uma forma menos condensada para maior clareza

#!/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" 
        }

    

Answer 2

Eu sei que não é bastante o que você procura, mas pessoalmente - eu não me dou bem com awk e, portanto, sugiro ter uma chance em perl.

Algo parecido com isto:

#!/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";
}

Qual dado dados da sua amostra fornecem:

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

Você deve ser capaz de estender isso para fazer uma verificação / validação de faixa mais complicada também!