This is a test version of Biostars. For the public version, visit https://www.biostars.org.
How to parse a multiple sequence alignment to quantify the number of SNPs in different regions

Hello All,

I am new to programming and would like to figure out a way to count the number of SNPs inside and outside of a given probe-to-read match, separately. There are two main files I am working with: 1) a FASTA alignment for >2000 loci across 16 individuals and 2) a BLAST output table. Columns G and H from the BLAST table (see below) indicate where the probe matches the read from start to end (in base pairs). Any insight as to how I could tackle this problem would be greatly appreciated.

Snippet from FASTA alignment:

>locus1_individual-1
------------------------------------AGAAAGACAACTCTAACACAGTGATATCCATTTCCATTAGGCAAATAAACAACTACAAATCATTCAGATTATGTAAATGCTCAATGATATTTGCCTCCCGGTATCGCAGAGTTAGGACTGTCCAAAGTTTTATAAGCATATTTTCATCAGAGAACCTACTAGAGCATTGATACATTTGATTTTCAGTTTTTTCTAATGACAAATCACCGTCTTGGAAGCAAAAAGGCAAATAACTGATTATAAGTGCAATGAAAGTCTCTTGAATTTACATATACATTTACTTGTACCCTGTCACATTTTCTTAACCTCTAAACCTCCACACAACATATAAATGAACTTGTCTGGAATTGTTTTTTTCCCAAATTGCAGCCTATATTTTCTGCAGCACTGATGCAATTATGCCACCTCAATTGCCTCGCAGTGGAAAAAAGCACAATGCAAAACTTAGTGTATCAAGAAAATAAAGCCCAGAAAATTGCCTGTGACACACACCAACACAAAGACAGCTACGCTTCCTCAGAAAGTGGACAGAATAAAAGTAGACAGTTTCTATTAACAATAACAAGACAAAATACAGAGGCCCTTTATCTCTTAAAATGTAAATTTATACAAGGATTCATGCTGCCCTGCTGTAAATATCAAATACCCCTTAAGACCATAGAACAAGAATGAATTGATTTTCCTCCTAGTCTTCGCCATCTGCAGATTCTCATAGCTTGGCATTATATC
>locus1_individual-2
TAGGTAGTGACTGATGAGAAGTGAACTGCATCCCACAGAAAGACAACTCTAACACAGTGATATCCATTTCCATTAGGCAAATAAACAACTACAAATCATTCAGATTATGTAAATGCTCAATGATATTTGCCTCCCGGTATCGCAGAGTTAGGACTGTCCAAAGTTTTATAAGCATATTTTCATCAGAGAACCTACTAGAGCATTGATACATTTGATTTTCAGTTTTTTCTAATGACAAATCACCGTCTTGGAAGCAAAAAGGCAAATAACTGATTATAAGTGCAATGAAAGTCTCTTGAATTTACATATACATTTACTTGTACCCTGTCACATTTTCTTAACCTCTAAACCTCCACACAACATATAAATGAACTTGTCTGGAATTGTTTTTTTCCCAAATTGCAGCCTATATTTTCTGCAGCACTGATGCAATTATGCCACCTCAATTGCCTCGCAGTGGAAAAAAGCACAATGCAAAACTTAGTGTATCAAGAAAATAAAGCCCAGAAAATTGCCTGTGACACACACCAACACAAAGACAGCTACGCTTCCTCAGAAAGTGGACAGAATAAAAGTAGACAGT--------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------
.....
>locus1_individual-16
.....
>locus2000_individual-1
.....

Snippet from BLAST output:

A                   B               C   D   E   F   G   H   I   J   K           L
locus1_individual-1 locus1_probe1   100 120 0   0   240 359 1   120 3.00E-59    222
locus1_individual-2 locus1_probe1   100 120 0   0   276 395 1   120 3.00E-59    222

Probe matching locus 1:

>locus1_probe1
ATAAGTGCAATGAAAGTCTCTTGAATTTACATATACATTTACTTGTACCCTGTCACATTTTCTTAACCTCTAAACCTCCACACAACATATAAATGAACTTGTCTGGAATTGTTTTTTTCC
snp sequence next-gen blast alignment
• 3,416 views • link
updated by Ram Ram • written by modernsynthesis • modernsynthesis

Earlier today I did the galaxy 101 (@ https://usegalaxy.org/) tutorial and it seems they go over a simplified version of your problem using their web application.

Here's the link. https://usegalaxy.org/u/aun1/p/galaxy101

I apologize if I am missing a sublety in your problem, but either way this will not take long to investigate.

Good Hunting

• 0 views • link
updated by Ram Ram • written by dmccoll • dmccoll

Is this something I could do with sed and awk, or some other command-line tool? I will check out Galaxy, thanks.

• 1 views • link
written by modernsynthesis • modernsynthesis

0 answers

No answers yet.

Log in to answer this question.

More posts like this