How do you merge SNP data with a reference genome?

Question

My Data

I have a 23andMe file listing SNPs in the form:

rsid chromosome position genotype rsXXXXX 1 PPPPPP CT rsXXXXX 1 PPPPPP GG

Fields are TAB-separated and each line corresponds to a single SNP. For each SNP, four fields of data are supplied.

1. An identifier (an rsid or an internal id)
2. Its location on the reference genome.
   • The chromosome it is located on.
   • The position within the chromosome is is located on.
3. The genotype call oriented with respect to the plus strand on the human reference sequence.

The reference genome is the human assembly build 37 (also known as Annotation Release 104).

My Question

How do I merge the SNPs into the reference genome?

For example, take the first line in my SNP file:

rsXXXXX 1 PPPPPP CT

Part 1

I can see that I need to replace the nucleotide at position PPPPPP on chomosome 1 of the reference genome with a nucleotide from the genotype field, but which nucleotide am I supposed to use? C or T? And why?

Part 2

Where am I supposed to start counting from on the reference genome? Looking at chromosome 1 of the human assembly build 37, the first ~10,000 characters (excluding the first line description) are N. Is the first N number 1? eg. If PPPPPP was 100,000 would I replace the 100,000th character in the reference genome with the correct nucleotide from Part 1 of this question? Or should I start counting from the first non-N character in the fasta file?

Answer 1

First, you need to know which genome sequence does the SNP file refer to. They must have mentioned the reference sequence that they used.

As others mentioned the case of CT is heterozygosity. If you just want to mark the changes then discard the residue that is already present in the reference genome and use the other allele. However, you want to keep a track of the haplotype then you have make sure that a set of SNPs come from the same chromatid. This is difficult- you might still be able to know it for SNPs that are close enough to be mapped by a single read but it is almost impossible for SNPs that are separated well enough.

As Endre said, you have to start from the first nucleotide. However, it seems dubious that you are getting $\ (NNNN)_n$ in the beginning of chromosome 1. Complete assembled chromosomes do not have such stretches. Below are the first 10 lines of the chromosome 1 fasta file. See for yourself.

>gi|568815364|ref|NT_077402.3| Homo sapiens chromosome 1 genomic scaffold, GRCh38 Primary Assembly HSCHR1_CTG1
TAACCCTAACCCTAACCCTAACCCTAACCCTAACCCTAACCCTAACCCTAACCCTAACCCTAACCCTAAC
CCTAACCCTAACCCTAACCCTAACCCTAACCCTAACCCAACCCTAACCCTAACCCTAACCCTAACCCTAA
CCCTAACCCCTAACCCTAACCCTAACCCTAACCCTAACCTAACCCTAACCCTAACCCTAACCCTAACCCT
AACCCTAACCCTAACCCTAACCCCTAACCCTAACCCTAAACCCTAAACCCTAACCCTAACCCTAACCCTA
ACCCTAACCCCAACCCCAACCCCAACCCCAACCCCAACCCCAACCCTAACCCCTAACCCTAACCCTAACC
CTACCCTAACCCTAACCCTAACCCTAACCCTAACCCTAACCCCTAACCCCTAACCCTAACCCTAACCCTA
ACCCTAACCCTAACCCTAACCCCTAACCCTAACCCTAACCCTAACCCTCGCGGTACCCTCAGCCGGCCCG
CCCGCCCGGGTCTGACCTGAGGAGAACTGTGCTCCGCCTTCAGAGTACCACCGAAATCTGTGCAGAGGAC
AACGCAGCTCCGCCCTCGCGGTGCTCTCCGGGTCTGTGCTGAGGAGAACGCAACTCCGCCGTTGCAAAGG

How to replace $N^{th}$ residue is quite a straightforward task. But that is a programming question and not the scope of this forum. Assuming that you have solved the part 1 problem and have a tab separated sorted file like this:

chromosome  position    residue
 1           79989           G
 1           100232          T
 3           341342          A

This script may not be the best but would work in a linux/*nix/Cygwin terminal, to replace the residues (ensure that you have gawk version >=4.0):

gawk -F "\t" '(FNR==1){x++} (x==1){a[$1][$2]=$3;next} (x==2){if($0~/>/){h=$0;sub(/^.*chromosome /,"",h);sub(/ .*/,"",h)} else{seq[h]=seq[h]$0}} END{for(i in a){s=0; for(j in a[i]){m=m substr(seq[i],s,j-1) a[i][j];s=j+1} m=m substr(seq[i],s); print ">Chr"i"\n"m}}' SNP_file Genome.fa | fold -w 60

Answer 2

Genetics 101, you have 2 copies of all your DNA at every position, one copy from your mother, one from your father. So for the "CT" one, you have one copy with a C, and one with a T.

And yea, it's normal for the first several thousand, or million letters to be N's. The genome is repetitive and yucky there, but it's counted for numbering purposes anyway.

Honestly, I wouldn't do this with a giant text file of the genome. Just search for your SNP in ensembl.org using the rs number, and you'll get the SNP, and some flanking sequence, and some context. Search it in PubMed if you want to see if it's ever turned up in any publication

Answer 3

Part 1:

According to Lior Pachter, 23andme data are not phased. Which means that for each entry in the genotype field, you do not know which chromosome copy it came from. This happens since modern microarray platforms are not able to tell which of the two copies of a chromosome a snp came from.

You can solve this problem for most snps by comparing your alleles to the reference genome, but this would take some programming chops. You could use https://github.com/endrebak/qc_gwas as an example, which does the same thing, but for plink files.

Part 2:

I'm assuming that you would like to do this programatically, and not by copying and pasting the snps into the reference genome.

The short answer is that the first N is the first nucleotide. But, you should rather use a package such as Biopython to do the counting for you, it might be knottier than you think (you need to adjust for line endings in the fasta file, for example).