most eukaryotic transcription factors ( tfs ) are members of protein families that share a common deoxyribonucleic acid ( dna ) binding domain and have highly similar dna binding preferences .
however , individual tf family members ( i.e. paralogous tfs ) often have different functions and bind to different genomic regions in vivo , as observed from chromatin immunoprecipitation assays followed by microarray analysis or high - throughput sequencing ( chip - chip or chip - seq ) ( 1,2 ) . despite the large amount of in vivo chip - seq data currently available , especially through the encode project ( 3 ) , computational tools for analyzing differences between the genomic binding profiles of paralogous tfs are still lacking .
first , some pairs of paralogous tfs exhibit subtle differences in dna binding specificity  either for the core binding site ( 4 ) or for the binding site flanks ( 1)and such differences can explain , at least in part , how each tf selects its unique targets .
second , paralogous tfs may interact with different protein co - factors that modulate their dna binding specificity ( 5 ) , or they may respond differently to certain chromatin environments .
third , some paralogous tfs are expressed in different cells or at different stages during cellular differentiation or during the cell cycle ; in such cases , the precise chromatin environment in the cell where each paralogous tf is expressed will dictate where the tf binds in the genome . here , we focus on paralogous tfs that are present in the cell at the same time , have highly similar dna binding specificities , but still show significant differences in their in vivo genomic binding profiles , as measured by chip - seq . for such paralogous tfs , interactions with different sets of protein co -
we present an extensive web implementation of our recently published algorithm couger ( co - factors associated with uniquely - bound genomic regions ) ( 6 ) , a classification - based framework for identifying protein co - factors that might provide specificity to paralogous tfs .
couger can be applied to any two sets of genomic regions bound by paralogous tfs ( e.g. regions derived from chip - seq experiments ) .
the framework uses state - of - the - art classification algorithms ( support vector machines and random forest ) with features that reflect the dna - binding specificities of putative co - factors .
a custom feature selection procedure is used to obtain a small subset of non - redundant putative co - factors that are most important for distinguishing between genomic regions bound by the considered pair of paralogous tfs .
the identified co - factors are presented in a user - friendly output page , together with information about the importance of each co - factor feature , and the classification accuracy .
users can run couger through a web interface : http://couger.oit.duke.edu , or as a stand - alone python software tool ( available for download on the couger website ) .
couger uses support vector machine ( svm ) ( 7 ) and random forest ( rf ) ( 8) , two state - of - the - art classification algorithms with free software packages : libsvm ( 9 ) and random jungle ( 10 ) .
both algorithms are highly accurate , can successfully handle high - dimensional data and are robust on data with highly correlated features .
svm is a non - probabilistic binary linear classifier with great performance on both linear and nonlinear classification problems .
rf is an ensemble of multiple classification trees , which explicitly computes a measure of the importance of each variable for the classification task .
we trained svms with both linear and radial basis function kernels ( \documentclass[12pt]{minimal } 
 \usepackage{amsmath } 
 \usepackage{wasysym } 
 \usepackage{amsfonts } 
 \usepackage{amssymb } 
 \usepackage{amsbsy } 
 \usepackage{upgreek } 
 \usepackage{mathrsfs } 
 \setlength{\oddsidemargin}{-69pt } 
 \begin{document } 
 } { } $ \mathrm { svm } _
{ \mathrm lin } $ \end{document } and \documentclass[12pt]{minimal } 
 \usepackage{amsmath } 
 \usepackage{wasysym } 
 \usepackage{amsfonts } 
 \usepackage{amssymb } 
 \usepackage{amsbsy } 
 \usepackage{upgreek } 
 \usepackage{mathrsfs } 
 \setlength{\oddsidemargin}{-69pt } 
 \begin{document } 
 } { } $ \mathrm { svm}_\mathrm { \mathrm { rbf } } $ \end{document } , respectively ) ( 9 ) , and rf with the unscaled permutation importance ( \documentclass[12pt]{minimal } 
 \usepackage{amsmath } 
 \usepackage{wasysym } 
 \usepackage{amsfonts } 
 \usepackage{amssymb } 
 \usepackage{amsbsy } 
 \usepackage{upgreek } 
 \usepackage{mathrsfs } 
 \setlength{\oddsidemargin}{-69pt } 
 \begin{document } 
 } { } $ \mathrm { rf}_\mathrm{\mathrm { pi}}$\end{document } ) .
