This is a test version of Biostars. For the public version, visit https://www.biostars.org.
Finding SNPs per gene for SNP Density in UNIX

I am attempting find out which gene harbors the most variation by using SNP density in a particular isolate. I am at a loss to figure out SNPs per gene, however. There's 445 SNPs present in the isolate with 9 genes. I used VarScan for that. I know their exact lengths but figuring out SNPs per gene is tripping me up. I can supposedly grep for patterns in UNIX (ssh shell) and then use an option to count the number of occurrences, but so far I'm not seeing what command to use/what pattern to grep for. Can anyone make sense of what grep pattern I'm looking for or lead me to an alternative way of figuring out SNPs per gene in unix...

grep snp varscan unix
• 3,688 views • link
updated by Ram Ram • written by crawlthrutacos • crawlthrutacos

2 answers

Via a combination of BEDOPS vcf2bed conversion and bedmap --count operation:

$ vcf2bed < snps.vcf | bedmap --echo --count --delim '\t' genes.bed - > genes_with_number_of_overlapping_snps.bed

You'll need a file called snps.vcf containing your SNPs, and a sorted file called genes.bed containing gene annotations of interest.

To calculate SNPs per gene, you sum all the counts and divide by the number of genes:

$ awk 'BEGIN { s = 0; } \
             { s += $NF; } \
         END { print (s / NR); }' genes_with_number_of_overlapping_snps.bed
• 1 views • link
written by Alex Reynolds Alex Reynolds
I guess you have a table/text file containing the information you need but not in the snp per gene format you are seekong. Use data aggregation. If you are not that into command line approaches, use pivot tables from microsoft excel. Learning pivot tables (or any other way to do data aggregation) will make your life easier
• 0 views • link
written by Irsan Irsan

Log in to answer this question.

More posts like this