BACKGROUND
nanomaterials are being manufactured on a commercial scale for use in medical, diagnostic, energy, component and communications industries  <cit> .
engineered nanomaterials range considerably in their physicochemical properties making them more desirable than their micro- and macro-counterparts due to, for example, their increased surface area, tensile strength, tunability, etc.
<cit> .
from limited early reports, concerns over the safety of engineered nanomaterials have surfaced  <cit> .
humans can be exposed to nanomaterials in different ways such as inhalation or exposure through the integumentary system.
however, the skin is a unique organ in the body not only because it gives the body such a large surface area for exposure but also because of the avascular property of epidermis, in which particles can reside without being removed by phagocytosis  <cit> .
gene expression microarrays have become a tool to investigate the interactions of biological systems by observing the simultaneous activities of tens of thousands of genes.
over recent years, this tool has been applied to toxicology forming a new discipline, toxicogenomics  <cit> .
microarrays have most recently been a tool used by pharmaceutical drug discovery and development to screen for efficacy and adverse effects thereby prioritizing drug candidates and redeveloping ones which show off-target toxicities  <cit> .
the approach described here combines a global screening technology, gene expression microarray profiling, with systems biology, to investigate the interactions of engineered nanomaterials with primary human cells.
the biological and cellular system is perturbed and reiteratively sampled over both time and dose to compile a more comprehensive picture of nanomaterial-cellular interactions.
from over  <dig> papers which were reviewed by the authors in  <cit> , only  <dig> papers have dealt with the effect of concentration plus time while the remaining papers dealt only with the parameter of time.
initial studies that were published previously focused on reporting significantly-expressed genes and using clustering methods to identify similarities and differences between expression profiles  <cit> .
in addition, the study cited in  <cit>  was the only previously combined study investigating time-course and dosage-effect simultaneously, while the initial  <dig> cited studies investigated time-effect and dosage effect separately.
in the present study, we propose different approaches to this kind of analysis.
we considered the gene expression of primary human epidermal keratinocytes, under exposure to the following low-micron to nano-scale materials: carbonyl iron , carbon black  silica  and single-walled carbon nanotubes , at noncytotoxic and cytotoxic doses for each.
the nanomaterials used, except for the carbonyl iron  and swnt, are not intended for medical use.
these materials are currently being used in construction materials, consumer goods, and communications and it applications.
the fc nanoparticles were used as a negative control compound and have been approved for use by the fda as a pharmaceutical carrier formulation.
the single-walled carbon nanotubes are being developed for medical applications  only after being functionalized with other components.
we remark here that the goal of the experimental design was not to study particle size or penetration effects.
it was to study whether there was an overall interaction with the nanomaterials.
in particular, the cytotoxic dose  used with the carbonyl iron was due to its toxic effect of an overload of iron on the cells.
in the approach discussed here, we first used heat maps and multi-dimensional scaling  plots to visualize the dose response of nanomaterials over time.
then, in order to find out the most common patterns in gene-expression profiles, we used self-organizing maps  combined with two different criteria to determine the number of clusters.
the consistency of the som results is discussed in context of the results from the mds plots.
finally, in order to identify the genes that have significantly different responses among the different dose levels while simultaneously accounting for the effect of time, we have used a two-way anova model, in connection with tukey's additivity test and the box-cox transformation.
these analysis results are discussed in the context of the biological and cellular interactions of the engineered nanomaterials.
methods
biological experimental design
primary human epidermal keratinocytes  were cultured with serum-free defined media at 5% co <dig> and 37°c.
the hek were exposed to the following low-micron to nano-scale materials: carbonyl iron , carbon black  silica  and single-walled carbon nanotubes .
fc was supplied by isp technologies with a mean particle size of  <dig>  μm, cb was supplied by degussa with a mean particle size of  <dig> nm, and sio <dig> was supplied by us silica with a mean particle size of  <dig>  μm.
swnt was prepared by southwest nanotechnologies with a mean diameter of  <dig>  nm.
the purity of swnt was checked, and it was found that all heavy metal contamination was very low, less than 1%.
the nanotubes were acid-purified from the raw manufactured intermediate.
the average length of the swnt nanotubes was  <dig> nm; this is the only nanomaterial among the ones considered here that is not spherical.
from previous viability assays  <cit> , two different doses were selected for the exposure of each nanomaterial – a noncytotoxic dose and a cytotoxic dose.
the cytotoxic dose was extrapolated from the viability curve of each substance where 50% of the cells were still viable.
the noncytotoxic and cytotoxic doses for each nanomaterial were found to be, respectively,  <dig>  and  <dig> ,  <dig>  and  <dig>  ,  <dig>  and  <dig> , and  <dig>  and  <dig> , all concentrations being expressed in mg/ml.
the cells were exposed to fc, cb, sio <dig>  and swnt and harvested at  <dig>   <dig>   <dig>   <dig>   <dig>   <dig>  and  <dig> hours after exposure.
cells cultured under the same conditions and not exposed to any of the nanomaterials were harvested at the same time points for a time-matched control baseline.
the cells were trysinized and cell counts taken.
the cells were then collected by centrifugation, snap frozen in liquid nitrogen, and stored at -80°c.
the frozen pellets were shipped to genus biosystems for isolation of the total rna and processing of the gene expression microarrays.
total rna was isolated from the cell pellets, reverse-transcribed to biotinylated crnas and hybridized onto  <dig> k human gene expression microarrays .
the corresponding crna for each biological sample was hybridized to triplicate microarrays.
the arrays were rinsed, dried, scanned and image analyzed and the flat files returned to the houston advanced research center.
data analysis and preprocessing
the microarray flat files contained the quantitative expression values for all probes  and all discovery probe values were assessed with a quality flag.
only those probes that have "discovery" and "g"  tags across all time points and dose levels of all treatments were considered for further analysis to guarantee maximal reproducibility of the results.
by using this criterion, the number of genes was reduced from  <dig>  to  <dig> .
then the average of the three replicate measurements was considered as the actual expression level of each gene at that time in its corresponding treatment.
heatmaps
the heatmaps are visualizations of hierarchical clustering  <cit>  using "average" linkage and pearson correlation as the distance metric between expression levels of each gene after taking the log ratio of the data.
it is worthy to mention that pearson correlation is in fact the normalized correlation between values of two random variables that have had their mean subtracted from them.
the dendrogram added to the left side of each heat map is obtained by pearson correlation.
the heatmaps were drawn by first subtracting the mean of each row from each data cell and then normalizing each row to obtain a variance of  <dig>  note that normalization was performed after taking the log ratio of the expression of each gene between the given treatment and the control baseline.
heatmaps obtained in this way reveal the underlying shape of expression pattern of genes better and are less affected by some large values which may be the result of noise.
for more information on applying hierarchical clustering to gene expression patterns, see  <cit>  and  <cit>  and the references mentioned there.
multi-dimensional scaling
multidimensional scaling  is a method to project high dimensional data onto a lower dimension while maintaining the approximate distances between data points  <cit> .
the accuracy of the representation is measured by the "stress."
a stress of 10% is usually considered good.
first, the data are collected by taking the log of the expression values of genes in each treatment and the baseline.
then at each time point, a matrix of distances between the gene-expression profiles corresponding to each treatment and the baseline is created, where the distance is calculated as  <dig> – the pearson correlation between the profiles.
finally, mds is performed based on this matrix, in order to obtain a low-dimensional representation of the gene-expression profiles and the distance among them.
self-organizing maps
som  <cit>  is a powerful method for clustering.
after log-normalizing the data as explained above, this method was employed in two ways.
first, the underlying hidden number of clusters of the data was estimated using two different criteria the silhouette criterion  <cit>  and the clest criterion  <cit> .
in  <cit> , it has been shown that the clest criterion is the stronger tool to find number of clusters.
however, from the comparison of the clest and the silhouette criteria in  <cit> , it is inferred that whenever the true number of clusters has been equal to  <dig>  the silhouette criterion is the same as the clest or has done a better job than the clest criterion, while for greater number of clusters, the clest performs better.
therefore, the criterion that we have used to determine the number of clusters is a combination of the two criteria.
first, the silhouette criterion is implemented to determine the number of clusters.
if the number of clusters for any treatment is equal to  <dig>  and if the silhouette width for this number of clusters is fairly close to  <dig>  then we assigned  <dig> as the optimal number of clusters.
otherwise, the clest criterion is implemented to determine the number of clusters.
in the second way of implementing som, larger numbers of clusters than the previous approach were tried, as a means of obtaining clusters that show tighter expression patterns than in the first approach.
anova model
a two-way anova model was used to identify genes for which there are significant differences in the response to the various dose levels, and have them ranked by the p-values of dose main effect.
the general procedure was as follows:
<dig> – for each gene, compile the data in a two-dimensional array consisting of the expression values from the three levels of dose  and the eight time points.
use tukey's additivity test for two-way anova  to see if there is interaction between time and dose for the given gene.
if there is no interaction go to step  <dig>
<dig> – if there is significant interaction between these two factors, try to remove the interaction by using the box-cox transformation  <cit> .
the removal of interaction can be tested by tukey's additivity test for transformed values corresponding to each parameter of the box-cox transformation.
we have chosen this parameter to vary between - <dig> to  <dig> with  <dig>  spacing.
if the interaction is removed, use these as the new data for the given gene.
if the interaction could not be removed, set the p-value of dose main effect to na and go to step  <dig> to start the process for the next gene.
<dig> – use two-way anova to find whether there are significant differences among the levels of the dose and time factor.
if there is not any significant difference on any factor, go to step  <dig>
<dig> – use tukey's hsd post-hoc test to determine between which levels of either factor there are significant differences and record them.
since tukey's hsd test is conservative, in rare cases the test cannot find any significant difference between levels.
<dig> – save the p-value for the difference in dose main effect for this gene and go to step  <dig> to start the process for the next gene.
at the end of this process, sort all genes with respect to p-values on dose main effect.
the lower the p-value is for a gene, the stronger the evidence is that there are significant differences in response among the various levels of dose for that gene.
RESULTS
a. heatmap plots
for the heatmaps, we consider the log ratio between the gene expression in one of the treatment time-courses and the corresponding gene expression in the baseline time-course, for each gene and each time point, and plot the log ratio for each gene  versus time point .
there are  <dig> heatmaps, corresponding to each of the  <dig> treatments , which are displayed in figures  <dig> and  <dig>  in these plots, the closer the color of a cell is to white  and red , the more difference there is between the expression value of that gene between the given treatment and the baseline at the particular time point; by contrast, orange indicates an absence of differential expression with respect to the baseline.
considering, for instance, the sio <dig> ld  heatmap, we can see that differential expression tends to occur early in the time course  and tends to disappear with time.
it can also be observed in this heatmap that the genes fall in basically two clusters with respect to their time-course expression: a majority of genes that are up-regulated early, and a minority that is down-regulated early and slightly up-regulated later in the time course.
b. multi-dimensional scaling  plots
for the mds plots, we consider the log of the absolute gene expression for each of the treatments and the baseline .
there are  <dig> plots, one for each time point, which are displayed in figures  <dig> and  <dig>  in these plots, the closer two points are to each other, the more similar their gene expression profiles are, according to the pearson correlation .
considering, for instance, the plot at time  <dig> hrs, it can be easily seen that most of the treatments show similar gene expression to the baseline at this time , while the cb hd  and fc hd treatments present a stronger response.
this is clear by the fact that all the treatments are concentrated in one point except cb hd and fc hd.
in fact, one observes that all treatments tend to become closer to the baseline as time proceeds, except for cb hd and fc hd.
this shows that these two treatments, generally, have the most effect on changing gene expression.
the same fact can also be observed in the heatmaps of figures  <dig> and  <dig>
c. self-organizing maps
here, we again consider the log ratio between the gene expression in one of the treatment time-courses and the corresponding gene expression in the baseline time-course, for each gene and each time point.
table  <dig> shows the optimal number of clusters for each treatment that were obtained by using a combination of the silhouette and clest criteria, as described in the methods section.
it is interesting to compare these numbers with what is deduced from visualizing the heatmaps in figures  <dig> and  <dig>  for example, table  <dig> says that the number of clusters in sio <dig> ld and cb hd is  <dig> and  <dig>  respectively.
this fact is consistent with what one may infer from figures  <dig> and  <dig>
d. anova model
the two-way anova model that we consider has as factors time and dose, where the former has  <dig> levels  and the latter has  <dig> levels, corresponding to the seven time points in the experimental design.
since we are interested in the dose main effect on a given gene, we attempt to remove interaction of dose with time.
after determining those genes affected significantly by the dose main effect, it is worthy to investigate between which specific levels there is a difference, which is accomplished by the tukey's hsd post-hoc test .
discussion
we comment here on the results obtained for the various treatments.
the analysis resulted in several gene families or genes associated with the same cellular processes which responded significantly to the treatments.
very little toxicity overall was seen with any nanomaterial at the noncytotoxic dose except for silica, the positive control substance.
pathways for inflammation and irritation were identified; e.g., these pathways were observed for the silica at a noncytotoxic dose.
at the cytotoxic doses, the most toxicity was seen with carbonyl iron – the explanation for this being the iron overload experienced by the cells.
on the other hand, the majority of the genes found show significant differences at both the low and high dose levels.
we observed that gene expression level becomes closer to basal level in later time points.
this trend is seen with all gene expression profiling where the initial hours show the most activity  <cit> .
next we provide next detailed comments on the genes displayed in table  <dig>  which are associated with significant response to the cb treatment.
many of these genes are notable.
three different chromodomain-helicase-dna binding protein genes are seen: chd <dig>  chd <dig> and chd <dig>  these genes code for proteins which form a histone deacetylase complex involved in chromatin remodeling, an important transcriptional process  <cit> .
two different genes coding for proteins in the protein tyrosine phosphatase or protein tyrosine phosphatase-like families, ptp1b and ptplb, are found to be highly significant.
in addition, four genes associated with apoptosis are found.
il1a and tradd are both involved in the nf-kappab signaling pathways  <cit> .
nr4a <dig> is an early response gene which codes for a nuclear orphan receptor and has an active role in cellular stress  <cit> .
taf <dig> is a tata box-binding protein associated factor binds histones and is important in transcriptional regulation and cell viability  <cit> .
additionally, another protease, psme <dig>  which is related to the ubiquitin pathway, was found  <cit> .
stx1a is a gene that may be involved parkinson's disease as it translates to proteases located in brain tissue  <cit> .
the gene paip <dig> is involved in the regulation of translation.
paip <dig> is very similar to eif4a  which has a role in binding mrna molecules to the ribosome.
paip <dig> is active in mrna turnover and their deadenylation and decay  <cit> .
one gene, frag <dig>  appears twice as two different frag <dig> probes were found to be significant within the top  <dig> genes.
genes that show differential expression at both dose levels include mapk8ip <dig>  trio, atp6v1f, and hs.
<dig> .
trio regulates the reorganization of the cytoskeleton and is essential for cell migration and cell growth  <cit> .
atp6v1f is active in transport mechanisms; it mediates the acidification of intracellular compartments essential to ion transport across cellular membranes  <cit> .
mapk8ip <dig> is involved in signal transduction by coding for a scaffold protein in the jnk signaling pathway  <cit> .
it is activated by cytokines and environmental stresses.
one gene  is unannotated, being a transcribed locus that translates to a protein with no known functions.
several genes show differential expression only with the high dose.
these genes include il1a, lamb <dig>  impa <dig>  arid5b, chd <dig> and arhgap <dig>  the genes arhgap <dig>  arid5b and chd <dig> are involved in the regulation of transcription.
arhgap <dig> codes for a rac gptase active in signal transduction and seems to play an indirect role in cytoskeletal reorganization and cell motility  <cit> .
arid5b regulates transcription by binding to dna and has important roles in development and differentiation  <cit> .
chd <dig> represses transcription by modulating the amount of protein at the centrosome during mitosis  <cit> .
il1a is involved in apoptosis as well as many roles in immune and inflammatory responses and hematopoiesis  <cit> .
the last two genes, lamb <dig> and impa <dig> are both involved in development.
impa <dig> has an essential role in phosphatidylinositol pathways and lamb <dig> codes for a laminin isoform which is involved in embryonic development and the organization of cells into tissues  <cit> .
CONCLUSIONS
the analysis presented here lead to interesting and complementary conclusions about the response across time of human epidermal keratinocytes after exposure to nanomaterials.
for example, we observed  that gene expression for each treatment become closer to the expression of the baseline cultures as time proceeds, except for cb hd and fc hd, which generally have the most changing expression levels.
the results from the mds analysis mostly confirmed the results from the heatmaps and som clustering.
the genes found to be differentially-expressed by the anova analysis are involved in a number of cellular processes, including regulation of transcription and translation, protein localization, transport, cell cycle progression, cell migration, cytoskeletal reorganization, signal transduction, and development.
the majority of these genes show significant differences at both the low and high dose levels.
supplementary information
http://www.ece.tamu.edu/~ulisses/harc/index.html
competing interests
the authors declare that they have no competing interests.
authors' contributions
az: bioinformatics, drafted the manuscript; for mjc: cell culture, coordinated rna isolation and microarray processing, conceived the study, helped draft the manuscript; ubn: bioinformatics, helped draft the manuscript; erd: participated in design and coordination and helped draft the manuscript.

 nanomaterials are being manufactured on a commercial scale for use in medical, diagnostic, energy, component and communications industries. engineered nanomaterials range considerably in their physicochemical properties making them more desirable than their micro- and macro-counterparts due to, for example, their increased surface area, tensile strength, tunability, etc.. from limited early reports, concerns over the safety of engineered nanomaterials have surfaced. humans can be exposed to nanomaterials in different ways such as inhalation or exposure through the integumentary system. however, the skin is a unique organ in the body not only because it gives the body such a large surface area for exposure but also because of the avascular property of epidermis, in which particles can reside without being removed by phagocytosis. gene expression microarrays have become a tool to investigate the interactions of biological systems by observing the simultaneous activities of tens of thousands of genes. over recent years, this tool has been applied to toxicology forming a new discipline, toxicogenomics. microarrays have most recently been a tool used by pharmaceutical drug discovery and development to screen for efficacy and adverse effects thereby prioritizing drug candidates and redeveloping ones which show off-target toxicities. the approach described here combines a global screening technology, gene expression microarray profiling, with systems biology, to investigate the interactions of engineered nanomaterials with primary human cells. the biological and cellular system is perturbed and reiteratively sampled over both time and dose to compile a more comprehensive picture of nanomaterial-cellular interactions. from over papers which were reviewed by the authors in, only papers have dealt with the effect of concentration plus time while the remaining papers dealt only with the parameter of time. initial studies that were published previously focused on reporting significantly-expressed genes and using clustering methods to identify similarities and differences between expression profiles. in addition, the study cited in was the only previously combined study investigating time-course and dosage-effect simultaneously, while the initial cited studies investigated time-effect and dosage effect separately. in the present study, we propose different approaches to this kind of analysis. we considered the gene expression of primary human epidermal keratinocytes, under exposure to the following low-micron to nano-scale materials: carbonyl iron, carbon black silica and single-walled carbon nanotubes, at noncytotoxic and cytotoxic doses for each. the nanomaterials used, except for the carbonyl iron and swnt, are not intended for medical use. these materials are currently being used in construction materials, consumer goods, and communications and it applications. the fc nanoparticles were used as a negative control compound and have been approved for use by the fda as a pharmaceutical carrier formulation. the single-walled carbon nanotubes are being developed for medical applications only after being functionalized with other components. we remark here that the goal of the experimental design was not to study particle size or penetration effects. it was to study whether there was an overall interaction with the nanomaterials. in particular, the cytotoxic dose used with the carbonyl iron was due to its toxic effect of an overload of iron on the cells. in the approach discussed here, we first used heat maps and multi-dimensional scaling plots to visualize the dose response of nanomaterials over time. then, in order to find out the most common patterns in gene-expression profiles, we used self-organizing maps combined with two different criteria to determine the number of clusters. the consistency of the som results is discussed in context of the results from the mds plots. finally, in order to identify the genes that have significantly different responses among the different dose levels while simultaneously accounting for the effect of time, we have used a two-way anova model, in connection with tukey's additivity test and the box-cox transformation. these analysis results are discussed in the context of the biological and cellular interactions of the engineered nanomaterials. methods biological experimental design primary human epidermal keratinocytes were cultured with serum-free defined media at 5% co and 37°c. the hek were exposed to the following low-micron to nano-scale materials: carbonyl iron, carbon black silica and single-walled carbon nanotubes. fc was supplied by isp technologies with a mean particle size of μm, cb was supplied by degussa with a mean particle size of nm, and sio was supplied by us silica with a mean particle size of μm. swnt was prepared by southwest nanotechnologies with a mean diameter of nm. the purity of swnt was checked, and it was found that all heavy metal contamination was very low, less than 1%. the nanotubes were acid-purified from the raw manufactured intermediate. the average length of the swnt nanotubes was nm; this is the only nanomaterial among the ones considered here that is not spherical. from previous viability assays, two different doses were selected for the exposure of each nanomaterial – a noncytotoxic dose and a cytotoxic dose. the cytotoxic dose was extrapolated from the viability curve of each substance where 50% of the cells were still viable. the noncytotoxic and cytotoxic doses for each nanomaterial were found to be, respectively, and, and, and, and and, all concentrations being expressed in mg/ml. the cells were exposed to fc, cb, sio and swnt and harvested at and hours after exposure. cells cultured under the same conditions and not exposed to any of the nanomaterials were harvested at the same time points for a time-matched control baseline. the cells were trysinized and cell counts taken. the cells were then collected by centrifugation, snap frozen in liquid nitrogen, and stored at -80°c. the frozen pellets were shipped to genus biosystems for isolation of the total rna and processing of the gene expression microarrays. total rna was isolated from the cell pellets, reverse-transcribed to biotinylated crnas and hybridized onto k human gene expression microarrays. the corresponding crna for each biological sample was hybridized to triplicate microarrays. the arrays were rinsed, dried, scanned and image analyzed and the flat files returned to the houston advanced research center. data analysis and preprocessing the microarray flat files contained the quantitative expression values for all probes and all discovery probe values were assessed with a quality flag. only those probes that have "discovery" and "g" tags across all time points and dose levels of all treatments were considered for further analysis to guarantee maximal reproducibility of the results. by using this criterion, the number of genes was reduced from to. then the average of the three replicate measurements was considered as the actual expression level of each gene at that time in its corresponding treatment. heatmaps the heatmaps are visualizations of hierarchical clustering using "average" linkage and pearson correlation as the distance metric between expression levels of each gene after taking the log ratio of the data. it is worthy to mention that pearson correlation is in fact the normalized correlation between values of two random variables that have had their mean subtracted from them. the dendrogram added to the left side of each heat map is obtained by pearson correlation. the heatmaps were drawn by first subtracting the mean of each row from each data cell and then normalizing each row to obtain a variance of note that normalization was performed after taking the log ratio of the expression of each gene between the given treatment and the control baseline. heatmaps obtained in this way reveal the underlying shape of expression pattern of genes better and are less affected by some large values which may be the result of noise. for more information on applying hierarchical clustering to gene expression patterns, see and and the references mentioned there. multi-dimensional scaling multidimensional scaling is a method to project high dimensional data onto a lower dimension while maintaining the approximate distances between data points. the accuracy of the representation is measured by the "stress." a stress of 10% is usually considered good. first, the data are collected by taking the log of the expression values of genes in each treatment and the baseline. then at each time point, a matrix of distances between the gene-expression profiles corresponding to each treatment and the baseline is created, where the distance is calculated as – the pearson correlation between the profiles. finally, mds is performed based on this matrix, in order to obtain a low-dimensional representation of the gene-expression profiles and the distance among them. self-organizing maps som is a powerful method for clustering. after log-normalizing the data as explained above, this method was employed in two ways. first, the underlying hidden number of clusters of the data was estimated using two different criteria the silhouette criterion and the clest criterion. in, it has been shown that the clest criterion is the stronger tool to find number of clusters. however, from the comparison of the clest and the silhouette criteria in, it is inferred that whenever the true number of clusters has been equal to the silhouette criterion is the same as the clest or has done a better job than the clest criterion, while for greater number of clusters, the clest performs better. therefore, the criterion that we have used to determine the number of clusters is a combination of the two criteria. first, the silhouette criterion is implemented to determine the number of clusters. if the number of clusters for any treatment is equal to and if the silhouette width for this number of clusters is fairly close to then we assigned as the optimal number of clusters. otherwise, the clest criterion is implemented to determine the number of clusters. in the second way of implementing som, larger numbers of clusters than the previous approach were tried, as a means of obtaining clusters that show tighter expression patterns than in the first approach. anova model a two-way anova model was used to identify genes for which there are significant differences in the response to the various dose levels, and have them ranked by the p-values of dose main effect. the general procedure was as follows: – for each gene, compile the data in a two-dimensional array consisting of the expression values from the three levels of dose and the eight time points. use tukey's additivity test for two-way anova to see if there is interaction between time and dose for the given gene. if there is no interaction go to step – if there is significant interaction between these two factors, try to remove the interaction by using the box-cox transformation. the removal of interaction can be tested by tukey's additivity test for transformed values corresponding to each parameter of the box-cox transformation. we have chosen this parameter to vary between - to with spacing. if the interaction is removed, use these as the new data for the given gene. if the interaction could not be removed, set the p-value of dose main effect to na and go to step to start the process for the next gene. – use two-way anova to find whether there are significant differences among the levels of the dose and time factor. if there is not any significant difference on any factor, go to step – use tukey's hsd post-hoc test to determine between which levels of either factor there are significant differences and record them. since tukey's hsd test is conservative, in rare cases the test cannot find any significant difference between levels. – save the p-value for the difference in dose main effect for this gene and go to step to start the process for the next gene. at the end of this process, sort all genes with respect to p-values on dose main effect. the lower the p-value is for a gene, the stronger the evidence is that there are significant differences in response among the various levels of dose for that gene.  a. heatmap plots for the heatmaps, we consider the log ratio between the gene expression in one of the treatment time-courses and the corresponding gene expression in the baseline time-course, for each gene and each time point, and plot the log ratio for each gene versus time point. there are heatmaps, corresponding to each of the treatments, which are displayed in figures and in these plots, the closer the color of a cell is to white and red, the more difference there is between the expression value of that gene between the given treatment and the baseline at the particular time point; by contrast, orange indicates an absence of differential expression with respect to the baseline. considering, for instance, the sio ld heatmap, we can see that differential expression tends to occur early in the time course and tends to disappear with time. it can also be observed in this heatmap that the genes fall in basically two clusters with respect to their time-course expression: a majority of genes that are up-regulated early, and a minority that is down-regulated early and slightly up-regulated later in the time course. b. multi-dimensional scaling plots for the mds plots, we consider the log of the absolute gene expression for each of the treatments and the baseline. there are plots, one for each time point, which are displayed in figures and in these plots, the closer two points are to each other, the more similar their gene expression profiles are, according to the pearson correlation. considering, for instance, the plot at time hrs, it can be easily seen that most of the treatments show similar gene expression to the baseline at this time, while the cb hd and fc hd treatments present a stronger response. this is clear by the fact that all the treatments are concentrated in one point except cb hd and fc hd. in fact, one observes that all treatments tend to become closer to the baseline as time proceeds, except for cb hd and fc hd. this shows that these two treatments, generally, have the most effect on changing gene expression. the same fact can also be observed in the heatmaps of figures and c. self-organizing maps here, we again consider the log ratio between the gene expression in one of the treatment time-courses and the corresponding gene expression in the baseline time-course, for each gene and each time point. table shows the optimal number of clusters for each treatment that were obtained by using a combination of the silhouette and clest criteria, as described in the methods section. it is interesting to compare these numbers with what is deduced from visualizing the heatmaps in figures and for example, table says that the number of clusters in sio ld and cb hd is and respectively. this fact is consistent with what one may infer from figures and d. anova model the two-way anova model that we consider has as factors time and dose, where the former has levels and the latter has levels, corresponding to the seven time points in the experimental design. since we are interested in the dose main effect on a given gene, we attempt to remove interaction of dose with time. after determining those genes affected significantly by the dose main effect, it is worthy to investigate between which specific levels there is a difference, which is accomplished by the tukey's hsd post-hoc test. discussion we comment here on the results obtained for the various treatments. the analysis resulted in several gene families or genes associated with the same cellular processes which responded significantly to the treatments. very little toxicity overall was seen with any nanomaterial at the noncytotoxic dose except for silica, the positive control substance. pathways for inflammation and irritation were identified; e.g., these pathways were observed for the silica at a noncytotoxic dose. at the cytotoxic doses, the most toxicity was seen with carbonyl iron – the explanation for this being the iron overload experienced by the cells. on the other hand, the majority of the genes found show significant differences at both the low and high dose levels. we observed that gene expression level becomes closer to basal level in later time points. this trend is seen with all gene expression profiling where the initial hours show the most activity. next we provide next detailed comments on the genes displayed in table which are associated with significant response to the cb treatment. many of these genes are notable. three different chromodomain-helicase-dna binding protein genes are seen: chd chd and chd these genes code for proteins which form a histone deacetylase complex involved in chromatin remodeling, an important transcriptional process. two different genes coding for proteins in the protein tyrosine phosphatase or protein tyrosine phosphatase-like families, ptp1b and ptplb, are found to be highly significant. in addition, four genes associated with apoptosis are found. il1a and tradd are both involved in the nf-kappab signaling pathways. nr4a is an early response gene which codes for a nuclear orphan receptor and has an active role in cellular stress. taf is a tata box-binding protein associated factor binds histones and is important in transcriptional regulation and cell viability. additionally, another protease, psme which is related to the ubiquitin pathway, was found. stx1a is a gene that may be involved parkinson's disease as it translates to proteases located in brain tissue. the gene paip is involved in the regulation of translation. paip is very similar to eif4a which has a role in binding mrna molecules to the ribosome. paip is active in mrna turnover and their deadenylation and decay. one gene, frag appears twice as two different frag probes were found to be significant within the top genes. genes that show differential expression at both dose levels include mapk8ip trio, atp6v1f, and hs.. trio regulates the reorganization of the cytoskeleton and is essential for cell migration and cell growth. atp6v1f is active in transport mechanisms; it mediates the acidification of intracellular compartments essential to ion transport across cellular membranes. mapk8ip is involved in signal transduction by coding for a scaffold protein in the jnk signaling pathway. it is activated by cytokines and environmental stresses. one gene is unannotated, being a transcribed locus that translates to a protein with no known functions. several genes show differential expression only with the high dose. these genes include il1a, lamb impa arid5b, chd and arhgap the genes arhgap arid5b and chd are involved in the regulation of transcription. arhgap codes for a rac gptase active in signal transduction and seems to play an indirect role in cytoskeletal reorganization and cell motility. arid5b regulates transcription by binding to dna and has important roles in development and differentiation. chd represses transcription by modulating the amount of protein at the centrosome during mitosis. il1a is involved in apoptosis as well as many roles in immune and inflammatory responses and hematopoiesis. the last two genes, lamb and impa are both involved in development. impa has an essential role in phosphatidylinositol pathways and lamb codes for a laminin isoform which is involved in embryonic development and the organization of cells into tissues.  the analysis presented here lead to interesting and complementary conclusions about the response across time of human epidermal keratinocytes after exposure to nanomaterials. for example, we observed that gene expression for each treatment become closer to the expression of the baseline cultures as time proceeds, except for cb hd and fc hd, which generally have the most changing expression levels. the results from the mds analysis mostly confirmed the results from the heatmaps and som clustering. the genes found to be differentially-expressed by the anova analysis are involved in a number of cellular processes, including regulation of transcription and translation, protein localization, transport, cell cycle progression, cell migration, cytoskeletal reorganization, signal transduction, and development. the majority of these genes show significant differences at both the low and high dose levels. supplementary information http://www.ece.tamu.edu/~ulisses/harc/index.html competing interests the authors declare that they have no competing interests. authors' contributions az: bioinformatics, drafted the manuscript; for mjc: cell culture, coordinated rna isolation and microarray processing, conceived the study, helped draft the manuscript; ubn: bioinformatics, helped draft the manuscript; erd: participated in design and coordination and helped draft the manuscript. 
