Incorrect sequence retrieved by coordinates

Question

[Igor Filippov]

I want to fetch some sequences from reference genome by coordinates. I've downloaded Homo_sapiens.GRCh38.dna.primary_assembly.fa from the Ensembl's FTP server. I've tried both pysam and Biopython to get the sequence from the reference genome like so:

import pysam
ref_genome = pysam.FastaFile("data/Homo_sapiens.GRCh38.dna.primary_assembly.fa")
ref_genome.fetch("1", 37566812, 37595985)

from Bio import SeqIO
chrs = {}

for seq in SeqIO.parse("data/Homo_sapiens.GRCh38.dna.primary_assembly.fa", "fasta"):
    chrs[seq.id] = seq.seq
str(chrs['1'][37566812:37595985])

Both fetch the same sequence. However, when I look up the sequence at NCBI https://www.ncbi.nlm.nih.gov/nuccore/NC_000001.11?report=fasta&from=37566812&to=37595985&strand=true I get different data, while I assumed they're supposed to match. Am I doing something wrong?

[GenoMax]

Your link above (with strand=true) brings back sequence that is this (abbreviated):

>NC_000001.11:c37595985-37566812 Homo sapiens chromosome 1, GRCh38.p13 Primary Assembly
AGTCGCTCTGAGGCCGAGAGGGACGCGAGCGGAAGTGACGTACGTCGTGCACACGTGGTCCGGCGTGGTT
CAGGCGGGTGTCTTCGGCCGGGCTTGGGAACATAAAAGTTTGTTTCACCACGTAAGCCGGACCTCGCACT
CCGGTCCCGGTCTCGTCGCCAAGATGGTGAAGCCCAAGTACAAAGGACGGAGCACCATCAACCCGTCCAA
GGCCAGCACAAACCCAGGTACAGGCGGCAGGGCTGATCGAGGGTGGAGCTTGGAGGAAGTATAGGGAAGG
GCTCCGAGGAGCACTGGAGGGGAAGGAGGGTTACTCGAGTGGCATTGCTTAGAATGACTACGAGCGGGGT
GGCTGTTGGGGATGAAAAGGGAAAAAGGAAAGCAATATTTGTGGTGGTTTCTTCTCAGGGGGGAACTGAG
GAGGGCTCAGTGCACGCCTGGCAGGCCAAAGTGTAGTCAACTGGTCTGAAGTCAAGGAAAGTGAACACTC
AGGCTGACAGCCTTGGGTCCTCTCCAGTTGAGAGATGCAGAGAAAGTCAGATTAAAAATGGTGTAGGTTG
ACCCTTCGGGCCTTTTACACGTGATCATACTTGCCTGCCACCTTAGCTTTGACCTTTAAAATGAGAAAGG
GCTTATAGTTTTGGCCCTTGGAATTATTCTATCTGGTAAAGTGCTTAGAGAACCCTACAAAGAAGTCAGT

If I take a chunk of that sequence and blat it at UCSC then the top hit is (truncated)

   ACTIONS      QUERY                           SCORE START   END QSIZE IDENTITY  CHROM                 STRAND  START       END   SPAN
--------------------------------------------------------------------------------------------------------------------------------------
browser details NC_000001.11:c37595985-37566812  2310     1  2310  2310   100.0%  chr1                  -    37593676  37595985   2310

The alignment is on the - strand (truncated)

00000001 agtcgctctgaggccgagagggacgcgagcggaagtgacgtacgtcgtgc 00000050
<<<<<<<< |||||||||||||||||||||||||||||||||||||||||||||||||| <<<<<<<<
37595985 agtcgctctgaggccgagagggacgcgagcggaagtgacgtacgtcgtgc 37595936

00000051 acacgtggtccggcgtggttcaggcgggtgtcttcggccgggcttgggaa 00000100
<<<<<<<< |||||||||||||||||||||||||||||||||||||||||||||||||| <<<<<<<<
37595935 acacgtggtccggcgtggttcaggcgggtgtcttcggccgggcttgggaa 37595886

