Variability in synovial inflammation in rheumatoid arthritis investigated by microarray technology

doi:10.1186/ar1903

Arthritis Res Ther. 2006; 8(2): R47.

Published online 2006 February 16. doi: 10.1186/ar1903

PMCID: PMC1526587

Variability in synovial inflammation in rheumatoid arthritis investigated by microarray technology

Johan Lindberg,¹ Erik af Klint,² Ann-Kristin Ulfgren,² André Stark,³ Tove Andersson,¹ Peter Nilsson,¹ Lars Klareskog,² and Joakim Lundeberg¹

¹Department of Biotechnology, AlbaNova University Center, Royal Institute of Technology, S-106 91 Stockholm, Sweden

²Department of Rheumatology, Karolinska Institutet, Karolinska University Hospital, 171 76 Stockholm, Sweden

³Department of Orthopedics, Karolinska University Hospital, 171 76 Stockholm, Sweden

Corresponding author.

Johan Lindberg: johanl/at/biotech.kth.se; Erik af Klint: erik.af.klint/at/medks.ki.se; Ann-Kristin Ulfgren: ann-Kristin.Ulfgren/at/cmm.ki.se; André Stark: andreas.stark/at/karolinska.se; Tove Andersson: tove/at/jimmy.harvard.edu; Peter Nilsson: nipe/at/biotech.kth.se; Lars Klareskog: lars.klareskog/at/medks.ki.se; Joakim Lundeberg: joakim.lundeberg/at/biotech.kth.se

Received July 3, 2005; Revisions requested September 9, 2005; Revised November 18, 2005; Accepted January 23, 2006.

This is an open access article distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/2.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

This article has been cited by other articles in PMC.

Abstract

In recent years microarray technology has been used increasingly to acquire knowledge about the pathogenic processes involved in rheumatoid arthritis. The present study investigated variations in gene expression in synovial tissues within and between patients with rheumatoid arthritis. This was done by applying microarray technology on multiple synovial biopsies obtained from the same knee joints. In this way the relative levels of intra-patient and inter-patient variation could be assessed. The biopsies were obtained from 13 different patients: 7 by orthopedic surgery and 6 by rheumatic arthroscopy. The data show that levels of heterogeneity varied substantially between the biopsies, because the number of genes found to be differentially expressed between pairs of biopsies from the same knee ranged from 6 to 2,133. Both arthroscopic and orthopedic biopsies were examined, allowing us to compare the two sampling methods. We found that the average number of differentially expressed genes between biopsies from the same patient was about three times larger in orthopedic than in arthroscopic biopsies. Using a parallel analysis of the tissues by immunohistochemistry, we also identified orthopedic biopsies that were unsuitable for gene expression analysis of synovial inflammation due to sampling of non-inflamed parts of the tissue. Removing these biopsies reduced the average number of differentially expressed genes between the orthopedic biopsies from 455 to 171, in comparison with 143 for the arthroscopic biopsies. Hierarchical clustering analysis showed that the remaining orthopedic and arthroscopic biopsies had gene expression signatures that were unique for each patient, apparently reflecting patient variation rather than tissue heterogeneity. Subsets of genes found to vary between biopsies were investigated for overrepresentation of biological processes by using gene ontology. This revealed representative 'themes' likely to vary between synovial biopsies affected by inflammatory disease.

Introduction

Rheumatoid arthritis (RA) is a common chronic inflammatory disease, so far defined by a set of criteria [1] rather than by a knowledge of the underlying molecular pathogenesis. Substantial efforts have been made to characterize the synovial inflammation in RA, and during these studies it has become evident that there is a large variability in cell content and in protein expression, both within single joints and between patients with RA [2-7]. This variation also exists at the gene expression level [8]. Microarray (MA) technology allows the expression of thousands of genes to be monitored simultaneously and can thus increase the understanding of the complicated molecular processes of joint inflammation in more detail than has been possible with immunohistochemistry and related techniques [9-14]. Recently reviewed [15], MA has been used to acquire knowledge about RA in various experimental systems with the use of both cell cultures [16-22] and biopsies [23-30] obtained from the synovium. So far, MA has been used to investigate tissue heterogeneity between synovial biopsies obtained from different patients in both juvenile RA [23] and long-standing RA [25,30]. Tsubaki and colleagues [23] used laser capture microdissection on biopsies retrieved by rheumatic arthroscopy from patients with juvenile RA to characterize proliferative lesions in the synovial lining. Two subgroups were discovered; one had a gene expression profile similar to that of long-standing RA. Van der Pouw Kraan and colleagues [25,30] used MA to investigate heterogeneity between synovial biopsies obtained by orthopedic surgery from different patients. In both of these studies, at least two different gene expression profiles were observed, which were suggested to correspond to high and low inflammatory status.

These and other studies therefore suggest that the MA technique might indeed be able to discern variable molecular features of the joint inflammation that would be both biologically and clinically meaningful. However, further investigation of the potential of these gene expression patterns to predict disease course as well as the response to various therapies is hampered by an incomplete knowledge of the natural variability of gene expression within the inflamed joints of single patients and between different patients with RA. In this study we therefore compared variation in gene expression patterns at the biopsy site, between different sites, and between patients. We used inflamed synovial tissues of patients with RA obtained during open surgery and during rheumatic arthroscopy, which were our methods of choice for synovial tissue retrieval.

Materials and methods

Patients

Thirteen patients, all fulfilling the American College of Rheumatology classification criteria for RA [1], were included in this study. Synovial tissues were taken from seven of these patients with erosive, end-stage disease during knee joint replacement surgery at the Department of Orthopedic Surgery, Karolinska University Hospital, Sweden. No further data on the characteristics of this subgroup of patients were available. Synovial tissue was obtained from the other six patients by rheumatic arthroscopy solely for research purposes. The clinical characteristics of these patients (five women and one man) are shown in Table Table1.1. All six arthroscopic patients were recruited from the outpatient clinic of the Karolinska University Hospital Rheumatology Unit, and all except one (patient 13) had clinical arthritis with effusion in at least one knee joint at the time of the investigation. All patients except one (patient 11) were using the disease-modifying anti-rheumatic drug methotrexate, four in conjunction with low-dose corticosteroids, and all except one were using nonsteroidal anti-inflammatory drugs. Patient 13 had been taking methotrexate for two months; the others had been doing so for more than six months. Patients 10, 11 and 13 had erosive disease. The Ethical Committee at the Karolinska Institute approved the study protocol and all patients gave informed consent.

Table 1

Clinical data for patients 8 to 13

Synovial tissue, sampling and handling

