Python: Find single-nucleotide polymorphism from aligned sequences

This short blog deal with the detection of SNPs from two aligned sequences

Brainstorming

I was working on comparative chloroplast genome analysis. My purpose was to find SNPs between two chloroplast genomes in order to appreciate sequence divergence among two species. Thanks to Biostar platform for the inspiration.

Solution

Assuming that I have two chloroplast genome assemblies cp1.fasta and cp2.fasta.

First step: Concatenate the two sequences

$ cat cp1.fasta cp2.fasta > cp1_cp2.fasta 

Second step: Multiple sequence alignment with MAFFT or other program

$ mafft --auto cp1_cp2.fata > cp1_cp2.mafft

Third step: Create the python script

$ vi find_snp.py

Fourth step:Scripting

from Bio import AlignIO

y=0
alignment = AlignIO.read("fasta.fas", "fasta")
for r in range(0,len(alignment[1].seq)):
    if alignment[0,r] != alignment[1,r]:
        if alignment[0,r] != "-" and alignment[1,r] != "-":
            y=y+1
            print (r, alignment[0,r], alignment[1,r], y)
        else:
            y=0

Fifth step: Run the script

$ python find_snp.py > snp.output.txt

A view of the result

$ more snp.output.txt

10 a t 1
11 a g 2
92 t c 3
113 a t 4
161 a c 1
234 c t 2
243 c t 1
246 a t 2

To quickly count the number of detected SNP:

$ cat snp.output.txt | wc -l

In summary, this returns position of SNP, nucleotide in cp1, nucleotide in cp2, running tally of the number of SNPs upstream of each indel.