{
"cells": [
{
"cell_type": "markdown",
"metadata": {},
"source": [
"**Source of the materials**: Biopython cookbook (adapted)\n",
"Status: Draft"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Sequence motif analysis using Bio.motifs\n",
"========================================\n",
"\n",
"This chapter gives an overview of the functionality of the `Bio.motifs`\n",
"package included in Biopython. It is intended for people who are\n",
"involved in the analysis of sequence motifs, so I’ll assume that you are\n",
"familiar with basic notions of motif analysis. In case something is\n",
"unclear, please look at Section \\[sec:links\\] for some relevant links.\n",
"\n",
"Most of this chapter describes the new `Bio.motifs` package included in\n",
"Biopython 1.61 onwards, which is replacing the older `Bio.Motif` package\n",
"introduced with Biopython 1.50, which was in turn based on two older\n",
"former Biopython modules, `Bio.AlignAce` and `Bio.MEME`. It provides\n",
"most of their functionality with a unified motif object implementation.\n",
"\n",
"Speaking of other libraries, if you are reading this you might be\n",
"interested in [TAMO](http://fraenkel.mit.edu/TAMO/), another python\n",
"library designed to deal with sequence motifs. It supports more\n",
"*de-novo* motif finders, but it is not a part of Biopython and has some\n",
"restrictions on commercial use.\n",
"\n",
"Motif objects {#sec:object}\n",
"-------------\n",
"\n",
"Since we are interested in motif analysis, we need to take a look at\n",
"`Motif` objects in the first place. For that we need to import the\n",
"Bio.motifs library:\n",
"\n"
]
},
{
"cell_type": "code",
"execution_count": 1,
"metadata": {
"collapsed": true
},
"outputs": [],
"source": [
"from Bio import motifs"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"\n",
"and we can start creating our first motif objects. We can either create\n",
"a `Motif` object from a list of instances of the motif, or we can obtain\n",
"a `Motif` object by parsing a file from a motif database or motif\n",
"finding software.\n",
"\n",
"### Creating a motif from instances\n",
"\n",
"Suppose we have these instances of a DNA motif:\n",
"\n"
]
},
{
"cell_type": "code",
"execution_count": 2,
"metadata": {
"collapsed": true
},
"outputs": [],
"source": [
"from Bio.Seq import Seq\n",
"instances = [Seq(\"TACAA\"),\n",
" Seq(\"TACGC\"),\n",
" Seq(\"TACAC\"),\n",
" Seq(\"TACCC\"),\n",
" Seq(\"AACCC\"),\n",
" Seq(\"AATGC\"),\n",
" Seq(\"AATGC\")]"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"\n",
"then we can create a Motif object as follows:\n",
"\n"
]
},
{
"cell_type": "code",
"execution_count": 3,
"metadata": {
"collapsed": true
},
"outputs": [],
"source": [
"m = motifs.create(instances)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"\n",
"The instances are saved in an attribute `m.instances`, which is\n",
"essentially a Python list with some added functionality, as described\n",
"below. Printing out the Motif object shows the instances from which it\n",
"was constructed:\n",
"\n"
]
},
{
"cell_type": "code",
"execution_count": 4,
"metadata": {
"collapsed": false
},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"TACAA\n",
"TACGC\n",
"TACAC\n",
"TACCC\n",
"AACCC\n",
"AATGC\n",
"AATGC\n",
"\n"
]
}
],
"source": [
"print(m)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"\n",
"The length of the motif is defined as the sequence length, which should\n",
"be the same for all instances:\n",
"\n"
]
},
{
"cell_type": "code",
"execution_count": 5,
"metadata": {
"collapsed": false
},
"outputs": [
{
"data": {
"text/plain": [
"5"
]
},
"execution_count": 5,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"len(m)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"\n",
"The Motif object has an attribute `.counts` containing the counts of\n",
"each nucleotide at each position. Printing this counts matrix shows it\n",
"in an easily readable format:\n",
"\n"
]
},
{
"cell_type": "code",
"execution_count": 6,
"metadata": {
"collapsed": false
},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
" 0 1 2 3 4\n",
"A: 3.00 7.00 0.00 2.00 1.00\n",
"C: 0.00 0.00 5.00 2.00 6.00\n",
"G: 0.00 0.00 0.00 3.00 0.00\n",
"T: 4.00 0.00 2.00 0.00 0.00\n",
"\n"
]
}
],
"source": [
"print(m.counts)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"\n",
"You can access these counts as a dictionary:\n",
"\n"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": true
},
"outputs": [],
"source": [
"m.counts['A']"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"\n",
"but you can also think of it as a 2D array with the nucleotide as the\n",
"first dimension and the position as the second dimension:\n",
"\n"
]
},
{
"cell_type": "code",
"execution_count": 7,
"metadata": {
"collapsed": false
},
"outputs": [
{
"data": {
"text/plain": [
"4"
]
},
"execution_count": 7,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"m.counts['T', 0]"
]
},
{
"cell_type": "code",
"execution_count": 8,
"metadata": {
"collapsed": false
},
"outputs": [
{
"data": {
"text/plain": [
"2"
]
},
"execution_count": 8,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"m.counts['T', 2]"
]
},
{
"cell_type": "code",
"execution_count": 9,
"metadata": {
"collapsed": false
},
"outputs": [
{
"data": {
"text/plain": [
"0"
]
},
"execution_count": 9,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"m.counts['T', 3]"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"\n",
"You can also directly access columns of the counts matrix\n",
"\n"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": true
},
"outputs": [],
"source": [
"m.counts[:, 3]"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"\n",
"Instead of the nucleotide itself, you can also use the index of the\n",
"nucleotide in the sorted letters in the alphabet of the motif:\n",
"\n"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": true
},
"outputs": [],
"source": [
"m.alphabet"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": true
},
"outputs": [],
"source": [
"m.alphabet.letters"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": true
},
"outputs": [],
"source": [
"sorted(m.alphabet.letters)"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": true
},
"outputs": [],
"source": [
"m.counts['A',:]"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": true
},
"outputs": [],
"source": [
"m.counts[0,:]"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"\n",
"The motif has an associated consensus sequence, defined as the sequence\n",
"of letters along the positions of the motif for which the largest value\n",
"in the corresponding columns of the `.counts` matrix is obtained:\n",
"\n"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": true
},
"outputs": [],
"source": [
"m.consensus"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"\n",
"as well as an anticonsensus sequence, corresponding to the smallest\n",
"values in the columns of the `.counts` matrix:\n",
"\n"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": true
},
"outputs": [],
"source": [
"m.anticonsensus"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"\n",
"You can also ask for a degenerate consensus sequence, in which ambiguous\n",
"nucleotides are used for positions where there are multiple nucleotides\n",
"with high counts:\n",
"\n"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": true
},
"outputs": [],
"source": [
"m.degenerate_consensus"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"\n",
"Here, W and R follow the IUPAC nucleotide ambiguity codes: W is either A\n",
"or T, and V is A, C, or G @cornish1985. The degenerate consensus\n",
"sequence is constructed following the rules specified by Cavener\n",
"@cavener1987.\n",
"\n",
"We can also get the reverse complement of a motif:\n",
"\n"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": true
},
"outputs": [],
"source": [
"r = m.reverse_complement()\n",
"r.consensus"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": true
},
"outputs": [],
"source": [
"r.degenerate_consensus"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": true
},
"outputs": [],
"source": [
"print(r)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"\n",
"The reverse complement and the degenerate consensus sequence are only\n",
"defined for DNA motifs.\n",
"\n",
"### Creating a sequence logo\n",
"\n",
"If we have internet access, we can create a\n",
"[weblogo](http://weblogo.berkeley.edu):\n",
"\n"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": true
},
"outputs": [],
"source": [
"m.weblogo(\"mymotif.png\")"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"\n",
"We should get our logo saved as a PNG in the specified file.\n",
"\n",
"Reading motifs {#sec:io}\n",
"--------------\n",
"\n",
"Creating motifs from instances by hand is a bit boring, so it’s useful\n",
"to have some I/O functions for reading and writing motifs. There are not\n",
"any really well established standards for storing motifs, but there are\n",
"a couple of formats that are more used than others.\n",
"\n",
"### JASPAR\n",
"\n",
"One of the most popular motif databases is\n",
"[JASPAR](http://jaspar.genereg.net). In addition to the motif sequence\n",
"information, the JASPAR database stores a lot of meta-information for\n",
"each motif. The module `Bio.motifs` contains a specialized class\n",
"`jaspar.Motif` in which this meta-information is represented as\n",
"attributes:\n",
"\n",
"- `matrix_id` - the unique JASPAR motif ID, e.g. ’MA0004.1’\n",
"\n",
"- `name` - the name of the TF, e.g. ’Arnt’\n",
"\n",
"- `collection` - the JASPAR collection to which the motif\n",
" belongs, e.g. ’CORE’\n",
"\n",
"- `tf_class` - the structual class of this TF, e.g. ’Zipper-Type’\n",
"\n",
"- `tf_family` - the family to which this TF belongs, e.g.\n",
" ’Helix-Loop-Helix’\n",
"\n",
"- `species` - the species to which this TF belongs, may have multiple\n",
" values, these are specified as taxonomy IDs, e.g. 10090\n",
"\n",
"- `tax_group` - the taxonomic supergroup to which this motif\n",
" belongs, e.g. ’vertebrates’\n",
"\n",
"- `acc` - the accession number of the TF protein, e.g. ’P53762’\n",
"\n",
"- `data_type` - the type of data used to construct this motif, e.g.\n",
" ’SELEX’\n",
"\n",
"- `medline` - the Pubmed ID of literature supporting this motif, may\n",
" be multiple values, e.g. 7592839\n",
"\n",
"- `pazar_id` - external reference to the TF in the\n",
" [PAZAR](http://pazar.info) database, e.g. ’TF0000003’\n",
"\n",
"- `comment` - free form text containing notes about the construction\n",
" of the motif\n",
"\n",
"The `jaspar.Motif` class inherits from the generic `Motif` class and\n",
"therefore provides all the facilities of any of the motif formats —\n",
"reading motifs, writing motifs, scanning sequences for motif instances\n",
"etc.\n",
"\n",
"JASPAR stores motifs in several different ways including three different\n",
"flat file formats and as an SQL database. All of these formats\n",
"facilitate the construction of a counts matrix. However, the amount of\n",
"meta information described above that is available varies with the\n",
"format.\n",
"\n",
"#### The JASPAR `sites` format {#the-jaspar-sites-format .unnumbered}\n",
"\n",
"The first of the three flat file formats contains a list of instances.\n",
"As an example, these are the beginning and ending lines of the JASPAR\n",
"`Arnt.sites` file showing known binding sites of the mouse\n",
"helix-loop-helix transcription factor Arnt.\n",
"\n"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"```\n",
">MA0004 ARNT 1\n",
"CACGTGatgtcctc\n",
">MA0004 ARNT 2\n",
"CACGTGggaggtac\n",
">MA0004 ARNT 3\n",
"CACGTGccgcgcgc\n",
"...\n",
">MA0004 ARNT 18\n",
"AACGTGacagccctcc\n",
">MA0004 ARNT 19\n",
"AACGTGcacatcgtcc\n",
">MA0004 ARNT 20\n",
"aggaatCGCGTGc\n",
"\n",
"```"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"\n",
"The parts of the sequence in capital letters are the motif instances\n",
"that were found to align to each other.\n",
"\n",
"We can create a `Motif` object from these instances as follows:\n",
"\n"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": true
},
"outputs": [],
"source": [
"from Bio import motifs\n",
"with open(\"Arnt.sites\") as handle:\n",
" arnt = motifs.read(handle, \"sites\")"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"\n",
"The instances from which this motif was created is stored in the\n",
"`.instances` property:\n",
"\n"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": true
},
"outputs": [],
"source": [
"print(arnt.instances[:3])"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": true
},
"outputs": [],
"source": [
"for instance in arnt.instances:\n",
" print(instance)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"\n",
"The counts matrix of this motif is automatically calculated from the\n",
"instances:\n",
"\n"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": true
},
"outputs": [],
"source": [
"print(arnt.counts)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"\n",
"This format does not store any meta information.\n",
"\n",
"#### The JASPAR `pfm` format {#the-jaspar-pfm-format .unnumbered}\n",
"\n",
"JASPAR also makes motifs available directly as a count matrix, without\n",
"the instances from which it was created. This `pfm` format only stores\n",
"the counts matrix for a single motif. For example, this is the JASPAR\n",
"file `SRF.pfm` containing the counts matrix for the human SRF\n",
"transcription factor:\n",
"\n"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"```\n",
"2 9 0 1 32 3 46 1 43 15 2 2\n",
" 1 33 45 45 1 1 0 0 0 1 0 1\n",
"39 2 1 0 0 0 0 0 0 0 44 43\n",
" 4 2 0 0 13 42 0 45 3 30 0 0\n",
"\n",
"```"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"\n",
"We can create a motif for this count matrix as follows:\n",
"\n"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": true
},
"outputs": [],
"source": [
"with open(\"SRF.pfm\") as handle:\n",
"srf = motifs.read(handle, \"pfm\")"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": true
},
"outputs": [],
"source": [
"print(srf.counts)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"\n",
"As this motif was created from the counts matrix directly, it has no\n",
"instances associated with it:\n",
"\n"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": true
},
"outputs": [],
"source": [
"print(srf.instances)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"\n",
"We can now ask for the consensus sequence of these two motifs:\n",
"\n"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": true
},
"outputs": [],
"source": [
"print(arnt.counts.consensus)"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": true
},
"outputs": [],
"source": [
"print(srf.counts.consensus)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"\n",
"As with the instances file, no meta information is stored in this\n",
"format.\n",
"\n",
"#### The JASPAR format `jaspar` {#the-jaspar-format-jaspar .unnumbered}\n",
"\n",
"The `jaspar` file format allows multiple motifs to be specified in a\n",
"single file. In this format each of the motif records consist of a\n",
"header line followed by four lines defining the counts matrix. The\n",
"header line begins with a `>` character (similar to the Fasta file\n",
"format) and is followed by the unique JASPAR matrix ID and the TF name.\n",
"The following example shows a `jaspar` formatted file containing the\n",
"three motifs Arnt, RUNX1 and MEF2A:\n",
"\n"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"```\n",
">MA0004.1 Arnt\n",
"A [ 4 19 0 0 0 0 ]\n",
"C [16 0 20 0 0 0 ]\n",
"G [ 0 1 0 20 0 20 ]\n",
"T [ 0 0 0 0 20 0 ]\n",
">MA0002.1 RUNX1\n",
"A [10 12 4 1 2 2 0 0 0 8 13 ]\n",
"C [ 2 2 7 1 0 8 0 0 1 2 2 ]\n",
"G [ 3 1 1 0 23 0 26 26 0 0 4 ]\n",
"T [11 11 14 24 1 16 0 0 25 16 7 ]\n",
">MA0052.1 MEF2A\n",
"A [ 1 0 57 2 9 6 37 2 56 6 ]\n",
"C [50 0 1 1 0 0 0 0 0 0 ]\n",
"G [ 0 0 0 0 0 0 0 0 2 50 ]\n",
"T [ 7 58 0 55 49 52 21 56 0 2 ]\n",
"\n",
"```"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"\n",
"The motifs are read as follows:\n",
"\n"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": true
},
"outputs": [],
"source": [
"fh = open(\"jaspar_motifs.txt\")\n",
"for m in motifs.parse(fh, \"jaspar\"))\n",
" print(m)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"\n",
"Note that printing a JASPAR motif yields both the counts data and the\n",
"available meta-information.\n",
"\n",
"#### Accessing the JASPAR database {#accessing-the-jaspar-database .unnumbered}\n",
"\n",
"In addition to parsing these flat file formats, we can also retrieve\n",
"motifs from a JASPAR SQL database. Unlike the flat file formats, a\n",
"JASPAR database allows storing of all possible meta information defined\n",
"in the JASPAR `Motif` class. It is beyond the scope of this document to\n",
"describe how to set up a JASPAR database (please see the main\n",
"[JASPAR](http://jaspar.genereg.net) website). Motifs are read from a\n",
"JASPAR database using the `Bio.motifs.jaspar.db` module. First connect\n",
"to the JASPAR database using the JASPAR5 class which models the the\n",
"latest JASPAR schema:\n",
"\n"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": true
},
"outputs": [],
"source": [
"from Bio.motifs.jaspar.db import JASPAR5"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": true
},
"outputs": [],
"source": [
"JASPAR_DB_HOST = \n",
"JASPAR_DB_NAME = \n",
"JASPAR_DB_USER = \n",
"JASPAR_DB_PASS = \n"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": true
},
"outputs": [],
"source": [
"jdb = JASPAR5(\n",
"host=JASPAR_DB_HOST,\n",
"name=JASPAR_DB_NAME,\n",
"user=JASPAR_DB_USER,\n",
"password=JASPAR_DB_PASS\n",
")"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"\n",
"Now we can fetch a single motif by its unique JASPAR ID with the\n",
"`fetch_motif_by_id` method. Note that a JASPAR ID conists of a base ID\n",
"and a version number seperated by a decimal point, e.g. ’MA0004.1’. The\n",
"`fetch_motif_by_id` method allows you to use either the fully specified\n",
"ID or just the base ID. If only the base ID is provided, the latest\n",
"version of the motif is returned.\n",
"\n"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": true
},
"outputs": [],
"source": [
"arnt = jdb.fetch_motif_by_id(\"MA0004\")"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"\n",
"Printing the motif reveals that the JASPAR SQL database stores much more\n",
"meta-information than the flat files:\n",
"\n"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": true
},
"outputs": [],
"source": [
"print(arnt)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"\n",
"We can also fetch motifs by name. The name must be an exact match\n",
"(partial matches or database wildcards are not currently supported).\n",
"Note that as the name is not guaranteed to be unique, the\n",
"`fetch_motifs_by_name` method actually returns a list.\n",
"\n"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": true
},
"outputs": [],
"source": [
"motifs = jdb.fetch_motifs_by_name(\"Arnt\")\n",
"print(motifs[0])"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"\n",
"The `fetch_motifs` method allows you to fetch motifs which match a\n",
"specified set of criteria. These criteria include any of the above\n",
"described meta information as well as certain matrix properties such as\n",
"the minimum information content (`min_ic` in the example below), the\n",
"minimum length of the matrix or the minimum number of sites used to\n",
"construct the matrix. Only motifs which pass ALL the specified criteria\n",
"are returned. Note that selection criteria which correspond to meta\n",
"information which allow for multiple values may be specified as either a\n",
"single value or a list of values, e.g. `tax_group` and `tf_family` in\n",
"the example below.\n",
"\n"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": true
},
"outputs": [],
"source": [
"motifs = jdb.fetch_motifs(\n",
"collection = 'CORE',\n",
"tax_group = ['vertebrates', 'insects'],\n",
"tf_class = 'Winged Helix-Turn-Helix',\n",
"tf_family = ['Forkhead', 'Ets'],\n",
"min_ic = 12\n",
")\n",
"for motif in motifs:\n",
" pass # do something with the motif\n"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"\n",
"#### Compatibility with Perl TFBS modules {#compatibility-with-perl-tfbs-modules .unnumbered}\n",
"\n",
"An important thing to note is that the JASPAR `Motif` class was designed\n",
"to be compatible with the popular [Perl TFBS\n",
"modules](http://tfbs.genereg.net/). Therefore some specifics about the\n",
"choice of defaults for background and pseudocounts as well as how\n",
"information content is computed and sequences searched for instances is\n",
"based on this compatibility criteria. These choices are noted in the\n",
"specific subsections below.\n",
"\n",
"- **Choice of background:**\\\n",
" The Perl `TFBS` modules appear to allow a choice of custom\n",
" background probabilities (although the documentation states that\n",
" uniform background is assumed). However the default is to use a\n",
" uniform background. Therefore it is recommended that you use a\n",
" uniform background for computing the position-specific scoring\n",
" matrix (PSSM). This is the default when using the Biopython\n",
" `motifs` module.\n",
"\n",
"- **Choice of pseudocounts:**\\\n",
" By default, the Perl `TFBS` modules use a pseudocount equal to\n",
" $\\sqrt{N} * \\textrm{bg}[\\textrm{nucleotide}]$, where $N$ represents\n",
" the total number of sequences used to construct the matrix. To apply\n",
" this same pseudocount formula, set the motif `pseudocounts`\n",
" attribute using the `jaspar.calculate\\_pseudcounts()` function:\n",
"\n"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": true
},
"outputs": [],
"source": [
"motif.pseudocounts = motifs.jaspar.calculate_pseudocounts(motif)\n"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"\n",
" Note that it is possible for the counts matrix to have an unequal\n",
" number of sequences making up the columns. The pseudocount\n",
" computation uses the average number of sequences making up\n",
" the matrix. However, when `normalize` is called on the counts\n",
" matrix, each count value in a column is divided by the total number\n",
" of sequences making up that specific column, not by the average\n",
" number of sequences. This differs from the Perl `TFBS` modules\n",
" because the normalization is not done as a separate step and so the\n",
" average number of sequences is used throughout the computation of\n",
" the pssm. Therefore, for matrices with unequal column counts, the\n",
" PSSM computed by the `motifs` module will differ somewhat from the\n",
" pssm computed by the Perl `TFBS` modules.\n",
"\n",
"- **Computation of matrix information content:**\\\n",
" The information content (IC) or specificity of a matrix is computed\n",
" using the `mean` method of the\n",
" `PositionSpecificScoringMatrix` class. However of note, in the Perl\n",
" `TFBS` modules the default behaviour is to compute the IC without\n",
" first applying pseudocounts, even though by default the PSSMs are\n",
" computed using pseudocounts as described above.\n",
"\n",
"- **Searching for instances:**\\\n",
" Searching for instances with the Perl `TFBS` motifs was usually\n",
" performed using a relative score threshold, i.e. a score in the\n",
" range 0 to 1. In order to compute the absolute PSSM score\n",
" corresponding to a relative score one can use the equation:\n",
"\n"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": true
},
"outputs": [],
"source": [
"abs_score = (pssm.max - pssm.min) * rel_score + pssm.min\n"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"\n",
" To convert the absolute score of an instance back to a relative\n",
" score, one can use the equation:\n",
"\n"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": true
},
"outputs": [],
"source": [
"rel_score = (abs_score - pssm.min) / (pssm.max - pssm.min)\n"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"\n",
" For example, using the Arnt motif before, let’s search a sequence\n",
" with a relative score threshold of 0.8.\n",
"\n"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": true
},
"outputs": [],
"source": [
"test_seq=Seq(\"TAAGCGTGCACGCGCAACACGTGCATTA\", unambiguous_dna)\n",
"arnt.pseudocounts = motifs.jaspar.calculate_pseudocounts(arnt)\n",
"pssm = arnt.pssm\n",
"max_score = pssm.max\n",
"min_score = pssm.min\n",
"abs_score_threshold = (max_score - min_score) * 0.8 + min_score\n",
"for position, score in pssm.search(test_seq,\n"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": true
},
"outputs": [],
"source": [
"rel_score = (score - min_score) / (max_score - min_score)\n",
"print(\"Position %d: score = %5.3f, rel. score = %5.3f\" % (\n"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"\n",
"### MEME\n",
"\n",
"MEME @bailey1994 is a tool for discovering motifs in a group of related\n",
"DNA or protein sequences. It takes as input a group of DNA or protein\n",
"sequences and outputs as many motifs as requested. Therefore, in\n",
"contrast to JASPAR files, MEME output files typically contain multiple\n",
"motifs. This is an example.\n",
"\n",
"At the top of an output file generated by MEME shows some background\n",
"information about the MEME and the version of MEME used:\n",
"\n"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"```\n",
"********************************************************************************\n",
"MEME - Motif discovery tool\n",
"********************************************************************************\n",
"MEME version 3.0 (Release date: 2004/08/18 09:07:01)\n",
"...\n",
"\n",
"```"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"\n",
"Further down, the input set of training sequences is recapitulated:\n",
"\n"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"```\n",
"********************************************************************************\n",
"TRAINING SET\n",
"********************************************************************************\n",
"DATAFILE= INO_up800.s\n",
"ALPHABET= ACGT\n",
"Sequence name Weight Length Sequence name Weight Length\n",
"------------- ------ ------ ------------- ------ ------\n",
"CHO1 1.0000 800 CHO2 1.0000 800\n",
"FAS1 1.0000 800 FAS2 1.0000 800\n",
"ACC1 1.0000 800 INO1 1.0000 800\n",
"OPI3 1.0000 800\n",
"********************************************************************************\n",
"\n",
"```"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"\n",
"and the exact command line that was used:\n",
"\n"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"```\n",
"********************************************************************************\n",
"COMMAND LINE SUMMARY\n",
"********************************************************************************\n",
"This information can also be useful in the event you wish to report a\n",
"problem with the MEME software.\n",
"\n",
"command: meme -mod oops -dna -revcomp -nmotifs 2 -bfile yeast.nc.6.freq INO_up800.s\n",
"...\n",
"\n",
"```"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"\n",
"Next is detailed information on each motif that was found:\n",
"\n"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"```\n",
"********************************************************************************\n",
"MOTIF 1 width = 12 sites = 7 llr = 95 E-value = 2.0e-001\n",
"********************************************************************************\n",
"--------------------------------------------------------------------------------\n",
" Motif 1 Description\n",
"--------------------------------------------------------------------------------\n",
"Simplified A :::9:a::::3:\n",
"pos.-specific C ::a:9:11691a\n",
"probability G ::::1::94:4:\n",
"matrix T aa:1::9::11:\n",
"\n",
"```"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"\n",
"To parse this file (stored as `meme.dna.oops.txt`), use\n",
"\n"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": true
},
"outputs": [],
"source": [
"handle = open(\"meme.dna.oops.txt\")\n",
"record = motifs.parse(handle, \"meme\")\n",
"handle.close()\n"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"\n",
"The `motifs.parse` command reads the complete file directly, so you can\n",
"close the file after calling `motifs.parse`. The header information is\n",
"stored in attributes:\n",
"\n"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": true
},
"outputs": [],
"source": [
"record.version\n"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": true
},
"outputs": [],
"source": [
"record.datafile\n"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": true
},
"outputs": [],
"source": [
"record.command\n"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": true
},
"outputs": [],
"source": [
"record.alphabet\n"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": true
},
"outputs": [],
"source": [
"record.sequences\n"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"\n",
"The record is an object of the `Bio.motifs.meme.Record` class. The class\n",
"inherits from list, and you can think of `record` as a list of Motif\n",
"objects:\n",
"\n"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": true
},
"outputs": [],
"source": [
"len(record)\n"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": true
},
"outputs": [],
"source": [
"motif = record[0]\n",
"print(motif.consensus)\n"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": true
},
"outputs": [],
"source": [
"print(motif.degenerate_consensus)\n"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"\n",
"In addition to these generic motif attributes, each motif also stores\n",
"its specific information as calculated by MEME. For example,\n",
"\n"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": true
},
"outputs": [],
"source": [
"motif.num_occurrences\n"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": true
},
"outputs": [],
"source": [
"motif.length\n"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": true
},
"outputs": [],
"source": [
"evalue = motif.evalue\n",
"print(\"%3.1g\" % evalue)\n"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": true
},
"outputs": [],
"source": [
"motif.name\n"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"\n",
"In addition to using an index into the record, as we did above, you can\n",
"also find it by its name:\n",
"\n"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": true
},
"outputs": [],
"source": [
"motif = record['Motif 1']\n"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"\n",
"Each motif has an attribute `.instances` with the sequence instances in\n",
"which the motif was found, providing some information on each instance:\n",
"\n"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": true
},
"outputs": [],
"source": [
"len(motif.instances)\n"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": true
},
"outputs": [],
"source": [
"motif.instances[0]\n"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": true
},
"outputs": [],
"source": [
"motif.instances[0].motif_name\n"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": true
},
"outputs": [],
"source": [
"motif.instances[0].sequence_name\n"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": true
},
"outputs": [],
"source": [
"motif.instances[0].start\n"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": true
},
"outputs": [],
"source": [
"motif.instances[0].strand\n"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": true
},
"outputs": [],
"source": [
"motif.instances[0].length\n"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": true
},
"outputs": [],
"source": [
"pvalue = motif.instances[0].pvalue\n",
"print(\"%5.3g\" % pvalue)\n"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"\n",
"#### MAST {#mast .unnumbered}\n",
"\n",
"### TRANSFAC\n",
"\n",
"TRANSFAC is a manually curated database of transcription factors,\n",
"together with their genomic binding sites and DNA binding profiles\n",
"@matys2003. While the file format used in the TRANSFAC database is\n",
"nowadays also used by others, we will refer to it as the TRANSFAC file\n",
"format.\n",
"\n",
"A minimal file in the TRANSFAC format looks as follows:\n",
"\n"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"```\n",
"ID motif1\n",
"P0 A C G T\n",
"01 1 2 2 0 S\n",
"02 2 1 2 0 R\n",
"03 3 0 1 1 A\n",
"04 0 5 0 0 C\n",
"05 5 0 0 0 A\n",
"06 0 0 4 1 G\n",
"07 0 1 4 0 G\n",
"08 0 0 0 5 T\n",
"09 0 0 5 0 G\n",
"10 0 1 2 2 K\n",
"11 0 2 0 3 Y\n",
"12 1 0 3 1 G\n",
"//\n",
"\n",
"```"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"\n",
"This file shows the frequency matrix of motif `motif1` of 12\n",
"nucleotides. In general, one file in the TRANSFAC format can contain\n",
"multiple motifs. For example, this is the contents of the example\n",
"TRANSFAC file `transfac.dat`:\n",
"\n"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"```\n",
"VV EXAMPLE January 15, 2013\n",
"XX\n",
"//\n",
"ID motif1\n",
"P0 A C G T\n",
"01 1 2 2 0 S\n",
"02 2 1 2 0 R\n",
"03 3 0 1 1 A\n",
"...\n",
"11 0 2 0 3 Y\n",
"12 1 0 3 1 G\n",
"//\n",
"ID motif2\n",
"P0 A C G T\n",
"01 2 1 2 0 R\n",
"02 1 2 2 0 S\n",
"...\n",
"09 0 0 0 5 T\n",
"10 0 2 0 3 Y\n",
"//\n",
"\n",
"```"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"\n",
"To parse a TRANSFAC file, use\n",
"\n"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": true
},
"outputs": [],
"source": [
"handle = open(\"transfac.dat\")\n",
"record = motifs.parse(handle, \"TRANSFAC\")\n",
"handle.close()\n"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"\n",
"The overall version number, if available, is stored as `record.version`:\n",
"\n"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": true
},
"outputs": [],
"source": [
"record.version\n"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"\n",
"Each motif in `record` is in instance of the `Bio.motifs.transfac.Motif`\n",
"class, which inherits both from the `Bio.motifs.Motif` class and from a\n",
"Python dictionary. The dictionary uses the two-letter keys to store any\n",
"additional information about the motif:\n",
"\n"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": true
},
"outputs": [],
"source": [
"motif = record[0]\n",
"motif.degenerate_consensus # Using the Bio.motifs.Motif method\n"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": true
},
"outputs": [],
"source": [
"motif['ID'] # Using motif as a dictionary\n"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"\n",
"TRANSFAC files are typically much more elaborate than this example,\n",
"containing lots of additional information about the motif. Table\n",
"\\[table:transfaccodes\\] lists the two-letter field codes that are\n",
"commonly found in TRANSFAC files:\n",
"\n",
"\\[table:transfaccodes\\]\n",
"\n",
" ------ -------------------------------------------------\n",
" `AC` Accession number\n",
" `AS` Accession numbers, secondary\n",
" `BA` Statistical basis\n",
" `BF` Binding factors\n",
" `BS` Factor binding sites underlying the matrix\n",
" `CC` Comments\n",
" `CO` Copyright notice\n",
" `DE` Short factor description\n",
" `DR` External databases\n",
" `DT` Date created/updated\n",
" `HC` Subfamilies\n",
" `HP` Superfamilies\n",
" `ID` Identifier\n",
" `NA` Name of the binding factor\n",
" `OC` Taxonomic classification\n",
" `OS` Species/Taxon\n",
" `OV` Older version\n",
" `PV` Preferred version\n",
" `TY` Type\n",
" `XX` Empty line; these are not stored in the Record.\n",
" ------ -------------------------------------------------\n",
"\n",
" : Fields commonly found in TRANSFAC files\n",
"\n",
"Each motif also has an attribute `.references` containing the references\n",
"associated with the motif, using these two-letter keys:\n",
"\n",
" ------ -------------------\n",
" `RN` Reference number\n",
" `RA` Reference authors\n",
" `RL` Reference data\n",
" `RT` Reference title\n",
" `RX` PubMed ID\n",
" ------ -------------------\n",
"\n",
" : Fields used to store references in TRANSFAC files\n",
"\n",
"Printing the motifs writes them out in their native TRANSFAC format:\n",
"\n"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": true
},
"outputs": [],
"source": [
"print(record)\n"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"\n",
"You can export the motifs in the TRANSFAC format by capturing this\n",
"output in a string and saving it in a file:\n",
"\n"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": true
},
"outputs": [],
"source": [
"text = str(record)\n",
"handle = open(\"mytransfacfile.dat\", 'w')\n",
"handle.write(text)\n",
"handle.close()\n"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"\n",
"Writing motifs\n",
"--------------\n",
"\n",
"Speaking of exporting, let’s look at export functions in general. We can\n",
"use the `format` method to write the motif in the simple JASPAR `pfm`\n",
"format:\n",
"\n"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": true
},
"outputs": [],
"source": [
"print(arnt.format(\"pfm\"))\n"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"\n",
"Similarly, we can use `format` to write the motif in the JASPAR `jaspar`\n",
"format:\n",
"\n"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": true
},
"outputs": [],
"source": [
"print(arnt.format(\"jaspar\"))\n"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"\n",
"To write the motif in a TRANSFAC-like matrix format, use\n",
"\n"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": true
},
"outputs": [],
"source": [
"print(m.format(\"transfac\"))\n"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"\n",
"To write out multiple motifs, you can use `motifs.write`. This function\n",
"can be used regardless of whether the motifs originated from a TRANSFAC\n",
"file. For example,\n",
"\n"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": true
},
"outputs": [],
"source": [
"two_motifs = [arnt, srf]\n",
"print(motifs.write(two_motifs, 'transfac'))\n"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"\n",
"Or, to write multiple motifs in the `jaspar` format:\n",
"\n"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": true
},
"outputs": [],
"source": [
"two_motifs = [arnt, mef2a]\n",
"print(motifs.write(two_motifs, \"jaspar\"))\n"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"\n",
"Position-Weight Matrices\n",
"------------------------\n",
"\n",
"The `.counts` attribute of a Motif object shows how often each\n",
"nucleotide appeared at each position along the alignment. We can\n",
"normalize this matrix by dividing by the number of instances in the\n",
"alignment, resulting in the probability of each nucleotide at each\n",
"position along the alignment. We refer to these probabilities as the\n",
"position-weight matrix. However, beware that in the literature this term\n",
"may also be used to refer to the position-specific scoring matrix, which\n",
"we discuss below.\n",
"\n",
"Usually, pseudocounts are added to each position before normalizing.\n",
"This avoids overfitting of the position-weight matrix to the limited\n",
"number of motif instances in the alignment, and can also prevent\n",
"probabilities from becoming zero. To add a fixed pseudocount to all\n",
"nucleotides at all positions, specify a number for the `pseudocounts`\n",
"argument:\n",
"\n"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": true
},
"outputs": [],
"source": [
"pwm = m.counts.normalize(pseudocounts=0.5)\n",
"print(pwm)\n"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"\n",
"Alternatively, `pseudocounts` can be a dictionary specifying the\n",
"pseudocounts for each nucleotide. For example, as the GC content of the\n",
"human genome is about 40%, you may want to choose the pseudocounts\n",
"accordingly:\n",
"\n"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": true
},
"outputs": [],
"source": [
"pwm = m.counts.normalize(pseudocounts={'A':0.6, 'C': 0.4, 'G': 0.4, 'T': 0.6})\n",
"print(pwm)\n"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"\n",
"The position-weight matrix has its own methods to calculate the\n",
"consensus, anticonsensus, and degenerate consensus sequences:\n",
"\n"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": true
},
"outputs": [],
"source": [
"pwm.consensus\n"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": true
},
"outputs": [],
"source": [
"pwm.anticonsensus\n"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": true
},
"outputs": [],
"source": [
"pwm.degenerate_consensus\n"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"\n",
"Note that due to the pseudocounts, the degenerate consensus sequence\n",
"calculated from the position-weight matrix is slightly different from\n",
"the degenerate consensus sequence calculated from the instances in the\n",
"motif:\n",
"\n"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": true
},
"outputs": [],
"source": [
"m.degenerate_consensus\n"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"\n",
"The reverse complement of the position-weight matrix can be calculated\n",
"directly from the `pwm`:\n",
"\n"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": true
},
"outputs": [],
"source": [
"rpwm = pwm.reverse_complement()\n",
"print(rpwm)\n"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"\n",
"Position-Specific Scoring Matrices\n",
"----------------------------------\n",
"\n",
"Using the background distribution and PWM with pseudo-counts added, it’s\n",
"easy to compute the log-odds ratios, telling us what are the log odds of\n",
"a particular symbol to be coming from a motif against the background. We\n",
"can use the `.log_odds()` method on the position-weight matrix:\n",
"\n"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": true
},
"outputs": [],
"source": [
"pssm = pwm.log_odds()\n",
"print(pssm)\n"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"\n",
"Here we can see positive values for symbols more frequent in the motif\n",
"than in the background and negative for symbols more frequent in the\n",
"background. $0.0$ means that it’s equally likely to see a symbol in the\n",
"background and in the motif.\n",
"\n",
"This assumes that A, C, G, and T are equally likely in the background.\n",
"To calculate the position-specific scoring matrix against a background\n",
"with unequal probabilities for A, C, G, T, use the `background`\n",
"argument. For example, against a background with a 40% GC content, use\n",
"\n"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": true
},
"outputs": [],
"source": [
"background = {'A':0.3,'C':0.2,'G':0.2,'T':0.3}\n",
"pssm = pwm.log_odds(background)\n",
"print(pssm)\n"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"\n",
"The maximum and minimum score obtainable from the PSSM are stored in the\n",
"`.max` and `.min` properties:\n",
"\n"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": true
},
"outputs": [],
"source": [
"print(\"%4.2f\" % pssm.max)\n"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": true
},
"outputs": [],
"source": [
"print(\"%4.2f\" % pssm.min)\n"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"\n",
"The mean and standard deviation of the PSSM scores with respect to a\n",
"specific background are calculated by the `.mean` and `.std` methods.\n",
"\n"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": true
},
"outputs": [],
"source": [
"mean = pssm.mean(background)\n",
"std = pssm.std(background)\n",
"print(\"mean = %0.2f, standard deviation = %0.2f\" % (mean, std))\n"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"\n",
"A uniform background is used if `background` is not specified. The mean\n",
"is particularly important, as its value is equal to the Kullback-Leibler\n",
"divergence or relative entropy, and is a measure for the information\n",
"content of the motif compared to the background. As in Biopython the\n",
"base-2 logarithm is used in the calculation of the log-odds scores, the\n",
"information content has units of bits.\n",
"\n",
"The `.reverse_complement`, `.consensus`, `.anticonsensus`, and\n",
"`.degenerate_consensus` methods can be applied directly to PSSM objects.\n",
"\n",
"Searching for instances {#sec:search}\n",
"-----------------------\n",
"\n",
"The most frequent use for a motif is to find its instances in some\n",
"sequence. For the sake of this section, we will use an artificial\n",
"sequence like this:\n",
"\n"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": true
},
"outputs": [],
"source": [
"test_seq=Seq(\"TACACTGCATTACAACCCAAGCATTA\", m.alphabet)\n",
"len(test_seq)\n"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"\n",
"### Searching for exact matches\n",
"\n",
"The simplest way to find instances, is to look for exact matches of the\n",
"true instances of the motif:\n",
"\n"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": true
},
"outputs": [],
"source": [
"for pos, seq in m.instances.search(test_seq):\n",
"print(\"%i %s\" % (pos, seq))\n"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"\n",
"We can do the same with the reverse complement (to find instances on the\n",
"complementary strand):\n",
"\n"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": true
},
"outputs": [],
"source": [
"for pos, seq in r.instances.search(test_seq):\n",
"print(\"%i %s\" % (pos, seq))\n"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"\n",
"### Searching for matches using the PSSM score\n",
"\n",
"It’s just as easy to look for positions, giving rise to high log-odds\n",
"scores against our motif:\n",
"\n"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": true
},
"outputs": [],
"source": [
"for position, score in pssm.search(test_seq, threshold=3.0):\n",
"print(\"Position %d: score = %5.3f\" % (position, score))\n"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"\n",
"The negative positions refer to instances of the motif found on the\n",
"reverse strand of the test sequence, and follow the Python convention on\n",
"negative indices. Therefore, the instance of the motif at `pos` is\n",
"located at `test_seq[pos:pos+len(m)]` both for positive and for negative\n",
"values of `pos`.\n",
"\n",
"You may notice the threshold parameter, here set arbitrarily to $3.0$.\n",
"This is in $log_2$, so we are now looking only for words, which are\n",
"eight times more likely to occur under the motif model than in the\n",
"background. The default threshold is $0.0$, which selects everything\n",
"that looks more like the motif than the background.\n",
"\n",
"You can also calculate the scores at all positions along the sequence:\n",
"\n"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": true
},
"outputs": [],
"source": [
"pssm.calculate(test_seq)\n"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"\n",
"In general, this is the fastest way to calculate PSSM scores. The scores\n",
"returned by `pssm.calculate` are for the forward strand only. To obtain\n",
"the scores on the reverse strand, you can take the reverse complement of\n",
"the PSSM:\n",
"\n"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": true
},
"outputs": [],
"source": [
"rpssm = pssm.reverse_complement()\n",
"rpssm.calculate(test_seq)\n"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"\n",
"### Selecting a score threshold\n",
"\n",
"If you want to use a less arbitrary way of selecting thresholds, you can\n",
"explore the distribution of PSSM scores. Since the space for a score\n",
"distribution grows exponentially with motif length, we are using an\n",
"approximation with a given precision to keep computation cost\n",
"manageable:\n",
"\n"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": true
},
"outputs": [],
"source": [
"distribution = pssm.distribution(background=background, precision=10**4)\n"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"\n",
"The `distribution` object can be used to determine a number of different\n",
"thresholds. We can specify the requested false-positive rate\n",
"(probability of “finding” a motif instance in background generated\n",
"sequence):\n",
"\n"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": true
},
"outputs": [],
"source": [
"threshold = distribution.threshold_fpr(0.01)\n",
"print(\"%5.3f\" % threshold)\n"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"\n",
"or the false-negative rate (probability of “not finding” an instance\n",
"generated from the motif):\n",
"\n"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": true
},
"outputs": [],
"source": [
"threshold = distribution.threshold_fnr(0.1)\n",
"print(\"%5.3f\" % threshold)\n"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"\n",
"or a threshold (approximately) satisfying some relation between the\n",
"false-positive rate and the false-negative rate\n",
"($\\frac{\\textrm{fnr}}{\\textrm{fpr}}\\simeq t$):\n",
"\n"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": true
},
"outputs": [],
"source": [
"threshold = distribution.threshold_balanced(1000)\n",
"print(\"%5.3f\" % threshold)\n"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"\n",
"or a threshold satisfying (roughly) the equality between the\n",
"false-positive rate and the $-log$ of the information content (as used\n",
"in patser software by Hertz and Stormo):\n",
"\n"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": true
},
"outputs": [],
"source": [
"threshold = distribution.threshold_patser()\n",
"print(\"%5.3f\" % threshold)\n"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"\n",
"For example, in case of our motif, you can get the threshold giving you\n",
"exactly the same results (for this sequence) as searching for instances\n",
"with balanced threshold with rate of $1000$.\n",
"\n"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": true
},
"outputs": [],
"source": [
"threshold = distribution.threshold_fpr(0.01)\n",
"print(\"%5.3f\" % threshold)\n"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": true
},
"outputs": [],
"source": [
"for position, score in pssm.search(test_seq, threshold=threshold):\n",
"print(\"Position %d: score = %5.3f\" % (position, score))\n"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"\n",
"Each motif object has an associated Position-Specific Scoring Matrix\n",
"--------------------------------------------------------------------\n",
"\n",
"To facilitate searching for potential TFBSs using PSSMs, both the\n",
"position-weight matrix and the position-specific scoring matrix are\n",
"associated with each motif. Using the Arnt motif as an example:\n",
"\n"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": true
},
"outputs": [],
"source": [
"from Bio import motifs\n",
"with open(\"Arnt.sites\") as handle:\n",
"motif = motifs.read(handle, 'sites')\n"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": true
},
"outputs": [],
"source": [
"print(motif.counts)\n"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": true
},
"outputs": [],
"source": [
"print(motif.pwm)\n"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"\n"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": true
},
"outputs": [],
"source": [
"print(motif.pssm)\n"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"\n",
"The negative infinities appear here because the corresponding entry in\n",
"the frequency matrix is 0, and we are using zero pseudocounts by\n",
"default:\n",
"\n"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": true
},
"outputs": [],
"source": [
"for letter in \"ACGT\":\n",
"print(\"%s: %4.2f\" % (letter, motif.pseudocounts[letter]))\n"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"\n",
"If you change the `.pseudocounts` attribute, the position-frequency\n",
"matrix and the position-specific scoring matrix are recalculated\n",
"automatically:\n",
"\n"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": true
},
"outputs": [],
"source": [
"motif.pseudocounts = 3.0\n",
"for letter in \"ACGT\":\n",
"print(\"%s: %4.2f\" % (letter, motif.pseudocounts[letter]))\n"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"\n"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": true
},
"outputs": [],
"source": [
"print(motif.pwm)\n"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"\n"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": true
},
"outputs": [],
"source": [
"print(motif.pssm)\n"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"\n",
"You can also set the `.pseudocounts` to a dictionary over the four\n",
"nucleotides if you want to use different pseudocounts for them. Setting\n",
"`motif.pseudocounts` to `None` resets it to its default value of zero.\n",
"\n",
"The position-specific scoring matrix depends on the background\n",
"distribution, which is uniform by default:\n",
"\n"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": true
},
"outputs": [],
"source": [
"for letter in \"ACGT\":\n",
"print(\"%s: %4.2f\" % (letter, motif.background[letter]))\n"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"\n",
"Again, if you modify the background distribution, the position-specific\n",
"scoring matrix is recalculated:\n",
"\n"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": true
},
"outputs": [],
"source": [
"motif.background = {'A': 0.2, 'C': 0.3, 'G': 0.3, 'T': 0.2}\n",
"print(motif.pssm)\n"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"\n",
"Setting `motif.background` to `None` resets it to a uniform\n",
"distribution:\n",
"\n"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": true
},
"outputs": [],
"source": [
"motif.background = None\n",
"for letter in \"ACGT\":\n",
"print(\"%s: %4.2f\" % (letter, motif.background[letter]))\n"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"\n",
"If you set `motif.background` equal to a single value, it will be\n",
"interpreted as the GC content:\n",
"\n"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": true
},
"outputs": [],
"source": [
"motif.background = 0.8\n",
"for letter in \"ACGT\":\n",
"print(\"%s: %4.2f\" % (letter, motif.background[letter]))\n"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"\n",
"Note that you can now calculate the mean of the PSSM scores over the\n",
"background against which it was computed:\n",
"\n"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": true
},
"outputs": [],
"source": [
"print(\"%f\" % motif.pssm.mean(motif.background))\n"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"\n",
"as well as its standard deviation:\n",
"\n"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": true
},
"outputs": [],
"source": [
"print(\"%f\" % motif.pssm.std(motif.background))\n"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"\n",
"and its distribution:\n",
"\n"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": true
},
"outputs": [],
"source": [
"distribution = motif.pssm.distribution(background=motif.background)\n",
"threshold = distribution.threshold_fpr(0.01)\n",
"print(\"%f\" % threshold)\n"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"\n",
"Note that the position-weight matrix and the position-specific scoring\n",
"matrix are recalculated each time you call `motif.pwm` or `motif.pssm`,\n",
"respectively. If speed is an issue and you want to use the PWM or PSSM\n",
"repeatedly, you can save them as a variable, as in\n",
"\n"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": true
},
"outputs": [],
"source": [
"pssm = motif.pssm\n"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"\n",
"Comparing motifs {#sec:comp}\n",
"----------------\n",
"\n",
"Once we have more than one motif, we might want to compare them.\n",
"\n",
"Before we start comparing motifs, I should point out that motif\n",
"boundaries are usually quite arbitrary. This means we often need to\n",
"compare motifs of different lengths, so comparison needs to involve some\n",
"kind of alignment. This means we have to take into account two things:\n",
"\n",
"- alignment of motifs\n",
"\n",
"- some function to compare aligned motifs\n",
"\n",
"To align the motifs, we use ungapped alignment of PSSMs and substitute\n",
"zeros for any missing columns at the beginning and end of the matrices.\n",
"This means that effectively we are using the background distribution for\n",
"columns missing from the PSSM. The distance function then returns the\n",
"minimal distance between motifs, as well as the corresponding offset in\n",
"their alignment.\n",
"\n",
"To give an example, let us first load another motif, which is similar to\n",
"our test motif `m`:\n",
"\n"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": true
},
"outputs": [],
"source": [
"with open(\"REB1.pfm\") as handle:\n",
"m_reb1 = motifs.read(handle, \"pfm\")\n"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": true
},
"outputs": [],
"source": [
"m_reb1.consensus\n"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": true
},
"outputs": [],
"source": [
"print(m_reb1.counts)\n"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"\n",
"To make the motifs comparable, we choose the same values for the\n",
"pseudocounts and the background distribution as our motif `m`:\n",
"\n"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": true
},
"outputs": [],
"source": [
"m_reb1.pseudocounts = {'A':0.6, 'C': 0.4, 'G': 0.4, 'T': 0.6}\n",
"m_reb1.background = {'A':0.3,'C':0.2,'G':0.2,'T':0.3}\n",
"pssm_reb1 = m_reb1.pssm\n",
"print(pssm_reb1)\n"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"\n",
"We’ll compare these motifs using the Pearson correlation. Since we want\n",
"it to resemble a distance measure, we actually take $1-r$, where $r$ is\n",
"the Pearson correlation coefficient (PCC):\n",
"\n"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": true
},
"outputs": [],
"source": [
"distance, offset = pssm.dist_pearson(pssm_reb1)\n",
"print(\"distance = %5.3g\" % distance)\n"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": true
},
"outputs": [],
"source": [
"print(offset)\n"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"\n",
"This means that the best PCC between motif `m` and `m_reb1` is obtained\n",
"with the following alignment:\n",
"\n"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"```\n",
"m: bbTACGCbb\n",
"m_reb1: GTTACCCGG\n",
"\n",
"```"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"\n",
"where `b` stands for background distribution. The PCC itself is roughly\n",
"$1-0.239=0.761$.\n",
"\n",
"*De novo* motif finding {#sec:find}\n",
"-----------------------\n",
"\n",
"Currently, Biopython has only limited support for *de novo* motif\n",
"finding. Namely, we support running and parsing of AlignAce and MEME.\n",
"Since the number of motif finding tools is growing rapidly,\n",
"contributions of new parsers are welcome.\n",
"\n",
"### MEME {#sec:meme}\n",
"\n",
"Let’s assume, you have run MEME on sequences of your choice with your\n",
"favorite parameters and saved the output in the file `meme.out`. You can\n",
"retrieve the motifs reported by MEME by running the following piece of\n",
"code:\n",
"\n"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": true
},
"outputs": [],
"source": [
"from Bio import motifs\n",
"with open(\"meme.out\") as handle:\n",
"motifsM = motifs.parse(handle, \"meme\")\n"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"\n"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": true
},
"outputs": [],
"source": [
"motifsM\n"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"\n",
"Besides the most wanted list of motifs, the result object contains more\n",
"useful information, accessible through properties with self-explanatory\n",
"names:\n",
"\n",
"- `.alphabet`\n",
"\n",
"- `.datafile`\n",
"\n",
"- `.sequence_names`\n",
"\n",
"- `.version`\n",
"\n",
"- `.command`\n",
"\n",
"The motifs returned by the MEME Parser can be treated exactly like\n",
"regular Motif objects (with instances), they also provide some extra\n",
"functionality, by adding additional information about the instances.\n",
"\n"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": true
},
"outputs": [],
"source": [
"motifsM[0].consensus\n"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": true
},
"outputs": [],
"source": [
"motifsM[0].instances[0].sequence_name\n"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": true
},
"outputs": [],
"source": [
"motifsM[0].instances[0].start\n"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": true
},
"outputs": [],
"source": [
"motifsM[0].instances[0].strand\n"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"\n"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": true
},
"outputs": [],
"source": [
"motifsM[0].instances[0].pvalue\n"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"\n",
"### AlignAce {#sec:alignace}\n",
"\n",
"We can do very similar things with the AlignACE program. Assume, you\n",
"have your output in the file `alignace.out`. You can parse your output\n",
"with the following code:\n",
"\n"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": true
},
"outputs": [],
"source": [
"from Bio import motifs\n",
"with open(\"alignace.out\") as handle:\n",
"motifsA = motifs.parse(handle, \"alignace\")\n"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"\n",
"Again, your motifs behave as they should:\n",
"\n"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": true
},
"outputs": [],
"source": [
"motifsA[0].consensus\n"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"\n",
"In fact you can even see, that AlignAce found a very similar motif as\n",
"MEME. It is just a longer version of a reverse complement of the MEME\n",
"motif:\n",
"\n"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": true
},
"outputs": [],
"source": [
"motifsM[0].reverse_complement().consensus\n"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"\n",
"If you have AlignAce installed on the same machine, you can also run it\n",
"directly from Biopython. A short example of how this can be done is\n",
"shown below (other parameters can be specified as keyword parameters):\n",
"\n"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": true
},
"outputs": [],
"source": [
"command=\"/opt/bin/AlignACE\"\n",
"input_file=\"test.fa\"\n",
"from Bio.motifs.applications import AlignAceCommandline\n",
"cmd = AlignAceCommandline(cmd=command, input=input_file, gcback=0.6, numcols=10)\n",
"stdout, stderr= cmd()\n"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"\n",
"Since AlignAce prints all of its output to standard output, you can get\n",
"to your motifs by parsing the first part of the result:\n",
"\n"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": true
},
"outputs": [],
"source": [
"motifs = motifs.parse(stdout, \"alignace\")\n"
]
}
],
"metadata": {
"kernelspec": {
"display_name": "Python 3",
"language": "python",
"name": "python3"
},
"language_info": {
"codemirror_mode": {
"name": "ipython",
"version": 3
},
"file_extension": ".py",
"mimetype": "text/x-python",
"name": "python",
"nbconvert_exporter": "python",
"pygments_lexer": "ipython3",
"version": "3.5.1"
}
},
"nbformat": 4,
"nbformat_minor": 0
}