00000101 cataaaagtttgtttcaccacgtaagccggacctcgcactccggtcccgg 00000150
<<<<<<<< |||||||||||||||||||||||||||||||||||||||||||||||||| <<<<<<<<
37595885 cataaaagtttgtttcaccacgtaagccggacctcgcactccggtcccgg 37595836

00000151 tctcgtcgccaagatggtgaagcccaagtacaaaggacggagcaccatca 00000200
<<<<<<<< |||||||||||||||||||||||||||||||||||||||||||||||||| <<<<<<<<
37595835 tctcgtcgccaagatggtgaagcccaagtacaaaggacggagcaccatca 37595786

00000201 acccgtccaaggccagcacaaacccaggtacaggcggcagggctgatcga 00000250
<<<<<<<< |||||||||||||||||||||||||||||||||||||||||||||||||| <<<<<<<<
37595785 acccgtccaaggccagcacaaacccaggtacaggcggcagggctgatcga 37595736

00000251 gggtggagcttggaggaagtatagggaagggctccgaggagcactggagg 00000300
<<<<<<<< |||||||||||||||||||||||||||||||||||||||||||||||||| <<<<<<<<
37595735 gggtggagcttggaggaagtatagggaagggctccgaggagcactggagg 37595686

00000301 ggaaggagggttactcgagtggcattgcttagaatgactacgagcggggt 00000350
<<<<<<<< |||||||||||||||||||||||||||||||||||||||||||||||||| <<<<<<<<
37595685 ggaaggagggttactcgagtggcattgcttagaatgactacgagcggggt 37595636

[GenoMax]

If you remove strand=true in the URL then you get this sequence (truncated):

>NC_000001.11:37566812-37595985 Homo sapiens chromosome 1, GRCh38.p13 Primary Assembly
TAACTGCCTGTTTTTGTATAATTTAATAAAAACCTTTTAAACATTACTGCTTTTGTCTGAATTTTTTGCG
TTTGTGTTTCTGTCCCTCTGAGTCATTGGTCTTCTTTTTGTTCCTGTTCCTATTTTTCACGTTGTGTGTT
TCATAGTAGCGCACACCAACTTTTTTCGGCCGTTGCTGTCGTACTGCTCGCCTCCGCTGTAAAAAATGGC
AAGAAAAGATGAAGAAAAAAAACAACAAAACCCTGAAAGTCATTCACAAAACAGGAGTCTCTGCTTCTCA
TCTCTGTGTTCTATTTCCCGTGCTGCTTCCACAGCTACTGGAAGCAGTGCTAATGAGGGAGGGCTCCCCC
GTAGACCTAAGGTTACTCTTCCATTGTGACTTATGCAACTCCAAGTTACAAATCCACATGTGATCCAACC
CTTAGCCAGTGCAAGAATCTGTTCTGCATCACTCCTGACAAGTCCCCTGCAGCAAGGCTGAGTGCTGCTA
ACCACTCCCCTCCCCAGAGGGCCCACTCTAGCTCTGAACAGCTCAGCTCCAAAGCTGTCCTTCACAGAAG

Which when blatted at UCSC produces a hit that starts at right base pair (truncated query)

   ACTIONS      QUERY                          SCORE START   END QSIZE IDENTITY  CHROM                 STRAND  START       END   SPAN
-------------------------------------------------------------------------------------------------------------------------------------
browser details NC_000001.11:37566812-37595985  2800     1  2800  2800   100.0%  chr1                  +    37566812  37569611   2800

And results in the alignment on + strand:

00000001 taactgcctgtttttgtataatttaataaaaaccttttaaacattactgc 00000050
>>>>>>>> |||||||||||||||||||||||||||||||||||||||||||||||||| >>>>>>>>
37566812 taactgcctgtttttgtataatttaataaaaaccttttaaacattactgc 37566861

00000051 ttttgtctgaattttttgcgtttgtgtttctgtccctctgagtcattggt 00000100
>>>>>>>> |||||||||||||||||||||||||||||||||||||||||||||||||| >>>>>>>>
37566862 ttttgtctgaattttttgcgtttgtgtttctgtccctctgagtcattggt 37566911

