Incorrect sequence retrieved by coordinates
2
0
Entering edit mode
9 months ago

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?

python fasta biopython pysam • 273 views
ADD COMMENT
1
Entering edit mode

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
ADD REPLY
1
Entering edit mode

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
ADD REPLY
0
Entering edit mode

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.

ADD REPLY
0
Entering edit mode

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.

ADD REPLY
2
Entering edit mode
9 months ago

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
ADD COMMENT
2
Entering edit mode
9 months ago
GenoMax 99k

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
ADD COMMENT

Login before adding your answer.

Traffic: 1382 users visited in the last hour
Help About
FAQ
Access RSS
API
Stats

Use of this site constitutes acceptance of our User Agreement and Privacy Policy.

Powered by the version 2.3.6