Highlighting SNPs for better presentation of your sequence alignment

Imagine the following FASTA alignment:

>header1
ATGGATTTATCTGCTCTTCGCGTTGAAGAAGTACAAAATGTCATTAATGCTATGCAGAAAATCTTA
>header2
ATGGATTTATCTGCTCTTCGCGTTGAAGAAGTACAAAATGTCATTAATGCTATGCAGAAAATCTTA
>header3
ATGGATTTATCTGCTCTTCGCGTTGAAGGAGTACAAAATGTCATTAATGCTATGCAGAAAATCTTA
>header4
ATGGATTTATCTGTTCTTCGCGTTGAAGGAGTACAAAATGTCACTAATGCTATGCAGAAAATCTTA

Based only in the above alignment, is difficult to instantly identify where the SNPs are. Set dots where we have alignment matches is a very intuitive way to spotlight the SNPs.

header4 ATGGATTTATCTGTTCTTCGCGTTGAAGGAGTACAAAATGTCACTAATGCTATGCAGAAAATCTTA
header2 .............C..............A..............T......................
header3 .............C.............................T......................
header1 .............C..............A..............T......................

Changing similar nucleotides or amino acids by dots is the most simple way of arranging your sequence alignment in order to highlight the regions of inequality.

Here’s how it could be done in Python.

#!/usr/bin/python
#usage: python script.py file.fasta
import sys

filefasta = sys.argv[1]

try:
    f = open(filefasta)
except IOError:
    print ("File doesn't exist!")

seqs={}

for line in f:
    line = line.rstrip()
    if line[0] == '>':
        words=line.split()
        name=words[0][1:]
        seqs[name]=''
    else:
        seqs[name] = seqs[name] + line

first_header = seqs.keys()[-1]
first_sequence = list(seqs.values()[-1])

print first_header,''.join(first_sequence)

comparison = {}

for i in seqs.items():
    header= i[0]
    seq = i[1]
    if header != first_header:
        comparison[header] = ''
        for i in range(0,len(seq)):
            letter=seq[i]
            base=first_sequence[i]
            if letter == base:
                comparison[header] = comparison[header] + '.'
            else:
                comparison[header] = comparison[header] + letter

for i in comparison.items():
    header=i[0]
    seq=i[1]
    print header,seq

Advertisements
This entry was posted in Blog. Bookmark the permalink.

Leave a Reply

Fill in your details below or click an icon to log in:

WordPress.com Logo

You are commenting using your WordPress.com account. Log Out /  Change )

Google photo

You are commenting using your Google account. Log Out /  Change )

Twitter picture

You are commenting using your Twitter account. Log Out /  Change )

Facebook photo

You are commenting using your Facebook account. Log Out /  Change )

Connecting to %s