Orthopedic samples (patients 1 to 7)

Knee joint replacement surgery was performed in accordance with standard procedures, during which three synovial tissue specimens were obtained from random sites and immediately handled by research personnel. Each biopsy was visually inspected to minimize non-inflammatory synovial tissue contamination. Each orthopedic biopsy was then split into two parts, one of which was used for the MA experiment; the other was saved for histochemical analysis (except for biopsy 3 of patient 2 which was used only for MA). After dividing the biospies they were snap frozen (within two minutes) in precooled isopentane and stored at -80°C until further use, to ensure high RNA quality. For patients 1 to 3 each half of a biopsy that was to be used in the MA experiment was further divided into three parts, hereafter referred to as sub-biopsies. In total this resulted in 39 specimens from the orthopedic patients (nine sub-biopsies from patients 1 to 3 and three biopsies from each of patients 4 to 7). The average weight of the biopsies from patients 1 to 3 before division into three parts was 99 mg, yielding an average sub-biopsy weight of 33 mg. The average weight of the biopsies from patients 4 to 7 was 29 mg.

Arthroscopic samples (patients 8 to 13)

Rheumatic arthroscopy was performed by a technique previously described [31], including biopsy site-scoring for signs of inflammation (vascularity and proliferation), photography and mapping according to local standards (not shown). Biopsies were taken at the site of inflammation close to cartilage or not close to cartilage, defined as less than or more than 1.5 cm away from cartilage, respectively. Multiple biopsies were taken from two sites in patients 8, 11, 12 and 13 and from four sites in patients 9 and 10. The samples were frozen, as above, within two minutes. The average weight of the arthroscopic samples was 19 mg.

Histochemistry

Frozen biopsies were embedded in Optimal Cutting Temperature (OCT; Tissue-Tek, SAKURA Finetek, Zoeterwoude, Netherlands) and cut with a cryostat into 7 μm thick. Sections were placed on SuperFrost^®Plus slides (Menzel-Gläser, Braunschweig, Germany) and air-dried for 30 minutes, then stained with Mayer's hematoxylin and eosin to confirm the histopathology of each biopsy.

RNA extraction

RNA was successfully extracted from all biopsies except biopsy 3 of patient 7, in which the RNA was degraded. For both types of biopsy (arthroscopic and orthopedic) one biopsy yielded enough RNA to perform MA experiments. To extract the RNA the biopsies were placed in steel-bead matrix tubes (Lysing Matrix D; Qbiogene, Irvine, CA, USA) containing buffer (600 μl of phenol, 600 μl RLT of buffer from an RNeasy kit (Qiagen, Hilden, Germany) and 0.6 μl of 2-mercaptoethanol) and homogenized with a tabletop FastPrep homogenizer (Qbiogene). The tubes were shaken for 30 seconds at speed setting 6 and then put on ice for 30 seconds. This procedure was repeated four times to ensure thorough homogenization. The tubes were then centrifuged for 5 minutes at 12,000 r.p.m. All steps up to this point were performed at 4°C. The water phase was collected and transferred to Qiashredder (Qiagen) columns and centrifuged (13,000 r.p.m. at room temperature) for two minutes to ensure complete homogenization. To the flow-through was added 600 μl of 70% ethanol.

The mixture was loaded onto RNeasy spin columns (Qiagen) and centrifuged for 15 seconds at 13,000 r.p.m. An RNeasy kit from Qiagen was used to wash and elute the extracted RNA (in 30 μl of RNase-free water). Before eluting the RNA the columns were treated with 2 units of DNAse H (Omega Biotech, Victoria, Canada) for 15 minutes at room temperature to remove residual DNA contamination. For further details see 'Preparation of RNA from tissues' at the KTH microarray core facility web site [32], under 'Protocols'. The average biopsy weight used for RNA extraction was 28.7 mg and the average RNA yield was 411.4 ng/μl in 30 μl. All concentration measurements were made with the Nanodrop (Nanodrop Technologies, Wilmington, DE USA). RNA quality was ensured with the RNA 6000 Nano LabChip kit of the Bioanalyzer system (Agilent Technologies, Palo Alto, CA, USA) where pass or fail judgments were based on an evaluation of Bioanalyzer electropherograms [33]. Two samples showed signs of partial degradation (patient 3, biopsy 3, and patient 6, biopsy 1). The average ratio of 28S to 18S rRNA among the remaining samples was 1.6.

RNA amplification

Because of the small amounts of RNA extracted, the RNA was amplified with a RiboAmp RNA amplification kit (Arcturus, Mountain View, CA, USA). RiboAmp uses T7-based in vitro transcription to generate amplified RNA (aRNA), the bulk of which consists of sequences 250 to 1,800 base pairs long. Total RNA (300 ng to 1 μg) was used in each RNA amplification, and the average yield was 503 ng/μl in 11 μl of water.

RNA reference

Universal Human Reference RNA from Stratagene (La Jolla, CA, USA) was used as reference RNA and was amplified in the same manner as the sample RNA. The reference RNA was pooled before use for hybridization.

Labeling and cDNA synthesis

To prime the reaction, 1 μl of random hexamer primer (5 μg/μl; Operon, Alameda CA, USA) was added to 1 μg of amplified aRNA. The volume was adjusted to 18.4 μl with RNase-free water. The sample was mixed and incubated for 10 minutes at 70°C to denature the aRNA, then incubated for a further 5 minutes on ice and centrifuged briefly. A cDNA synthesis mixture (11.6 μl) consisting of 6 μl of 5× first-strand buffer, 3 μl of 0.1 M dithiothreitol, 2 μl of Superscript III (Invitrogen, San Diego, CA, USA) and 0.6 μl of 50× aa-dUTP+dNTP mix (Sigma-Aldrich, St. Louis, MO, USA) was added to each sample. The whole mixture was gently mixed by pipetting and incubated at 25°C. After 10 minutes at this temperature the mixture was incubated at 46°C for a further 2 hours. To terminate the reaction and to hydrolyze the RNA strand, 3 μl of 0.2 M EDTA pH 8.0 and 4.5 μl of 1 M NaOH were added. The sample was vortex-mixed briefly, incubated for 15 minutes at 70°C, cooled to room temperature and centrifuged briefly. Then 4.5 μl of 1 M HCl was added to restore the pH to neutrality. The sample was vortex-mixed and centrifuged briefly.

