mjhjmtu

I have this code that compares genes to a large list of snps on the same chromosome. For this I only want to compare genes and snps that are within +/- 1000000 bases of each other, but when I try to filter with awk its not working.

My file that I'm extracting from looks like this

CHR# SNP_ID POS samp_1 samp_2 ...
chr1 rs1212 174654646 0 2 ...
chr1 rs1331 321311111 1 1 ...
... ... ... ... ... ...

my filtering process looks like this

upper_bound=$(expr $gene_stop + 1000000)
lower_bound=$(expr $gene_start - 1000000)
zcat chr1.genotypes.txt.gz | tail -n +2 | awk 'if ($3 >= $lower_bound && $3 <= $upper_bound) print $0' > tmp_filtered

It is currently outputting empty files. When I change the awk conditional to only ($3 >= $lower_bound) nothing is printed, and when I change the conditional to be ($3 <= $upper) it prints but doesn't filter anything. I've tried checking that the lower and upper bound variables are reasonable. 1st, manually checking the positions of my snps I see that there are snps that lie in between the two thresholds. 2nd by printing out the length of the variable with $#foo prints out the correct length, so we can assume there are no hidden characters making it act as a string.

Can any one advise me?

TL;DR trying to grab items with a position between a given range, awk is not working as I expect

Shell variables are single quoted. In single quotes, variables are not expanded.

$ start=100
$ echo '$start'
$start

The same happens to awk:

$ start=100
$ echo awk '$3>=$start'
awk $3>=$start

The usual solution is to set the values with -v:

awk -vvar1=$lower -vvar2=$upper 'if ($3 >= var1 && $3 <= $var2) print $0'

So, your script should work with:

up_b=$(expr $gene_stop + 1000000)
lo_b=$(expr $gene_start - 1000000)
zcat chr1.genotypes.txt.gz | tail -n +2 | 
awk -vlo=$lo_b -vup=$up_b 'if ($3 >= lo && $3 <= up) print $0' > tmp_filtered