then, in order to find out the most common patterns in gene-expression profiles, we used self-organizing maps combined with two different criteria to determine the number of clusters. the genes found to be differentially-expressed are involved in a number of cellular processes, including regulation of transcription and translation, protein localization, transport, cell cycle progression, cell migration, cytoskeletal reorganization, signal transduction, and development. finally, in order to identify the genes that have significantly different responses among different levels of dose of each treatment while accounting for the effect of time at the same time, we used a two-way anova model, in connection with tukey's additivity test and the box-cox transformation. the consistency of som results is discussed in context of the information derived from the mds plots. the results are discussed in the context of the cellular responses of engineered nanomaterials. for example, we observed that gene expression for most treatments become closer to the expression of the baseline cultures as time proceeds. the analysis presented here lead to interesting and complementary conclusions about the response across time of human epidermal keratinocytes after exposure to nanomaterials. the interactions of engineered nanomaterials with primary human cells was investigated, using a systems biology approach combining gene expression microarray profiling with dynamic experimental parameters. first, we used heat maps and multi-dimensional scaling plots to visualize the dose response of nanomaterials over time. here, we propose different approaches to this kind of analysis. in this experiment, primary human epidermal keratinocytes cells were exposed to several low-micron to nano-scale materials, and gene expression was profiled over both time and dose to compile a comprehensive picture of nanomaterial-cellular interactions. nanomaterials are being manufactured on a commercial scale for use in medical, diagnostic, energy, component and communications industries. very few gene-expression studies so far have dealt with both time and dose response simultaneously. humans can be exposed to nanomaterials in different ways such as inhalation or exposure through the integumentary system.
 nanomaterials are being manufactured on a commercial scale for use in medical, diagnostic, energy, component and communications industries. however, concerns over the safety of engineered nanomaterials have surfaced. humans can be exposed to nanomaterials in different ways such as inhalation or exposure through the integumentary system.  the interactions of engineered nanomaterials with primary human cells was investigated, using a systems biology approach combining gene expression microarray profiling with dynamic experimental parameters. in this experiment, primary human epidermal keratinocytes cells were exposed to several low-micron to nano-scale materials, and gene expression was profiled over both time and dose to compile a comprehensive picture of nanomaterial-cellular interactions. very few gene-expression studies so far have dealt with both time and dose response simultaneously. here, we propose different approaches to this kind of analysis. first, we used heat maps and multi-dimensional scaling plots to visualize the dose response of nanomaterials over time. then, in order to find out the most common patterns in gene-expression profiles, we used self-organizing maps combined with two different criteria to determine the number of clusters. the consistency of som results is discussed in context of the information derived from the mds plots. finally, in order to identify the genes that have significantly different responses among different levels of dose of each treatment while accounting for the effect of time at the same time, we used a two-way anova model, in connection with tukey's additivity test and the box-cox transformation. the results are discussed in the context of the cellular responses of engineered nanomaterials.  the analysis presented here lead to interesting and complementary conclusions about the response across time of human epidermal keratinocytes after exposure to nanomaterials. for example, we observed that gene expression for most treatments become closer to the expression of the baseline cultures as time proceeds. the genes found to be differentially-expressed are involved in a number of cellular processes, including regulation of transcription and translation, protein localization, transport, cell cycle progression, cell migration, cytoskeletal reorganization, signal transduction, and development. 20- february sixth annual mcbios conference. transformational bioinformatics: delivering value from genomes starkville, ms, usa 