the latter measure represents the average decrease in classification accuracy when the values of the respective variable are randomly permuted ( 10 ) .
we use different classifiers in order to assess the reliability of the results and their independence of particular techniques .
in addition , each method has specific strengths and weaknesses ( \documentclass[12pt]{minimal } 
 \usepackage{amsmath } 
 \usepackage{wasysym } 
 \usepackage{amsfonts } 
 \usepackage{amssymb } 
 \usepackage{amsbsy } 
 \usepackage{upgreek } 
 \usepackage{mathrsfs } 
 \setlength{\oddsidemargin}{-69pt } 
 \begin{document } 
 } { } $ \mathrm { svm}_\mathrm { \mathrm { rbf } } $ \end{document } usually yields better performance than \documentclass[12pt]{minimal } 
 \usepackage{amsmath } 
 \usepackage{wasysym } 
 \usepackage{amsfonts } 
 \usepackage{amssymb } 
 \usepackage{amsbsy } 
 \usepackage{upgreek } 
 \usepackage{mathrsfs } 
 \setlength{\oddsidemargin}{-69pt } 
 \begin{document } 
 } { } $ \mathrm { { \rm svm}}_\mathrm { \mathrm { { \rm lin } } } $ \end{document } , while results obtained with \documentclass[12pt]{minimal } 
 \usepackage{amsmath } 
 \usepackage{wasysym } 
 \usepackage{amsfonts } 
 \usepackage{amssymb } 
 \usepackage{amsbsy } 
 \usepackage{upgreek } 
 \usepackage{mathrsfs } 
 \setlength{\oddsidemargin}{-69pt } 
 \begin{document } 
 } { } $ \mathrm { { \rm svm}}_\mathrm { \mathrm { { \rm lin } } } $ \end{document } are more interpretable ) .
the two classes are the dna sequences under the chip - seq peaks for two paralogous tfs ( tf1 and tf2 ) , which the user can specify either in fasta format or , for convenience , directly with chip - seq peak coordinates in encode narrowpeak format , encode broadpeak format , or even bed format .
tf1- and tf2-specific sequences are defined by excluding the chip - seq peaks that overlap any peak of the other tf ( figure 1 , step 1 ) . in order to avoid a potential classification bias toward one of the two classes , an equal number of dna sequences from each set is selected .
then , in the case of narrowpeak input files , which we strongly recommend , each sequence is trimmed to 100 bp on each side of the chip - seq peak summit .
if possible , couger considers only the close vicinity of the tf1 and tf2 chip - seq peak summits because our goal is to identify co - factors that bind together with tf1 or tf2 .
we note that for high - quality chip - seq data , the tf - dna binding events are thought to occur , in general , within 50 bp of the peak summit . the couger framework .
step 1 represents the derivation of the two classes ( tf1- and tf2-specific ) and is omitted in the case of fasta input files . in step 2 , all features are computed for the two classes , from pbm or pwm data .
the circles represent the number of features that are considered in each iteration ; the darker circles correspond to the sets of features that are used in classification . in step 4
, classifiers are learned on the training set and then predictions are made on the  test set .
steps 3 and 4 are repeated five times , according to the 5-fold cv setting .
features reflecting the binding specificity of putative co - factors are computed from : ( i ) high - throughput in vitro tf - dna binding data from universal protein - binding microarray ( pbm ) assays ( 11,12 ) , or ( ii ) from large collections of dna binding motifs ( i.e. position weight matrices , pwms ) ( 13,14 ) ( figure 1 , step 2 ) . from each universal pbm data set we use the enrichment scores ( e - scores ) for all possible 8-mers .
mann whitney statistic and ranges from 0.5 ( least favored sequence ) to + 0.5 ( most favored sequence ) . for a given pbm data set or pwm , and a given dna sequence , couger uses two features :
 max ( the maximum score over all the k - mers in that sequence ) and  top3avg ( the average score over the top 3 highest - scoring k - mers in that sequence  this takes into account the fact that tf binding sites may occur in clusters ) ; see ( 6 ) and supplementary figure s1 for more details on computing the features .