The cDNA was purified and eluted by the following procedure. First, 60 μl of water and 500 μl of PB buffer (MinElute Reaction Cleanup Kit; Qiagen) were added. The mixture was thoroughly mixed and transferred to a MinElute Reaction Cleanup Kit spin column and centrifuged for 30 seconds at 13,000 r.p.m. The flow-through was reapplied to the column and the centrifugation step was repeated. The flow-through was discarded and 650 μl of 80% ethanol was added to the column. The column was centrifuged for 30 seconds at 13,000 r.p.m. and the flow-through was again discarded. The ethanol wash step was repeated and then the membrane was dried by centrifugation for 1 minute at 13,000 r.p.m. The column was transferred to a new tube, and 10 μl of 100 mM NaHCO₃pH 9.0 was added. The column was then incubated for 1 minute at room temperature; the sample was eluted by centrifugation for 30 seconds at 13,000 r.p.m. The elution step was repeated after a further addition of 10 μl of 100 mM NaHCO₃pH 9.0 to ensure high yield.

To couple fluorophores, the eluate was mixed with a dried aliquot of either Cy3 or Cy5 mono-reactive esters (Amersham-Biosciences, Little Chalfont, Bucks., UK) and incubated for 30 minutes at room temperature in a dark container, after which 70 μl of water and 500 μl of PB buffer were added. The mixture was thoroughly mixed and transferred to a MinElute Reaction Cleanup Kit spin column, which was centrifuged for 30 seconds at 13,000 r.p.m. The flow-through was reapplied to the column and the centrifugation step was repeated. The flow-through was discarded and 650 μl of PE buffer (MinElute Reaction Cleanup Kit) was added. The column was centrifuged for 30 seconds at 13,000 r.p.m. and the flow-through was discarded. The wash step was repeated and then the membrane was washed by centrifugation for 1 minute at 13,000 r.p.m. The column was transferred to a new tube, 10 μl of EB buffer (MinElute Reaction Cleanup Kit) was added, and the column was incubated for 1 minute at room temperature. The sample was eluted by centrifugation for 30 seconds at 13,000 r.p.m. The elution step was repeated after a further addition of 10 μl of EB buffer, to ensure high yield. The concentrations of the incorporated fluorophore and cDNA were measured with the Nanodrop to confirm success in the labeling reaction. The sample was then ready for hybridization. For further information about the preparation of N-hydroxysuccinimide-ester fluorophores and indirect labeling of cDNA see SOP 001 and SOP 002 at the KTH microarray core facility web site [32] under 'Protocols'.

cDNA microarray

The cDNA arrays used in this study were produced at the KTH microarray core facility. The clones on the array originate from the first 310 96-well plates of a commercial clone collection containing 46,000 sequence-verified human cDNA clones (Research Genetics; now Invitrogen). The clones have been prepared by cell culture, plasmid preparation, PCR amplification and purification with PCR filter plates (Millipore, Bedford MA, USA). The cDNA was spotted in 30% dimethylsulphoxide onto UltraGAPS slides (Corning, NY, USA) with a QArray spotter (Genetix, Hampshire, UK) with 24 SMP2.5 pins (Telechem, Sunnyvale, CA, USA) in 48 blocks, each of which contained 25 × 25 clones, spotted with a center-to-center distance of 170 μm, non-specifically attached to the surface by UV crosslinking. According to a UniGene mapping performed in September 2004 based on GenBank accession numbers (29,717 on the whole chip), 25,087 of the 30,000 spots on the chip have a UniGene ID and 16,164 of those were unique. For more information about the chip see the KTH microarray core facility web site [32] under 'HUM 30k cDNA array'.

Hybridization

After prehybridizing of the slides for 30 minutes at 42°C in prehybridization buffer consisting of 1% BSA (Sigma-Aldrich), 5 × SSC (where SSC consists of 0.15 M NaCl and 0.015 M sodium citrate) and 0.1% SDS, the arrays were washed, first in a trough containing water and then in a trough containing propan-2-ol, and centrifuged dry. After that the samples (one labeled with Cy5 and the other labeled with Cy3) were pooled and dried to a volume of 13.6 μl. Hybridization buffer consisting of 5 × SSC, 50% formamide (Sigma-Aldrich), 0.1% SDS and 0.2 μg/μl Cot-1 DNA (Invitrogen) was added to the pooled samples to a final volume of 64.5 μl. The hybridization mixture was then denatured for 3 minutes at 95°C and cooled for 2 minutes on ice before being applied to the array. Lifter-slips (Erie Scientific Company, Shelton, CT, USA) were used to contain the hybridization mixture on the array during hybridization. The arrays were then placed, in hybridization chambers (Corning), in a water bath at 42°C for 14 to 18 hours. After hybridization the slides were washed once for 5 minutes at 42°C in wash buffer 1 (2 × SSC, 0.1% SDS), once for 5 minutes in wash buffer 2 (0.1 × SSC, 0.1% SDS) and five times, for 1 minute each, in wash buffer 3 (0.1 × SSC). The slides were then centrifuged dry and scanned. For further information about the hybridization see SOP 003 at the KTH microarray core facility web site [32] under 'Protocols'.

Scanning and image processing

An Agilent G2565BA scanner was used to scan the slides and acquire 50-megabyte TIFF images. The scanner resolution was set at 10 μm. GenePix 5.1.0.0 (Axon Instruments, Foster City, CA, USA) was used to extract the raw signals from the TIFF images and to assign each spot an ID. Spots defined as 'not found' by GenePix were flagged with a negative flag (- 50) and removed downstream in the analysis. Spots with clearly abnormal morphology due to dust particles or other factors were manually flagged as bad (- 100) and were also removed in downstream analysis. No further processing of the slides was performed in GenePix. The data are available at ArrayExpress, a public repository for MA data (accession number E-MEXP-367) [34].

Data analysis

The data were analyzed mainly with the help of packages in R [35] except for the Expression Analysis Systematic Explorer (EASE) analysis [36], which was performed in MEV [37]. EASE uses the hierarchical structure of gene ontology (GO) [38] to find biological themes among sets of differentially expressed (DE) genes. Each GO category is given an EASE score, which is a conservative adjustment of Fisher's exact probability [39] in which Fisher's exact probability is jackknifed to weigh significance in favor of GO categories supported by many genes. R is a language and environment for statistical computing and graphics. The packages that were used in R were LIMMA [40], aroma [41], the KTH package [32] and bioconductor [42]. All operations performed on the data in R during analysis can be accomplished with these packages. After the result files (gpr) produced by GenePix had been imported into R, unreliable spots with abnormal physical properties were removed with four filters:

1. filterFlags, which removes spots flagged as not found or absent in GenePix.

2. filterSaturated, which removes spots saturated in both cy5 and cy3.

