This is a test version of Biostars. For the public version, visit https://www.biostars.org.
Is There A Tool To Extract The Reference Genomic Sequence From A Given Coordinate?

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
• 5,244 views • link
updated by Gabriel R. Gabriel R. • written by bibb77 bibb77

2 answers

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

  1. put all your fastas into one, samtools faidx your reference
  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

• 35 views • link
written by Gabriel R. Gabriel R.

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.

• 0 views • link
written by SRKR SRKR

Log in to answer this question.

More posts like this