0
$\begingroup$

I have two FASTA sequences. Each is from a different species' mitochondrial genome. The sequences are pasted below.

>speciesA_control_region_partial_sequence_mitochondrial
TTAAAGTGATGTAATCCTACATTAATGAAAAATCAAAAGTTATAAGAACTTAAATACATTATATCATCAAATAAATATGAAGGTAGACATAAACCACTAAACTTCAAACTCCATTAAATATGTTATAAAAATGACGATATTGAATTGTCCTATCATAACTCTCATCAGTCTAGATATACCAGGACTCACCACCTCTGCAAGTCAGAGTCAAATGTAGTAAGAGACCACCATCAGTTGATTTCTTAATGTACACGTTTATTGATGGTCAAGGACAAAAATCGTGGGGGTAGCACACAGTGAACTATTCCTGGCATTTGGTTCCTATTTCAGGTTCATCAATTGAATACATTCCCCTAACATTTCCTGACGCTTGCATAAGTTAATGGTGGAGTACATACTCCTCGTTACCCACCATGCCGAGCATTCACTCTATAGGGCTACTGGTATCTTTTTTCTATTTCCTTTCACTTGACATTTCAGAGTGCATACAGAAATGTTATATTAAGGTTGAACATTTCCTTGCCCGGCGGAAAATGTATGAATGTAATTAAACTTTTATTTATGAATTGCATAACTGATATCAAGAGCATAAAGTATTATTTTTTCCCCTGACTTTCCTGTCAAGCACCCCTCTCGGCTTTTACGGGCCAAAACCCCCCCTCCCCCCTAAACTCGTGAGATTCCATTGTGTCTGCAAACCCCCCGGAAACAGAGCAAATCTTACTAGTTTTATCCCTTCTCAAATTTTTGTTCACTTGTATTTTTA
>speciesB_control_region_partial_sequence_mitochondrial
TTAAAGTAATGTAACCCTACATTAATGAAAAATCAAAA-TTATAGGAACTTAAATACATTACACCATCAAATAAATATGAAGGTAGACATAAACCATTGAATTTTAAATTTCATTAAACATGTTAAAAAAATGACGATATTGAATTGCCCTATCACAACTCTCATCAGTCTAGATATACCAGGACTCACCACCTCTGCAAGTAAGAGTCAAATGCAGTAAGAGACCACCATCAGTTGATTTCTTAATGTACACGTTTATTGATGGTCACGGACAAAAATCGTGGGGGTCGCACACAGTGAACTATTCCTGGCATTTGGTTCCTATTTCAGGTTCATCAATTGAATACATTCCCCTAACATTTCTTGACGCTTGCATAAGTTAATGGTGGAGTACATACTCCTCGTTACCCACCATGCCGAGCCTTCACTCTATAGGGCTACTGGTATCTTTTTTCTATTTCCTTTCACTTGACATTTCAGAGTGCATACAGAAATGTTATATCAAGGTTGAACATTTCCTTGCCCGGCGGAAAATGTATGAATGTAATTAAACTTTTATTTATGAATTGCATAACTGATATCAAGAGCATAAAGTATCATTTTTTCCCCTGACTTTCCTATCAAGCACCCCTCTCGGCTTTTGCGGGCCAAAACCCCCCCTCCCCCCTAAACTCGTGAGATTCCATTGTGTCTGCAAACCCCCCGGAAACAGAGCAAATCTTACTAGTTTTATCCCTTCTCAAATTTTTGTTCACTTGTATTTTTA

My goal is to extract all SNPs in these sequences and format them in SNP FASTA format, as shown below:

>rs11540652
GAACC[A/G]GAGGCCC

I specifically need this format (rather than, say, VCF) so I can search for RFLP sites using the following website: http://bio.kuas.edu.tw/snp-rflping2/rflpUI.jsp

Does anyone know of an existing program which will accomplish this? Also, what additional materials do I need to detect these SNPs? I presume that I need assembled mitochondrial reference genomes for each species - anything else?

Improve this question
$\endgroup$
1
  • $\begingroup$ Do you have a VCF file? It should be possible to compute one using bcftools and it would give the position in the reference sequence so generating these “SNP” sequences could be done using this and the reference sequence. The full sequence for both samples will take longer to compute and is not necessary. $\endgroup$
    – Tom Kelly ケリー・トム
    Commented Mar 12, 2021 at 15:08

1 Answer 1

Reset to default
2
$\begingroup$

Your question is not very clear. Where does the flanking sequence in :