3. filterB2SD, which removes spots in which 70% of the pixels have below background intensity + 2 standard deviations.

4. filterSize, which removes spots that are enlarged due to spotting artefacts.

On average 4,030 spots were removed by the filters, leaving about 25,970 spots for downstream analysis. More information about the filters can be found at the KTH package web site [32]. After filtering, the slides were normalized with print-tip loess (local regression) normalization [43]. To identify DE genes a parametric empirical Bayes approach implemented in LIMMA [40,44] was used. This test statistic will assign a score (B-score) to each gene. The B-score was used to rank the genes so that the gene with the highest score has the highest probability of being DE. When differences were being investigated, two criteria had to be fulfilled for a gene to be regarded as DE: the genes had to have a B-score of more than 0 and an |M-value| of more than 1 (an M-value is the second logarithm of the fold change [43]). When the DE genes were used to cluster the biopsies a third criterion was added: a gene had to occur in at least two or more comparisons between biopsies (regardless of patient) to be used for clustering. This was done to remove noise, because no biological replication was possible. When we were identifying DE genes between tissues with an overrepresentation of adipose cells versus all the other biopsies, we defined a gene as DE if the gene had a B-score of more than 20. Here samples from all patients were used, allowing the approximation of the different parameters used for the test statistics, for example the standard error, to be improved, and thus the B-scores to be high. A cutoff was therefore set at a B-score of more than 20 to investigate a reasonable number of genes with the highest ranking in this comparison.

A moderated t test [40] was performed in parallel, with the use of a false discovery rate [45] correction for multiple testing. Technical replicates were dealt with in different ways depending on the comparison in question. When three different levels of replicates were available (for example, for patients 1 to 3 the levels were technical replicates (i), multiple samples from each biopsy (sub-biopsies) (ii) and the biopsy (iii) itself), instead of just taking the average of technical replicates, we used the duplicateCorrelation function [46] available in LIMMA [40] to acquire an approximation of gene-by-gene variance. This retains valuable information about the variance when fitting a linear model to the data so as to identify DE genes. When four levels were available, as when testing for DE genes between patients 1 and 3 (technical replicates, sub-biopsies, biopsies and patients), the first level was averaged and duplicateCorrelation was used for the sub-biopsies and biopsies. When only two levels of replicates were available (for example when testing for DE genes between biopsies in patient 4 there were technical replicates of each biopsy that was tested) the replicates were all treated in the same way.

Several hierarchical clusterings were performed [47], in which 1 minus the Pearson correlation was used as the distance measure. When creating the agglomerative dendrogram the average distance between each cluster was used. To evaluate the structure the clustering algorithm imposes on the data the cophenetic (coph) correlation coefficient [48] was determined. This measures how well the hierarchical structure from the dendrogram represents the actual distances; coph = 1 indicates perfect representation, whereas coph = 0 indicates no representation. To facilitate color representation in the hierarchical clustering the (log₂) expression value for each gene in each of the biopsies was adjusted by subtracting the respective mean log₂expression value across all biopsies.

Experimental design

Two series of hybridizations were performed in this study. The design chosen for both series was a common reference design that allows both the identification of DE genes in different contexts and unsupervised classification, such as hierarchical clustering. Each hybridization was performed with a technical replicate (the same amount of RNA taken from the same amplified RNA aliquot labeled in two separate reactions and hybridized onto two separate arrays). The average correlation between the M-values of technical replicates was 0.97. In the first series (Figure (Figure1a),1a), aRNA from the orthopedic sub-biopsies/biopsies (depending on the patient they were obtained from; see above) was labeled and hybridized in duplicate versus the reference (also amplified). In all, 76 hybridizations were performed in this series. In the second series (Figure (Figure1b),1b), aRNA from the arthroscopic biopsies was labeled and hybridized in duplicate versus the reference. A total of 32 hybridizations were performed in this series. The gene expression data were filtered and normalized as described above.

Figure 1

Experimental design. (a) Orthopedic samples. Three orthopedic biopsies (b1, b2, and b3) were taken from each of seven patients at random sites of the synovium. For patients 1 to 3 each biopsy was split into three parts. RNA from each sub-biopsy was hybridized (more ...)

Results

Heterogeneity between orthopedic biopsies

Two series of hybridizations were performed to investigate tissue heterogeneity, as indicated in Figure Figure1,1, covering variation both between orthopedic biopsies and between arthroscopic biopsies. In addition, the orthopedic biopsies of patients 1 to 3 were each divided into three parts (sub-biopsies) to study variation in gene expression between adjacent parts of the individual biopsies. First, an analysis of the heterogeneity between the orthopedic biopsies (b1, b2, and b3) from the same patients (Figure (Figure2a)2a) was performed. To measure heterogeneity we applied MA technology and estimated the differences between samples by determining the number of DE genes. Few, if any, genes should be DE between homogeneous tissues. As described in Materials and methods a gene was considered to be DE if it had a B-score of more than 0 and an |M-value| of more than 1 (that is, a more than twofold change). The B-score is derived from a statistical test (a parametric empirical Bayes approach, as described in Materials and methods) used to identify DE genes [44]. The higher the B-score, the more evidence there is of differential expression. The differences in gene expression seen in Figure Figure2a2a were caused by both variation in cellular composition and true changes in gene expression due to various factors. In all, 2,133 genes were identified as DE when comparing biopsies 1 and 2 from patient 4. In contrast, only six genes were found to be DE when comparing biopsies 2 and 3 of patient 2. This demonstrates that there was a high level of variation between biopsies originating from the same joint, indicating that there may be large differences in cellular composition between synovial biopsies obtained from patients with RA.

Figure 2

Differentially expressed genes between orthopedic biopsies and sub-biopsies. (a) As shown in the scheme at the left, comparisons were made patient by patient between all biopsies from all seven patients. Because of the degradation of RNA in the third (more ...)

Heterogeneity between sub-biopsies

Variation between adjacent biopsies was investigated by separately analyzing gene expression in each of the three sub-biopsies derived from the biopsies obtained from patients 1 to 3 (Figure (Figure2b).2b). Interestingly, many genes were DE (B-score>0 and an |M-value|>1) between sub-biopsies, despite apparent homogeneity between biopsies, as exemplified by the biopsies of patient 2 (Figure (Figure2a).2a). However, overall the sub-biopsies were quite similar, with the average number of DE genes being 180, 202, and 132 for patients 1, 2, and 3, respectively, although significant differences were found in these comparisons. For example, biopsies b1 and b3 from patient 1 seemed to be quite similar (Figure (Figure2a),2a), but further investigation found great variation between sub-biopsies of biopsy b1 (Figure (Figure2b).2b). Here, the sub-biopsies A and C seemed similar but there seemed to be large variations in cellular composition between sub-biopsies A and B. Furthermore, the hematoxylin/eosin staining patterns demonstrated that biopsies b1 and b3 of patient 1 contained an overrepresentation of adipose cells (Figure (Figure3),3), whereas b2 seemed to represent synovial hyperplasia. It should be noted that the differences between biopsies seemed to be larger than those between sub-biopsies.

