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?
•
written
by
Igor Filippov •
2 answers
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
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
Log in to answer this question.
-
Help with Biopython for Beginner
written by cput •
-
bio python, sequence retrieval
written by Balatheskulter •
-
look for sequences containing a specific motif
written by Florian BERNARD •
-
Error in pysamstats
written by PAn •
-
Best Way To Preprocess Barcoded Illumina Paired-End Data
written by Ying W
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.
The one I get when using python starts as:
And the one at NCBI is:
Doesn't seem to be off just by one, unless I'm missing something.
Your link above (with
strand=true) brings back sequence that is this (abbreviated):If I take a chunk of that sequence and blat it at UCSC then the top hit is (truncated)
The alignment is on the
-strand (truncated)If you remove
strand=truein the URL then you get this sequence (truncated):Which when blatted at UCSC produces a hit that starts at right base pair (truncated query)
And results in the alignment on
+strand: