Is There A Tool To Extract The Reference Genomic Sequence From A Given Coordinate?
2
0
Entering edit mode
8.9 years ago
bibb77 ▴ 90

Hello everyone, this is my problem:

I have a list of SNPs coordinates in this format (from chr 1 to chr Y, ~580k SNPs), they are under build GRCh37:

chr1:108681808
chr1:109440678
chr1:109479801
chr1:110655430
chr1:11193226
chr1:113933669
chr1:115258741
chr1:115527488
...


And I have the Reference Genome build37 in .fa format

>chr1
NNNNNNNNNNNNNNNNNNNNNNNNNNNNN
NNNNNNNNN....gctttatacaatctat
ttgttactttttattctattttgcatttt
gttcctttgcctgaataattcactttggt
ctgcaatggctaattgcaatagattt...


Is there tool to obtain the sequence corresponding to each given coordinate?

chr1:108681808  A
chr1:109440678  C
chr1:109479801  T
chr1:110655430  G
chr1:11193226    A
chr1:113933669  T
chr1:115258741  C
chr1:115527488  C


Hope you can help me

reference vcf coordinates • 4.2k views
2
Entering edit mode
1
Entering edit mode
8.9 years ago
Gabriel R. ★ 2.8k

You have a relatively small number of positions, I would just do

2. do a for loop in bash where you replace chr1:108681808 to chr1:108681808-108681808, a quick awk should wor, see below:
3. for each, do a faidx on the reference:

for i in cat file.pos |awk 'BEGIN{FS=":"}{print $1":"$2"-$2}' do echo -np$i"\t";

samtools faidx reference.fa \$i;

done

0
Entering edit mode
8.9 years ago
SRKR ▴ 180

What is the genome size that you are working with... I mean the sequence in the .fa file. If it's reasonably small. I can create an online tool for the purpose.