Figure 3

Hematoxylin/eosin staining of the synovial membrane of patients 1 to 7. Original magnification × 250. Biopsy 3 (b3) of patient 2 was available only for RNA extraction and not for staining. Biopsy 3 of patient 7 was not represented in the stainings (more ...)

Hierarchical cluster of orthopedic biopsies

Four levels of replicates were used in this study: technical replicates (the same amount of RNA taken from the same amplified RNA aliquot labeled in two separate reactions and hybridized onto two separate arrays); adjacent biopsies in the form of sub-biopsies; multiple biopsies obtained from different locations in the same joint and patient; and biological replicates in the form of biopsies from different patients. Hierarchical clustering (HCL) analyses (Figure (Figure4)4) were initially performed to obtain an overview of the differences between the technical replicates, resampled biological replicates (sub-biopsies/biopsies), and biological replicates (patients). The first HCL analysis was performed (Figure (Figure4a)4a) on the raw data after filtration and normalization. As can be seen in the resulting clustering, the biopsies separate patient by patient except for those from patient 1. This is consistent with data presented in Figure Figure2a2a for this sample.

Figure 4

Hybridizations of orthopedic samples (patients 1 to 7). (a) Results of hierarchical clustering (HCL) performed on the raw data after normalization. On average, 26,111 features were used per array. Cophenetic (coph) correlation coefficient = 0.9. (b) (more ...)

Hierarchical cluster of orthopedic biopsies using a subset of genes

The DE genes that were found to vary between biopsies within patients (4,175 genes; Figure Figure2a)2a) and were present in one or more comparisons (2,285 genes) were used to cluster the data (Figure (Figure4b).4b). As shown in the dendrogram, this subset of genes contains information that distinguishes between three biopsies (biopsies 1 of patient 1, biopsies 3 of patient 1, and biopsy 1 of patient 4). Overall the results of the HCL analysis in Figure Figure4b4b correlated well with the hematoxylin/eosin staining patterns (Figure (Figure3),3), which demonstrated that these three biopsies had a clear overrepresentation of adipose cells relative to the other biopsies. Note that biopsies 1 and 2 of patient 7 also clustered separately, which is again explained by the histological data although the difference was visually not as obvious. The remaining biopsies clustered according to patient. The three biopsies with an overrepresentation of adipose cells were further investigated by determining the DE genes between these biopsies and all other biopsies. The results can be seen in Figure Figure4c.4c. Genes with a B-score of more than 20 were chosen for further investigation as described in Materials and methods. False discovery rate [45] correction for multiple testing was applied on a moderated T-test, yielding a proportion of 1.03 × 10^-11false positives among the 418 chosen genes with a B-score of more than 20, 341 of which were upregulated in the biopsies with an overrepresentation of adipose cells relative to the rest. These genes were used to elucidate possible biological themes with the use of the GO tool EASE [36]. The results of the EASE analysis (Table (Table2)2) show that the overrepresentation of adipocytes detected in the staining analysis is confirmed by the gene expression data.

Table 2

EASE results

Heterogeneity between arthroscopic biopsies

