Plotting SNP density along a chromosome from VCF files

Question

[rc16955]

Hi all,

I have generated VCF files from a range of samples of my study organism using a published reference genome. I now want to use these to assess where diversity is along a given chromosome; in other words I want to know the density of SNPs along chromosome 2.

One way I had thought of doing this was to split the chromosome into 10Kb chunks, count the number of SNPs in each one, and plot this in a chart. So what I want is basically:

Chr2:1-10000, SNPs = X Chr2:10000-20000, SNPs = Y etc.

But I can't figure out an efficient way of actually extracting this information from my VCF files.I am able to use the command

tabix -h myfile.vcf.gz chr2:[10Kregion] > myoutput.vcf

to extract a given region of any VCF file and then use bcftools stats to count the number of SNPs therein, but the chromosome is 16Mb long; there must be a more efficient way to do this than for me to manually extract 1,600 VCF files and examine them with bcftools stats.

I wonder if there's some software that could help with this. I have Integrated Genome Viewer, MEGA and SeaView, but if any of them have this kind of functionality, I haven't found it yet.

Otherwise, any help or advice would be much appreciated. It must be doable because I see plenty of published papers with SNP density plotted against chromosome position, but I can't figure out how to do it myself.

Thanks in advance!

Answer 1

[sacha]

I did something similar with human genom here : https://github.com/dridk/snp_location
Briefly I create my bins ( 10000 pb) using bedtools makewindows :

bedtools makewindows -g hg19.sizes -w 10000 > windows.bed

Then I count SNP per bin using bedtools coverage

bedtools coverage -a windows.bed -b yourdata.vcf -counts > coverage.txt

Finally I plot my data using R with IdeoViz package. see https://github.com/dridk/snp_location/blob/master/plot_chrom.r .
I get the following plot showing SNP density per bin :

[spiral01]

I have tried using bedtools coverage and get errors at several lines of the file:

WARNING: line number 10210 of file chr21.vcf.gz has an imprecise variant, but no SVLEN specified. Defaulting to length 1.

The resulting output gives a count of 0 at every window. Can anyone help with this?

[nagarsaggi]

Hi Sacha, I also want to do something similar but I for each sample from a multisample vcf file. I tried extracting a single sample from the multi-vcf file and then followed your command but the number of variants counts for different samples remain the same (probably these are the sites where a variant has been being called across the samples). Could you guide me to extract unique variants for each individuate and then plot the variants in a bin size.

[Romeo.Dervishaj]

Greetings Sacha, just out of curiosity, I want to do something similar to your solution but for a non-model organism (in my field usually is used Circos but visually I find more pleasant your solution) is it possible to do so? And eventually how can I build the ideo file my self? Thank you in advance

Answer 2

[Pierre Lindenbaum]

using sqlite

 curl -s "ftp://ftp.1000genomes.ebi.ac.uk/vol1/ftp/phase3/integrated_sv_map/ALL.wgs.integrated_sv_map_v2.20130502.svs.genotypes.vcf.gz" |\
gunzip -c | grep -v "^#" |\
awk -F '\t' 'BEGIN{printf("create table T(c text,p int); BEGIN TRANSACTION;\n");} {printf("insert into T(c,p) values(\"%s\",%d);\n",$1,$2);} END {printf("COMMIT; SELECT c,CAST(p/1E7  AS INT)*1E7, count(*) from T group by c,CAST(p/1E7  AS INT)*1E7 ; drop table T;");}' |\
sqlite3 -separator ' ' tmp.sqlite  && rm tmp.sqlite

1 0.0 245
1 10000000.0 208
1 20000000.0 191
1 30000000.0 197
1 40000000.0 211
1 50000000.0 198
1 60000000.0 159
1 70000000.0 240
1 80000000.0 194
1 90000000.0 165
1 100000000.0 220
1 110000000.0 206
1 120000000.0 15
1 140000000.0 101
1 150000000.0 202
1 160000000.0 179
1 170000000.0 181
1 180000000.0 231
1 190000000.0 260
(...)

[m.m]

Dear Pierre Lindenbaum! I hope this message finds you well. How would I have to modify this code to have it print out the result into a file?

[Pierre Lindenbaum]

https://www.guru99.com/linux-redirection.html

[m.m]

Thank you for your response. I managed using:

grep -v "^#" chr25newids.vcf | awk -F '\t' 'BEGIN{printf("create table T(c text,p int); BEGIN TRANSACTION;\n");} {printf("insert into T(c,p) values(\"%s\",%d);\n",$1,$2);} END {printf("COMMIT; SELECT c,CAST(p/1E4 AS INT)1E4, count() from T group by c,CAST(p/1E4 AS INT)*1E4;");}' | sqlite3 -separator ' ' tmp.sqlite

sqlite3 -separator ' ' tmp.sqlite "SELECT c,CAST(p/1E4 AS INT)1E4, count() from T group by c,CAST(p/1E4 AS INT)*1E4;" > somefile.txt