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
2
-
1$\begingroup$ Please do not cross-post to BioStars and Stack Exchange biostars.org/p/273588 $\endgroup$– Emily_EnsemblSep 20, 2017 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, 2017 at 8:07
Add a comment
|
3 Answers
$\begingroup$
$\endgroup$
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
$\begingroup$
$\endgroup$
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.
$\begingroup$
$\endgroup$
If you are comfortable with Linux command line, try Entrez Direct and GNU parallel.
First, I recommend acquiring a NCBI API key to increase the maximum requests per second.
After that, if you already have the accession numbers that you want to download, just create a tab delimited file (eg arg.tsv) contains 3 columns (no header), including accession number, sequence start and sequence stop.
To fetch the desired sequence region and save it in separate file in parallel:
parallel -a arg.tsv -C '\t' -j8 "efetch -db nuccore -id {1} -seq_start {2} -seq_stop {3} -format fasta >{1}.fa"
Or you can use Biopython from Bio import Entrez. But I have not tried yet.