one of the most important steps of our classification approach is feature selection , because we expect only a small number of tfs to be potential co - factors and interact with the considered paralogous tfs .
prior to classification , the input set of sequences is randomly divided into a training set ( consisting of 80% of the data ) and a test set ( consisting of the remaining 20% of the data ) .
only sequences in the training set are used during feature selection , to ensure the complete independence of the test set , which is used only to evaluate the performance of our classification models .
couger performs feature selection using a combined procedure consisting of rf recursive feature elimination ( rf - fs ) and minimum redundancy maximum relevance feature selection ( mrmr ) ( 15 ) .
this procedure is illustrated in figure 1 , step 3 , which depicts the case of using 858 features derived from pbm data for 429 tfs .
after the feature selection step , classification is performed on five feature sets : ( i ) all features , ( ii ) under 100 features selected by rf - fs , ( iii ) under 10 features selected by rf - fs , ( iv ) the first five features selected by mrmr and ( v ) the first 10 features selected by mrmr .
rf - fs is an iterative process in which a random forest is grown at each step and a subset of variables are discarded .
we recursively eliminate half of the features ( with the smallest importance ) until < 100 or < 10 features remain .
then we apply a normalized variant of the mrmr algorithm , called nmifs ( normalized mutual information feature selection ) ( 16 ) , to the set of features from the rf - fs iteration in which their number was below 100 .
the nmifs technique ranks the features by considering both the relevance / importance for distinguishing between classes , and the redundancy between pairs of features .
the method is based on mutual information ( mi ) and works best with discrete - valued features ( 17 ) .
for this reason , we derived a discretization approach that computes the 20-quantiles ( separately for  max features and  top3avg features ) and maps the features to 20 integer values .
this procedure also allows us to avoid the mi 's bias toward features with larger sets of possible values .
we note that we tested several discretization approaches , as well as several feature selection algorithms ( 15,16,1820 ) , and we found that nmifs with discretized features achieved the highest classification accuracy ( data not shown ) .
we evaluate the performance of our classification models using a 5-fold cross - validation ( cv ) approach . in each of the five runs
, couger performs three procedures on the training data : feature selection , grid search over the parameter space and training of the classifier(s ) .
for each set of features and each type of classifier , couger computes median values for accuracy , sensitivity , specificity and precision , which are reported in separate files .
in addition , it shows the classification accuracy ( the fraction of true positives and true negatives in the test set ) in a user - friendly format ( figure 2a ) .
the results correspond to tfs c - fos ( henceforth referred to as fos ) and jund , with pbm - derived features .
( a ) median classification accuracies for fos versus jund ( left ) , and options for interactive selection ( right ) . the accuracies are presented in a heatmap manner , where green corresponds to the minimum value and red to the maximum .
each column represents a selected feature from  rf - fs under 100 ( i.e. random forest feature selection run to select < 100 features ) .
( c ) a set of selected features ( sorted by their score ) , together with their correlation ( i.e. the pearson correlation coefficient ) with the class label . the first class , in this example fos , is considered class  0. the second class , in this example jund , is considered class  1.
thus , a negative correlation for a particular feature suggests that the feature is important for tf1 , while a negative correlation suggests that the feature is important for tf2 .
the name of each feature contains the name of pbm file used to generate that feature , as well as  max or 
top3avg , which specifies whether the feature represents the score of the best site in each sequence or an average over the top three sites , respectively .
couger web runs under apache2 ( httpd.apache.org ) with mod_wsgi on debian 7.0  wheezy ( www.debian.org ) .
it was developed using the django web framework ( www.djangoproject.com ) , which is written in python ( www.python.org ) and thus allowed direct integration with the couger source code .
also , part of the web functions are implemented with javascript and jquery ( jquery.com ) .
couger takes as input two sets of genomic regions bound by ( paralogous ) tfs .
users can submit these files either in encode narrowpeak format , encode broadpeak format , bed format or in fasta format .
if the user provides fasta files , couger will consider that the classes are predetermined and will start with step 2 ( figure 1 ) , so only the nucleotide sequences will be read ( we note that in the absence of genomic coordinates , it is much harder and time - consuming to identify the overlapping sequences between the two sets ) . for all other input formats , couger will start with step 1 ( figure 1 ) , removing from each set the sequences that have an overlap with any sequence from the other set .
the resulting sets of sequences are reported by couger in fasta format , and can be used to rerun the algorithms using different settings , such as different feature sets .
moreover , if narrowpeak files are uploaded , couger will focus the search for putative co - factor binding sites by trimming the peaks to 100 bases centered at the peak summit , which will reduce the running time and improve the results .
thus , narrowpeak is the recommended format for the input sets of sequence . in the case of a bed - like format , the user may choose the reference genome among five different versions of the human genome , four versions of the mouse genome and three versions of the fly genome .
we note that the user can specify other genomes or genome versions in the downloadable version of our framework .
the web version of couger also enforces a restriction of maximum n sequences per class , due to time and resources constraints ( n is a threshold set by the user , and can vary between 300 and 1000 ) .
if the number of unique targets exceeds the threshold , then couger will run using only the top n sequences .
couger also requires a set of features that reflect the dna binding specificities of putative co - factors .
the current version of the web server offers six choices of such features for human and mouse data :  pbm data from uniprobe ( i.e. pbm 8-mer e - scores for 429 mammalian tfs from the uniprobe database ( 12 ) ) ,  pwms from uniprobe ( i.e. pwms derived from pbm data in uniprobe ( 12 ) ) ,  pwms from transfac ( i.e. 1226 pwms from the transfac database ( 13 ) ) ,  pwms from ht - selex data ( i.e. 239 pwms for human or mouse tfs , derived from ht - selex data of jolma et al .
( 14 ) ) ,  pwms from jaspar core vertebrata ( i.e. 205 pwms from the jaspar database ( 21 ) ) and  pwms from uniprobe & ht - selex & jaspar core vertebrata ( all 876 pwms from the three databases ) .
couger also offers three choices for drosophila melanogaster data : pwms from transfac ( i.e. 1226 pwms from transfac database ( 13 ) ) ,  pwms from jaspar core insecta ( i.e. 131 pwms from jaspar database ( 21 ) ) and  pwms from transfac and jaspar core insecta ( all 1357 pwms from the two databases ) . after submission and validation , the user is redirected to a status page , where job details and the running log are provided .
when the job is completed , the results replace the status page , and an e - mail is sent to the user if an email address was provided .
the results can be viewed online , or can be downloaded as a zip file ( which is recommended , because the results may be deleted after 48 h ) .
the median classification accuracies ( before and after feature selection ) are displayed in a heatmap - like color coded table with values for each type of algorithm and each set of features ( figure 2 a ) .
the user can interact with the results page to view more detailed information : the variation in classification accuracy over the cv runs ( the information is displayed in three boxplots , one for each classifier : svm - lin , svm - rbf and rf - pi ) ; the accuracies for individual cv runs ; a heatmap showing all the features selected by rf - fs and their values for the sequences in each of the two classes ( figure 2 b ) ; the ranking of all the features in main feature sets ( rf - fs ( under 100 ) , rf - fs ( under 10 ) and nmifs ) ( figure 2 b and c ) .
we tested couger on 20 pairs of paralogous tfs ( table 1 ) with chip - seq data from encode ( 3 ) in the k562 cell line , processed using a uniform pipeline .
briefly , we applied the idr framework ( 22 ) together with the macs ( 23 ) peak caller ( version 2.0.10 ) , for which the size of the shift was previously estimated by spp ( 24 ) ( see supporting materials and methods ) .
next , we analyzed the results and filtered out the tfs with idr scores that did not follow the restrictions recommended by encode ( 25 ) ( see supporting material for details )
. using the peak calling pipeline , we identified high - quality chip - seq data ( in terms of data reproducibility ) for 20 pairs of paralogous tf ( table 1 ) . 

