Imbriqué 'awk' dans une boucle 'while', parsingr deux files ligne par ligne et comparer les valeurs des colonnes

J'ai besoin d'aide avec une combinaison de awk & while loop. J'ai deux files simples avec des colonnes (les normales sont très grandes), l'une représentant des intervalles simples pour un ID = 10 (des régions codantes (exons), pour le chromosome 10 ici):

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

et l'autre représentant les lectures séquencées (= juste des intervalles à nouveau mais plus petit) avec une autre valeur comme dernière colonne, dont j'ai besoin plus tard:

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

Je voudrais donc effectuer une search rapide et efficace et find quels intervalles de lecture (dont ligne dans le file) et combien, tombent dans une région de encoding:

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

afin que je puisse get la valeur de la quasortingème colonne de ces lignes plus tard, pour chaque exon.

J'ai écrit un code, qui a probablement besoin de corrections sur la boucle while, puisque je ne peux pas faire parsingr les lignes une par une pour chaque awk. C'est ici:

 while read chr ab 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 

Pour le moment je peux faire fonctionner la ligne awk quand je donne manuellement a, b, mais je veux la faire fonctionner automatiquement pour chaque paire a, b par file.

Toute suggestion sur la modification de la syntaxe, ou la façon de le faire, est très appréciée!

SUIVRE

Finalement, j'ai travaillé avec ce 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 

SORTIE:

  Exon 2 : Reads with their COV: 2|40 3|22 4|35 5|41 Sum= 138 Mean= 34.5 etc. 

Le premier problème est que vous ne pouvez pas utiliser les variables bash à l'intérieur de awk comme ça. $a dans awk évalue le champ a mais a est vide puisqu'il n'est pas défini dans awk , mais dans bash . Une façon de contourner ceci est d'utiliser l'option -v awk pour définir la 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. 

Donc, vous pourriez faire:

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

Vous avez une autre erreur là-bas cependant. Pour qu'une lecture tombe dans un exon, la position de départ de la lecture doit être supérieure à la position de départ de l'exon et sa position de fin plus petite que la position finale de l'exon. Vous utilisez $2<=a && b <= $3 qui sélectionnera les lectures dont le début est en dehors des limites de l'exon. Ce que vous voulez est $2>=a && $3<=b .

En tout cas, exécuter ce type de chose dans une boucle bash est très inefficace puisqu'il faut lire le file d'input une fois pour chaque paire de a et b . Pourquoi ne pas tout faire en 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 

Le script ci-dessus produit la sortie suivante s'il est exécuté sur vos files d'exemple:

 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 

Voici la même chose sous une forme less condensée pour plus de clarté

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

Je sais que ce n'est pas tout à fait ce que vous cherchez, mais personnellement – je ne m'entends pas avec awk et je suggère donc d'avoir une perl.

Quelque chose comme ça:

 #!/usr/bin/perl #REALLY GOOD IDEA at the start of any perl code use ssortingct; 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"; } 

Ce que donne votre échantillon de données donne:

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

Vous devriez pouvoir étendre ceci pour faire un contrôle / validation de gamme plus compliqué aussi!