5
$\begingroup$

So I have a list of start and stop positions along chromosomes in different species, and I'd like to get the corresponding DNA sequence for each set of coordinates. In the past, I've just download the genome as a fasta file and then use pyfaidx to extract the sequences at the given positions. But now that I'm working with several species at once, I was wondering if there's any kind of tool in Python or R that can fetch your sequences of interest without downloading a bunch of large files. Thanks

$\endgroup$
  • 1
    $\begingroup$ Please do not cross-post to BioStars and Stack Exchange biostars.org/p/273588 $\endgroup$ – Emily_Ensembl Sep 20 '17 at 7:42
  • 2
    $\begingroup$ @Emily_Ensembl it would be great if you could flesh out your answer from BioStars so it could be posted as an answer here! Maybe with a simple example like Pierre has done below? $\endgroup$ – terdon Sep 20 '17 at 8:07
8
$\begingroup$

using a http request.

if there is a DAS server, you can always use this protocol to download the xml -> fasta. see https://www.biostars.org/p/56/

$ curl -s "http://genome.ucsc.edu/cgi-bin/das/hg19/dna?segment=chrM:100,200"  | xmllint --xpath '/DASDNA/SEQUENCE/DNA/text()' - | tr -d '\n'
ggagccggagcaccctatgtcgcagtatctgtctttgattcctgcctcattctattatttatcgcacctacgttcaatattacaggcgaacatacctacta

or use the UCSC utility twoBitToFa which works with remote files.

If you have an NCBI accession number you can use the E-utilities with seq_start and seq_end.

$ wget -q  -O -  "https://eutils.ncbi.nlm.nih.gov/entrez/eutils/efetch.fcgi?db=nuccore&id=AE014134.1&seq_start=100&seq_stop=300&rettype=fasta"
>AE014134.1:100-300 Drosophila melanogaster chromosome 2L complete sequence
TGCCAACATATTGTGCTCTTTGATTTTTTGGCAACCCAAAATGGTGGCGGATGAACGAGATGATAATATA
TTCAAGTTGCCGCTAATCAGAAATAAATTCATTGCAACGTTAAATACAGCACAATATATGATCGCGTATG
CGAGAGTAGTGCCAACATATTGTGCTAATGAGTGCCTCTCGTTCTCTGTCTTATATTACCG
$\endgroup$
1
$\begingroup$

I am not aware of any tools that offer this functionality out-of-the-box, although that doesn't necessarily mean it doesn't exist.

The tool you describe would have to be some kind of web service, where you submit a request with something like the following:

  • species
  • assembly version
  • list of intervals, each with:
    • chromosome ID
    • start position
    • end position

The response would be one or more subsequences in Fasta format that you could save to a file.

This approach pre-supposes that there will be a relatively small number of subsequences, that those subsequences won't be too large, and you won't be requesting them very frequently. Otherwise, this approach will offer very little practical benefit over your current approach.

$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.