00000101 cttctttttgttcctgttcctatttttcacgttgtgtgtttcatagtagc 00000150
>>>>>>>> |||||||||||||||||||||||||||||||||||||||||||||||||| >>>>>>>>
37566912 cttctttttgttcctgttcctatttttcacgttgtgtgtttcatagtagc 37566961

00000151 gcacaccaacttttttcggccgttgctgtcgtactgctcgcctccgctgt 00000200
>>>>>>>> |||||||||||||||||||||||||||||||||||||||||||||||||| >>>>>>>>
37566962 gcacaccaacttttttcggccgttgctgtcgtactgctcgcctccgctgt 37567011

00000201 aaaaaatggcaagaaaagatgaagaaaaaaaacaacaaaaccctgaaagt 00000250
>>>>>>>> |||||||||||||||||||||||||||||||||||||||||||||||||| >>>>>>>>
37567012 aaaaaatggcaagaaaagatgaagaaaaaaaacaacaaaaccctgaaagt 37567061

00000251 cattcacaaaacaggagtctctgcttctcatctctgtgttctatttcccg 00000300
>>>>>>>> |||||||||||||||||||||||||||||||||||||||||||||||||| >>>>>>>>
37567062 cattcacaaaacaggagtctctgcttctcatctctgtgttctatttcccg 37567111

[Istvan Albert]

My guess would be that the python code is zero based coordinate system, the NCBI server is 1 based, thus you are off by one.

[Igor Filippov]

The one I get when using python starts as:

AACTGCCTGTTTTTGTATAATTTAATAAAAACCTTTTAAACAT

And the one at NCBI is:

AGTCGCTCTGAGGCCGAGAGGGACGC

Doesn't seem to be off just by one, unless I'm missing something.

Answer 1

[Istvan Albert]

You are off by one in the Python representation, but in addition you are also producing the sequence as the reverse complement, a double whammy. Here is what the sequences may look depending on coordinate systems and orientaton:

1. Zero based coordinates, forward strand.

efetch -db nuccore -chr_start 37566812 -chr_stop 37595985  -format fasta -id NC_000001.11 | head -2

>NC_000001.11:37566813-37595986 Homo sapiens chromosome 1, GRCh38.p13 Primary Assembly
AACTGCCTGTTTTTGTATAATTTAATAAAAACCTTTTAAACATTACTGCTTTTGTCTGAATTTTTTGCGT

2. One based coordinates, forward strand.

efetch -db nuccore -seq_start 37566812 -seq_stop 37595985  -format fasta -id NC_000001.11 | head -2

>NC_000001.11:37566812-37595985 Homo sapiens chromosome 1, GRCh38.p13 Primary Assembly
TAACTGCCTGTTTTTGTATAATTTAATAAAAACCTTTTAAACATTACTGCTTTTGTCTGAATTTTTTGCG

3. Zero based coordinates, reverse complement.

efetch -db nuccore -chr_start 37566812 -chr_stop 37595985  -format fasta -id NC_000001.11 -strand 2 | head -2

>NC_000001.11:c37595986-37566813 Homo sapiens chromosome 1, GRCh38.p13 Primary Assembly
AAGTCGCTCTGAGGCCGAGAGGGACGCGAGCGGAAGTGACGTACGTCGTGCACACGTGGTCCGGCGTGGT

4. One based coordinates, reverse complement.

efetch -db nuccore -seq_start 37566812 -seq_stop 37595985  -format fasta -id NC_000001.11 -strand 2 | head -2

>NC_000001.11:c37595985-37566812 Homo sapiens chromosome 1, GRCh38.p13 Primary Assembly
AGTCGCTCTGAGGCCGAGAGGGACGCGAGCGGAAGTGACGTACGTCGTGCACACGTGGTCCGGCGTGGTT

Answer 2

[GenoMax]

I think the real answer is on this page, excerpted below. My comments above illustrate this.

flip=true - Flips the strand. Default is false. (Example: flip=true or flip=false) . Optional parameter.
strand=true - Synonym for flip parameter. Optional parameter