Trying to get binned GC content in hg19

[mlodato]

Hi,

I have a large (~2 million) set of SNPs and I would like to know the GC content +/- 500 bp from the variants. I have some Python and Unix skills, but am not an expert. I have download the hg19 gc5Base from the UCSC table browser, which has the GC content in 5 bp bins across the genome:

$ more hg19.gc5Base.txt
variableStep chrom=chr1 span=5
10001 40
10006 40
10011 40
10016 60
10021 60
10026 60


But I am not sure of the best way to get specific regions from this file based on the SNPs I am interested in. Any help would be appreciated!

[vivek_]

One simple way I can think of with minimal coding work is to create a bed file with regions +/- 500 bp for each SNP. You could likely do it with a simple unix command. If the SNPs are in a VCF, something like this could work:

awk '! /\#/' SNPs.vcf | awk '{print $1"\t"($2-500)"\t"($2+500)}' > SNP_Regions.bed

You can then use bedtools nuc module to get the GC content for each region.

bedtools nuc -fi hg19.fa -bed SNP_Regions.bed

as detailed in this post:

https://www.biostars.org/p/47047/

[mlodato]

This is a great suggestion. I will give it a try today!