I have a long list of chromosome positions in the format below, except there are 15,000 of them. I would like to get the sequence for each of the spans (or alternatively 150bp up from the start position). Can anyone tell me how to do this for this number of positions?
chr15:60212080-60213230
chr10:60242850-60243000
chr11:60469240-60469390
chr19:60954240-60954390
chr12:61260820-61260970
chr21:61576770-61576920
chr1:61586420-61586570
chr1:61927840-61927990
Thanks in advance
•
written
by
sebastiz •
4 answers
K hopefully this should be easy enough.
Follow these steps.
samtools faidx <reference genome>and a little bash loop.
while read region do samtools faidx <reference genome> $region >> output_file.txt done < FILE_REGIONS
•
written
by
smilefreak
Using bedtools getfasta. Input is: chromosome positions & reference genome ( https://www.biostars.org/p/1796/ ).
Using Bioconductor/R:
- create GRanges object containing your ranges
- use the appropriate BSgenome package (or FaFile) and
getSeq()
Using DAS after replacing "-" with ",". This URL gives an example result:
http://genome.ucsc.edu/cgi-bin/das/hg19/dna?segment=chr1:100000,200000
Many other possiblities!
•
written
by
Sean Davis
Using Python, you can use the pyfaidx module:
pip install --user pyfaidx
faidx genome.fasta --bed regions.bed
Log in to answer this question.