table 1.pairs of paralogous tf with high - quality chip - seq data for the k562 cell line in encodetf1tf2 # seqstf1tf2 # seqsnfyanfyb70fosjunb1842cjunjund90elf1ets11962maxmxi1382elf1gabpa1976maffmafk506fosjund2032cmycmxi1538junbjund2194cmycmax*586cmycmax**2598sp1sp2736ets1gabpa3516gata1gata2802cfoscjun5186e2f4e2f61002cebpbcebpd7000atf1atf31682bhlhe40tal17152the pairs of factors are sorted by the  # seqs column , which contains the number of sequences selected by couger for both classes ( tf1- and tf2-specific ) .
the pair marked with * corresponds to a chip - seq data set using an igg control .
the pair marked with * * corresponds to a chip - seq data set using a standard control .
pairs of paralogous tf with high - quality chip - seq data for the k562 cell line in encode the pairs of factors are sorted by the  # seqs column , which contains the number of sequences selected by couger for both classes ( tf1- and tf2-specific ) .
the pair marked with * corresponds to a chip - seq data set using an igg control .
the pair marked with * * corresponds to a chip - seq data set using a standard control .
although we report here results for data processed with the idr pipeline , which is the current standard for the encode project ( 25 ) , couger has no requirements or limitations regarding the peak finding methodology .
we recommend , however , a uniform and careful processing of the input data because it is important that the chip data for both tf1 and tf2 is comparable and of high quality .
large differences in data quality between tfs will likely result in one of the two paralogous tfs being the dominant signal when comparing the unique peaks of the two factors .
we ran couger on 20 pairs of tfs ( table 1 ) , using features derived either from pbm data from uniprobe , or the joint set of pwms from uniprobe , ht - selex and jaspar core vertebrata .
after feature selection , the classification accuracy varies between 71.9 and 99.5% depending on the pair of paralogous tfs , with no correlation between classification performance and the number of sequences in the training set ( detailed results are presented in supplementary tables s1 and s2 ) .
couger performed well with both pbm - derived and pwm - derived features ( figure 4 , accuracy and precision , and supplementary figure s2 , sensitivity and specificity ) and returned similar sets of putative co - factors , which is not surprising given the high overlap between the tfs represented in the two data sets .
in general , we recommend running couger with both options , as some tfs have data in only one of the two types of feature sets .
couger classification performance for 20 pairs of tfs , with features derived from pwm data ( from uniprobe , ht - selex and jaspar core vertebrata ) .
the values correspond to the best result from all three classifiers ( \documentclass[12pt]{minimal } 
 \usepackage{amsmath } 
 \usepackage{wasysym } 
 \usepackage{amsfonts } 
 \usepackage{amssymb } 
 \usepackage{amsbsy } 
 \usepackage{upgreek } 
 \usepackage{mathrsfs } 
 \setlength{\oddsidemargin}{-69pt } 
 \begin{document } 
 } { } $ \mathrm { svm}_\mathrm { \mathrm { lin } } $ \end{document}\documentclass[12pt]{minimal } 
 \usepackage{amsmath } 
 \usepackage{wasysym } 
 \usepackage{amsfonts } 
 \usepackage{amssymb } 
 \usepackage{amsbsy } 
 \usepackage{upgreek } 
 \usepackage{mathrsfs } 
 \setlength{\oddsidemargin}{-69pt } 
 \begin{document } 
 } { } $ \mathrm { svm}_\mathrm { \mathrm { lin } } $ \end{document } , \documentclass[12pt]{minimal } 
 \usepackage{amsmath } 
 \usepackage{wasysym } 
 \usepackage{amsfonts } 
 \usepackage{amssymb } 
 \usepackage{amsbsy } 
 \usepackage{upgreek } 
 \usepackage{mathrsfs } 
 \setlength{\oddsidemargin}{-69pt } 
 \begin{document } 
 } { } $ \mathrm { svm}_\mathrm { \mathrm { rbf } } $ \end{document}\documentclass[12pt]{minimal } 
 \usepackage{amsmath } 
 \usepackage{wasysym } 
 \usepackage{amsfonts } 
 \usepackage{amssymb } 
 \usepackage{amsbsy } 
 \usepackage{upgreek } 
 \usepackage{mathrsfs } 
 \setlength{\oddsidemargin}{-69pt } 
 \begin{document } 
 } { } $ \mathrm { svm}_\mathrm { \mathrm { rbf } } $ \end{document } and \documentclass[12pt]{minimal } 
 \usepackage{amsmath } 
 \usepackage{wasysym } 
 \usepackage{amsfonts } 
 \usepackage{amssymb } 
 \usepackage{amsbsy } 
 \usepackage{upgreek } 
 \usepackage{mathrsfs } 
 \setlength{\oddsidemargin}{-69pt } 
 \begin{document } 
 } { } $ \mathrm { rf}_\mathrm{\mathrm { pi}}$\end{document}\documentclass[12pt]{minimal } 
 \usepackage{amsmath } 
 \usepackage{wasysym } 
 \usepackage{amsfonts } 
 \usepackage{amssymb } 
 \usepackage{amsbsy } 
 \usepackage{upgreek } 
 \usepackage{mathrsfs } 
 \setlength{\oddsidemargin}{-69pt } 
 \begin{document } 
 } { } $ \mathrm { rf}_\mathrm{\mathrm { pi}}$\end{document } ) and all four sets of features derived by the fs procedure ( under 100 and under 10 features selected by rf - fs , and first 5 and 10 features selected by nmifs ) .
