BACKGROUND
gene expression data from blood genomic studies are widely used for investigation of human disease  <cit> .
predictive gene signatures have been developed to carry out differential diagnosis of infectious diseases  <cit> , identify specific disease states  <cit>  and characterize the immune response to vaccination  <cit> .
however, is some cases gene expression signatures from blood can be weakly expressed and highly variable  <cit> .
the identification of these signatures is complicated by the fact that blood is a mixed tissue, composed of multiple cell subsets, so that differential expression profiles can reflect changes in cell subset proportions, changes in subset-specific gene expression or both.
in cases where the relevant disease genes are expressed in a subset-specific manner, as has been shown for sle  <cit> , analyses based on mixed cell expression data are inherently limited since differential expression of genes in one cell subset  will be diluted by rna from other cells.
experimental studies that isolate specific cell subsets before expression profiling can provide important biological insight by demonstrating subset-specific gene expression as well as increased predictive signal  <cit> .
given the large number of potential cell subsets that can be defined, the ability to identify the most informative subset to isolate would be a great aid to these studies.
most genome-wide expression studies are based on analysis of total peripheral blood mononuclear cells .
pbmcs are composed of over a dozen cell subsets that are derived from a common progenitor in the bone marrow .
these cells are commonly divided into myeloid and lymphoid cells.
myeloid cells include monocytes and their descendants, as well as granulocytes like neutrophils and basophils.
lymphoid cells are primarily composed of b cells, t cells and nk cells.
proportions of these cells can vary widely between individuals, but t cells and b cells together usually make up ~75% of pbmcs, while nk cells and monocytes make up around 10-15% each.
the remaining cell types, such as dendritic cells, are much more rare, and account for <1% of total pbmcs  <cit> .
neutrophils, which normally compose the majority of cells in a blood sample , are normally excluded by the methods used to isolate pbmcs, but may account for up to 20% of a pbmc sample due to contamination  <cit> .
genome-wide expression measurements based on total pbmcs reflect both condition-specific gene expression as well as the proportion of different cell subsets in the sample.
microarray deconvolution methods have been developed to take advantage of this latter dependence in order to quantify the relative proportion of different cell subsets  <cit> .
in this approach, the expression level of each gene is modeled as a linear function of the expression from each cell subset comprising the sample.
deconvolution thus depends on prior knowledge of quantitative expression levels for each subset.
in cases where the subset proportions have been measured experimentally, methods have been developed to compute subset-specific differential-expression for individual genes  <cit> .
however, in addition to requiring measurements that are not available for many existing datasets, these methods do not directly address the problem of predicting the cell subset source of a pre-defined gene signature .
in some cases, clues to the cellular source of a gene expression signature may be found in the signature itself.
this would be true when the signature is enriched for genes known to be expressed in a particular cell type  <cit> , and may occur when the proportion of cells in the blood changes in the condition of interest.
for example, berry et al.
found that changes in t-cell related genes included in a transcriptional signature of active tuberculosis resulted from changes in t cell number in the blood  <cit> .
in cases where the predictive signature does not include a significant number of cell subset-specific genes, the problem may be addressed experimentally through the isolation of specific cell subsets followed by transcriptional profiling  <cit> .
however, it is often not obvious which cell subset will be the most informative.
the time and expense of isolating and analyzing all possibly relevant cell subsets provides the motivation for the method developed here: a computational approach to predict the most likely cell subset source of a gene expression signature.
the development of this method was also motivated by the clinical problem of predicting the outcome of antiviral therapy in patients chronically infected with hepatitis c virus .
in this case, gene expression changes that are predictive when measured in hepatocytes are significantly blunted in pbmcs, making it difficult to develop highly predictive models from this cell population  <cit> .
the isolation of specific cell subsets prior to expression profiling may allow for the identification of predictive signatures based on a blood sample, with significant advantages for the patient.
we have developed a computational method to predict the most likely cellular source for a pre-defined gene expression signature using transcriptional profiling data from total pbmcs.
the general approach we propose is to:  estimate the relative proportion of each pbmc subset from individual patients using genome-wide transcriptional profiling data, and  correlate these proportions with the gene expression signature that predicts therapy.
we theorize that the pbmc subset proportion that correlates most closely with the predicted gene expression signature over a large set of patients is the most likely source of the signal.
as a specific example, let's consider a study where we find that interferon stimulated genes  are up-regulated in a subset of patients that fail to respond to a particular therapy.
now, we focus on the non-responding patients.
if isg expression tends to be higher in patients where the proportion of dendritic cells  is predicted to be high, then we would predict that dcs are the most likely source of the observed gene expression signature.
note that there can be two interpretations of this result with respect to the original study.
first, the dcs may always express isgs, but the proportion of dcs in the blood is higher for non-responders.
second, the dcs undergo specific gene expression changes in non-responders which cause the level of isg expression to increase in the blood sample.
in either case, dcs are the relevant population to follow-up experimentally.
here, we develop this general idea into a computational method called spec , which we then validate using both healthy and disease data.
while the development of spec was motivated by the problem of predicting the response to therapy in hcv patients, the method is generally applicable to link gene expression signatures with specific pbmc subsets.
since only gene expression profiles from total pbmcs are required, the technique can easily be used to mine the massive amount of existing microarray data.
<dig>
RESULTS
subset prediction from enrichment correlation  is a method that takes advantage of the fluctuations of cell proportions in a mixed cell sample to find the most likely source for a gene expression signal.
the general approach can be broken down into four steps:
<dig>  estimate the relative proportion of each pbmc subset  across a population.
we propose to accomplish this by calculating the enrichment score  associated with a set of genes expressed specifically in each subset s, called a subset signature .
this is summarized by the following pseudo-code:
code:
n_p = number of individuals
n_s = number of pbmc subsets
exprdata = gene expression profile for individual i
subsetsignature = list of genes representing pbmc subset s
querysignature = a list of genes input by the user
for each
for each
e = enrichment of subsetsignature in exprdata
<dig>  calculate the enrichment score of the query signature for each individual i in the population:
for each
q = enrichment of querysignature in exprdata
<dig>  determine the pearson correlation  between the enrichment scores of the query signature and each pbmc subset s:
for each
c = correlation between q and e across i
<dig>  predict the subset with the highest correlation to the query signature:
predicted subset = s that gives rise to the max
the calculation of enrichment scores is explained in methods and figure  <dig>
<dig>  estimating relative pbmc subset proportions using enrichment scores
the first step in our proposed method is to estimate the relative proportion of each pbmc subset.
while a number of methods exist that can estimate these values, we chose to use enrichment scores calculated from cell subset-specific gene signatures for our initial implementation of spec due to its relative simplicity.
as described in the methods, the enrichment score calculates the extent to which a gene signature is concentrated among the most highly expressed genes in a single genome-wide expression experiment.
the cell subset-specific signatures are defined as groups of genes which have been found to be up-regulated in only a specific subset of cells in a mixed cell population.
the resulting lists of up-regulated genes represent cell markers for that blood cell type, and include well known cell surface markers such as cd <dig> for t cells and cd <dig> for b cells.
lists of genes describing specific pbmc subsets have previously been defined in  <cit>  and in  <cit> .
here, we choose to use signatures representing the b cell, t cell, nk cell, and monocyte subsets, we also included a signature for neutrophils despite the fact that they are not considered pbmcs, since neutrophil contamination can account for up to 20% of a pbmc sample  <cit> .
in addition, we also use a more general set of signatures for the lymphoid and myeloid subsets .
the lymphoid and myeloid signatures are composed of genes which are expressed in more than one cell type from those lineages, thus these subsets should vary according to the sum of all the cells in that lineage.
we run spec separately on the detailed  and general  signatures, as described below.
we hypothesize that the enrichment scores of the subset-specific gene signatures will be related to the relative cell subset proportions across patients.
although microarray deconvolution can be used to accomplish this task  <cit> , these methods require quantitative gene expression levels for each subset to be estimated, which is not always available.
in addition, it is important to be clear that we do not require the enrichment scores of different cell subsets to predict their proportions in a single individual .
rather, we propose that the enrichment score for a single cell subset will vary along with the relative size of the cell subset across patients.
for example, if subject b has a higher fraction of nk cells than subject a then the nk signature enrichment score should also be higher, but both of these enrichment scores may be higher than the t cell signature enrichment scores even if nk cells constitute a lower fraction of cells in both subjects.
in order to validate this enrichment score approach, we compared it with predictions from microarray deconvolution  <cit>  using data from healthy subjects where we expect deconvolution should perform relatively well.
we predicted the relative proportions of cell subsets for  <dig> healthy gene expression profiles using both microarray deconvolution and the cell signature enrichment scores.
figure  <dig> is an example showing the correlation between the nk cell signature enrichment score and the predicted nk cell fraction estimated using deconvolution.
the correlation of  <dig>  implies that the enrichment score is positively varying with the relative proportion of nk cells in these healthy pbmc samples.
we found similarly high correlations for b cells, neutrophils and monocytes .
in contrast, t cells  had a much weaker correlation, perhaps reflecting the quality of this particular gene signature.
overall, these results suggest that enrichment scores can be used to track the fluctuations in the levels of each cell type in a mixed cell expression profile .
<dig>  validation of spec using split subset signatures
as an initial validation of spec, we sought to predict the source of a gene expression signature where the correct cell subset is known.
this can be done by using our pre-defined cell subset expression signatures.
for example, consider the b cell signature defined in  <cit> .
the genes that are part of this signature can be split into two disjoint sets.
one of these sets can be labeled as the cell subset signature and used to estimate the changes in the b cell population, while the other can be retained as the query gene signature being tested.
a successful method would be able to identify this query signature as coming from b cells or from lymphoid cells.
when running spec, we separately analyze the set of detailed signatures  and general signatures , which represent cell types that are a superset of those covered by detailed signatures .
for each test, one of the subset signatures was randomly divided into two halves, with one half used as the query signature and the other half retained as the subset signature.
this was repeated  <dig> times for each of the pbmc subset signatures.
in order to determine the significance of the correlation value for each run we used a monte carlo permutation test .
the following two subsections describe the performance of spec on healthy and disease data.
<dig> .
<dig> performance using healthy data
a set of  <dig> pbmc gene expression profiles from healthy subjects was collected from six published studies  <cit>  .
figure  <dig> summarizes the results of running spec on these data using the split subset signature validation test described above.
results were binned into significant vs. non-significant groups using a p-value cutoff of  <dig> .
overall, spec very accurately predicts the source of the query signal for all split signatures.
nearly all significant predictions correctly linked the query to its proper subset and, even including the non-significant results, there were relatively few false positives.
for most of the subset signatures, approximately 75% of the runs were significant and correct.
as an important negative control, we found that most results from the neutrophil subset did not reach statistical significance reflecting that fact that granulocytes  are mostly absent in pbmcs  <cit> .
interestingly, a large portion of the non-significant results were still correct, most likely because contamination from neutrophils still provides a small amount of signal in the pbmc samples.
<dig> .
<dig> performance in patients with disease
the previous validation experiments show that spec works for pbmc expression profiles from healthy subjects, but it is possible that the altered expression patterns in disease patients could interfere with the accuracy of spec.
to verify that spec works in a disease setting, we repeated the validation steps described in section  <dig> .
<dig> using pbmc expression data from  <dig> patients with hcv  <cit>  , as well as  <dig> patients with systemic lupus erythematosus   <cit>  .
for the hcv data, nearly all significant results correctly predicted the source of the query signature, although overall there were fewer results that were considered significant than in the healthy data.
the sle data, on the other hand, actually performed better on the b and t signatures with 100% of the runs being correct and significant but, as expected, spec was often unable to predict the correct source of a neutrophil signature.
the monocyte signature also was incorrectly associated with neutrophils about half the time, suggesting that performance could be improved by removing cell subset signatures that should be absent from a sample.
we conclude that, despite some quantitative differences in the performance measures compared with data from healthy subjects, spec is not overly affected by potential differential gene expression induced in sle or hcv, and can be effectively applied to data from individuals in a disease setting.
values are the percentage of time the query signature  had the highest correlation with each subset signature .
values are of the form: % significant .
cells marked in bold indicate the highest value for each query signature.
rows may not add up to 100% due to rounding.
values are the percentage of time the query signature  had the highest correlation with each subset signature .
values are of the form: % significant .
cells marked in bold indicate the highest value for each query signature.
rows may not add up to 100% due to rounding.
<dig>  validation of spec using independently-derived signatures
<dig> .
<dig> performance using independent subset signatures
an alternative method for validating spec is to use signatures for the same cell subset that were generated from independent sources.
to carry this out, we rely on studies by abbas et al.
<cit>  and palmer et al.
<cit> , which both defined signatures for b cells and t cells, as well as a signature for either granulocytes or neutrophils .
while these signatures have surprisingly little overlap with each other, we nevertheless removed the small number of common genes from the query signatures before running spec to avoid biasing the performance measurements.
as shown in table  <dig>  spec correctly predicts the source of the b cell signature, and the granulocyte signature was appropriately identified with the neutrophil signature.
the granulocyte signature is also correctly predicted as coming from the myeloid cell subset .
the t cell signature on the other hand does not show a significant correlation with any subset signature.
this is likely caused by problems with one of the t cell signatures, as we already found that this signature had the poorest performance in the validation of enrichment scores using deconvolution values .
values indicate the correlation between the enrichment scores calculated by the primary signatures  and the independently derived signatures .
highlighted cells are correlations with p <  <dig>
<dig> .
<dig> performance using experimentally-derived subset-specific disease signatures
our final validation approach involved testing the ability of spec to accurately predict the source of experimentally-determined query signatures derived from purified cell subsets from sle patients.
in this case transcriptional profiling data was measured specifically from purified b cells, t cells and myeloid cells  <cit> .
we defined cell subset-signatures as the set of genes that were up-regulated in sle patients  and that did not overlap with differentially-expressed genes from other cell subsets .
sle signatures for the three cell subsets were used as query signatures and the subset source was predicted with spec.
as expected, query signatures created using a particular cell subset were linked by spec to the correct subset using only the gene expression data from total pbmcs.
the b cell specific disease signature has the highest correlation with the b cell subset signature .
similarly, the t-cell specific disease signature has the highest correlation with the t cell subset signature .
the myeloid signature had high correlation with myeloid cells in the myeloid vs. lymphoid comparison , and also had the highest correlations with all of the myeloid-derived subset signatures in the more detailed comparisons.
thus, spec is capable of predicting the source of a disease-associated subset-specific gene expression pattern using only gene expression from total pbmcs.
<dig> .
<dig> effect of query signature size on accuracy
the size of the sle subset-specific signatures used in the previous section was fixed as the top- <dig> up-regulated genes.
we used a fixed cutoff so that we could easily compare between results on the different subsets.
in practice, drawing a cutoff or choosing the number of genes to include in a gene expression signature can be difficult, and the number of genes that reach statistical significance for differential expression may be small.
thus, we sought to investigate the effect of query gene signature size on the accuracy of spec predictions using the sle disease signatures from section  <dig> .
<dig>  the size of the query signature was successively decreased, and spec was applied to each of the query signatures .
as expected, p-values generally increased with smaller signatures and, in this case, the results were no longer significant when gene signatures contained less than  <dig> genes.
interestingly, this is also the minimum number of genes recommended for calculating robust enrichment scores in gsea  <cit> .
all significant results were correct and, even with the smaller subset sizes that failed to reach statistical significance, spec often predicted the correct subset source.
these results demonstrate the importance of robust subset signatures in the application of spec, but also highlight the fact that spec can be effective even when using relatively small query gene sets.
<dig>  case study: identifying the isg signature source for predicting therapy response in chronic hcv patients
the current standard therapy for chronic hcv consists of pegylated  interferon -α and ribavirin.
however, this treatment is costly, poorly tolerated due to adverse side effects, and ineffective in many patients  <cit> .
the majority of hcv infections in the united states are with genotype  <dig>  which responds poorly to therapy   <cit> .
methods to identify patients who are likely to be responsive to therapy would be a great help in clinical decision making.
in addition, the ability to predict non-responsiveness could allow alternate treatments to be explored.
previous studies have found that patients who fail to achieve a sustained virologic response  after treatment display high baseline levels of ifn-stimulated gene  expression in hepatocytes prior to the initiation of therapy  <cit> , and the expression of these isgs is not increased by ifn treatment.
however, gene expression changes in total pbmcs in hcv patients appear blunted compared with looking directly at hepatocytes, making it difficult to develop highly predictive models from this cell population  <cit> .
because there are significant advantages to using pbmcs for interrogating disease processes and performing clinical tests, it is important to investigate ways to improve the predictive power of pbmc gene expression profiling.
we hypothesize that the accuracy of existing models based on pbmcs is limited because differential expression of isgs is likely to be cell type-specific, and thus their signal will be diluted in total pbmc by rna from other cells.
this logic applies to non-isg signatures as well.
therefore, using specific pbmc subsets for transcriptional profiling should improve our ability to extract knowledge from these data.
unfortunately, it is not obvious which cell subset are the most informative, so we have applied spec to predict the most likely cellular source for predictive signatures in hcv.
we first generated a list of genes that were significantly up-regulated in patients classified as non-responders  based on gene expression profiling data from total pbmcs in  <dig> chronic hcv patients generated by sarasin-filipowicz et al.
<cit> .
using this gene list as the query signature, spec was then applied to the pbmc expression profiles from only the  <dig> non-responders in this study  in order to predict the cellular source of the up-regulated genes.
spec identified the myeloid subset , as the most likely source for this query signature.
none of the associations with the detailed signatures reached statistical significance.
this could be due in part to the relatively small number of samples used in the analysis.
it is also possible that the therapy-response signature arises from multiple cell types in the myeloid family, or possibly from dendritic cells , which were not included in our subsets for reasons previously described.
overall, these results suggest that gene expression analysis on the myeloid subset should improve our ability to predict therapy response in hcv patients.
<dig> discussion
blood is a mixture of many different cell types.
the ability to associate a predictive gene expression signature with a specific cell subset would allow better classifiers to be developed through the experimental isolation of the most relevant cells.
in addition, knowledge of which cells are being affected would provide important biological insights.
no other computational methods exist to directly address this problem.
in this study we developed a computational method, subset prediction from enrichment correlation , to predict the most likely cellular source of a predictive gene expression signature.
our method requires only gene expression data from total pbmcs, and does not depend on knowing the subset proportions in the blood sample.
spec uses enrichment scores to estimate the strength of gene expression signatures across individuals in a population.
the pbmc subset signature showing the highest correlation with the query gene signature is predicted to be the source.
we show that spec has good performance on both healthy and disease data, and correctly predicts the source of experimentally-derived subset-specific genes in sle patients.
as a case study, we applied spec to transcriptional profiling data from chronic hcv patients in order to predict the potential source of a gene expression signature associated with non-responsiveness to standard therapy.
using spec, we predict that a stronger signal for non-responsiveness will be apparent when focusing on the myeloid subset.
we hope that looking at the cell subpopulations separately will allow us to develop a blood test that can predict the clinical outcome of antiviral therapy prior to treatment.
the quality of cell subset gene expression signatures is critical to the performance of spec.
it is clear from our results  that some signatures are simply better than others.
in this work, we use gene expression signatures for b cells, t cells, nk cells, neutrophils and monocytes available from  <cit>  and  <cit> .
we decided to exclude the dc subset signatures from  <cit>  since we found that it did not work well within the spec framework.
we speculate that this is due to the fact that dcs are such a small component of the total pbmcs .
more work is needed to define the qualities of a "good" signature for the spec framework.
other groups have proposed additional signatures that we plan to evaluate  <cit> , with the long-term goal of defining a set of non-overlapping signatures for a detailed set of pbmc subsets.
as an alternative to applying spec simultaneously to all available subsets, it is possible spec may be applied in a hierarchical manner: first predicting whether the query comes from the lymphoid vs. myeloid subsets then, if lymphoid, predict b cell vs. t cell vs. nk cell then, if b cell, predict naïve b cell vs. memory b cell vs... and so on.
an additional area for future investigation concerns alternatives to the enrichment score, which is currently used for estimating both the relative cell subset proportions and expression of the query signature.
one possibility that we have already investigated is using microarray deconvolution, described in  <cit>  and in section  <dig> , to measure cell subset fractions.
however, we found this method actually reduced the accuracy of spec .
future work could explore other methods for estimating cell proportions from gene expression data  <cit> .
one thing to keep in mind is that spec does not require estimating the absolute proportion of all cell types within a single individual as provided by deconvolution methods, but rather assumes that we can estimate the relative proportion for each individual cell subset across a population.
by trying to solve a harder problem, deconvolution methods may actually introduce more error compared with enrichment scores.
another advantage of enrichment scores is that they are likely to be less dependent on sample source .
<dig>
CONCLUSIONS
this work describes the implementation and validation of a computational methodology to support blood genomics studies with general applicability to basic and translational research.
spec predicts the cellular source of predictive gene expression signatures.
when applied to data from chronic hcv patients, spec predicts that focusing on myeloid cells may enrich for an interferon signature that has been observed in hepatocytes from patients that do not respond to standard therapy.
the approach requires only transcriptional profiling data from total pbmcs, allowing wide applicability, and should improve the predictive power of blood gene expression profiling studies by allowing experiments to focus on the most informative cell subsets.
it is also easy to envision how this framework could be combined with classification methods to define predictive signatures that are more likely to be cell subset-specific, and therefore amenable to experimental signal enrichment through sorting.
<dig> methods
microarray data and normalization
raw affymetrix microarray data were downloaded from geo.
transcriptional profiles from pbms were obtained for healthy controls , sle patients , and chronic hcv patients .
cell specific transcriptional profiles for b cells, t cells and myeloid cells from sle patients were obtained from gse <dig>  these data were normalized with gcrma  <cit>  using the bioconductor software package in r  <cit> .
generation of sle cell specific disease signatures
differentially-expressed genes were determined for each cell subset  by comparing sle and healthy individuals using the limma package in bioconductor  <cit> .
genes were separated into groups of upregulated and downregulated genes, ranked by fdr, and the top  <dig> genes were pulled out.
cell specific disease signatures were created by using only the genes unique to each signature.
generation of hcv response signature
gene expression data was obtained from the study by sarasin-filipowicz et al.
<cit> .
differentially-expressed genes were determined by comparing non-responders vs. responders using the limma package in bioconductor  <cit> .
significance was defined using a false discovery rate cutoff of  <dig>  using the method of benjamini and hochberg.
enrichment score calculation
the enrichment score is calculated as described in  <cit>  with p =  <dig> and weighted by the gene rank.
genes are ranked based on their expression values from a single transcriptional profile .
the enrichment score is the maximum deviation from zero of a running sum statistic which is calculated by walking down the ranked list of genes.
the sum is incremented whenever a gene in the signature is encountered and decremented if the gene is not in the signature.
increments and decrements are weighted so the statistic sums to  <dig> over all the genes.
microarray deconvolution
deconvolution was performed to estimate the fraction of each cell subset from a single pbmc transcriptional profile using the methods described in  <cit> .
the cell subset expression data used for deconvolution was downloaded from the supplementary material.
the predicted cell proportions of related subsets were added together to generate the final cell fraction .
monte carlo permutation test for p value estimation
to calculate a p value for the maximum correlation coefficient found by spec, we estimated the distribution of these values under the null hypothesis by randomly permuting the assignment of genes on the microarray to the query signature, while maintaining the size of the query signature.
we then ran spec to determine the maximum correlation between the enrichment scores for these new queries and each subset signature.
this process was repeated  <dig> times to generate a distribution of maximum correlations.
p-values were then estimated by fitting a normal distribution to these data and calculating the cumulative density for each spec result.
implementation and availability
spec was implemented in the r statistical programming language, and is available at http://clip.med.yale.edu/spec.
additional requests can be made by contacting crb or shk.
authors' contributions
crb helped design the study, programmed the spec algorithm and carried out the validation and analysis.
mu helped with programming and the validation using the sle data.
shk conceived of the study and participated in its design and coordination.
crb and shk drafted the manuscript.
all authors have read and approved the final manuscript.
competing interests
the authors declare that they have no competing interests.

 gene expression data from blood genomic studies are widely used for investigation of human disease. predictive gene signatures have been developed to carry out differential diagnosis of infectious diseases, identify specific disease states and characterize the immune response to vaccination. however, is some cases gene expression signatures from blood can be weakly expressed and highly variable. the identification of these signatures is complicated by the fact that blood is a mixed tissue, composed of multiple cell subsets, so that differential expression profiles can reflect changes in cell subset proportions, changes in subset-specific gene expression or both. in cases where the relevant disease genes are expressed in a subset-specific manner, as has been shown for sle, analyses based on mixed cell expression data are inherently limited since differential expression of genes in one cell subset will be diluted by rna from other cells. experimental studies that isolate specific cell subsets before expression profiling can provide important biological insight by demonstrating subset-specific gene expression as well as increased predictive signal. given the large number of potential cell subsets that can be defined, the ability to identify the most informative subset to isolate would be a great aid to these studies. most genome-wide expression studies are based on analysis of total peripheral blood mononuclear cells. pbmcs are composed of over a dozen cell subsets that are derived from a common progenitor in the bone marrow. these cells are commonly divided into myeloid and lymphoid cells. myeloid cells include monocytes and their descendants, as well as granulocytes like neutrophils and basophils. lymphoid cells are primarily composed of b cells, t cells and nk cells. proportions of these cells can vary widely between individuals, but t cells and b cells together usually make up ~75% of pbmcs, while nk cells and monocytes make up around 10-15% each. the remaining cell types, such as dendritic cells, are much more rare, and account for <1% of total pbmcs. neutrophils, which normally compose the majority of cells in a blood sample, are normally excluded by the methods used to isolate pbmcs, but may account for up to 20% of a pbmc sample due to contamination. genome-wide expression measurements based on total pbmcs reflect both condition-specific gene expression as well as the proportion of different cell subsets in the sample. microarray deconvolution methods have been developed to take advantage of this latter dependence in order to quantify the relative proportion of different cell subsets. in this approach, the expression level of each gene is modeled as a linear function of the expression from each cell subset comprising the sample. deconvolution thus depends on prior knowledge of quantitative expression levels for each subset. in cases where the subset proportions have been measured experimentally, methods have been developed to compute subset-specific differential-expression for individual genes. however, in addition to requiring measurements that are not available for many existing datasets, these methods do not directly address the problem of predicting the cell subset source of a pre-defined gene signature. in some cases, clues to the cellular source of a gene expression signature may be found in the signature itself. this would be true when the signature is enriched for genes known to be expressed in a particular cell type, and may occur when the proportion of cells in the blood changes in the condition of interest. for example, berry et al. found that changes in t-cell related genes included in a transcriptional signature of active tuberculosis resulted from changes in t cell number in the blood. in cases where the predictive signature does not include a significant number of cell subset-specific genes, the problem may be addressed experimentally through the isolation of specific cell subsets followed by transcriptional profiling. however, it is often not obvious which cell subset will be the most informative. the time and expense of isolating and analyzing all possibly relevant cell subsets provides the motivation for the method developed here: a computational approach to predict the most likely cell subset source of a gene expression signature. the development of this method was also motivated by the clinical problem of predicting the outcome of antiviral therapy in patients chronically infected with hepatitis c virus. in this case, gene expression changes that are predictive when measured in hepatocytes are significantly blunted in pbmcs, making it difficult to develop highly predictive models from this cell population. the isolation of specific cell subsets prior to expression profiling may allow for the identification of predictive signatures based on a blood sample, with significant advantages for the patient. we have developed a computational method to predict the most likely cellular source for a pre-defined gene expression signature using transcriptional profiling data from total pbmcs. the general approach we propose is to: estimate the relative proportion of each pbmc subset from individual patients using genome-wide transcriptional profiling data, and correlate these proportions with the gene expression signature that predicts therapy. we theorize that the pbmc subset proportion that correlates most closely with the predicted gene expression signature over a large set of patients is the most likely source of the signal. as a specific example, let's consider a study where we find that interferon stimulated genes are up-regulated in a subset of patients that fail to respond to a particular therapy. now, we focus on the non-responding patients. if isg expression tends to be higher in patients where the proportion of dendritic cells is predicted to be high, then we would predict that dcs are the most likely source of the observed gene expression signature. note that there can be two interpretations of this result with respect to the original study. first, the dcs may always express isgs, but the proportion of dcs in the blood is higher for non-responders. second, the dcs undergo specific gene expression changes in non-responders which cause the level of isg expression to increase in the blood sample. in either case, dcs are the relevant population to follow-up experimentally. here, we develop this general idea into a computational method called spec, which we then validate using both healthy and disease data. while the development of spec was motivated by the problem of predicting the response to therapy in hcv patients, the method is generally applicable to link gene expression signatures with specific pbmc subsets. since only gene expression profiles from total pbmcs are required, the technique can easily be used to mine the massive amount of existing microarray data.  subset prediction from enrichment correlation is a method that takes advantage of the fluctuations of cell proportions in a mixed cell sample to find the most likely source for a gene expression signal. the general approach can be broken down into four steps: estimate the relative proportion of each pbmc subset across a population. we propose to accomplish this by calculating the enrichment score associated with a set of genes expressed specifically in each subset s, called a subset signature. this is summarized by the following pseudo-code: code: n_p = number of individuals n_s = number of pbmc subsets exprdata = gene expression profile for individual i subsetsignature = list of genes representing pbmc subset s querysignature = a list of genes input by the user for each for each e = enrichment of subsetsignature in exprdata calculate the enrichment score of the query signature for each individual i in the population: for each q = enrichment of querysignature in exprdata determine the pearson correlation between the enrichment scores of the query signature and each pbmc subset s: for each c = correlation between q and e across i predict the subset with the highest correlation to the query signature: predicted subset = s that gives rise to the max the calculation of enrichment scores is explained in methods and figure estimating relative pbmc subset proportions using enrichment scores the first step in our proposed method is to estimate the relative proportion of each pbmc subset. while a number of methods exist that can estimate these values, we chose to use enrichment scores calculated from cell subset-specific gene signatures for our initial implementation of spec due to its relative simplicity. as described in the methods, the enrichment score calculates the extent to which a gene signature is concentrated among the most highly expressed genes in a single genome-wide expression experiment. the cell subset-specific signatures are defined as groups of genes which have been found to be up-regulated in only a specific subset of cells in a mixed cell population. the resulting lists of up-regulated genes represent cell markers for that blood cell type, and include well known cell surface markers such as cd for t cells and cd for b cells. lists of genes describing specific pbmc subsets have previously been defined in and in. here, we choose to use signatures representing the b cell, t cell, nk cell, and monocyte subsets, we also included a signature for neutrophils despite the fact that they are not considered pbmcs, since neutrophil contamination can account for up to 20% of a pbmc sample. in addition, we also use a more general set of signatures for the lymphoid and myeloid subsets. the lymphoid and myeloid signatures are composed of genes which are expressed in more than one cell type from those lineages, thus these subsets should vary according to the sum of all the cells in that lineage. we run spec separately on the detailed and general signatures, as described below. we hypothesize that the enrichment scores of the subset-specific gene signatures will be related to the relative cell subset proportions across patients. although microarray deconvolution can be used to accomplish this task, these methods require quantitative gene expression levels for each subset to be estimated, which is not always available. in addition, it is important to be clear that we do not require the enrichment scores of different cell subsets to predict their proportions in a single individual. rather, we propose that the enrichment score for a single cell subset will vary along with the relative size of the cell subset across patients. for example, if subject b has a higher fraction of nk cells than subject a then the nk signature enrichment score should also be higher, but both of these enrichment scores may be higher than the t cell signature enrichment scores even if nk cells constitute a lower fraction of cells in both subjects. in order to validate this enrichment score approach, we compared it with predictions from microarray deconvolution using data from healthy subjects where we expect deconvolution should perform relatively well. we predicted the relative proportions of cell subsets for healthy gene expression profiles using both microarray deconvolution and the cell signature enrichment scores. figure is an example showing the correlation between the nk cell signature enrichment score and the predicted nk cell fraction estimated using deconvolution. the correlation of implies that the enrichment score is positively varying with the relative proportion of nk cells in these healthy pbmc samples. we found similarly high correlations for b cells, neutrophils and monocytes. in contrast, t cells had a much weaker correlation, perhaps reflecting the quality of this particular gene signature. overall, these results suggest that enrichment scores can be used to track the fluctuations in the levels of each cell type in a mixed cell expression profile. validation of spec using split subset signatures as an initial validation of spec, we sought to predict the source of a gene expression signature where the correct cell subset is known. this can be done by using our pre-defined cell subset expression signatures. for example, consider the b cell signature defined in. the genes that are part of this signature can be split into two disjoint sets. one of these sets can be labeled as the cell subset signature and used to estimate the changes in the b cell population, while the other can be retained as the query gene signature being tested. a successful method would be able to identify this query signature as coming from b cells or from lymphoid cells. when running spec, we separately analyze the set of detailed signatures and general signatures, which represent cell types that are a superset of those covered by detailed signatures. for each test, one of the subset signatures was randomly divided into two halves, with one half used as the query signature and the other half retained as the subset signature. this was repeated times for each of the pbmc subset signatures. in order to determine the significance of the correlation value for each run we used a monte carlo permutation test. the following two subsections describe the performance of spec on healthy and disease data.. performance using healthy data a set of pbmc gene expression profiles from healthy subjects was collected from six published studies. figure summarizes the results of running spec on these data using the split subset signature validation test described above. results were binned into significant vs. non-significant groups using a p-value cutoff of. overall, spec very accurately predicts the source of the query signal for all split signatures. nearly all significant predictions correctly linked the query to its proper subset and, even including the non-significant results, there were relatively few false positives. for most of the subset signatures, approximately 75% of the runs were significant and correct. as an important negative control, we found that most results from the neutrophil subset did not reach statistical significance reflecting that fact that granulocytes are mostly absent in pbmcs. interestingly, a large portion of the non-significant results were still correct, most likely because contamination from neutrophils still provides a small amount of signal in the pbmc samples.. performance in patients with disease the previous validation experiments show that spec works for pbmc expression profiles from healthy subjects, but it is possible that the altered expression patterns in disease patients could interfere with the accuracy of spec. to verify that spec works in a disease setting, we repeated the validation steps described in section. using pbmc expression data from patients with hcv, as well as patients with systemic lupus erythematosus. for the hcv data, nearly all significant results correctly predicted the source of the query signature, although overall there were fewer results that were considered significant than in the healthy data. the sle data, on the other hand, actually performed better on the b and t signatures with 100% of the runs being correct and significant but, as expected, spec was often unable to predict the correct source of a neutrophil signature. the monocyte signature also was incorrectly associated with neutrophils about half the time, suggesting that performance could be improved by removing cell subset signatures that should be absent from a sample. we conclude that, despite some quantitative differences in the performance measures compared with data from healthy subjects, spec is not overly affected by potential differential gene expression induced in sle or hcv, and can be effectively applied to data from individuals in a disease setting. values are the percentage of time the query signature had the highest correlation with each subset signature. values are of the form: % significant. cells marked in bold indicate the highest value for each query signature. rows may not add up to 100% due to rounding. values are the percentage of time the query signature had the highest correlation with each subset signature. values are of the form: % significant. cells marked in bold indicate the highest value for each query signature. rows may not add up to 100% due to rounding. validation of spec using independently-derived signatures. performance using independent subset signatures an alternative method for validating spec is to use signatures for the same cell subset that were generated from independent sources. to carry this out, we rely on studies by abbas et al. and palmer et al., which both defined signatures for b cells and t cells, as well as a signature for either granulocytes or neutrophils. while these signatures have surprisingly little overlap with each other, we nevertheless removed the small number of common genes from the query signatures before running spec to avoid biasing the performance measurements. as shown in table spec correctly predicts the source of the b cell signature, and the granulocyte signature was appropriately identified with the neutrophil signature. the granulocyte signature is also correctly predicted as coming from the myeloid cell subset. the t cell signature on the other hand does not show a significant correlation with any subset signature. this is likely caused by problems with one of the t cell signatures, as we already found that this signature had the poorest performance in the validation of enrichment scores using deconvolution values. values indicate the correlation between the enrichment scores calculated by the primary signatures and the independently derived signatures. highlighted cells are correlations with p <. performance using experimentally-derived subset-specific disease signatures our final validation approach involved testing the ability of spec to accurately predict the source of experimentally-determined query signatures derived from purified cell subsets from sle patients. in this case transcriptional profiling data was measured specifically from purified b cells, t cells and myeloid cells. we defined cell subset-signatures as the set of genes that were up-regulated in sle patients and that did not overlap with differentially-expressed genes from other cell subsets. sle signatures for the three cell subsets were used as query signatures and the subset source was predicted with spec. as expected, query signatures created using a particular cell subset were linked by spec to the correct subset using only the gene expression data from total pbmcs. the b cell specific disease signature has the highest correlation with the b cell subset signature. similarly, the t-cell specific disease signature has the highest correlation with the t cell subset signature. the myeloid signature had high correlation with myeloid cells in the myeloid vs. lymphoid comparison, and also had the highest correlations with all of the myeloid-derived subset signatures in the more detailed comparisons. thus, spec is capable of predicting the source of a disease-associated subset-specific gene expression pattern using only gene expression from total pbmcs.. effect of query signature size on accuracy the size of the sle subset-specific signatures used in the previous section was fixed as the top- up-regulated genes. we used a fixed cutoff so that we could easily compare between results on the different subsets. in practice, drawing a cutoff or choosing the number of genes to include in a gene expression signature can be difficult, and the number of genes that reach statistical significance for differential expression may be small. thus, we sought to investigate the effect of query gene signature size on the accuracy of spec predictions using the sle disease signatures from section. the size of the query signature was successively decreased, and spec was applied to each of the query signatures. as expected, p-values generally increased with smaller signatures and, in this case, the results were no longer significant when gene signatures contained less than genes. interestingly, this is also the minimum number of genes recommended for calculating robust enrichment scores in gsea. all significant results were correct and, even with the smaller subset sizes that failed to reach statistical significance, spec often predicted the correct subset source. these results demonstrate the importance of robust subset signatures in the application of spec, but also highlight the fact that spec can be effective even when using relatively small query gene sets. case study: identifying the isg signature source for predicting therapy response in chronic hcv patients the current standard therapy for chronic hcv consists of pegylated interferon -α and ribavirin. however, this treatment is costly, poorly tolerated due to adverse side effects, and ineffective in many patients. the majority of hcv infections in the united states are with genotype which responds poorly to therapy. methods to identify patients who are likely to be responsive to therapy would be a great help in clinical decision making. in addition, the ability to predict non-responsiveness could allow alternate treatments to be explored. previous studies have found that patients who fail to achieve a sustained virologic response after treatment display high baseline levels of ifn-stimulated gene expression in hepatocytes prior to the initiation of therapy, and the expression of these isgs is not increased by ifn treatment. however, gene expression changes in total pbmcs in hcv patients appear blunted compared with looking directly at hepatocytes, making it difficult to develop highly predictive models from this cell population. because there are significant advantages to using pbmcs for interrogating disease processes and performing clinical tests, it is important to investigate ways to improve the predictive power of pbmc gene expression profiling. we hypothesize that the accuracy of existing models based on pbmcs is limited because differential expression of isgs is likely to be cell type-specific, and thus their signal will be diluted in total pbmc by rna from other cells. this logic applies to non-isg signatures as well. therefore, using specific pbmc subsets for transcriptional profiling should improve our ability to extract knowledge from these data. unfortunately, it is not obvious which cell subset are the most informative, so we have applied spec to predict the most likely cellular source for predictive signatures in hcv. we first generated a list of genes that were significantly up-regulated in patients classified as non-responders based on gene expression profiling data from total pbmcs in chronic hcv patients generated by sarasin-filipowicz et al.. using this gene list as the query signature, spec was then applied to the pbmc expression profiles from only the non-responders in this study in order to predict the cellular source of the up-regulated genes. spec identified the myeloid subset, as the most likely source for this query signature. none of the associations with the detailed signatures reached statistical significance. this could be due in part to the relatively small number of samples used in the analysis. it is also possible that the therapy-response signature arises from multiple cell types in the myeloid family, or possibly from dendritic cells, which were not included in our subsets for reasons previously described. overall, these results suggest that gene expression analysis on the myeloid subset should improve our ability to predict therapy response in hcv patients. discussion blood is a mixture of many different cell types. the ability to associate a predictive gene expression signature with a specific cell subset would allow better classifiers to be developed through the experimental isolation of the most relevant cells. in addition, knowledge of which cells are being affected would provide important biological insights. no other computational methods exist to directly address this problem. in this study we developed a computational method, subset prediction from enrichment correlation, to predict the most likely cellular source of a predictive gene expression signature. our method requires only gene expression data from total pbmcs, and does not depend on knowing the subset proportions in the blood sample. spec uses enrichment scores to estimate the strength of gene expression signatures across individuals in a population. the pbmc subset signature showing the highest correlation with the query gene signature is predicted to be the source. we show that spec has good performance on both healthy and disease data, and correctly predicts the source of experimentally-derived subset-specific genes in sle patients. as a case study, we applied spec to transcriptional profiling data from chronic hcv patients in order to predict the potential source of a gene expression signature associated with non-responsiveness to standard therapy. using spec, we predict that a stronger signal for non-responsiveness will be apparent when focusing on the myeloid subset. we hope that looking at the cell subpopulations separately will allow us to develop a blood test that can predict the clinical outcome of antiviral therapy prior to treatment. the quality of cell subset gene expression signatures is critical to the performance of spec. it is clear from our results that some signatures are simply better than others. in this work, we use gene expression signatures for b cells, t cells, nk cells, neutrophils and monocytes available from and. we decided to exclude the dc subset signatures from since we found that it did not work well within the spec framework. we speculate that this is due to the fact that dcs are such a small component of the total pbmcs. more work is needed to define the qualities of a "good" signature for the spec framework. other groups have proposed additional signatures that we plan to evaluate, with the long-term goal of defining a set of non-overlapping signatures for a detailed set of pbmc subsets. as an alternative to applying spec simultaneously to all available subsets, it is possible spec may be applied in a hierarchical manner: first predicting whether the query comes from the lymphoid vs. myeloid subsets then, if lymphoid, predict b cell vs. t cell vs. nk cell then, if b cell, predict naïve b cell vs. memory b cell vs... and so on. an additional area for future investigation concerns alternatives to the enrichment score, which is currently used for estimating both the relative cell subset proportions and expression of the query signature. one possibility that we have already investigated is using microarray deconvolution, described in and in section, to measure cell subset fractions. however, we found this method actually reduced the accuracy of spec. future work could explore other methods for estimating cell proportions from gene expression data. one thing to keep in mind is that spec does not require estimating the absolute proportion of all cell types within a single individual as provided by deconvolution methods, but rather assumes that we can estimate the relative proportion for each individual cell subset across a population. by trying to solve a harder problem, deconvolution methods may actually introduce more error compared with enrichment scores. another advantage of enrichment scores is that they are likely to be less dependent on sample source.  this work describes the implementation and validation of a computational methodology to support blood genomics studies with general applicability to basic and translational research. spec predicts the cellular source of predictive gene expression signatures. when applied to data from chronic hcv patients, spec predicts that focusing on myeloid cells may enrich for an interferon signature that has been observed in hepatocytes from patients that do not respond to standard therapy. the approach requires only transcriptional profiling data from total pbmcs, allowing wide applicability, and should improve the predictive power of blood gene expression profiling studies by allowing experiments to focus on the most informative cell subsets. it is also easy to envision how this framework could be combined with classification methods to define predictive signatures that are more likely to be cell subset-specific, and therefore amenable to experimental signal enrichment through sorting. methods microarray data and normalization raw affymetrix microarray data were downloaded from geo. transcriptional profiles from pbms were obtained for healthy controls, sle patients, and chronic hcv patients. cell specific transcriptional profiles for b cells, t cells and myeloid cells from sle patients were obtained from gse these data were normalized with gcrma using the bioconductor software package in r. generation of sle cell specific disease signatures differentially-expressed genes were determined for each cell subset by comparing sle and healthy individuals using the limma package in bioconductor. genes were separated into groups of upregulated and downregulated genes, ranked by fdr, and the top genes were pulled out. cell specific disease signatures were created by using only the genes unique to each signature. generation of hcv response signature gene expression data was obtained from the study by sarasin-filipowicz et al.. differentially-expressed genes were determined by comparing non-responders vs. responders using the limma package in bioconductor. significance was defined using a false discovery rate cutoff of using the method of benjamini and hochberg. enrichment score calculation the enrichment score is calculated as described in with p = and weighted by the gene rank. genes are ranked based on their expression values from a single transcriptional profile. the enrichment score is the maximum deviation from zero of a running sum statistic which is calculated by walking down the ranked list of genes. the sum is incremented whenever a gene in the signature is encountered and decremented if the gene is not in the signature. increments and decrements are weighted so the statistic sums to over all the genes. microarray deconvolution deconvolution was performed to estimate the fraction of each cell subset from a single pbmc transcriptional profile using the methods described in. the cell subset expression data used for deconvolution was downloaded from the supplementary material. the predicted cell proportions of related subsets were added together to generate the final cell fraction. monte carlo permutation test for p value estimation to calculate a p value for the maximum correlation coefficient found by spec, we estimated the distribution of these values under the null hypothesis by randomly permuting the assignment of genes on the microarray to the query signature, while maintaining the size of the query signature. we then ran spec to determine the maximum correlation between the enrichment scores for these new queries and each subset signature. this process was repeated times to generate a distribution of maximum correlations. p-values were then estimated by fitting a normal distribution to these data and calculating the cumulative density for each spec result. implementation and availability spec was implemented in the r statistical programming language, and is available at http://clip.med.yale.edu/spec. additional requests can be made by contacting crb or shk. authors' contributions crb helped design the study, programmed the spec algorithm and carried out the validation and analysis. mu helped with programming and the validation using the sle data. shk conceived of the study and participated in its design and coordination. crb and shk drafted the manuscript. all authors have read and approved the final manuscript. competing interests the authors declare that they have no competing interests. 
