Comparative Genomics

• Brent Pedersen
• Freeling Lab / Plant Biology

• plants undergo whole genome duplication (WGD).
• 3 genomes in 2 organisms. following polyploidy, lose genes where 2 copies arent needed.
• allow for comparative analyses.

Arabidopsis Thaliana/Lyrata

• BLAST: define, dot: define, box: define, axis define.
• clusters of dots are centromere == lots of crap / repetitive sequence.
• Arabidopsis is simpler fewer chrs. both have 2 genome dups.
• lose genes from 1 or other homeolog.
• data stored in Rtree, rendered with MPL.
• too much data for sqlite.
• number of blast points for sorghum vs. maize
• use other things for data storage like pytables.
• tend to do spatial queries on n-dimensional data.
• use in memory or mmap'ed npy array when possible.

Rtree

>>> from rtree import Rtree
>>> from biostuff import BlastLine
>>> r = Rtree('ex.rtree')
>>> for i, line in enumerate(open('data/tabd.blast')):
...     b = BlastLine(line)
...     smin, smax = sorted([b.sstart, b.sstop])
...     #     id, (xmin, ymin, xmax, ymax), object
...     r.add(i, (b.qstart, smin, b.qstop, smax), b)
>>> del r # flush for next block.

# somewhere else.
>>> from rtree import Rtree
>>> r = Rtree('ex.rtree')
>>> r.intersection((44384798, 9693297, 44384798, 9693297), objects="raw")
[BlastLine(chr6[44384798:44387498]-10[9693297:9696007], ptcid=93.540, eval=0.000)]

• Rtree can be in memory or disk based. _very_ fast.
• handles huge amounts of data (where sqlite failed)
• takes any object that can be pickled. (or define your own serializer).
• c++ lib marios h, hobu for c api, ctypes wrapper.
• well-tested, documented, on pypi. (version 0.6)
• N-dimensions, R-Tree/TPR/MVR/etc.
• can adjust properties leafsize, pagesize, pool_capacity etc.
• used to find features around gene-pair of interest.
• since n-dimensional, it can be used for n-wise pairs.
• can also do n-nearest neighbors

GenomeTools

• CLI for GFF feature extraction and filtering, LTR, suffixerator, etc.
• C API with bindings for Lua, Python, Ruby.

• genometools.org
• combines c api and CLI (python, ruby, lua bindings)
• nicely written "object-oriented" C with api.h files.
• cairo bindings with bumping for rendering.
• gff3 index. in memory.
• lots of bioinformatics frameworks in python but gt CLI very nice for shell scripting pipelines. pipe to various tools, visualize.

pyfasta: pythonic access to fasta sequence files

>>> from pyfasta import Fasta
>>> f = Fasta('data/three_chrs.fasta')
>>> sorted(f.keys())
['chr1', 'chr2', 'chr3']

>>> f['chr1'][:5] # normal python slicing.
'ACTGA'

>>> f['chr1'][2::3] # last basepair in every codon.
'TCAGTCAGTCAGTCAGTCAGTCAGTC'

>>> import numpy as np
>>> a = np.array(f['chr1'])
>>> a.dtype, a.shape
(dtype('|S1'), (80,))

>>> a[a == 'C'] = np.array('N', dtype='S1') # mask cytosines
>>> a[:10].tostring()
'ANTGANTGAN'

• what is 'fasta' (>headernsequence) can have newlines in the sequence. (so key=> value matches dict protocol)
• pythonic access. getitem. very simple, do one thing well.
• supports numpy array interface. useful for broadcasting mask repetitve sequence or coding sequences or non-coding.
• CLI for sequence extraction.
• on pypi