couger classification performance ( accuracy and precision ) for 20 pairs of tfs , with pbm and pwm features .
the values correspond to the best result from all three classifiers ( \documentclass[12pt]{minimal } 
 \usepackage{amsmath } 
 \usepackage{wasysym } 
 \usepackage{amsfonts } 
 \usepackage{amssymb } 
 \usepackage{amsbsy } 
 \usepackage{upgreek } 
 \usepackage{mathrsfs } 
 \setlength{\oddsidemargin}{-69pt } 
 \begin{document } 
 } { } $ \mathrm { svm}_\mathrm { \mathrm { lin } } $ \end{document}\documentclass[12pt]{minimal } 
 \usepackage{amsmath } 
 \usepackage{wasysym } 
 \usepackage{amsfonts } 
 \usepackage{amssymb } 
 \usepackage{amsbsy } 
 \usepackage{upgreek } 
 \usepackage{mathrsfs } 
 \setlength{\oddsidemargin}{-69pt } 
 \begin{document } 
 } { } $ \mathrm { svm}_\mathrm { \mathrm { lin } } $ \end{document } , \documentclass[12pt]{minimal } 
 \usepackage{amsmath } 
 \usepackage{wasysym } 
 \usepackage{amsfonts } 
 \usepackage{amssymb } 
 \usepackage{amsbsy } 
 \usepackage{upgreek } 
 \usepackage{mathrsfs } 
 \setlength{\oddsidemargin}{-69pt } 
 \begin{document } 
 } { } $ \mathrm { svm}_\mathrm { \mathrm { rbf } } $ \end{document}\documentclass[12pt]{minimal } 
 \usepackage{amsmath } 
 \usepackage{wasysym } 
 \usepackage{amsfonts } 
 \usepackage{amssymb } 
 \usepackage{amsbsy } 
 \usepackage{upgreek } 
 \usepackage{mathrsfs } 
 \setlength{\oddsidemargin}{-69pt } 
 \begin{document } 
 } { } $ \mathrm { svm}_\mathrm { \mathrm { rbf } } $ \end{document } and \documentclass[12pt]{minimal } 
 \usepackage{amsmath } 
 \usepackage{wasysym } 
 \usepackage{amsfonts } 
 \usepackage{amssymb } 
 \usepackage{amsbsy } 
 \usepackage{upgreek } 
 \usepackage{mathrsfs } 
 \setlength{\oddsidemargin}{-69pt } 
 \begin{document } 
 } { } $ \mathrm { rf}_\mathrm{\mathrm { pi}}$\end{document}\documentclass[12pt]{minimal } 
 \usepackage{amsmath } 
 \usepackage{wasysym } 
 \usepackage{amsfonts } 
 \usepackage{amssymb } 
 \usepackage{amsbsy } 
 \usepackage{upgreek } 
 \usepackage{mathrsfs } 
 \setlength{\oddsidemargin}{-69pt } 
 \begin{document } 
 } { } $ \mathrm { rf}_\mathrm{\mathrm { pi}}$\end{document } ) and all four sets of features derived by the fs procedure .
