0

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

  • 1
    the shell variables are not accessible inside awk. Try something like zcat your_file | awk -vstart=$start -v stop=$stop '$3 >= start - 1000000 && $3 <= stop + 1000000'. Notice that you don't need an explicit if and print to awk; simply use your condition as the "pattern". – mosvy Dec 1 '18 at 0:08
  • Please post always input file and expected output – Praveen Kumar BS Dec 1 '18 at 6:29
0

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

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.