Pharmacological Effect of Caulophyllum robustum on Collagen-Induced Arthritis and Regulation of Nitric Oxide, NF-κB, and Proinflammatory Cytokines In Vivo and In Vitro

Caulophyllum robustum Maxim (C. robustum) has commonly been used as traditional Chinese medicine for the treatment of rheumatic pain and rheumatoid arthritis (RA) in China. This paper first investigated the anti-inflammation effect of C. robustum extraction (CRME) on RAW264.7 cells stimulated by lipopolysaccharide (LPS) and gene expression levels of inflammatory factors. Moreover, we first evaluated the anti-RA effects of CRME using collagen-induced arthritis (CIA) in DBA/1J mice, and the incidence, clinical score, and joint histopathology were evaluated. The levels of IL-1, IL-6, TNF-α, and PGE2 inflammatory factors in sera of mice were detected by enzyme-linked immunosorbent assay. The expression of NF-κB p65 in the joint was tested by immune histochemical technique. The results showed that, compared with the model group, CRME significantly improved symptoms of the arthritis index, limb swelling, and histological findings by decreasing synovial membrane damage, the extent of inflammatory cell infiltration, and the expansion of capillaries in CIA mice. The results also showed that CRME can reduce the levels of IL-1, IL-6, TNF-α, and PGE2 and inhibit the expression of NF-κB p65. All these results indicated the anti-inflammatory efficacy of CRME as a novel botanical extraction for the treatment of RA.

Introduction
As an autoimmune disease, human rheumatoid arthritis (RA) is a major cause of disability characterized by chronic inflammation and joint damage. The pathogenesis of RA is a complex process involving synovial cell proliferation and fibrosis, pannus formation, and cartilage and bone erosion [1,2]. To date, it has been found that the joint destruction is caused by several mechanisms including the overexpression of inflammatory response cytokines and transcription factors relevant to RA inflammation during the occurrence of RA [3]. In RA, the transcription factor nuclear factor-kappa B (NF-B) is overexpressed in the inflamed synovium [4,5]. It participates in the transactivation of various genes and inducible nitric oxide (iNOS) [6], which mediate immune and inflammatory responses [7]. It can activate and be activated by proinflammatory cytokines, such as Tumor Necrosis Factoralpha (TNF-), Interleukin-(IL-) 1, Prostaglandin E (PGE), and IL-6, resulting in a positive regulatory circle that causes amplification of local inflammatory responses [8]. TNFstimulates synovial fibroblasts to express intercellular adhesion molecules, which increases the migration of leukocytes into the RA joints. It is a central factor in the pathogenesis of RA as it amplifies inflammation and causes joint damage. IL-6 is also known as another proinflammatory cytokine that closely contributes to vascular endothelial growth and that stimulates the growth of blood vessels, including those in hypertrophy of joints. IL-1 is another key mediator involved in bone resorption and cartilage destruction that can modulate the production of nitric oxide and PGE2. PGE2 can act to mediate hyperalgesia-induced sensitizing pain receptors and induce fever [9]. Relieving inflammation 2 Evidence-Based Complementary and Alternative Medicine and reducing proinflammatory cytokine overexpression are therefore regarded as one of the major approaches to help prevent the progressive damage to RA joints.
Historically, C. robustum is a traditional Chinese herb used in the treatment of RA [10]. It is acrid and bitter in taste and warm in nature. The Dictionary of Chinese Medicinal Plants records that C. robustum roots and rhizomes have the ability to expel wind, promote blood circulation, and regulate the flow of vital energy to stop pain. It is mainly used for external injuries, rheumatism, and RA treatment. Pharmacological study showed that its decoction inhibits tissue edema caused by croton oil, suppresses capillary permeability induced by acetic acid, and suppresses cotton pellet granuloma [11]. Modern studies have demonstrated that alkaloids and triterpenoid saponins seem to contribute most of the physiological anti-inflammatory activity [12]. Among them, Caulosides A∼D could dramatically suppress the expression of proinflammatory factors, for instance, TNF-, IL-1 , and IL-6 [13][14][15]. Magnoflorine has shown activity against FCAinduced pyrexia [16]. However, some potential toxic and side effects of C. robustum have been found. For its convenient clinical use, we improved the extraction technique and obtained CRME with lower toxicity (LD50 is 4.6 g).
We set out to investigate the anti-inflammation effect of CRME on RAW264.7 cells stimulated by LPS. Then, collageninduced arthritis (CIA) by the immunization of DAB/1 mice with type II collagen in complete Freund's adjuvant was used for activity evaluation in vivo. The levels of IL-1 , IL-6, TNF-, and PGE2 in sera of the CIA rats and of NF-B protein were determined in synovial tissue. All the above were helpful to our investigation into the therapeutic effects of CRME on RA and its possible mechanism.

Sample Extraction and
Isolation of C. robustum. C. robustum was purchased from Heilongjiang province of China and authenticated by Professor Shaowa Lv (College of Pharmacy, Heilongjiang University of Chinese Medicine). A specimen (20150301) of it was deposited in the College of Pharmacy, Heilongjiang University of Chinese Medicine. Its roots and rhizomes were extracted twice with 70% ethanol (v/v). The filter was combined and reclaimed and then loaded onto an AB-8 macropore resin column and eluted with pure water and 70% ethanol (v/v). The 70% ethanol eluate was concentrated to dryness (the powder was named CRME). We determined its main chemical components by high-performance liquid chromatography (HPLC). The result showed that there were six components in CREM, including Cauloside H, Leonticin D, Cauloside G, Cauloside D, Cauloside B, and Cauloside C and their content was 6.21%, 5.14%, 28.88%, 10.19%, 3.68%, and 2.26%, respectively (Figures 1(a) and 1(b)).  , MO). IL-1 , IL-6, TNF-, PGE2, and NO were purchased from R&D system in the USA. NF-Bp65 was purchased from Wuhan Boster Biological Engineering Co., Ltd. The rabbit-anti-mouse second antibody and DAB were purchased from Beijing Chinese Fir Golden Bridge Biotechnology Co., Ltd.

Drugs and
Triptolide-(Tri-) tablets were purchased from Broad medicine Yellowstone Feiyun Pharmaceutical Co., Ltd., lot number: 20130601, standard: 10 mg. Methotrexate (MTX) tablets were purchased from Shanghai Sym. Pharma Co., Ltd., lot number: 03614040, standard: 2.5 mg. To determine the therapeutic efficiency of CRME for the treatment of RA, we chose MTX and Tri as positive medicines. MTX is the first-line therapy for this disease, and Tri is isolated from Tripterygium Wilfordii Hook f., which is a kind of Chinese medicinal plant used to ameliorate the symptoms of rheumatic diseases [17].

Cell Cultures and Their
Use to Test CRME. RAW264.7 cells were maintained in DMEM containing 10% fetal bovine serum (FBS) and 1% penicillin streptomycin (PS) at 37 ∘ C in a 5% CO2 incubator. Before starting any experiment, RAW264.7 cells were seeded in culture plates and maintained for 24 h until approximately 70%-80% confluence. RAW264.7 cells were stimulated with LPS of Escherichia coli (serotype, 055: B5) at a concentration of 1 g/ml for establishing the inflammatory cell model. RAW264.7 cells were incubated with CRME by the addition of different concentrations of the former (0.250 mg/mL, 0.125 mg/mL, and 0.063 mg/mL) for 1 h and then treated with 1 g/mL LPS for an additional 24 h. After 24 h, the cell culture supernatant was collected to detect NO (nitric oxide) by Griess reagent assay. Total RNA was isolated immediately after harvesting the RAW264.7 cells for RT-PCR as described below. Unstimulated RAW264.7 cells were treated as a negative control.

RNA Extraction and RT-PCR Analysis.
Total RNA was separated from RAW264.7 cells by Trizol reagent according to the protocol. For cDNA synthesis, total RNA (2 g/sample) was reverse transcribed with ReverTra Ace qPCR RT Master Mix with the gDNA Remove. The cDNA was then amplified with THUNDERBIRD SYBR qPCR Mix. The primer sequences were as follows: for IL-1 , forward primer (FP) =

Animals, Induction of CIA, and Compound Admini-
Stration. DBA/1 mice (half male and female, 6∼8 weeks old) were purchased from HFK Bioscience Co., Ltd. (Beijing, China, Certificate: 2010A050). All mice were maintained in a room with a temperature of 25 ± 2 ∘ C and humidity of 60%. Mice drank water freely and were given 12 h of light each day. Bedding was kept dry. All animal care and treatment procedures were approved by the committee of Experimental Animal Administration of Heilongjiang University of Chinese Medicine.
Forty-two DBA/1 mice were divided into seven groups with identical numbers of mice ( = 6). One of the groups served as "normal," whereas the other 6 groups were subjected to CIA induction. Bovine type II collagen was dissolved in 0.1 M acetic acid, mixed to 2 mg/ml mixed liquor, and incubated overnight at 4 ∘ C. The mixture and an equal volume of CFA were emulsified at low temperature. The DBA/1 mice were subcutaneously immunized at 1-2 cm from the base of the tail with 100 L hybrid emulsion (containing 100 g of bovine type II collagen), and this was treated as the first immunization on day 0. Mice were boosted intraperitoneally with 100 g bovine type II collagen emulsified in an equal volume of IFA on day 21, and this was treated as the second immunization. CRME, Tri tablets, and MTX tablets were separately dissolved in distilled water and administered orally once per day for 30 days, from day 22 to day 52. The groups are as follows: normal group (Normal), CIA model group (Vehicle), Tri group 11.3 mg/kg, MTX group 0.9 mg/kg, CRME group (L group) 24.8 mg/(kg⋅day), CRME group (M group) 49.7 mg/(kg⋅day), and CRME group (H group) 99.4 mg/(kg⋅day). The dosage selection for the present study was according to the LD50 (4.6 g/kg) for CRME from our previous study and reference to classical Chinese medicine books and literature. For example, Guizhou folk medicine records that 5 g C. robustum powder taken with wine has the most significant effect. In addition, experiments show that C. robustum ethanol extract (69.23, 34.61, and 17.31 mg/kg) 4 Evidence-Based Complementary and Alternative Medicine inhibits the swelling of the joint and decreases the spleen index and the level of inflammatory factors [18].

Morphologic and Immune Studies.
Mice were observed once every 3 days after primary immunization. Arthritis severity was evaluated by arthritis index, arthritis incidence, and percentage of arthritic limbs as assessed by two independent, blinded observers. All four limbs of the mice were evaluated and scored from 0 to 4 as follows: 0, normal; 1, 1∼2 toe joints swelling; 2, the whole ball with mild swelling and redness; 3, severe swelling and redness of the whole ball; 4, overall deformation and/or joint stiffness. The arthritis index of all four limbs was the total score, and the highest score for each mouse was 16 points. When the score of a paw was more than 2, the mice were considered to suffer from RA. The volume of arthritic limbs was detected with a digital plethysmometer (LE 7500, Panlab, Spain). According to the digit volume of mice and the following formula, the limb swelling of mice was calculated: There is a close connection between the aggressive progression of RA and the loss of body mass. After the secondary immunization, we detected the changes in body weight once every 3 days in order to evaluate CRME for the treatment of RA. The weight of individual CIA mice was calculated with the following formula: Change of body weight (%) = body weight (day 31 of arthritis) body weight (day 1 of arthritis) × 100%. ( Mice were sacrificed by cervical dislocation after a month of oral administration. The spleen and thymus were obtained, fur and blood on the surface of the viscera were rinsed clean with normal saline, and they were weighed with an analytical balance. We measured the wet weight of the immune organs to preliminarily assess CRME for treatment of RA. The computation formula is as follows: 2.7. Enzyme-Linked Immunosorbent (ELISA) Assay. On the 52nd day, 1 mL of blood was collected from each mouse's eye socket vein after more than 12 h fasting. Then, the blood samples stood at room temperature for 20 min. Serum was separated from the rest of the blood using a high-speed centrifuge (H2050R, China) at 10000 rpm for 10 min, and IL-1 , IL-6, TNF-, and PGE2 levels were measured with ELISA kits. According to the manufacturer's instructions, the processed samples were measured for OD value at 450 nm. There is a certain linear relationship between cytokine concentrations in mice and the OD value, which can be determined by acknowledgement of standard curve sample concentrations of cytokines in mice, and, then, the data were analyzed to sum up the influence of CRME on the inflammatory factors.

Histological.
On the 52nd day, the mice were sacrificed. The hind limbs and immune organs, such as spleen and thymus, were fixed immediately in 4% buffered formalin solution, decalcified, and embedded in paraffin, after which they were divided into 6 m pieces. The left hind limbs and visceral organs were treated with Hematoxylin-Eosin (HE) staining for morphological examination, while the right hind limbs were used for immunohistochemistry with NF-Bp65. A Motic 3000 micrograph system served as radiography at 400 times, and Image-pro plus 6.0 pathological Image was used to analyze positive expression quantitatively. The expression of protein was represented by integral optical density (IOD), namely, integral optical density (IOD) = positive area × average optical density. Each case was analyzed randomly in two vision fields at high magnification, and, then, the relative expression was averaged, for which there were 3 cases in each group.
2.9. Statistics. The software SPSS version 19.0 for Windows was applied for statistical analysis. Measurement data were expressed as the mean ± SD, and the two-way ANOVA test followed by least-significant differences (LSD) test was used to compare between groups; count data were expressed as percentages and compared with the chi-square test between groups. Statistical significance was set at * < 0.05, and extreme significance was set at # < 0.01.

Levels of Nitric Oxide (NO) in Culture Supernatant of RAW 264.7 Cells.
The CRME showed dose-dependent inhibition of the generation of NO in LPS-stimulated RAW264.7 cells (Figure 2(a)). This finding suggested that the CRME may suppress the LPS-induced inflammatory response through inhibition of NO generation. Therefore, this result showed that CRME has anti-inflammatory response.

Gene Expression of IL-1 , IL-6, and TNF-in Cell
Culture. CRME exerted an anti-inflammatory effect on LPSinduced responses accompanied by the gene expression of proinflammatory cytokines. RAW264.7 cells were treated with CRME and LPS for 24 h. The gene expression levels of TNF-, IL-1 , and IL-6 were reduced by treatment with CRME in a dose-dependent manner (Figures 2(b), 2(c), and 2(d)). These results indicated that CRME has an antiinflammatory effect on the gene expression of LPS-induced proinflammatory cytokines in RAW264.7 cells. with the normal group, joints of mice in the vehicle group appeared obviously swollen and red, which indicated the success of the modeling. Starting from the 10th day, arthritis indexes were significantly lower for the CIA mice in each administered group than for the model group (Figure 3(a)). As shown in Figure 3(b), the arthritis index of the H group had a significant difference from vehicle group ( < 0.01) starting from the 13th day, and with the passage of time, differences between the H group and positive medicine group increased gradually with statistical significance (All < 0.05). As shown in Figure 3(c), the limb swelling degree of CIA mice in the H group was significantly lower than that in the positive medicine group (All < 0.05). On day 51, the CIA swelling degree of the H group had dropped to 29.8%.

Morphologic Study of Survival
Remarkably, CRME could effectively suppress the loss of body weight of CIA mice (for the L, M, and H groups versus Effects of CRME, MTX, and Tri treatment on the severity of disease progression as evidenced by arthritic index (b) and swelling degree (c) are presented here. There was a significant suppression of both the arthritic index and swelling degree by CRME and positive drug treatment within days 27∼51. CRME on severity of arthritis of CIA mice was evaluated through body weight (d) and visceral index (e). Data are expressed as the mean ± SD of = 6 animals per group. Data were analyzed by two-way ANOVA and Student-Newman-Keuls test as post hoc test ( # < 0.01, compared with the normal group. * < 0.05 and * * < 0.01, compared with vehicle group. △ < 0.05, △△ < 0.01, and △△△ < 0.001, compared with MTX group. f < 0.05, ff < 0.01, and fff < 0.001, compared with Tri group). vehicle: < 0.05). The weight loss of the H group was obviously lower than that of MTX and Tri, and the differences were statistically significant (All < 0.05) (Figure 3(d)). In addition, the spleen index and thymus index of mice were lower in the CRME groups and positive drug groups than in the vehicle group, with significant differences (All < 0.01) (Figure 3(e)). (Figure 4) showed that synovial membrane hyperplasia in the normal group had not been found, its cell nuclear and cytoplasmic staining were normal, and loose underlying connective tissue was seen clearly; vascular arteries and veins showed no expansion and no blood stasis or inflammatory cell infiltration. Compared with the normal group, the synovial membrane of the vehicle group was thick, partly falling off, and uneven. It had a large number of inflammatory cells; at the same time, it was accompanied by capillary expansion. While compared with vehicle group, the synovium of CIA mice had thinned, with neat edges. The number of inflammatory cells was reduced, and capillary expansion was not obvious. Expansion of vascular arteries and veins and blood stasis were less than in the MTX (0.9 mg/kg) and Tri groups (11.3 mg/kg). Results of histology of the left hind limbs and immune organs. At the end of the experiment, the hind limbs and immune organs of CIA mice were examined for histology. Mice in H group were orally given 99.4 mg/kg CRME. Pathological sections of joints (a part of the synovial cells in the red box, with arrows showing vascular artery and vein expansion) are shown in Figure 4.

Analysis of Proinflammatory
Cytokines. TNF-, IL-1 , IL-6, and PGE2 levels are indicated in Figure 5. As shown in the results, the levels of IL-1 , IL-6, TNF-, and PGE2 in sera of CIA mice were apparently higher in the vehicle group than the normal group, differences that were statistically significant (All < 0.01). A significant decrease in TNF-, IL-1 , IL-6, and PGE2 levels was observed in all administered groups (H, L, and M group) when compared to the vehicle group. In addition, the levels of IL-1 , IL-6, TNF-, and PGE2 of CIA mice were distinctly lower in the H group than in the MTX group or Tri group ( < 0.05), while the differences among the L group, M group, and positive drug groups were not statistically significant ( > 0.05) ( Figure 5). (d) were monitored in the medium by using ELISA. Data represent the mean ± SD of three separate experiments. Data were analyzed by two-way ANOVA and Student-Newman-Keuls test as post hoc test ( # < 0.01, compared with the normal group. * < 0.05 and * * < 0.01, compared with vehicle group. △ < 0.05, △△ < 0.01, and △△△ < 0.001, compared with MTX group. f < 0.05, ff < 0.01, and fff < 0.001, compared with Tri group). Figure 6(a), the NF-B p65 positive expression parts are tan or brown particles. Data analysis indicated in Figure 6(b) shows that NF-B p65 is expressed more positively in vehicle group than in treatment group ( < 0.05). What is more, the IOD volume of NF-B p65 in CIA mice from the H group was 1325.18, which was lower than that from the positive drug group (Tri: 1765.43, MTX: 1751.42) ( < 0.001).

3.7.
Discussion. Thus, it may be functionally important to tightly regulate the degree of phlogistic pathway activation, which could be a good therapeutic target to resist RA. For the last decade, efforts have been made to develop these agents. However, the long-term administration of anti-inflammatory drugs is limited due to the side effects. Now, there are biological inhibitors such as TNF-inhibitors and infliximab [19,20]. These activators can also induce specific immunity to increase the self-protection ability of the organism. An excessive inflammatory reaction can also damage target cells and tissue. Therefore, novel anti-RA agents with fewer side effects are needed. TCM can offer a promising repertoire of potentially therapeutic agents for RA [21]. C. robustum has often been used for treating RA patients in China. In books, such as the compendium of Chinese traditional herbal medicine, it has been recorded that C. robustum tincture can be made with 15 g of roots and rhizomes soaked in 300 mL of alcohol for 7 days, and this was used to treat joint pain, bruises, and RA. In this paper, we found that CRME possessed inhibitory activity against NO, IL-1 , IL-6, and TNF-in RAW267.4 cells stimulated by LPS.
Our present data show that all the doses of CRME (24.8, 49.6, and 99.4 mg/kg) significantly inhibited edema compared to their untreated counterparts and with greater effect at 99.4 mg/kg than Tri at 11.3 mg/kg. In reducing the loss of body weight, CRME also has an advantage over MTX and Tri. Then, CRME improved histological findings by decreasing the extent of inflammatory cell infiltration and the capillary expansion of CIA mice more effectively than MTX and Tri. These results were also confirmed by the expression of inflammatory cytokines by ELISA, which showed that TNF-, IL-6, IL-1 , and PGE2 were distinctly lower than in the MTX group and Tri group. The results showed that oral administration of CRME may effectively preserve both synovial membrane damage and inflammatory cell infiltration of inflamed joints. This clearly showed the recuperative RA function and lower toxicity of CRME at the given dosages. However, it should be noted that long-term intake of high concentration of CRME (1.59 g/kg) will show gastrointestinal irritation and vomiting from our previous study. Furthermore, Kennelly et al. showed that the alkaloid contained in the C. robustum had a potential teratogenic effect in mice embryonic culture experiments in vitro [12]. It has also been reported that Cauloside B and Cauloside C had cytotoxicity on the growth of sea urchin embryos by altering cell permeability [13,14]. Kawai et al. studies have shown that Caulosides A∼C have a high hemolytic effect [22]. Therefore, patients with gastrointestinal diseases and pregnant women should not use CRME.
C. robustum consists of more than 60 constituents, including sterols, triterpenes, glycosides, and alkaloids. Triterpene saponins are believed to be the major bioactive compounds contributing to the anti-inflammatory effects on RA [10]. Triterpene saponins constitute an important class of potential pharmacological agents possessing a range of different physiological activities including anticancer, 10 Evidence-Based Complementary and Alternative Medicine anticoagulation, anti-inflammation, antibacterial, analgesic, and comparative immune-modulation activity [23]. Recently, triterpene saponins have attracted positive research interest. Of great interest is the possibility that this class of constituents could be a source of drugs for the treatment of several diseases. We investigated the triterpene saponin content in all CRME. Our findings suggest that CRME is mainly composed of six marker compounds including Cauloside H, Leonticin D, Cauloside G, Cauloside D, Cauloside B, and Cauloside C, which constitute more than 50 percent. In particular, the abovementioned compounds seem to act as inhibitors of inflammatory factors [24].
The inhibition of NF-B might result in a marked amelioration of chronic inflammation by suppressing the production of cytokines [25]. Measurable levels of cytokines such as TNF-, IL-6, IL-1 , and PGE2 are found locally in synovial tissue, the major site of inflammation in CIA as well as RA, where they facilitate the local and systemic inflammatory response [26,27]. Existing research showed that the NF-B in synovial cells induced by TNF-can mediate nuclear transfer leading to NF-B signaling pathway activation [28]. To be more specific, to activate the NF-B pathway, TNFhas a variety of signal molecules involved, including TNFR related factor, RIP (receptor interacting protein), MAP3K (mitogen activated protein kinase 3), and IKK compounds [29]. In turn, TNF-produces physiological reaction via NF-B acting as an important mediator, and activated NF-B can also induce the expression of TNF-and enhance TNFgene transcription. In some cells and tissues, NF-B p65 and TNF-both promote each other, which creates a vicious cycle that increases RA. In RA, IL-6 mainly induces liver cell synthesis of a variety of acute phase reaction proteins and promotes the synthesis of B cells to produce immunoglobulin and rheumatoid factor [30]. It has been confirmed that NF-B can promote gene transcription of inflammation factors such as TNF-, IL-6, and IL-1 [31]. In addition, it might enhance the inflammatory damage of TNF-and IL-1 by improving the generation and release of TNF-and IL-1 [32]. IL-1 is of importance in passing information and mediating immune cell activation, proliferation, and inflammatory response. IL-1 , acting as another important extracellular cue, can activate NF-B to amplify the inflammatory response. NF-B is promoted to transport into the nucleus by IL-1 , where it combines with the promoters of target genes to start target protein expression [33]. As one of the most important targets for the treatment of inflammation and tumor molecules, PGE2 can promote the division and proliferation of fibroblasts and induce microvascular formation and growth [34]. Like TNF-, IL-6, and IL-1 , PGE2 participates in synovium and joint inflammation of surrounding tissue, cartilage destruction, and hyperplasia lesions through the NF-B signaling pathway [35].
We confirmed that the pathogenesis of CIA was accompanied by a substantial increase in the levels of TNF-, IL-6, IL-1 , and PGE2 in the plasma of mice. Interestingly, the administration of CRME drastically reduced the levels of these four proinflammatory cytokines significantly. Our results indicated that CRME effectively improved arthritic symptoms including limb swelling by downregulating the production of proinflammatory cytokines. Further, we suspect that CRME, being an important negative regulator of these proinflammatory cytokines related to RA, might work by suppressing the expression of NF-B.

Conclusions
It has been thoroughly demonstrated that CRME is an antiinflammatory drug as observed under morphological study in CIA mice. This was evident from not only the significantly improved symptoms of the arthritis index and limb swelling but also the improved histological findings of decreased synovial membrane damage and inflammatory cell infiltration in CIA mice. We also revealed that the protective mechanisms of CRME treatment could be partially explained by a decrease in the proinflammatory cytokines TNF-, IL-6, IL-1 , and PGE2. Specifically, it has been found to be very effective at ameliorating the damage resulting from progression of arthritis by controlling the expression of NF-B. Additional studies are required with more physiological and pharmaceutical investigations to establish its clinical applicability.


Differential methylation analysis of reduced representation bisulfite sequencing experiments using edgeR

Cytosine methylation is an important DNA epigenetic modification. In vertebrates, methylation occurs at CpG sites, which are dinucleotides where a cytosine is immediately followed by a guanine in the DNA sequence from 5' to 3'. When located in the promoter region of a gene, DNA methylation is often associated with transcriptional silencing of the gene. Aberrant DNA methylation is associated with the development of various diseases such as cancer. Bisulfite sequencing (BS-seq) is the current "gold-standard" technology for high-resolution profiling of DNA methylation. Reduced representation bisulfite sequencing (RRBS) is an efficient form of BS-seq that targets CpG-rich DNA regions in order to save sequencing costs. A typical bioinformatics aim is to identify CpGs that are differentially methylated (DM) between experimental conditions. This workflow demonstrates that differential methylation analysis of RRBS data can be conducted using software and methodology originally developed for RNA-seq data. The RNA-seq pipeline is adapted to methylation by adding extra columns to the design matrix to account for read coverage at each CpG, after which the RRBS and RNA-seq pipelines are almost identical. This approach is statistically natural and gives analysts access to a rich collection of analysis tools including generalized linear models, gene set testing and pathway analysis. The article presents a complete start to finish case study analysis of RRBS profiles of different cell populations from the mouse mammary gland using the Bioconductor package edgeR. We show that lineage-committed cells are typically hyper-methylated compared to progenitor cells and this is true on all the autosomes but not the sex chromosomes. We demonstrate a strong negative correlation between methylation of promoter regions and gene expression as measured by RNA-seq for the same cell types, showing that methylation is a regulatory mechanism involved in epithelial linear commitment.

Introduction
Cytosine methylation is an important epigenetic DNA modification that is generally associated with transcriptional silencing 1 . In vertebrates, methylation occurs at CpG sites, which are dinucleotides where a cytosine (C) is immediately followed by a guanine (G) in the DNA sequence from 5' to 3'. CpG dinucleotides are relatively uncommon on the human genome but occur more frequently in gene promoters and exons 2 . About 72% of human gene promoters are enriched for CpGs 2 . CpGs in gene promoters tend to cluster in CpG islands (CGIs), which are regions of a few hundred to a couple of thousand base pairs with very strong enrichment of CpGs 3 .
The relationship of DNA methylation to transcription in vertebrates is complex 4 . Methylation of CGIs causes robust transcriptional repression and is required for long-term mono-allelic silencing including X inactivation and genomic imprinting 1 . Methylation of CpG-poor promoters is more weakly associated with gene expression, and the mechanisms by which this occurs are unclear 1 . Methylation of gene bodies on the other hand can be positively associated with gene expression 4,5 .
DNA methylation is relatively stable in that most CGIs do not change methylation state during normal cell development. Nevertheless DNA methylation is understood to play a regulatory role in differentiation and commitment in adult cell lineages 6,7 . Aberrant methylation patterns are also associated with the development of diseases such as cancer 8,9 .
Bisulfite sequencing (BS-seq) is increasingly used to profile DNA methylation 10,11 . Unmethylated cytosines (C) are converted to Uracils (U) by sodium bisulfite and then deaminated to thymines (T) during PCR amplification. Methylated Cs, on the other hand, remain intact after bisulfite treatment. This strategy produces whole genome bisulfite sequencing (WGBS) when combined with sequencing of the entire genome. WGBS is sometimes considered the "gold standard" for methylation profiling because it provides single-nucleotide resolution and whole-genome coverage 11 . However it requires large quantities of DNA and is expensive because of the amount of sequencing required.
The fact that CpG islands constitute only a small percentage of the genome makes the WGBS approach inefficient in terms of information content per sequenced read. To improve efficiency and reduce costs, enrichment strategies have been developed and combined with BS-seq to target a specific fraction of the genome. A common targeted approach is reduced representation bisulfite sequencing (RRBS) that targets CpG-rich regions 11,12 . Under the RRBS strategy, small fragments that compose only 1% of the genome are generated using MspI digestion, which means fewer reads need to be sequenced in total to provide reasonable coverage of the targeted regions. The RRBS approach can capture approximately 70% of gene promoters and 85% of CpG islands, while requiring only small quantities of input sample 13 . RRBS has great advantages in cost and efficiency, especially when CGI or geneorientated results are required. RRBS is also applicable for single cell studies 14 .
developed for detecting DMRs using the BS-seq technology. methylkit 23 and RnBeads 24 implement Fisher's Exact Test, which is a popular choice for two-group comparisons with no replicates. In the case of complex experimental designs, regression methods are widely used to model methylation levels or read counts. RnBeads offers a linear regression approach based on the moderated t-test and empirical Bayes method implemented in limma 25 . BSmooth 26 is another analysis pipeline that uses linear regression together with a local likelihood smoother. methylkit also has an option to apply logistic regression with overdispersion correction 23 . Some other methods have been developed based on beta-binomial distribution to achieve better variance modeling. For example, DSS fits a Bayesian hierarchical beta-binomial model to BS-seq data and uses Wald tests to detect DMRs 27 . Other software using beta-binomial model include BiSeq 28 , MOABS 29 and RADMeth 30 .
In this workflow, we demonstrate an edgeR pipeline for differential methylation analysis. edgeR is one of the most popular Bioconductor packages for assessing differential expression in RNA-seq data 31,32 . It is based on the negative binomial (NB) distribution and it models the variation between biological replicates through the NB dispersion parameter. Unlike other approaches to methylation sequencing data, the analysis explained in this workflow keeps the counts for methylated and unmethylated reads as separate observations. edgeR linear models are used to fit the total read count (methylated plus unmethylated) at each genomic locus, in such a way that the proportion of methylated reads at each locus is modeled indirectly as an over-dispersed binomial-like distribution. This approach has a number of advantages. First, it allows the differential methylation analysis to be undertaken using existing edgeR pipelines developed originally for RNA-seq differential expression analyses. The edgeR generalized linear model (GLM) framework offers great flexibility for analysing complex experimental designs while still accounting for the biological variability 33 . Second, keeping methylated and unmethylated read count as separate data observations allows the inherent variability of the data to be modeled more directly and perhaps more realistically. Differential methylation is assessed by likelihood ratio tests so we do not need to assume that the log-fold-changes or other coefficient estimators are normally distributed.
This article presents an analysis of an RRBS data set generated by the authors containing replicated RRBS profiles of basal and luminal cell populations from the mouse mammary epithelium. Our primary interest is in gene-orientated and pathway-orientated interpretations of the result. It is of particular importance to relate methylation changes to RNA-seq expression changes for the same genes. We show how to analyze differential methylation changes either for individual CpGs or for pre-specified genomic regions, such as chromosomes or genes, and we particularly focus on methylation changes by promoter regions.
As with other articles in the Bioconductor Gateway series, our aim is to provide an example analysis with complete start to finish code. As with other Bioconductor workflow articles, we illustrate one analysis strategy in detail rather than comparing different pipelines. The analysis approach illustrated in this article can in principle be applied to any BS-seq data but is especially appropriate for RRBS data. The approach is designed for experiments that included biological replication but can be used without replication if the NB dispersion is preset. The results shown in this article were generated using Bioconductor Release 3.7.
The next section gives an expository introduction to the edgeR approach to methylation data. The analysis of the mammary epithelial data starts afterwards.

Introducing the NB linear modeling approach to BS-seq data A very small example
To introduce the edgeR linear modeling approach to BS-seq data, consider a genomic locus that has mA methylated and u A unmethylated reads in condition A and m B methylated and u B unmethylated reads in condition B. Our approach is to model all four counts as NB distributed with the same dispersion but different means. Suppose the data is as given in Table 1. The counts can be entered into a matrix in R: > counts <-matrix(c(2,12,11,0),1,4) > dimnames(counts) <-list("Locus", c("A.Me","A.Un","B.Me","B.Un")) > counts A.Me A.Un B.Me B.Un Locus 2 12 11 0 The experimental design of the data can be summarized by indicating which counts correspond to each sample and each condition: In the design matrix, the first two columns are indicators for samples 1 and 2 while the third and fourth columns indicate the methylated count columns for conditions A and B.
If this were a complete data set, then it could be analyzed in edgeR as follows. First we fit a negative binomial generalized linear model to the counts: > library(edgeR) > fit <-glmFit(counts, design, lib.size=c(100,100,100,100), dispersion=0.0247) The first two coefficients in the linear model are used to model the total number of reads (methylated or unmethylated) for samples 1 and 2, respectively. Coefficient 3 estimates the log ratio of methylated to unmethylated reads for condition A, a quantity that can also be viewed as the logit proportion of methylated reads in condition A. In mathematical terms, we have βˆA= log 2 (2.125/12.125) = −2.51. (Note that edgeR adds 0.125 to each count when reporting the coefficients to avoid taking logarithms of zero.) Coefficient 4 estimates the log ratio of methylated to unmethylated reads for condition B, in mathematical terms βˆB = log 2 (11.125/0.125) = 6.48.
Then we conduct a likelihood ratio test to compare the methylation rate in condition B to that in condition A: > lrt <-glmLRT(fit, contrast = c(0,0,-1,1)) > topTags(lrt) The contrast vector indicates that we wish to compare the fourth to the third coefficient, ignoring the two sample coverage coefficients. The estimated value for the contrast is the difference in logit methylation rates, which is estimated as βˆB− βˆA= 6.48 − (−2.51) = 8.99, a quantity that is shown as the log2-fold-change (logFC) column of the results table. The chisquare likelihood ratio statistic (LR) is 20.7 and the P-value for differential methylation in condition B vs condition A is P = 5.27 × 10 −6 .
Note that the specific parametrization used to account for the sample effects is not important. The test for differential methylation would be unchanged if we represented the two sample effects using an intercept model. Indeed, the first two columns of the design matrix could be any two columns spanning the sample effects: > design <-cbind(Intercept = c(1,1,1,1), + Sample2 = c(0,0,1,1), + A = c(1,0,0,0), + B = c(0,0,1,0)) The dispersion parameter controls the degree of biological variability 33 . If we had set dispersion=0 in the above code, then the above analysis would be exactly equivalent to a logistic binomial regression, with the methylated counts as responses and the total counts as sizes, and with a likelihood ratio test for a difference in proportions between conditions A and B. Positive values for the dispersion produce over-dispersion relative to the binomial distribution. We have set the dispersion here equal to the value that is estimated below for the mammary epithelial data.
In the above code, the two library sizes for each sample should be equal. Otherwise, the library size values are arbitrary and any settings would lead to the same P-value.

Relationship to beta-binomial modeling
It is interesting to compare this approach with beta-binomial modeling. It is well known that if m and u are independent Poisson random variables with means µ m and µ u , then the conditional distribution of m given m + u is binomial with success probability p = µ m /(µ m + µ u ). If the Poisson means µ m and µ u themselves follow gamma distributions, then the marginal distributions of m and u are NB instead of Poisson. If the two NB distributions have different dispersions, and have expected values in inverse proportion to the dispersions, then the conditional distribution of m given m + u follows a beta-binomial distribution. The approach taken in this article is closely related to the beta-binomial approach but makes different and seemingly more natural assumptions about the NB distributions. We instead assume the two NB distributions to have the same dispersion but different means. The NB linear modeling approach allows the means and dispersions of the two NB distributions to be estimated separately, in concordance with the data instead of being artificially linked.

Automatic construction of the design matrix
In this expository introduction we have demonstrated the construction of the design matrix from basic principles. In practice, edgeR provides a function modelMatrixMeth to automate the construction of the design matrix. The user can simply specify the treatment conditions associated with each DNA sample (as in the second column of Table 2): > Condition <-factor(c("A","A","B","B")) Then the design matrix is constructed by: Alternatively, a user can construct any appropriate design matrix at the sample level (with 4 rows), exactly as one would do for a RNA-seq differential expression analysis, then the function will expand the sample-level design matrix to the appropriate observation-level matrix with 8 rows: This gives exactly the same result as above.
Another way to make the design matrix When there are just two treatment conditions to be compared, it can be slightly simpler to compute the log2-fold-change directly as a coefficient of the linear model, rather than forming it as a contrast. This can be done by including an intercept in the design matrix as follows: In this formulation, the second last coefficient still estimates the logit methylation level in condition A, but the last coefficient now estimates the log2-fold-change directly: > fit <-glmFit(counts, design, lib.size=100, dispersion=0.0247) > lrt <-glmLRT(fit, coef="ConditionB") > topTags(lrt) Again, this gives the same results as the previous analysis.
Analogy with paired-samples expression analyses The above design matrix might seem as if it is very special to methylation data, but in fact the same sort of design matrix would be used for any gene expression analysis when there are paired-samples and we wish to compare treatment effects between groups. To see this, consider the small RNA-seq experiment shown in Table 3. This experiment has four subjects, each of whom contributes both treated and untreated RNA samples. The subjects belong to two groups (A and B) and we wish to test whether the treatment effect differs between the groups.
Readers will note that Table 3 is the same as Table 2 except that the columns have been relabeled. Here, the four subjects are analogous to the four DNA samples in Table 2, the groups are analogous to conditions, and treatment is analogous to methylation. Exactly the same design matrix, as used above for the methylation experiment, would be appropriate here for the expression experiment. In edgeR, the only difference between the methylation and expression analyses would be in the normalization steps. The linear modeling and testing steps would be identical.
This shows that BS-seq data has a structure that is already familiar from RNA-seq experiments with pairedsamples. The theme of this article is that statistical analysis methods developed for RNA-seq can be beneficially applied to BS-seq data.

NB GLMs
In practice, the data will consist of many thousands of genomic loci and the counts matrix will have a row for each locus. The function glmFit fits a NB GLM to each row of read counts. If y gi is the ith read count for genomic locus g, then y gi is assumed to be NB distributed with expected value µ gi and dispersion φ g . The expected values are modeled by the design matrix, where the x ij are the ith row of the design matrix, β gj are the regression coefficients and N i are the library sizes.
In the small example above with four samples, the design matrix had 6 columns, so p = 6. In general, if we have an experiment comparing k treatment conditions using n samples, then the count matrix will have 2n columns and the design matrix will have p = n + k columns. The first n coefficients capture sample effects in read abundance, while the remaining k coefficients compare the methylation proportions between the conditions.

Dispersion estimation
In the small example above we did not estimate the dispersion but simply used a preset value. In practice, there will not only be replicates for each genomic locus but also many thousands of loci. The edgeR package estimates a dispersion for each locus from the replicate variability, using an empirical Bayes procedure to pool information across all the loci while estimating locus-specific dispersion estimates 33-35 . The full data example on mouse mammary epithelial cells will illustrate this.
The variance of the read counts is assumed to follow a quadratic mean-variance relationship 2 var( ) where φ g is the dispersion for locus g. The first term follows from sequencing variability and the second from biological variability 33 . The NB variance function captures the fact that larger counts have larger variances. The quadratic relationship captures the technical factors that affect variability, so that φ g itself reflects only the biological characteristics of each locus.

RRBS profiling of mammary epithelial cells
Aim of the study We now describe the main data set to be analyzed in this article. The epithelium of the mammary gland exists in a highly dynamic state, undergoing dramatic morphogenetic changes during puberty, pregnancy, lactation and regression 36 . Characterization of the lineage hierarchy of cells in the mammary epithelium is an important step toward understanding which cells are predisposed to oncogenesis. In this study, we profiled the methylation status of the two major functionally distinct epithelial compartments: basal and luminal cells.
The basal cells were further divided into those showing high or low expression of the surface marker Itga5 as part of our investigation of heterogeneity within the basal compartment. We carried out global RRBS DNA methylation assays on two biological replicates of each of the three cell populations to determine whether the epigenetic machinery played a potential role in (i) differentiation of luminal cells from basal and (ii) any compartmentalization of the basal cells associated with Itga5. Of particular interest is to relate transcriptional changes to methylation changes for the same genes.

Sample preparation
Inguinal mammary glands (minus lymph node) were harvested from FVB/N mice. All animal experiments were conducted using mice bred at and maintained in our animal facility according to the Walter and Eliza Hall Institute of Medical Research Animal Ethics Committee guidelines. Epithelial cells were suspended and fluorescence-activated cell sorting (FACS) was used to isolate basal and luminal cell populations 37 . Genomic DNA (gDNA) was extracted from freshly sorted cells using the Qiagen DNeasy kit. Around 25ng gDNA input was subjected to DNA methylation analysis by BS-seq using the Ovation RRBS Methyl-seq kit from NuGEN. The process includes MspI digestion of gDNA, sequencing adapter ligation, end repair, bisulfite conversion, and PCR amplification to produce the final sequencing library. The Qiagen EpiTect Bisulfite kit was used for bisulfite-mediated conversion of unmethylated cytosines.

Experimental design
There are three groups of samples: luminal population, Itga5-basal population and Itga5+ basal population. Two biological replicates were collected for each group. This experimental design is summarized in the table below.
The sequencing was carried out on the Illumina NextSeq 500 platform. About 30 million 75bp paired-end reads were generated for each sample.

Reading and annotating the methylation counts
Processing the BS-seq FASTQ files with Bismark The first step of the analysis is to map the sequence reads from the FASTQ files to the mouse genome and to perform methylation calls. This is the only step of the analysis that cannot currently be done in R.
Though many options are available, we used Bismark software (https://www.bioinformatics.babraham.ac.uk/ projects/bismark) to count the methylated and unmethylated reads at each genomic locus. Bismark was run using recommended default settings. We first ran trim_galore (https://www.bioinformatics.babraham. ac.uk/projects/trim_galore/) to remove adapters and to trim poor quality reads. Then Bismark version v0.13.0 was used to align the reads to the mouse mm10 genome using Bowtie2 22 . Finally, methylation calls were made using bismark_methylation_extractor.
The Bismark output consists of one coverage file for each sample. Readers wishing to reproduce the analysis presented in this article can download the coverage files produced by Bismark from http://bioinf.wehi. edu.au/ edgeR/F1000Research2017.

Reading in the data
The Bismark coverage files are just tab-delimited text files and so can be read into a dataframe in R using read. delim. Each of the files has the following format: In practice, we only need to read in columns V1, V2, V5 and V6. Column V3 is the same as V2, and the methylation proportion V4 is just a function of V5 and V6.
We now read in the Bismark coverage files for all the samples. The edgeR function readBismark2DGE reads all the files and collates the counts for all the sample into one data object: The edgeR package stores the counts and associated annotation in a simple list-based data object called a "DGEList". The count matrix is There is a row for each CpG locus found in any of the files. There columns of methylated and unmethylated counts for each sample. The chromosomes and genomic loci are stored in the annotation data.frame: Then we can subset the DGEList data object to remove the rows corresponding to the unwanted chromosomes: The option keep.lib.sizes=FALSE causes the library sizes to be recomputed.

Differential methylation analysis at CpG loci Filtering to remove low counts
We now turn to statistical analysis of differential methylation. Our first analysis will be for individual CpG loci.
CpG loci that have low coverage are removed prior to downstream analysis as they provide little information for assessing methylation levels. Filtering low-coverage CpGs also simplifies the subsequent analysis because of the reduction in data rows.
We sum up the counts of methylated and unmethylated reads to get the total read coverage at each CpG site for each sample: The analysis needs to be restricted to CpG sites that have enough coverage for the methylation level to be measurable in a meaningful way at that site. As a conservative rule of thumb, we require a CpG site to have a total count (both methylated and unmethylated) of at least 8 in every sample before it is considered in the study.
This filtering criterion could be relaxed somewhat in principle but the number of CpGs kept in the analysis is large enough for our purposes.
The DGEList object is subsetted to retain only the non-filtered loci: Again, the option keep.lib.sizes=FALSE causes the library sizes to be recomputed after the filtering. We generally recommend this, although the effect on the downstream analysis is usually small.

Normalization
A key difference between BS-seq and other sequencing data is that the pair of libraries holding the methylated and unmethylated reads for a particular sample are treated as a unit. To ensure that the methylated and unmethylated reads for the same sample are treated on the same scale, we need to set the library sizes to be equal for each pair of libraries. We set the library sizes for each sample to be the average of the total read counts for the methylated and unmethylated libraries: Other normalization methods developed for RNA-seq data, such as TMM 38 , are not required for BS-seq data.

Exploring differences between samples
In microarray methylation studies, a common measure of methylation level is the M-value, which is defined as M = log 2 {(Me + α)/(Un + α)} where Me and Un are the methylated and unmethylated intensities and α is some suitable offset to avoid taking logarithms of zero 39 . The M-value can be interpreted as the base2 logit transformation of the proportion of methylated signal at each locus.
We compute the corresponding methylation summary from the methylated and unmethylated counts.
Here M contains the empirical logit methylation level for each CpG site in each sample. We have used a prior count of 2 to avoid logarithms of zero. The exact value of the prior count is unimportant, but a value of 2 is common in other contexts such as RNA-seq.
Now we can generate a multi-dimensional scaling (MDS) plot to explore the overall differences between the methylation levels of the different samples ( Figure 1):

> plotMDS(M)
In this plot the distance between each pair of samples represents the average logit change between the samples for the top most differentially methylated CpG loci between that pair of samples. (We call this average the leading log-fold-change.) The two replicate samples from the luminal population (P6) are seen to be well separated from the four basal samples (populations P7 and P8), with a leading change of about 6 logit units.

Design matrix
One aim of this study is to identify differentially methylated (DM) loci between the different cell populations.
In edgeR, this can be done by fitting linear models under a specified design matrix and testing for corresponding coefficients or contrasts. A basic sample-level design matrix can be made as follows: The we expand this to the full design matrix modeling the sample and methylation effects: The first six columns represent the sample coverage effects. The last three columns represent the methylation levels (in logit units) in the three cell populations.

Dispersion estimation
With the design matrix specified, we can now proceed to the standard edgeR pipeline and analyze the data in the same way as for RNA-seq data. Similar to the RNA-seq data, the variability between biological replicates has also been observed in bisulfite sequencing data. This variability can be captured by the NB dispersion parameter under the generalized linear model (GLM) framework in edgeR.
The mean-dispersion relationship of BS-seq data has been studied in the past and no apparent mean-dispersion trend was observed 27 . This is also verified through our own practice. Therefore, we would not consider a meandependent dispersion trend as we normally would for RNA-seq data. A common dispersion estimate for all the loci, as well as an empirical Bayes moderated dispersion for each individual locus, can be obtained from the estimateDisp function in edgeR: > y <-estimateDisp(y, design, trend="none") > y$common.dispersion This returns a DGEList object with additional components (common.dispersion and tagwise.dispersion) added to hold the estimated dispersions. Here the estimation of trended dispersion has been turned off by setting trend="none". For this data, the estimated prior degrees of freedom (df) are infinite for all the loci, which implies all the CpG-wise dispersions are equal to the common dispersion. A BCV plot is often useful to visualize the dispersion estimates, but is not informative in this case.

Testing for differentially methylated CpG loci
We first fit NB GLMs for all the CpG loci using the glmFit function in edgeR.
> fit <-glmFit(y, design) Then we can proceed to testing for differentially methylated CpG sites between different populations. One of the most interesting comparisons is between the basal (P7 and P8) and luminal (P6) populations. The contrast corresponding to any specified comparison can be constructed conveniently using the makeContrasts function: The actual testing is performed using likelihood ratio tests (LRT) in edgeR: The top set of most differentially methylated (DM) CpG sites can be viewed with topTags: Here positive log-fold-changes represent CpG sites that have higher methylation level in the luminal population compared to the two basal populations. The Benjamini-Hochberg multiple testing correction is applied to control the false discovery rate (FDR).
The total number of DM CpG sites identified at an FDR of 5% can be shown with decideTests. There are in fact more than 50,000 differentially methylated CpGs in this comparison: > summary(decideTests(lrt)) 1*P6 -0.5*P7 -0.5*P8 Down 20302 NotSig 451589 Up 40111 The differential methylation results can be visualized with an MD plot (see Figure 2): The logFC of the methylation level for each CpG site is plotted against the average abundance in log2-CPM.
Significantly differentially methylated CpGs are highlighted.

Differential methylation by chromosome
We now explore overall methylation patterns by chromosome. To do this we make a list of index vectors, with each vector identifying the loci that belong to that chromosome:  This shows that methylation increases overall in committed luminal cells for every chromosome except for Chromosome X. This is consistent with the expectation that methylation is associated with silencing of genes not needed in a committed lineage. The column heading NGenes gives the number of CpG loci in each chromosome.

> ChrIndices <-list() > for (a in ChrNames) ChrIndices[[a]] <-which(y$genes$Chr==a)
The results for Chromosome X are interesting. The mixed P-value shows that the X chromosome is actually enriched for DM CpGs, but the direction of change is not consistent with hyper and hypomethylated CpGs balancing each other.

Differential methylation in gene promoters
Pre-defined gene promoters The majority of CpGs are methylated in mammals. On the other hand, unmethylated CpGs tend to group into clusters of CpG islands, which are often enriched in gene promoters. CpG methylation in promoter regions is often associated with silencing of transcription and gene expression 3 . Therefore it is of great biological interest to examine the methylation level within the gene promoter regions.
For simplicity, we define the promoter of a gene as the region from 2kb upstream to 1kb downstream of the transcription start site of that gene.

> InPromoter <-yall$genes$Distance >= -1000 & yall$genes$Distance <= 2000
We subset the CpGs to those contained in a promoter region: Summarizing counts in promoter regions One simple and effective way to conduct a gene-orientated analysis of the methylation changes is to collapse all the CpGs in each promotor into one locus, i.e., to compute the total number of methylated and unmethylated reads within each promoter. This strategy focuses on the aggregate methylation level within each promoter and allows us to test for an overall increase or decrease in methylation for each gene.
First we compute the total counts for each gene promoter: > ypr <-rowsum(yIP, yIP$genes$EntrezID, reorder=FALSE) > ypr$genes$EntrezID <-NULL The rowsum function here operated on the DGEList object yIP and produces a new DGEList object ypr with a row for each EntrezID. The integer matrix ypr$counts contains the total numbers of methylated and unmethylated CpGs observed within the promoter of each gene. Same as before, ypr$counts has 12 columns, two for each sample. The odd-numbered columns contain the numbers of methylated Cs, whereas the evennumbered columns contain the numbers of unmethylated Cs. The only difference is that each row of ypr$counts now represents a gene promoter instead of an individual CpG site.

Filtering to remove low counts
Filtering is performed in the same way as before. We sum up the read counts of both methylated and unmethylated Cs at each gene promoter within each sample. Since each row represents a 3,000-bps-wide promoter region that contains multiple CpG sites, we would expect less filtering than before.
> keeppr <-rowSums(Coverageprm >= 10) == 6 > Exploring differences between samples Same as before, we measure the methylation levels of gene promoter regions using M-values. A prior count of 2 is added to the calculation to avoid undefined values and to reduce the variability of M-values for gene promoters with low counts. Then a MDS plot is produced to examine the overall differences between the methylation levels of the different samples.

Dispersion estimation
We estimate the NB dispersions using the estimateDisp function in edgeR. For the same reason, we do not consider a mean-dependent dispersion trend as we normally would for RNA-seq data.

Testing for differential methylation in gene promoters
We first fit NB GLMs for all the gene promoters using glmFit.
> fitpr <-glmFit(ypr, design) Then we can proceed to testing for differential methylation in gene promoter regions between different populations. Suppose the comparison of interest is the same as before. The same contrast can be used for the testing.
> lrtpr <-glmLRT(fitpr, contrast=contr)  Here positive log-fold-changes represent gene promoters that have higher methylation level in the luminal population compared to the basal population. The Benjamini-Hochberg multiple testing correction is applied to control the false discovery rate (FDR). For future reference we make a dataframe of all the DM genes: > topME <-topTags(lrtpr, n=Inf, p=0.05)$table > dim(topME) [1] 1210 8 The differential methylation results can be visualized with an MD plot (see Figure 7):

RNA-seq profiles of mouse epithelium luminal and basal cells
To explore whether hypermethylation of promoter regions is associated with repressed gene expression, we relate the differential methylation results to differential expression results from RNA-seq for similar cell populations. The RNA-seq data used here is from a study of the epithelial cell lineage in the mouse mammary gland, in which the expression profiles were generated from basal stem-cell enriched cells and committed luminal cells in the mammary glands of virgin, pregnant and lactating mice 40 . The complete differential expression analysis of the data is described in Chen et al. 32 .
The RNA-seq data is stored as a DGEList object y_rna and saved in a RData file rna.RData. The object y_rna contains the count matrix, sample information, gene annotation, design matrix and dispersion estimates of the RNA-seq data. The gene filtering, normalization and dispersion estimation were performed in the same way as described in Chen et al. 32  We keep only the genes that are also included in our methylation analysis: > haveME <-row.names(y_rna) %in% row.names(ypr) > y_rna <-y_rna[haveME,] > dim(y_rna) [1] 13210 12 We assess differential expression between the luminal and basal virgin samples:

Gene set testing
The correlation in Figure 8 is convincing, but the above P-value computation assumes that genes are statistically independent of one another. We can perform a gene set test to get around this assumption.
First we make a dataframe of the logFCs of the DM genes, keeping only the genes for which we have RNA-seq results: > ME <-data.frame(ID=row.names(topME), logFC=topME$logFC, stringsAsFactors=FALSE) > inRNA <-ME$ID %in% row.names(y_rna) > ME <-ME[inRNA,] Then we test whether the DM genes are up or down-regulated using a "fry" gene set test. This analysis weights the DM genes by their methylation logFCs so that the test evaluates whether the expression changes are positively or negatively associated with the methylation changes: > fry(y_rna, index=ME, contrast=Contrastrna) NGenes Direction PValue PValue.Mixed set1 740 Down 1.28e-09 5e-11 The result Down in the Direction column indicates negative correlation between the methylation and expression changes. The small PValue confirms a highly significant result.

Discussion
This article has presented a complete start to finish analysis of an RRBS dataset from our own practice. The analysis demonstrates how BS-seq data can be analyzed using software designed for RNA-seq, thus benefiting from a large reservoir of highly-developed RNA-seq methodology. At first sight, BS-seq data is fundamentally different to RNAseq because of the need to focus on the proportion of methylated reads at each locus rather than on the absolute number of reads at each locus. The link is achieved by conducting the statistical inference conditional on the read coverage at each locus, and this in turn is achieved by including sample-specific locus effects in the linear model.
We have concentrated on analysing methylation changes for pre-defined genomic regions. At the highest resolution, conducted differential methylation tests for each distinct CpG site. Our main analysis was gene-level, whereby we aggregated counts over a putative promoter region around the TSS for each gene. This gives a "big picture" analysis and facilities a correlation of DM results with differential expression results from RNA-seq for the same cell populations. We also illustrated the use of gene set tests to examine whether there were overall methylation changes at the whole chromosome level.
The analysis presented here was designed for reduced representation BS-seq and was tuned to our own research interests. While the same analysis approach could be fruitfully applied to WGBS data, researchers with WGBS data may also want to discover DMRs de novo without the use of gene annotation or pre-specified genomic regions, something we haven't explored in this article. The edgeR approach presented here could in principle be extended to discover DMRs in a de novo fashion using similar methods to those developed for ChIP-seq data 41-43 .

Packages used
This workflow depends on various packages from version 3.7 of the Bioconductor project, running on R version 3.5.0 or higher. For all the code to work as presented, edgeR 3.22.2 or later is required. A complete list of the packages used for this workflow is shown below:

Data and software availability
All software packages used in this workflows are publicly available as part of Bioconductor 3.7. The data and analysis code used in this workflow are available from http://bioinf.wehi.edu.au/edgeR/F1000Research2017. The data and analysis code as at time of publication of this article has been archived at http://doi.org/10.5281/ zenodo.1052870 44 . software, to test for differential methylation from bisulfite-sequencing data, particularly edgeR reduced representation bisulfite-sequencing (RRBS). By leveraging the existing software, a edgeR popular tool in the analysis of diffential gene expression analysis from RNA-seq data, this method is immediately able to handle complex experimental designs and integrates with downstream analysis tools such as gene set tests provided by the software. The paper is well-written and I was able to limma reproduce the authors' analysis. It will be a useful workflow for people needing to develop an analysis of RRBS data.

Grant information
Like Simon Andrews, I initially struggled a little with some of the detail of the method itself. The method's elegance and power, like all those based on (generalized) linear models, is driven by careful formulation of the design matrix and the choice of contrasts. Necessarily, the design matrix for analysing bisulfite-sequencing data is more complex than that used to analyse RNA-seq data from an identical experimental design. As I know the authors are well aware, getting the design matrix and contrasts correct is 95% of the battle for most people analysing data with and . I will explain my concerns edgeR limma below (many are the same as raised by Simon in his review).

Main points
p4: The initial example has no replicates. Since the method is designed for "any BS-seq data that includes some replication", should this example include replicates? I appreciate the desire to keep the initial example simple (especially in light of my next comment). p4: I initially found the design matrix confusing. In fact, I had the same reaction/interpretation as Simon Andrews 4th comment: "In the small example description you say that A_MvsU estimates the log ratio for Sample1, but it wasn't clear to me why this would apply to only Sample 1 since the factor has a 1 against the meth count for both samples 1 and 2". I had to manually check a few quantities to convince myself, e.g., to rounding error, `coef(fit)[, "A_MvsU"]` is `logit((2 + prior.count) / (2 + prior.count + 12 + prior.count)`, where `prior.count = 0.125`. Because so much depends on constructing the appropriate design matrix, this description/section may warrant further explanation (e.g., comparing to some manually computed quantities). Like James MacDonald, although the code was clearly written, I was a little surprised that it didn't use more consistent integration with existing Bioconductor packages and data structures. To add to his example, almost all the work in the section 'Reading in the data' can be achieved with bsseq::read.bismark(fn)`, which will: read in an arbitrary number of Bismark `.cov.gz` files, appropriately combine samples with different sets of CpGs, and return a SummarizedExperiment -derived object (a instance) which could readily be used to construct the used in BSseq DGEList appropriately combine samples with different sets of CpGs, and return a SummarizedExperiment -derived object (a instance) which could readily be used to construct the used in BSseq DGEList the analysis. In my experience, loading the data and combining different sets of loci is a step fraught with danger of hard-to-track-down errors, so it may be better to advise workflow users to use a fairly well-tested function. Full disclosure: I am the author of `bsseq::read.bismark()`. p14: The aggregation of CpGs to promoters may lead to surprising results. An (extreme) example: the first half of a promoter is methylated in one condition and unmethylated in the other, and vice versa for the second half of the promoter. In aggregate over the promoter the proportion of methylated CpGs may be similar in both conditions, yet this promoter is clearly differentially methylated. I think a note encouraging workflow users to think carefully about their hypothesis when doing this form of aggregation is warranted.
Minor points p1: "The most commonly used technology of studying DNA methylation is bisulfite sequencing (BS-seq)". The Illumina 27k/450k/EPIC microarrays are the most commonly used 'genome-wide' assays for studying DNA methylation. However, (whole genome) BS-seq is arguably the gold standard genome-wide assay. p3: I think there's some confusion about CpGs and CpG islands (CGI). Approximately 0.9% of dinucleotides in the human genome (hg19) are CpGs, and approximately 0.7% of the genome is a CGI (using UCSC CGIs, which is not the only definition but perhaps the standard); see code below:``R library(BSgenome.Hsapiens.UCSC.hg19) hg19_size <-sum(as.numeric(seqlengths(BSgenome.Hsapiens.UCSC.hg19)[ paste0("chr", c(1:22, "X", "Y"))])) # CpGs on chr1-22,chrX,chrY in hg19 n_CpGs <-Reduce(sum, bsapply(BSParams = new("BSParams", X = BSgenome.Hsapiens.UCSC.hg19, FUN = countPattern, exclude = c("M", "_")), pattern = "CG")) 100 * n_CpGs / hg19_size # CGIs in hg19 library(rtracklayer) my_session <-browserSession("UCSC") genome(my_session) <-"hg19" cgi <-track(ucscTableQuery(my_session, track = "cpgIslandExt")) sum(width(cgi)) / hg19_sizè``p 3: Possible type, "with large genome" a p3: "WGBS is more suitable for studies where all CpG islands or promoters across the entire genome are of interest." Might also add 'distal regulatory elements' and CG-poor regions (RRBS targets CG-rich regions of the genome). p3: (implemented in ) doesn't use Empirical Bayes although it does use for BSmooth bsseq limma linear regression p4: Missing a `library(edgeR)` in order for the code to work p4: There's an extra parenthesis at the end of line 2 when constructing `dimnames(counts)p 4: The authors note that the method is "especially appropriate for RRBS data". Is the main challenge for running on WGBS data that of computational resources? p5: Typo, "mythlyated" should be "methylated" p5: Typo, "mythlyated" should be "methylated" Table 1: Condition should be 'A' or 'B' instead of '1' or '2' p6: Was Bowtie1 or Bowtie2 used as the Bismark backend for the mouse data? p8: The filtering step removes almost 90% of CpGs. Is this unavoidable, e.g., due to low sequencing coverage of these samples, or might the filtering be relaxed? Figure 1: Any thoughts for why the P8_6 sample is rather separated from the other Basal samples along dim2 of the MDS plot? Figure 2: What is the meaning of 'average abundance of each CpG site'? Is 'abundance' interpretable as 'sequencing depth'? p16: Possible typo, "Suppose the comparison of interest is same as before" the p22: In the DNA methylation literature, 'up-methylated' is typically called 'hypermethylated' and 'down-methylated' is typically called 'hypomethylated'. Thank you for your thoughtful and positive comments on our article. We particularly wish to thank you for your comments on the introductory material, which we have rewritten to take account of all your suggestions. We have specified the version of Bowtie, adopted the "hyper" and "hypo" terminology and made all the other corrections. In this revision we have added a lot of new material to the expository section, including an example with replicates, while retaining the first example without replicates. The purpose of the first example is to introduce the key mathematical idea that underlies our article, which is the use of linear model terms to account for total read coverage at each CpG, and to relate this to the more familiar beta-binomial modeling approach. This idea is already manifest in the very small dataset without replicates, and so that is the most fundamental context in which to exposit the idea. The need for replicates, and so that is the most fundamental context in which to exposit the idea. The need for replicates relates to dispersion estimation, which is a separate issue. For pedagogic purposes it is best to isolate these concepts separately. We now demonstrate in the first example how the logFC is computed in complete arithmetic detail. We have switched the toy examples to use the group-mean parametrization instead of the treatment-parametrization so that the computation of logFC is explicit rather than implicit. See also our reply to Simon Andrews. We have introduced a new function modelMatrixMeth() to automate the construction of the design matrix. From a user point of view, the design matrix is now made exactly as it would be for an RNA-seq analysis, with the sample effects added automatically. We have also added a new section "Analogy with paired-samples expression analyses" to show that the entire linear modeling approach is equivalent to an (already well understood) type of RNA-seq experiment. The purpose of our article is to show that methylation data can be analyzed exactly like RNA-seq and we use the same packages and data structures as our earlier RNA-seq workflow article (Chen et al, 2017). Adding extra packages and data structures would complicate rather than help our workflow. We are aware of your read.bismark() function in the bbseq package. When run on our data, read.bismark() produces a number of alarming warnings suggesting that incorrect genome annotation has been used for some of the samples. These warnings are false alarms, but would worry users, so we would have to either suppress them or explain them to readers, were we to adopt read.bismark. We are also aware of the readBismark() function in the BiSeq package and the methRead() function in the methylKit package. It seems that every Bioconductor package that analyses BS-seq has its own read function. The functions read the methylation data into a binomial format, whereas the whole point of our article is to show that the methylated and unmethylated counts can be analysed as they are instead of doing a binomial analysis. All produce data structures that are incompatible with edgeR. We could wrangle the data formats back to a DGEList, but the end result would be code that is no shorter or simpler than our original plain base R code. In fact the code would be harder to follow because of the extra packages and data structures involved. It seems unfortunate that the Bioconductor methylation packages all create their own read functions from scratch in order to create their own specialist S4 objects. Each defines a complex S4 class, unique to that one package. If Bioconductor provided a fast, light-weight, base-R style function to read Bismark files, then all the packages could use that as building block instead of re-writing from scratch. In this revision, we have decided to create a new function readBismark2DGE() to read the Bismark files directly into a DGEList object, substantially simplifying our workflow for users. The new function is faster than any existing options and over 3 times as fast as read.bismark() on the workflow data. The new function creates simple list output with minimal processing and so serves the purpose of a light-weight base-R style function. You comment on aggregation over promoters and guidance about hypotheses. We have added more text to clarify that aggregation of counts over gene promoters tests for overall changes only, but we think that users would have expected this already and shouldn't be surprised. The thought example you give is indeed extreme and does not match real that we have seen. Even if such a phenomenon could occur, our workflow includes a CpG-level analysis, so readers can easily check for CpGs changing in different directions within a promoter if that is of interest to them. We did this check ourselves and found that for 99% of genes whose promoters contain any DM CpGs, all the DM CpGs change in the same direction. In almost all the remaining cases, only one CpG changes in the "wrong" direction. The only exception is Mir1298 on Chromosome X that has two CpGs changing in the other direction to the majority.

Are sufficient details provided to allow replication of the method development and its use by others? Yes
Even if there were opposing effects within a few gene promoters, it wouldn't invalidate our Even if there were opposing effects within a few gene promoters, it wouldn't invalidate our aggregate analysis. It would simply mean that exceptional genes would have to analyzed and interpreted individually to understand the fully story, and we see no problem with that. We think there is a role for both aggregate analyses and high-resolution analyses of DNA methylation. However we think that aggregate analyses are under-used in the literature and should be encouraged more than at present. We think that genomic analysts sometimes try to accommodate too much complexity, leading to analyses that lack statistical power. This in turn leads to the danger of cherry picking, possibly over-interpreting results that might not be reproducible. Aggregate analyses ignore some of the complexity, but can lead to more reliable and interpretable "big picture" messages. In this article, our promoter region aggregated analysis leads to a clean and reliable correlation of methylation with RNA-seq results, a correlation that is more difficult to achieve by other means. Our remark that our method is especially appropriate for RRBS is prompted by the thought that researchers engaged in a WGBS study may have slightly different aims, and this goes back to your comment about hypotheses. In particular, it may be important in a WGBS analysis to discover DMRs in a fashion, pooling results for consecutive genomic loci without the use of gene de novo annotation. We have added remarks in the Introduction and the Discussion to clarify our aims and to briefly indicate possible extensions. There is no problem with computational resources -our pipeline is reasonably efficient. In this revision, we have relaxed the filtering threshold to 8 instead of 10 reads per CpG, although it makes no material difference to our analysis. While the filtering does indeed remove 85% of the CpGs, it keeps more than half the counts. The number of CpGs retained in the analysis is plenty for our purposes. In this revision, we form the MDS plot in a slightly better way and P8_6 is no longer an outlier. In Figure 2, the meaning of "average abundance" is explained by the axis label as log count per million, i.e., coverage for that CpG divided by sequencing depth. In effect, it is the coverage for each CpG relative to other CpGs. It is independent of sequencing depth. This is primarily a software pipeline article, showing how to use the Bioconductor edgeR package to analyze RRBS data, but to a certain extent is also a methods paper, as to my knowledge this is the first proposal for directly analyzing count data rather than converting to either ratios and using a beta-binomial, or to logits and using conventional linear modeling. This is an interesting idea, and should be explored further, but for this manuscript the main goal is to present the software pipeline.
The code is clearly written, and as straightforward as one could expect for a relatively complex analysis. However, I would prefer to see more consistent integration with other Bioconductor packages. In particular, when reading in the raw data, the authors use a clever trick to account for the fact that not all samples have reads for the same genomic positions. This step could just as easily be accomplished using the Bioconductor GenomicRanges package, which is intended for manipulating genomic data. In fact, the authors use GenomicRanges later in the pipeline to subset the methylation data to just gene promoter sites, so it would be more natural to start with a GRanges if you will need one later anyway.
Otherwise this is a good article that clearly shows how one could use an innovative method to analyze RRBS data using the edgeR package. Thank you for your positive report on our article. You're right, there is a novel methods proposal in our article, which we will seek to justify in more detail in a methodological article elsewhere, but the current article concentrates on the software pipeline. In this revision, we have introduced new functions readBismark2DGE() and nearestTSS() to simplify the workflow. These functions allow us to eliminate from the workflow any code chunks that involve substantial programming, making the workflow considerably easier for users. We felt the original code chunk to find CpGs in promoters in particular was too complex for the purposes of our article and assumed too much user knowledge.
I have used GenomicRanges for other projects ( ) but it https://f1000research.com/articles/4-1080 I have used GenomicRanges for other projects ( ) but it https://f1000research.com/articles/4-1080 doesn't help with the current workflow. The file reading code would be similar in length and speed with or without GenomicRanges, so we tend to give preference to R base. In any case, we have eliminated that code entirely, as explained in our discussion with Peter Hickey. The new nearestTSS() function is more informative than the old code, as well as being fast. It has the added advantage of using the Bioconductor organism package, which has more up-to-date TSS annotation than the transcript package.
No competing interests were disclosed. Competing Interests: 07  present an interesting re-application of the EdgeR analysis package to the analysis of et al bisulphite sequencing data. The method they propose would utilise the existing negative biomial models within EdgeR and would potentially provide the power which comes with the linear model framework to bisulphite data. The method described requires no changes to EdgeR itself, and merely describes a suitable formulation of design matrix to allow this to be applied to bisulphite data.
The article is generally well written and the authors go to great lengths to break down and describe the method. They also provide a page from which all of the underlying data and code can be obtained and I was able to reproduce the results, and independently verify them in a parallel analysis.
The main thing which I struggled with was some of the detail in the description of the method itself. There were some parts which I wasn't clear on, and some nomenclature which didn't help in understanding the explanation. I'll try to lay out my concerns below: 1) In the small example I completely understand that the authors wanted to keep this as simple as possible, but it might have helped to have had 2 samples per condition so that the full complexity of the method is visible.
2) There is a typo in the code for the small example so it doesn't run as is. The list function on line 2 has an extra bracket at the end.
3) The nomenclature in the small example is inconsistent. You have samples 1 and 2, but (in the table) also conditions 1 and 2, but in the code the conditions are A and B. If you had Samples 1,2,3,4 in conditions A and B this might help to alleviate some of the confusion. 4) In the small example description you say that A_MvsU estimates the log ratio for Sample1, but it wasn't clear to me why this would apply to only Sample 1 since the factor has a 1 against the meth count for both samples 1 and 2 In the expanded examples there were also some points on which I wasn't clear. 5) You calculate a single dispersion parameter for all data points and say that in contrast to RNA-Seq 5) You calculate a single dispersion parameter for all data points and say that in contrast to RNA-Seq there is no global trend to follow. It wasn't clear to be exactly why this is since read count and methlyation level would all affect the dispersion -is it simply because these factors are explicitly accounted for in the linear model? 6) In the design matrix for the RRBS it wasn't clear why the first column was all 1s, whereas the rest obviously matched the condition from which they came. This also contrasted with the simple example where the structure wasn't like this. Is this because you were comparing both P7 and P8 to P6? 7) I think this is possibly the same thing as point 4, but you say that the Me column represents the methylation level in P6, but again this highlights the methylated values in all samples, so why only P6?
For the final results obtained it would have been nice to show the general level of concordance with running the same analysis through one of the beta-distribution models to either show general agreement, or to generally explain any major differences.

Minor points:
In the introduction you say that "40% of mammalian genes and 70% of human genes have CpG islands enriched in their promoter regions". Enriched probably isn't the right word to use (or you need to say that CpGs are enriched rather than islands). The difference between 'humans' and 'mammals' is also somewhat contentious -non-human mammals certainly have weaker CpG islands which get missed by CpG island prediction tools, but for example in mouse Illingworth et al showed that if you use CpG binding protein ChIP that you can see about the same number of islands in both species.
It's also not really fair to say that CpG methylation in promoters is "generally" associated with repression of transcription. There is a categorical expression level shift associated with the presence/absense of CpG islands, but you can make a Dnmt1 knockout which removes pretty much all methlyation from the genome and for the vast majority of genes their transcription is completely unaffected.

Is the description of the method technically sound? Partly
Are sufficient details provided to allow replication of the method development and its use by others? Yes Thank you for your thoughtful comments on our article. We have rewritten the Introduction to take account of your comments and have fixed all typos. We have also expanded the expository section on the NB linear modeling approach to five pages instead of one. The purpose of the first toy example (Table 1) is to introduce the key mathematical idea that underlies our article, which is the use of linear model terms to account for read coverage at each CpG, and to relate this to the more familiar beta-binomial modeling approach. This idea is already manifest in the very small dataset without replicates, making the one-sample example the most fundamental and effective context in which to introduce the idea. We agree however that the condition labels in Table 1 should be A and B. Following your remarks, we have included a second example with 2 samples per condition ( Table  2). Keep in mind though that the edgeR approach to methylation however is very powerful and flexible and its full complexity is still far from being visible in a single locus example like this. To see more of the features in play, one has to work through the full data example. We have however taken the opportunity to preview more ideas in the expanded linear modeling introduction, including dispersion estimation and the use of contrasts. You asked about the A_MvsU coefficient computation. This is a standard linear modeling result in R but one that does confuse a lot of people when they see it for the first time. You will see analogous results explained in the limma and edgeR User's Guides. Anyway, we have rewritten the first example in a different way so that the coefficients are defined more explicitly and their calculation is explained in full arithmetic detail. The original implicit formulation of the logFC is now presented in the section "Another way to make the design matrix". Consider the last two columns of the design matrix in that section, which are now called "Intercept" and "ConditionB", and write β and β for the corresponding two coefficients. The interpretation of these coefficients becomes clear by multiplying the coefficients by the design matrix. For samples corresponding to condition A, the methylation level is represented by β alone because the last column of the design matrix is zero (rows 1 and 3). When condition B is applied, the methylation is modelled by β + β (rows 5 and 7). So it is apparent that β must represent the methylation status for condition A while β must be the difference in methylation between B and A. Users don't need to follow this mathematics, we just present it here for completeness. You also asked about the first column of the design matrix. It doesn't matter whether the first column is all 1s or not, as long as the leading columns of the design matrix span the sample effects. We have added a small example on pages 3-4 to demonstrate this. The first few columns of the design matrix model the sample coverages, and these columns are quite independent of the remaining columns which model the treatment conditions. We can compare the treatment remaining columns which model the treatment conditions. We can compare the treatment conditions in any way we wish regardless of how the design matrix is parametrized. We have added some material to the introduction to try to explain this but, in any case, it isn't something that a user needs to worry about. In this revision, we have introduced a new function modelMatrixMeth to create the design matrix automatically, so users now only need to focus on the treatment conditions, not on how to adjust for the sample effects.
There is no reason why the dispersion should be a function of read count or of methylation level. The purpose of the NB statistical model is to capture the technical mean-variance trends so that the dispersion parameter can be interpreted independently of these things. We have added a sub-section on dispersion estimation to page 7 to try to make this clear.
Regarding concordance with other methylation analysis software, we do not know of any beta-binomial modeling software that is able to conduct an equivalent analysis to that presented in our article. We have formed a contrast between luminal cells and the two basal populations, and no other software can do that as far as we know. This is an example of the extra flexibility provided by our approach. The DSS Bioconductor package is documented to be able to form general contrasts between treatment conditions, but this is not effective in practice because it does not distinguish hyper from hypo methylation for a contrast. Our experience is that edgeR gives better DM results than beta-binomial software but it would be unfair for us to claim that in our article without undertaking a full comparison study with rigorous benchmarking. Our plan is to publish such a comparison elsewhere.
No competing interests were disclosed.

Competing Interests:
The benefits of publishing with F1000Research: Your article is published within days, with no editorial bias You can publish traditional articles, null/negative results, case reports, data notes and more The peer review process is transparent and collaborative Your article is indexed in PubMed after passing peer review Dedicated customer support at every stage For pre-submission enquiries, contact research@f1000.com