Reference that matches vcf?

[brofallon]

I recently downloaded the reference (hg19) chromosome 10 fasta file from UCSC, and naively assumed that this reference would match the reference alleles in a vcf from 1000 genomes. For instance, if the vcf listed the reference allele as being a C at position 4, then the fasta file should have a C as the fourth base. This doesn't seem to be true - the first few million bases are all shifted by 9997 (yes, 9997). After about 5.7M bases or so, there is a long string of N's, and after that once again the vcf doesn't match the fasta. I tried trimming off all the N's in the fasta, but that didn't fix the problem.
Is there a fasta sequence out there that matches the vcf's provided by the 1000 Genomes? Is there a way of reliably mapping from one to the other?
Thanks!

[Richard Finney]

Position 4? or line 4?
I downloaded this ftp://ftp-trace.ncbi.nih.gov/1000gen...els_svs.vcf.gz

And check out the first entries:
zcat ALL.chr10.phase1_release_v2.20101123.snps_indels_svs.vcf.gz | grep -v "#" | cut -b-60 | head
10 60523 rs148087467 T G 100 PASS RSQ=0.8535;AN=2184;VT=SNP;
10 60969 rs187110906 C A 100 PASS AN=2184;RSQ=0.3757;VT=SNP;
10 61005 rs192025213 A G 100 PASS AN=2184;THETA=0.0077;LDAF=
10 61020 rs115033199 G C 100 PASS RSQ=0.3512;AN=2184;AC=8;ER
10 61334 rs183305313 G A 100 PASS AN=2184;VT=SNP;AA=.;THETA=
10 65978 rs186558141 A C 100 PASS AVGPOST=0.9986;RSQ=0.4133;
10 66040 . CAG C 58 PASS ERATE=0.0057;RSQ=0.3430;AN=2184;LDA
10 66269 rs190079063 C T 100 PASS THETA=0.0115;AN=2184;AVGPO
10 66326 rs12260013 A G 100 PASS THETA=0.0050;AN=2184;LDAF=0
10 66627 . TAAAC T 378 PASS AC=953;AVGPOST=0.7237;AN=2184;VT

[brofallon]

Position 4 was an arbitrary example... Using the data you listed above, the vcf indicates that the reference should be '"T" at position 60523, but this doesn't seem to be the case with the UCSC genome. Has anyone else noticed this? Maybe I got a corrupted download or something..

[Richard Finney]

Looks like an off by one issue:

Note the $2-1 in the awk statement in this command:

zcat ALL.chr10.phase1_release_v2.20101123.snps_indels_svs.vcf.gz | grep -v "#" | cut -b-60 | head -30 | awk '{print "./fetchdna chr"$1" "$2-1 " "$2" supposedtobe="$4".reallyis="}' | bash | sed 's/\./ /'
supposedtobe=T reallyis=: T
supposedtobe=C reallyis=: c
supposedtobe=A reallyis=: a
supposedtobe=G reallyis=: g
supposedtobe=G reallyis=: g
supposedtobe=A reallyis=: a
supposedtobe=CAG reallyis=: c
supposedtobe=C reallyis=: c
supposedtobe=A reallyis=: a
supposedtobe=TAAAC reallyis=: t
supposedtobe=C reallyis=: c
supposedtobe=GA reallyis=: g
supposedtobe=T reallyis=: t
supposedtobe=A reallyis=: a
supposedtobe=T reallyis=: t
supposedtobe=G reallyis=: g
supposedtobe=ATAACT reallyis=: a
supposedtobe=G reallyis=: G
supposedtobe=G reallyis=: g
supposedtobe=C reallyis=: c
supposedtobe=G reallyis=: G
supposedtobe=G reallyis=: G
supposedtobe=G reallyis=: g
supposedtobe=C reallyis=: c
supposedtobe=A reallyis=: a
supposedtobe=G reallyis=: g
supposedtobe=G reallyis=: g
supposedtobe=G reallyis=: g
supposedtobe=T reallyis=: t
supposedtobe=CAGACACA reallyis=: c

[brofallon]

As I mentioned in the OP, it looks to be an off-by-9997 issue, but only for the first 5.64M bases. After that, a long string of N's appears in the reference fasta, and after that shifting by 9997 doesn't work any longer, and I can't figure out what the correct shift is.
In your experience, if a vcf says the reference is X at position Y, then the (one-indexed) Xth base in the reference should be Y? If so, then I somehow managed to obtain a bad reference chromosome.
Thanks for the help,

[Bukowski]

http://www.1000genomes.org/category/reference