>rs11540652
GAACC[A/G]GAGGCCC

come from? I'll assume its meant to be the reference sequence and that speciesA is the reference and that the first nuc in the [] is from the reference... I don't know an existing existing program but here's a bespoke one.

speciesA='TTAAAGTGATGTAATCCTACATTAATGAAAAATCAAAAGTTATAAGAACTTAAATACATTATATCATCAAATAAATATGAAGGTAGACATAAACCACTAAACTTCAAACTCCATTAAATATGTTATAAAAATGACGATATTGAATTGTCCTATCATAACTCTCATCAGTCTAGATATACCAGGACTCACCACCTCTGCAAGTCAGAGTCAAATGTAGTAAGAGACCACCATCAGTTGATTTCTTAATGTACACGTTTATTGATGGTCAAGGACAAAAATCGTGGGGGTAGCACACAGTGAACTATTCCTGGCATTTGGTTCCTATTTCAGGTTCATCAATTGAATACATTCCCCTAACATTTCCTGACGCTTGCATAAGTTAATGGTGGAGTACATACTCCTCGTTACCCACCATGCCGAGCATTCACTCTATAGGGCTACTGGTATCTTTTTTCTATTTCCTTTCACTTGACATTTCAGAGTGCATACAGAAATGTTATATTAAGGTTGAACATTTCCTTGCCCGGCGGAAAATGTATGAATGTAATTAAACTTTTATTTATGAATTGCATAACTGATATCAAGAGCATAAAGTATTATTTTTTCCCCTGACTTTCCTGTCAAGCACCCCTCTCGGCTTTTACGGGCCAAAACCCCCCCTCCCCCCTAAACTCGTGAGATTCCATTGTGTCTGCAAACCCCCCGGAAACAGAGCAAATCTTACTAGTTTTATCCCTTCTCAAATTTTTGTTCACTTGTATTTTTA'
speciesB='TTAAAGTAATGTAACCCTACATTAATGAAAAATCAAAA-TTATAGGAACTTAAATACATTACACCATCAAATAAATATGAAGGTAGACATAAACCATTGAATTTTAAATTTCATTAAACATGTTAAAAAAATGACGATATTGAATTGCCCTATCACAACTCTCATCAGTCTAGATATACCAGGACTCACCACCTCTGCAAGTAAGAGTCAAATGCAGTAAGAGACCACCATCAGTTGATTTCTTAATGTACACGTTTATTGATGGTCACGGACAAAAATCGTGGGGGTCGCACACAGTGAACTATTCCTGGCATTTGGTTCCTATTTCAGGTTCATCAATTGAATACATTCCCCTAACATTTCTTGACGCTTGCATAAGTTAATGGTGGAGTACATACTCCTCGTTACCCACCATGCCGAGCCTTCACTCTATAGGGCTACTGGTATCTTTTTTCTATTTCCTTTCACTTGACATTTCAGAGTGCATACAGAAATGTTATATCAAGGTTGAACATTTCCTTGCCCGGCGGAAAATGTATGAATGTAATTAAACTTTTATTTATGAATTGCATAACTGATATCAAGAGCATAAAGTATCATTTTTTCCCCTGACTTTCCTATCAAGCACCCCTCTCGGCTTTTGCGGGCCAAAACCCCCCCTCCCCCCTAAACTCGTGAGATTCCATTGTGTCTGCAAACCCCCCGGAAACAGAGCAAATCTTACTAGTTTTATCCCTTCTCAAATTTTTGTTCACTTGTATTTTTA'
counter=0
for i, nuc in enumerate(speciesA):
    if nuc != speciesB[i]:
        counter+=1
        print ('>rs'+str(counter))
        print (f'{speciesA[i-5:i]}[{nuc}/{speciesB[i]}]{speciesA[i+1:i+6]}')

Gives:

>rs1
AAAGT[G/A]ATGTA
>rs2
TGTAA[T/C]CCTAC
>rs3
CAAAA[G/-]TTATA
>rs4
TTATA[A/G]GAACT
>rs5
CATTA[T/C]ATCAT
>rs6
TTATA[T/C]CATCA
>rs7
AACCA[C/T]TAAAC
>rs8
CCACT[A/G]AACTT
>rs9
CTAAA[C/T]TTCAA
>rs10
AACTT[C/T]AAACT
>rs11
TCAAA[C/T]TCCAT
...

As for 'I presume that I need assembled mitochondrial reference genomes for each species - anything else?' These sequences have clearly been aligned. So there's that.

I hope that helps.

Improve this answer
$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.