This is a test version of Biostars. For the public version, visit https://www.biostars.org.
Extract DNA sequence based on a single coordinate

Hello,

I know many have asked about extracting/fetching sequences from a range of coordinates. However, I would like to extract 10bp on the left and the right side from a specific position.

Let's say I have a genomic position at 50000. I would like to extract 10bp on each side, so from 49990-50010. Any suggestions how to do it?

Thank you,

TP

extract seq r python
• 5,342 views • link
written by ThePresident ThePresident

Do you need this just once or for multiple positions?

• 1 views • link
written by WouterDeCoster WouterDeCoster

Yup, multiple positions (coming from a screen)

• 1 views • link
written by ThePresident ThePresident

1 answer

I'm assuming you mean from a FASTA file. The easiest way to do this is to create a BED based on what you need (remember: 0-based indexing!).

Once you're set, you'll want to use bedtools getfasta.

This is also super quick to do for multiple files or multiple sequence alignments using Biostrings in R. Something like:

library(Biostrings)
a <- readDNAStringSet("your_file.fa")
b <- subseq(a, start=start, stop=stop)
writeXStringSet(b, file="new.fa")
• 1 views • link
written by Brice Sarver Brice Sarver

Update: I did it with bedtools. For those who might be asking the same question, here's my procedure:

  • I manually created a minimal BED file (tab-delimited txt file) with genome name (matching fasta file header) and genomic range I want to extract
  • Run the following: bedtools getfasta -fo BED_output.txt -tab -fi REF.fasta -bed BED_input.txt
• 1 views • link
written by ThePresident ThePresident

Although this is an easier method, one thing to take care is the co-ordinates. If you have given co-ordinates as 1 to 10 for extraction, your final sequence will have 9 nucleotides only. It will miss the first one (the same has been also shown in their example). It took me some time to figure this out and had to do complete reanalysis.

• 0 views • link
written by ankit hinsu • ankit hinsu

Log in to answer this question.

More posts like this