dbSNP

[yk-tanigawa]

• https://www.ncbi.nlm.nih.gov/variation/docs/human_variation_vcf/
• ftp://ftp.ncbi.nih.gov/snp/organisms/human_9606/VCF/

$ wget -O - ftp://ftp.ncbi.nih.gov/snp/organisms/human_9606/VCF/00-All.vcf.gz | zcat - > /home/ytanigaw/PI_HOME/data/dbsnp/dbsnp_all_20160527.vcf
$ awk '(! /^#/){ print "chr"$0 }(/^#/){ print $0 }' dbsnp_all_20160527.vcf > bgzip -c - > dbsnp_all_20160527.vcf.gz

[yk-tanigawa]

• I have installed tabix under /share/PI/mrivas/bin/ and add symbolic links for bgzip and tabix

[yk-tanigawa]

After indexing a vcf file,

awk '(! /^#/){ print "chr"$0 }(/^#/){ print $0 }' dbsnp_all_20160527.vcf | bgzip > dbsnp_all_20160527.100.vcf.gz
tabix -p vcf dbsnp_all_20160527.vcf.gz 

I should be able to query SNPs by index like

$tabix dbsnp_all_20160527.vcf.gz chr1:10000-10020

[yk-tanigawa]

• Indexing step has done
• I deleted the original vcf file
  • vcf.gz file and vcf.gz.tbi files are sufficient

[yk-tanigawa]

column names are

CHROM   POS ID  REF ALT QUAL    FILTER  INFO

[yk-tanigawa]

0-index/1-index issue

[yk-tanigawa]

[ytanigaw@sh-2-1 ~/data/nanopore-wgs-consortium]$ head nanopore-wgs.25000.sorted.10k.mapq50.ext.w0.dbsnps |awk '{print $1, $2, $3, $4, $5}'
chr10 36353487 rs147789727 G A
chr10 50479569 rs116159204 A G
chr10 50502624 rs568082579 T A
chr10 81365552 rs11189519 T A
chr11 4953946 rs115916471 G A
chr11 83638545 rs113399123 C T
chr11 87196993 rs777698009 T C
chr11 87199541 rs772763241 C G
chr11 89919898 rs75442183 A C
chr11 92048173 rs568985594 AAAG A

• When you look at chr10:36353487 on iGV, rs7893390 is present.
• I guess this is due to the difference between 1-index/0-index coordinate systems

[yk-tanigawa]

• Turns out that dbSNP database is based on GRCh38, not GRCh37(hg19)
  • https://www.ncbi.nlm.nih.gov/SNP/snp_ref.cgi?type=rs&rs=rs147789727

[yk-tanigawa]

• downloaded dbSNP: /share/PI/mrivas/data/dbsnp/dbsnp_all_20151104.vcf.gz

[yk-tanigawa]

• Looks like it works now.

[ytanigaw@sherlock-ln02 login_node ~/data/nanopore-wgs-consortium]$ cat nanopore-wgs.25000.sorted.10k.mapq50.ext.w0.snps|grep -v '>'|head|awk '{print $1, $2, $3, $4, $5}'
chr10 394714 rs4880610 G A
chr10 1227299 rs7894015 G A
chr10 1235018 rs1392831 G A
chr10 3461323 rs10795040 G A
chr10 3485920 rs7067778 C A
chr10 6978146 rs1417020 A G
chr10 8647329 rs10905422 A T
chr10 8653428 rs11255772 C T
chr10 12689914 rs12780223 T C
chr10 14445612 rs10796188 G A

• 1248 hits on dbSNP with base call Q >= 20

  [ytanigaw@sherlock-ln02 login_node ~/data/nanopore-wgs-consortium]$ cat nanopore-wgs.25000.sorted.10k.mapq50.ext.w0.snps|grep -v '>'|wc -l
  1248

[yk-tanigawa]

• Implemented in https://github.com/rivas-lab/nanopore/blob/20161217GIABextraction/src/dump_snps.py to check dbsnp.
  • you need to use -s option to pass a gzipped vcf file