most studies are based on analysis of total peripheral blood mononuclear cells, a mixed population. using spec, we predict that myeloid cells are the source of the interferon-therapy response gene signature associated with hcv patients who are non-responsive to standard therapy. spec does not rely on the occurrence of cell subset-specific genes in the signature, but rather takes advantage of correlations with subset-specific genes across a set of samples. we have developed a computational method to predict the cellular source for a pre-defined list of genes using only data from total pbmcs. in this case, accuracy is inherently limited since cell subset-specific differential expression of gene signatures will be diluted by rna from other cells. validation using multiple experimental datasets demonstrates that spec can accurately identify the source of a gene signature as myeloid or lymphoid, as well as differentiate between b cells, t cells, nk cells and monocytes. spec is widely applicable since only gene expression profiles from total pbmcs are required, and thus it can easily be used to mine the massive amount of existing microarray or rna-seq data. while using specific pbmc subsets for transcriptional profiling would improve our ability to extract knowledge from these data, it is rarely obvious which cell subset will be the most informative. genome-wide transcriptional profiling of patient blood samples offers a powerful tool to investigate underlying disease mechanisms and personalized treatment decisions. it can help identify specific cell subsets that are important for understanding disease and therapy response. spec is a powerful technique for blood genomic studies.
 genome-wide transcriptional profiling of patient blood samples offers a powerful tool to investigate underlying disease mechanisms and personalized treatment decisions. most studies are based on analysis of total peripheral blood mononuclear cells, a mixed population. in this case, accuracy is inherently limited since cell subset-specific differential expression of gene signatures will be diluted by rna from other cells. while using specific pbmc subsets for transcriptional profiling would improve our ability to extract knowledge from these data, it is rarely obvious which cell subset will be the most informative.  we have developed a computational method to predict the cellular source for a pre-defined list of genes using only data from total pbmcs. spec does not rely on the occurrence of cell subset-specific genes in the signature, but rather takes advantage of correlations with subset-specific genes across a set of samples. validation using multiple experimental datasets demonstrates that spec can accurately identify the source of a gene signature as myeloid or lymphoid, as well as differentiate between b cells, t cells, nk cells and monocytes. using spec, we predict that myeloid cells are the source of the interferon-therapy response gene signature associated with hcv patients who are non-responsive to standard therapy.  spec is a powerful technique for blood genomic studies. it can help identify specific cell subsets that are important for understanding disease and therapy response. spec is widely applicable since only gene expression profiles from total pbmcs are required, and thus it can easily be used to mine the massive amount of existing microarray or rna-seq data. 