# ACE contig to alignment

## Problem

Sometimes it is useful to be able represent a contig produced as part of genome or EST assembly as an alignment (eg to search for potential SNPs in runs from mixed samples or to be able to write a contig out in a way it can be viewed more easily). For assemblies that use the ACE file format we can use Biopython's ACE handling to add reads that make a contig to a generic alignment.'

## Solution

Let's represent contig in the ACE file that is used in biopython's testing framework: ~/biopython/Tests/Ace/contig1.ace as an example

```from Bio.Sequencing import Ace
from Bio.Align.Generic import Alignment
from Bio.Alphabet import IUPAC, Gapped

'''Replace residues on either end of a sequence with gaps.

In this case we want to cut out the sections of each read which the assembler has
decided are not good enough to include in the contig and replace them with gaps
'''
return (start-1) * '-' + read[start-1:end] + (end +1) * '-'

''' Pad out either end of a read so it fits into an alignment.

The start argument is the position of the first base of the reads sequence in
the contig it is part of. If the start value is negative (or 0 since ACE
files count from 1, not 0) we need to take some sequence off the start
otherwise each end is padded to the length of the consensus with gaps.
'''
if start < 1:
else:
seq = (start-1) * '-' + read
seq = seq + (conlength-len(seq)) * '-'
return seq[:conlength]

# We will use the Ace parser to read individual contigs from file. Be aware
# that using this iterator can miss WA, CT, RT and WR tags (which can be
# anywhere in the file, e.g. the end). Read the file specification here:
# If you need these tags you'll need to use  Ace.read() (and lots of RAM).

ace_gen = Ace.parse(open("contig1.ace", 'r'))
contig = ace_gen.next()
align = Alignment(Gapped(IUPAC.ambiguous_dna, "-"))

# Now we have started our alignment we can add sequences to it,
# we will loop through contig's reads and get quality clipping from

and when you print the alignment, or the sequences within it

```>>>print align
Gapped(IUPACAmbiguousDNA(), '-') alignment with 2 rows and 856 columns

tt*gtagagggaTGCTTCTGGGTAGCGACCGGTCTAAGTTCCTCGGAACAGGACGTCATAGAGGGTGAGAATCCCGTAT
TTTGTAGAGG*ATGCTTCTGGGTAGCGACCGGTCTAAGTTCCTCGGAACAGGACGTCATAGAGGGTGAGAATCCCGTAT```

You can also now write the alignment to any format you want to using AlignIO.

## Discussion

The details are given in the comments above, in broad strokes the ACE contig is read in with Ace.parse(), a generic alignment is started then the reads from the contig are added to the new alignment.