Cfrtjryk

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