In the second series of hybridizations (Figure (Figure1b),1b), comparisons were made between all of the arthroscopic biopsies from the respective patients. The barplot in Figure Figure5a5a shows the results of these comparisons. The largest difference (638 DE genes) was observed in patient 12 and the smallest between two biopsies from patient 9 (21 DE genes). The arthroscopic biopsies were clustered by using all of the features on the chip and also a subset of genes that were defined as DE in Figure Figure5a5a ((1,5361,536 genes) and were present in one or more comparisons (421 genes). The results of the cluster analysis can be seen in Figure 5b,c. In contrast to the orthopedic samples, the biopsies from the different patients clustered individually with both approaches.

Figure 5

Hybridizations of arthroscopic samples (patients 8 to 13). (a) Barplot of the differentially expressed genes identified in the patient-by-patient comparisons between biopsies. In total, 1,536 unique genes were found to vary between the biopsies. (b) (more ...)

Differences between biopsies and patients

The average number of DE genes was 143 between the arthroscopic biopsies (Figure (Figure5a)5a) and 455 between the orthopedic biopsies (Figure (Figure2a),2a), indicating that there was less heterogeneity among the arthroscopic biopsies than among the orthopedic samples. The distribution of differences between orthopedic and arthroscopic biopsies is illustrated in boxplots A and B of Figure Figure6,6, respectively. As can also be seen in boxplot C of Figure Figure6,6, when the biopsies containing an overrepresentation of adipose cells are removed from the orthopedic biopsies the distribution of DE genes becomes very similar to that of the arthroscopic biopsies, reducing the average number of DE genes between biopsies from 455 to 171. Interestingly, the heterogeneity between sub-biopsies (boxplot D of Figure Figure6;6; 173 DE genes on average) is similar to that of the reduced set of orthopedic biopsies. This indicates that the gene expression heterogeneity between adjacent biopsies is similar to that of biopsies farther apart if the biopsies consist of the same type of tissue. In addition, we compared sub-biopsy/biopsy variation with patient variation. As with the comparisons between sub-biopsies/biopsies, a gene was regarded as being DE if it had a B-score of more than 0 and an |M-value| of more than 1. The results are displayed in Figure Figure6,6, boxplots E to G. The same patterns that were observed between biopsies were also seen between patients, because the orthopedic patients were more heterogeneous than the arthroscopic patients. However, the differences were reduced when the biopsies that were not suitable for gene expression analysis of synovial inflammation were removed. EASE [36] was used to investigate the possible overrepresentation of biological processes by sorting the genes found to vary between biopsies into subsets based on GO categories. The numbers of 'hits' and the EASE scores for the 10 most overrepresented GO categories are shown in Table Table33 for both the orthopedic (4,175 unique genes) and the arthroscopic (1,536 unique genes) biopsies. No major differences in the relative representation of biological processes between the two types of samples were found in this analysis.

Figure 6

Boxplots of the differentially expressed (DE) genes in the different comparisons. Boxplot A shows the number of DE genes between the arthroscopic biopsies shown in Figure 5a. Boxplot B shows the number of DE genes between the orthopedic biopsies shown (more ...)

Table 3

EASE results for differentially expressed genes identified in hybridization series 1 and 2

Discussion

The purpose of this study was to investigate variations in gene expression in synovial tissues, within and between patients with RA, using MA technology. The study also allowed the sampling techniques used (orthopedic open surgery and rheumatic arthroscopy) to be compared. The results show large differences in the numbers of DE genes between biopsies, supporting previous studies on tissue heterogeneity in synovial joints affected by inflammatory disease [25,30]. The results also demonstrate that the differences between patients are larger than between biopsies obtained from the same joint (Figure (Figure6).6). For the arthroscopic biopsies the unique gene expression signature of each patient dominated the dendrograms obtained using both all genes and the subset of DE genes. This was not observed for the orthopedic biopsies. Most of the orthopedic biopsies clustered according to the patient, but some failed. Instead, they clustered according to their cellular composition; that is, mainly in relation to their content of adipose cells as detected by parallel histology. This has potential implications for how MA analysis can be used to classify and follow the course of patients with arthritis, and in particular for choosing the optimal sampling technique.

Rheumatological arthroscopy has the advantage that the investigator can choose a section of the synovium that is inflamed (vascular and proliferative changes) and avoid much fat or fibrous tissue. In contrast, orthopedic biopsies are most often taken during operation in an extremity where blood supply has been temporarily turned off and where it is difficult to separate inflamed from fibrous or fatty tissue at macroscopic sampling of biopsies. The increased precision in biopsy sampling from inflamed areas during arthroscopic compared with sampling during open joint surgery might therefore explain the patient-restricted clustering for arthroscopic but not for orthopedic biopsies. This observation also has implications for the choice between the use of arthroscopic or blind-needle biopsies, in which sampling of biopsies under direct vision during arthroscopy is more likely to focus on areas of inflammation than sampling with the blind-needle biopsy technique. In addition, the orthopedic samples were taken from patients undergoing end-stage arthroplastic surgery, whereas the samples taken during arthroscopy were from patients with a much shorter disease duration (2.3 years). When samples containing substantial fat tissue rather than inflammatory cells were removed from the analysis, biopsies clustered almost perfectly both at the level of biopsy site and individually.

It therefore seems that biopsies from actively inflamed synovial tissues of patients with RA show expression of unique patterns of mRNA, provided that the biopsy has been taken in such a way that the analysis is performed on cells from an inflamed site. Thus, both inter-individual and intra-individual variation must be taken into consideration when analysing gene expression in synovial tissue. One way is to take multiple biopsies, as suggested earlier [8]; however, this would reduce gene expression changes due to local differences in gene expression at each biopsy site. Instead we suggest taking one biopsy from a site of maximal macroscopic inflammation. If reproduced in still larger series of biopsy studies using the current as well as additional array methodologies, it is therefore feasible that the MA technology performed on single biopsies might provide information that allows us to investigate the value of this information in predicting disease course as well as response to therapy, and to follow the results of therapies over the longitudinal course of a chronic inflammatory disease such as RA.

Conclusion

The results in this paper demonstrate that levels of homogeneity vary between synovial biopsies obtained from the same knee. Nevertheless, the gene expression signature of each biopsy was patient-specific except for four orthopedic biopsies. These were identified by gene expression analysis and confirmed by histochemistry as not being suitable for the analysis of synovial inflammation because of tissue heterogeneity. The results also demonstrate differences between two different sampling methods: open surgery and arthroscopy. The number of DE genes between the orthopedic biopsies was on average almost threefold higher than in the arthroscopic biopsies.

Abbreviations

aRNA = amplified RNA; coph = cophenetic correlation coefficient; DE = differentially expressed; EASE = Expression Analysis Systematic Explorer; GO = gene ontology; HCL = hierarchical clustering; MA = microarray; PCR = polymerase chain reaction; RA = rheumatoid arthritis.

Competing interests

The authors declare that they have no competing interests.

Authors' contributions

J Lindberg performed parts of the MA-related laboratory work, contributed to the data analysis, and wrote parts of the article. EaK performed the arthroscopic surgery, wrote parts of the article, and contributed to the data analysis. A-KU contributed to the planning and design of the project, collected the orthopedic biopsies outside surgery, performed the hematoxylin/eosin staining analysis, and participated in both the data analysis and writing of the article. AS performed the joint replacement surgery and thereby provided the orthopedic biopsies. TA contributed to the planning and design of the project and performed parts of the MA-related laboratory work. PN was responsible for chip design and contributed to data analysis. LK was involved in planning the project, data analysis, and writing the article. J Lundeberg was involved in designing the project, the analysis of laboratory results, data analysis, and writing of the article. All authors read and approved the final manuscript.

Acknowledgements

Marianne Engström, Anna Westring, and Fredrik Lyngman are thanked for their technical assistance, and Anneli Walden is thanked for chip production. The work was supported by grants from the Knut och Alice Wallenbergs Foundation and the Knowledge Foundation.

References

Arnett FC. Revised criteria for the classification of rheumatoid arthritis. Bull Rheum Dis. 1989;38:1–6. [PubMed]
Dolhain RJ, Ter Haar NT, De Kuiper R, Nieuwenhuis IG, Zwinderman AH, Breedveld FC, Miltenburg AM. Distribution of T cells and signs of T-cell activation in the rheumatoid joint: implications for semiquantitative comparative histology. Br J Rheumatol. 1998;37:324–330. doi: 10.1093/rheumatology/37.3.324. [PubMed] [Cross Ref]
Lindblad S, Hedfors E. Intraarticular variation in synovitis. Local macroscopic and microscopic signs of inflammatory activity are significantly correlated. Arthritis Rheum. 1985;28:977–986. [PubMed]
Hutton CW, Hinton C, Dieppe PA. Intra-articular variation of synovial changes in knee arthritis: biopsy study comparing changes in patellofemoral synovium and the medial tibiofemoral synovium. Br J Rheumatol. 1987;26:5–8. [PubMed]
Kennedy TD, Plater-Zyberk C, Partridge TA, Woodrow DF, Maini RN. Representative sample of rheumatoid synovium: a morphometric study. J Clin Pathol. 1988;41:841–846. [PMC free article] [PubMed]
Kirkham B, Portek I, Lee CS, Stavros B, Lenarczyk A, Lassere M, Edmonds J. Intraarticular variability of synovial membrane histology, immunohistology, and cytokine mRNA expression in patients with rheumatoid arthritis. J Rheumatol. 1999;26:777–784. [PubMed]
Ulfgren AK, Grondal L, Lindblad S, Khademi M, Johnell O, Klareskog L, Andersson U. Interindividual and intra-articular variation of proinflammatory cytokines in patients with rheumatoid arthritis: potential implications for treatment. Ann Rheum Dis. 2000;59:439–447. doi: 10.1136/ard.59.6.439. [PMC free article] [PubMed] [Cross Ref]
Boyle DL, Rosengren S, Bugbee W, Kavanaugh A, Firestein GS. Quantitative biomarker analysis of synovial gene expression by real-time PCR. Arthritis Res Ther. 2003;5:R352–R360. doi: 10.1186/ar1004. [PMC free article] [PubMed] [Cross Ref]
Chung CH, Parker JS, Karaca G, Wu J, Funkhouser WK, Moore D, Butterfoss D, Xiang D, Zanation A, Yin X, et al. Molecular classification of head and neck squamous cell carcinomas using patterns of gene expression. Cancer Cell. 2004;5:489–500. doi: 10.1016/S1535-6108(04)00112-6. [PubMed] [Cross Ref]
Golub TR, Slonim DK, Tamayo P, Huard C, Gaasenbeek M, Mesirov JP, Coller H, Loh ML, Downing JR, Caligiuri MA, et al. Molecular classification of cancer: class discovery and class prediction by gene expression monitoring. Science. 1999;286:531–537. doi: 10.1126/science.286.5439.531. [PubMed] [Cross Ref]
Alizadeh AA, Eisen MB, Davis RE, Ma C, Lossos IS, Rosenwald A, Boldrick JC, Sabet H, Tran T, Yu X, et al. Distinct types of diffuse large B-cell lymphoma identified by gene expression profiling. Nature. 2000;403:503–511. doi: 10.1038/35000501. [PubMed] [Cross Ref]
Perou CM, Sorlie T, Eisen MB, van de Rijn M, Jeffrey SS, Rees CA, Pollack JR, Ross DT, Johnsen H, Akslen LA, et al. Molecular portraits of human breast tumours. Nature. 2000;406:747–752. doi: 10.1038/35021093. [PubMed] [Cross Ref]
Valk PJ, Verhaak RG, Beijen MA, Erpelinck CA, Barjesteh van Waalwijk van Doorn-Khosrovani S, Boer JM, Beverloo HB, Moorhouse MJ, van der Spek PJ, Lowenberg B, et al. Prognostically useful gene-expression profiles in acute myeloid leukemia. N Engl J Med. 2004;350:1617–1628. doi: 10.1056/NEJMoa040465. [PubMed] [Cross Ref]
van't Veer LJ, Dai H, van de Vijver MJ, He YD, Hart AA, Mao M, Peterse HL, van der Kooy K, Marton MJ, Witteveen AT, et al. Gene expression profiling predicts clinical outcome of breast cancer. Nature. 2002;415:530–536. doi: 10.1038/415530a. [PubMed] [Cross Ref]
Jarvis JN. Diagnostic and prognostic potential of gene microarrays in rheumatoid arthritis. Expert Rev Mol Diagn. 2005;5:655–659. doi: 10.1586/14737159.5.5.655. [PubMed] [Cross Ref]
Jarvis JN, Dozmorov I, Jiang K, Frank MB, Szodoray P, Alex P, Centola M. Novel approaches to gene expression analysis of active polyarticular juvenile rheumatoid arthritis. Arthritis Res Ther. 2004;6:R15–R32. doi: 10.1186/ar1018. [PMC free article] [PubMed] [Cross Ref]
Batliwalla FM, Baechler EC, Xiao X, Li W, Balasubramanian S, Khalili H, Damle A, Ortmann WA, Perrone A, Kantor AB, et al. Peripheral blood gene expression profiling in rheumatoid arthritis. Genes Immun. 2005;6:388–397. doi: 10.1038/sj.gene.6364209. [PubMed] [Cross Ref]
Maas K, Chan S, Parker J, Slater A, Moore J, Olsen N, Aune TM. Cutting edge: molecular portrait of human autoimmune disease. J Immunol. 2002;169:5–9. [PubMed]
Olsen N, Sokka T, Seehorn CL, Kraft B, Maas K, Moore J, Aune TM. A gene expression signature for recent onset rheumatoid arthritis in peripheral blood mononuclear cells. Ann Rheum Dis. 2004;63:1387–1392. doi: 10.1136/ard.2003.017194. [PMC free article] [PubMed] [Cross Ref]
Jeong JG, Kim JM, Cho H, Hahn W, Yu SS, Kim S. Effects of IL-1beta on gene expression in human rheumatoid synovial fibroblasts. Biochem Biophys Res Commun. 2004;324:3–7. doi: 10.1016/j.bbrc.2004.09.011. [PubMed] [Cross Ref]
Zhang HG, Hyde K, Page GP, Brand JP, Zhou J, Yu S, Allison DB, Hsu HC, Mountz JD. Novel tumor necrosis factor alpha-regulated genes in rheumatoid arthritis. Arthritis Rheum. 2004;50:420–431. doi: 10.1002/art.20037. [PubMed] [Cross Ref]
Gallagher J, Howlin J, McCarthy C, Murphy EP, Bresnihan B, FitzGerald O, Godson C, Brady HR, Martin F. Identification of Naf1/ABIN-1 among TNF-alpha-induced expressed genes in human synoviocytes using oligonucleotide microarrays. FEBS Lett. 2003;551:8–12. doi: 10.1016/S0014-5793(03)00823-8. [PubMed] [Cross Ref]
Tsubaki T, Arita N, Kawakami T, Shiratsuchi T, Yamamoto H, Takubo N, Yamada K, Nakata S, Yamamoto S, Nose M. Characterization of histopathology and gene-expression profiles of synovitis in early rheumatoid arthritis using targeted biopsy specimens. Arthritis Res Ther. 2005;7:R825–R836. doi: 10.1186/ar1751. [PMC free article] [PubMed] [Cross Ref]
Schmutz C, Hulme A, Burman A, Salmon M, Ashton B, Buckley C, Middleton J. Chemokine receptors in the rheumatoid synovium: upregulation of CXCR5. Arthritis Res Ther. 2005;7:R217–R229. doi: 10.1186/ar1475. [PMC free article] [PubMed] [Cross Ref]
van der Pouw Kraan TC, van Gaalen FA, Huizinga TW, Pieterman E, Breedveld FC, Verweij CL. Discovery of distinctive gene expression profiles in rheumatoid synovium using cDNA microarray technology: evidence for the existence of multiple pathways of tissue destruction and repair. Genes Immun. 2003;4:187–196. doi: 10.1038/sj.gene.6363975. [PubMed] [Cross Ref]
Devauchelle V, Marion S, Cagnard N, Mistou S, Falgarone G, Breban M, Letourneur F, Pitaval A, Alibert O, Lucchesi C, et al. DNA microarray allows molecular profiling of rheumatoid arthritis and identification of pathophysiological targets. Genes Immun. 2004;5:597–608. doi: 10.1038/sj.gene.6364132. [PubMed] [Cross Ref]
Kasperkovitz PV, Timmer TC, Smeets TJ, Verbeet NL, Tak PP, van Baarsen LG, Baltus B, Huizinga TW, Pieterman E, Fero M, et al. Fibroblast-like synoviocytes derived from patients with rheumatoid arthritis show the imprint of synovial tissue heterogeneity: evidence of a link between an increased myofibroblast-like phenotype and high-inflammation synovitis. Arthritis Rheum. 2005;52:430–441. doi: 10.1002/art.20811. [PubMed] [Cross Ref]
Lorenz P, Ruschpler P, Koczan D, Stiehl P, Thiesen HJ. From transcriptome to proteome: differentially expressed proteins identified in synovial tissue of patients suffering from rheumatoid arthritis and osteoarthritis by an initial screen with a panel of 791 antibodies. Proteomics. 2003;3:991–1002. doi: 10.1002/pmic.200300412. [PubMed] [Cross Ref]
Huh SJ, Paik DJ, Chung HS, Youn J. Regulation of GRB2 and FLICE2 expression by TNF-alpha in rheumatoid synovium. Immunol Lett. 2003;90:93–96. doi: 10.1016/j.imlet.2003.07.002. [PubMed] [Cross Ref]
van der Pouw Kraan TC, van Gaalen FA, Kasperkovitz PV, Verbeet NL, Smeets TJ, Kraan MC, Fero M, Tak PP, Huizinga TW, Pieterman E, et al. Rheumatoid arthritis is a heterogeneous disease: evidence for differences in the activation of the STAT-1 pathway between rheumatoid tissues. Arthritis Rheum. 2003;48:2132–2145. doi: 10.1002/art.11096. [PubMed] [Cross Ref]
Baeten D, Van den Bosch F, Elewaut D, Stuer A, Veys EM, De Keyser F. Needle arthroscopy of the knee with synovial biopsy sampling: technical experience in 150 patients. Clin Rheumatol. 1999;18:434–441. doi: 10.1007/s100670050134. [PubMed] [Cross Ref]
KTH microarray core facility homepage http://www.ktharray.se/
Imbeaud S, Graudens E, Boulanger V, Barlet X, Zaborski P, Eveno E, Mueller O, Schroeder A, Auffray C. Towards standardization of RNA quality assessment using user-independent classifiers of microcapillary electrophoresis traces. Nucleic Acids Res. 2005;33:e56. doi: 10.1093/nar/gni054. [PMC free article] [PubMed] [Cross Ref]
Parkinson H, Sarkans U, Shojatalab M, Abeygunawardena N, Contrino S, Coulson R, Farne A, Lara GG, Holloway E, Kapushesky M, et al. ArrayExpress – a public repository for microarray gene expression data at the EBI. Nucleic Acids Res. 2005:D553–D555. [PMC free article] [PubMed]
Team RDC. A Language and Environment for Statistical Computing. Vienna: R Foundation for Statistical Computing; 2004.
Hosack DA, Dennis G, Jr, Sherman BT, Lane HC, Lempicki RA. Identifying biological themes within lists of genes with EASE. Genome Biol. 2003;4:R70. doi: 10.1186/gb-2003-4-10-r70. [PMC free article] [PubMed] [Cross Ref]
Saeed A, Sharov V, White J, Li J, Liang W, Bhagabati N, Braisted J, Klapa M, Currier T, Thiagarajan M, et al. TM4: a free, open-source system for microarray data management and analysis. Biotechniques. 2003;34:374–378. [PubMed]
Ashburner M, Ball CA, Blake JA, Botstein D, Butler H, Cherry JM, Davis AP, Dolinski K, Dwight SS, Eppig JT, et al. Gene ontology: tool for the unification of biology. The Gene Ontology Consortium. Nat Genet. 2000;25:25–29. doi: 10.1038/75556. [PubMed] [Cross Ref]
Fleiss JL. Statistical Methods for Rates and Proportions. New York: John Wiley; 1981.
Smyth GK. Linear models and empirical Bayes methods for assessing differential expression in microarray experiments. Stat Appl Genet Mol Biol. 2004;3:Article 3. [PubMed]
Bengtsson H. aroma – An R Object-oriented Microarray Analysis Environment [Preprints in Mathematical Sciences, no 18] Lund: Centre for Mathematical Sciences, Lund Institute of Technology; 2004.
Gentleman RC, Carey VJ, Bates DM, Bolstad B, Dettling M, Dudoit S, Ellis B, Gautier L, Ge Y, Gentry J, et al. Bioconductor: open software development for computational biology and bioinformatics. Genome Biol. 2004;5:R80. doi: 10.1186/gb-2004-5-10-r80. [PMC free article] [PubMed] [Cross Ref]
Smyth GK, Speed TP. Normalization of cDNA microarray data. Methods. 2003;31:265–273. doi: 10.1016/S1046-2023(03)00155-5. [PubMed] [Cross Ref]
Lonnstedt I, Speed TP. Replicated microarray data. Stat Sin. 2002;12:31–46.
Reiner A, Yekutieli D, Benjamini Y. Identifying differentially expressed genes using false discovery rate controlling procedures. Bioinformatics. 2003;19:368–375. doi: 10.1093/bioinformatics/btf877. [PubMed] [Cross Ref]
Smyth GK, Michaud J, Scott H. The use of within-array replicate spots for assessing differential expression in microarray experiments. Bioinformatics. 2005;21:2067–2075. doi: 10.1093/bioinformatics/bti270. [PubMed] [Cross Ref]
Eisen MB, Spellman PT, Brown PO, Botstein D. Cluster analysis and display of genome-wide expression patterns. Proc Natl Acad Sci USA. 1998;95:14863–14868. doi: 10.1073/pnas.95.25.14863. [PubMed] [Cross Ref]
Sneath PHA, Sokal RR. Numerical Taxonomy; The Principles and Practice of Numerical Classification. San Francisco: WH Freeman; 1973.

Articles from Arthritis Research & Therapy are provided here courtesy of
BioMed Central