the horizontal lines represent the median value over all pairs of tfs , for pbm or pwm features .
we note that for c - myc and mxi1 , we obtained a classification accuracy of up to 93.5% ( supplementary table s1 ) , compared to our previously reported accuracy of 88.4% ( 6 ) .
the increased performance is due to improvements in the couger framework , mainly in the feature selection steps .
to determine whether the accuracy on the tested tf pairs is significant , we randomized the classes for all pairs of tfs and ran couger on the randomized sets ( see supporting materials and methods ) . as expected , the classification accuracy for randomized data varied between 44.44 and 53.94% , with a median of 49.75% ( supplementary tables s3 and s4 ) , which demonstrates that the high accuracies on real tf binding data are not due to chance .
we note that the median variation of randomized accuracies for five different class - label shuffling is 3.7% ( supplementary table s5 ) .
importantly , the couger output page is highly interactive and allows users to visualize details of the classification results , as well as details regarding the co - factor features that enabled a successful classification .
this is in contrast to many implementations of classification algorithms , which provide an accuracy measure but no indication of what features drive the classification . through its user - friendly and interactive design ,
the couger web tool makes it easy for users to understand and explore the contribution of individual features .
importantly , for several ( tf1 , tf2 ) pairs we have found evidence in the literature supporting our hypothesis that factors identified by couger interact with tf1 or tf2 and thus may contribute to their in vivo dna binding differences .
for example , in the case of fos and jund with pbm features ( figure 2 ) , we found support for several putative co - factors . both fos and jund are basic leucine zipper ( bzip ) proteins from the ap-1 subfamily .
fos binds dna as a heterodimer with c - jun and other members of the ap-1 subfamily , while jund can homodimerize and , interestingly , it can interact with tfs from the atf / creb subfamily , another branch of the bzip family . as shown in figure 2
, couger found that features reflecting the specificity of the atf / creb subfamily ( such as  atf1_3026_contig8mers_v1_max ) are strongly associated with jund - unique sequences , consistent with a direct interaction between jund and atf / creb factors .
gata factors , also found to be important for distinguishing fos- from jund - unique targets , have been previously reported to interact with proteins from the ap-1 subfamily ( 26,27 ) ; the exact identity of the gata co - factor that might interact specifically with jund remains to be determined .
tf sfpi1 , associated with jund - unique targets , is known to interact directly with jund ( 28 ) .
however , it has very similar sequence preferences to tgif2 , a factor that has been shown to interact specifically with jund ( 29 ) . in the fos - unique sequences ,
the most important features reflect the general specificity of fos : jun complexes , which might indicate that at those targets fos binds together with ap-1 proteins other than jund . in a comparison of c
we found e2f factors associated with c - fos targets , consistent with the previously reported roles of c - fos and e2fs in the same signaling cascade that links ras activity to cyclin a transcription ( 30 ) . in the same comparison , tf mitf was associated with c - jun , consistent with their direct interaction and synergistic effects on gene regulation ( 31,32 ) . in a comparison of atf1 versus atf3 ,
myc / max / mad tfs were found enriched in atf3-unique targets ; a potential interaction between atf3 and max has been reported previously in the encode project and is under further investigation . among the tested tf pairs ,
this could indicate that one tf is present at a higher concentration , binds with higher affinity overall , or has higher quality data .
importantly , couger allows the user to remove these tfs from the set of putative co - factors ( using an option in the input page ) and re - run the classification framework .
we note that even after eliminating tf1 and tf2 from the set of features , the classification accuracy remains very high .
to illustrate this , in supplementary figures s3 and s4 we show our results for tfs c - myc versus mxi1 , before and after removing features reflecting the specificity of these two factors .
importantly , we see that the predictions accuracy , as well as the selected putative co - factors ( namely , rfx proteins ) remained the same ( supplementary figures s3 and s4 ) .
our peak calling preprocessing pipeline allowed us to determine high confidence peaks for each tf ( which allowed us to compare paralogous tfs ) , as well as self - consistent peaks for each replicate ( which allowed us to compare replicate experiments performed for the same tf ) . therefore , as a control , we also ran couger on the 31 pairs of replicate experiments .
the classification accuracy for pairs of replicates ranged between 52.3 and 94.2% , with a median of 76.9% ( supplementary table s6 and supplementary figure s5 ) .
this was far from the expected random classification result , which was obtained for the randomized classes of pairs of paralogous tfs .
randomizing the data for the replicate pairs results in an expected accuracy level ( 47.0553.84% , with a median of 50% ; supplementary table s7 ) . by investigating further the behavior for pairs of replicates , we found that the putative co - factors selected by couger are different between replicate experiments versus paralogous tfs ( see supplementary table s8 ) .
the classification between replicate experiments is driven mostly by gata factors and by tfs that bind gc repeat regions ( such as zfp161 and e2f ) .
indeed , the nucleotide frequencies , in particular the gc content , play a major role in distinguishing between sequences unique to only one replicate .
supplementary figure s6 shows the difference in gc% for all 20 pairs of paralogous tfs and 31 pairs of tf replicates .
this difference correlates very well with the classification accuracy in the case of replicate experiments ( pearson correlation coefficient 0.706 ) , but not in the case of tf pairs ( pearson correlation coefficient 0.217 ) ( supplementary tables s9 and s10 ) .
the high classification accuracy for replicate data sets could be due to experimental bias : recent studies have found strong biases in chip - seq data for regions close to the tss of highly expressed genes ( 33,34 ) , which are oftentimes enriched in cg dinucleotides .
( indeed , many of the peaks unique to only one of the chip - seq replicates are close to tsss . ) given that control experiments and replicate experiments reported in encode were not performed at the same time , it is currently not possible to correct for this bias in the chip - seq data .
we note , however , that for the tf pairs in our analysis , even when the tfs shows relatively large differences in gc - content , the identified putative co - factors were not tfs that bound cg - repeats or high gc - content sites , which suggests that the potential bias driving the classification between replicates is not influencing couger 's ability to identify co - factors . in cases where there is uncertainty regarding the quality of the replicate experiments , we recommend running couger for tf1tf1 and tf2tf2 , in addition to tf1tf2 , and comparing the sets of putative co - factors obtained for paralogous tfs versus replicates .
the goal of the couger framework is to help users generate hypotheses regarding potential co - factors that might provide in vivo specificity to paralogous tfs .
these putative co - factors are selected from the set of tfs with known binding specificities provided as features .
couger can use either pbm scores or pwm motifs for extensive sets of proteins , both types of data reflecting protein - dna interactions .
we compared the use of pbm scores from uniprobe with that of pwms from uniprobe , and although the results were similar , we note that each approach has its own advantages and disadvantages .
using the pwms derived from pbm data has the advantage that one can easily combine these pwms with additional pwms from other databases ( such as jaspar and/or ht - selex ) , thus expanding the number of tested putative co - factors .
however , pwms represent only a summary of the pbm data , and they make the assumption that individual positions within tf binding sites contribute independently to the binding affinity .
this assumption is not always true , so in some cases the pwms may not accurately reflect tf binding specificity
. the pbm 8-mer data does not suffer from this drawback , because it simply represents the binding preferences of a tf for all possible 8-mers .
but its disadvantage is that the number of eukaryotic factors with pbm data available is relatively small compared to the number of tfs with known pwms , although we expect this situation to change as more and more pbm data sets are being generated . after a set of putative co
- factors are identified , follow - up studies should be performed to test : ( i ) whether the putative co - factors ( or factors with very similar specificities ) are expressed in the cell type of interest ; ( ii ) whether the identified putative co - factors are bound in vivo to tf1- or tf2-unique sequences ; and ( iii ) whether the putative co - factors interact physically with tf1 and tf2 .
such follow - up studies are necessary to test whether identified co - factors contribute to the regulatory specificity of the paralogous tfs of interest .
first , the current version of our web server allows users to run the classification with either pbm data from uniprobe ( 12 ) , or pwms from transfac ( 13 ) , ht - selex data ( 14 ) , jaspar core ( 21 ) and/or uniprobe ( 12 ) .
we will include an option for the users to upload a custom set of dna motifs or pbm datasets .
future work also includes adding features derived from chip - seq data for putative co - factors , extending the web server to other organisms , giving users the option to run only one of the three classification algorithms included in our tool , and adding more interactive options for feature analysis ( such as direct links from the co - factor features to databases with more information on each co - factor ) .
this work was supported by the phrma foundation through a research starter grant ( to r.g . ) .
 most transcription factors ( tfs ) belong to protein families that share a common dna binding domain and have very similar dna binding preferences 
 . however , many paralogous tfs ( i.e. members of the same tf family ) perform different regulatory functions and interact with different genomic regions in the cell . 
 a potential mechanism for achieving this differential in vivo specificity is through interactions with protein co - factors . 
 computational tools for studying the genomic binding profiles of paralogous tfs and identifying their putative co - factors are currently lacking . here , we present an interactive web implementation of couger , a classification - based framework for identifying protein co - factors that might provide specificity to paralogous tfs . 
 couger takes as input two sets of genomic regions bound by paralogous tfs , and it identifies a small set of putative co - factors that best distinguish the two sets of sequences . to achieve this task 
 , couger uses a classification approach , with features that reflect the dna - binding specificities of the putative co - factors . 
 the identified co - factors are presented in a user - friendly output page , together with information that allows the user to understand and to explore the contributions of individual co - factor features . 
 couger can be run as a stand - alone tool or through a web interface : http://couger.oit.duke.edu . 
INTRODUCTION
MATERIALS AND METHODS
RESULTS
DISCUSSION
SUPPLEMENTARY DATA
FUNDING