Transcriptomic dynamics of breast cancer progression in the MMTV-PyMT mouse model
© The Author(s). 2017
Received: 18 June 2016
Accepted: 7 February 2017
Published: 17 February 2017
Malignant breast cancer with complex molecular mechanisms of progression and metastasis remains a leading cause of death in women. To improve diagnosis and drug development, it is critical to identify panels of genes and molecular pathways involved in tumor progression and malignant transition. Using the PyMT mouse, a genetically engineered mouse model that has been widely used to study human breast cancer, we profiled and analyzed gene expression from four distinct stages of tumor progression (hyperplasia, adenoma/MIN, early carcinoma and late carcinoma) during which malignant transition occurs.
We found remarkable expression similarity among the four stages, meaning genes altered in the later stages showed trace in the beginning of tumor progression. We identified a large number of differentially expressed genes in PyMT samples of all stages compared with normal mammary glands, enriched in cancer-related pathways. Using co-expression networks, we found panels of genes as signature modules with some hub genes that predict metastatic risk. Time-course analysis revealed genes with expression transition when shifting to malignant stages. These may provide additional insight into the molecular mechanisms beyond pathways.
Thus, in this study, our various analyses with the PyMT mouse model shed new light on transcriptomic dynamics during breast cancer malignant progression.
KeywordsRNA-sequencing PyMT mouse model Breast cancer
Breast cancer is a globally prevalent disease and a leading cause of cancer-related mortality among women of all ages . In 2008 there were 1.38 million new cases worldwide, and the disease caused ~458,400 deaths . The number increased to 1.67 million in 2012, accounting for 11.9% of all cancer new cases . A malignancy with a multistep pathological process, breast cancer in humans starts with the premalignant atypical ductal hyperplasia (ADH), followed by ductal carcinoma in situ (DCIS) and subsequent malignant invasive ductal carcinoma (IDC) . Patients surviving the primary tumors often die of carcinoma-culminated metastasis . Despite widely recognized evidence that ADH and DCIS are precursors of IDC, few biomarkers identified from the early stages can explain and predict tumor progression. Many genes have been shown to contribute to breast cancer development , but the molecular mechanisms of its progression remains largely unknown, which greatly limits our abilities for early diagnosis and treatment of breast cancer patients with metastasis risk .
Transgenic mouse models have been widely used to study breast cancer, and the PyMT mouse model is one of them [7–9]. Expression of the oncoprotein, polyoma middle T (PyMT) antigen from mouse polyoma virus, is under the control of the mouse mammary tumor virus (MMTV) long terminal repeats (LTR) and is restricted to mammary epithelia . By stimulating multiple signaling including Shc and PI3-kinase, the membrane scaffold protein PyMT activates MAPK and PI3K pathways that function in cell proliferation and survival [11, 12]. Sharing both morphological and transcriptional features with the human disease  and resembling the human luminal B subtype of breast cancer on gene expression profiles , the MMTV-PyMT transgenic mice provide us a reliable animal model for breast cancer progression. The primary tumors developed in this mouse model go through four stereotypical stages of cancer progression – hyperplasia, adenoma/mammary intraepithelial neoplasia (MIN), early and late carcinoma – while progress from pre-malignancy to malignancy.
Because most previous studies of PyMT mice focused only on the carcinoma stages, little is known about gene expression alterations in the early stages as well as their impacts to the later stages. In this study, we examined gene expression dynamics in the full range of breast tumor development from hyperplasia to late carcinoma. Using RNA sequencing instead of microarray, our data had a wider dynamic range and a higher sensitivity to better detect differentially or lowly expressed genes. Moreover, to go beyond single gene inspected by earlier studies, we explored biological networks to learn the connections and interplays among genes. Networks are powerful in interpreting the underlying mechanisms of diseases by revealing disease modules, which are groups of highly connected genes or gene products . In this study, we carried out differential gene expression profiling, time-course analysis and network-based gene screening to identify candidate genes that may contribute to breast cancer progression. We found that many genes differentially expressed in the late carcinoma stage initiated the expression alteration at the hyperplasia stage. We also identified genes with disrupted expression during the transition from premalignance to malignance. Last, we found gene modules that co-expressed in tumors with hub genes predict future metastasis. Thus, by proposing novel candidate oncogenes that may promote tumor progression and malignant transition, our study helps to find genes as potential biomarkers and drug targets for breast cancer treatment.
Animals and tissue collection
This study of cancer in mice was approved by the Institutional Animal Care and Use Committee (IACUC) of Albert Einstein College of Medicine. All procedures involving mice were conducted in accordance with the National Institutes of Health guidelines concerning the use and care of experimental animals. Male PyMT mice (FVB/N-Tg(MMTV-PyVT)634Mul/J mice, Stock Number: 002374, the Jackson Laboratory) were randomly bred with homozygous FVB females to obtain F1 female mice (PyMT mice hereafter) heterozygous for the PyMT transgene; they developed breast cancer and were used as cases. Homozygous FVB females were used as controls.
We selected four time points corresponding to different tumor progression stages: hyperplasia at week 6, adenoma/MIN at week 8, early carcinoma at week 10, and late carcinoma at week 12 . At each time point, mammary tumors and normal mammary glands were collected from three PyMT mice and three age-matched FVB controls, respectively. We snap froze all samples and kept them at –80°C. We had a pathologist examined the morphology of our PyMT tumor samples and observed typical features of various tumor developmental stages (Additional file 1: Figure S1A). We also estimated the percentage of infiltrating cells as the sample purity. By H&E staining, the pathology report showed less than 5% stromal and muscle cells, with rare inflammatory cells in carcinoma samples. In general, we had greater than 90% tumor cells in our carcinoma samples, which surpassed the Cancer Genome Atlas (TCGA) criteria (over 60% tumor cells for human tumor samples) [14, 15].
Directional RNA sequencing
Total RNAs were extracted from frozen samples using the miRNeasy mini kit (Qiagen) according to the manufacturer’s protocol. Agilent 2100 Bioanalyzer was used to check RNA quality. The total RNA was treated with DNaseI, depleted of ribosomal RNA with Ribo-Zero Magnetic Gold Kit (Epicentre), followed by ethanol precipitation. Next, RNA was converted to cDNA using SuperScript III First-Strand Synthesis Kit (Invitrogen) with 80ng random hexamers and 50mM oligo-dT and subsequently ethanol precipitated. Single-stranded cDNA was converted to dsDNA by DNA polymerase I while incorporating dU/VTPs (10mM). Samples were fragmented to 200–300bp using Covaris. After fragmentation, samples were purified using the MinElute PCR purification kit (Qiagen). Fragmented samples underwent standard end-repair, dA-tailing and adapter ligation using Illumina TruSeq adapters for multiplexing. Adapter-ligated cDNA was treated with uracil-DNA glycosylase followed by enrichment PCR using Kapa reagents for 14 cycles. Libraries were size selected for 150–600 bp on a 2% low-melt ultra low-range agarose gel stained with SYBR Gold (Invitrogen) to eliminate adaptor dimers. Purified libraries were used to sequence on Hiseq2500 according to standard protocols. The PyMT RNA-sequencing data is available in the Gene Expression Omnibus (GEO) database as GSE76772.
Sequence data were preprocessed by WASP 3.0 , an in-house pipeline. FastQC  was used for reads quality control. The raw FASTQ files were trimmed for adapter sequences using quart. Then GSNAP v2012-07-02  was used to align reads to mm9 reference genome with default settings. Gene counts were given by HTSeq v0.5.3p3 .
All statistical analyses were carried out with R v3.0.1 . After removing transcriptionally inactive genes (read count per million < 1 in more than half of the samples) from raw RNA-sequencing gene counts, we got high confident gene counts, which were then normalized by DESeq . The PCA analysis used method adopted from DESeq. The R package edgeR v3.4.2  was used to perform statistical analysis on gene counts and to detect differentially expressed genes. Differentially expressed genes (DEGs) at each stage were analyzed for enrichment on REACTOME  and KEGG  pathways using the R package GOseq , which corrects bias owing to gene length and expression variability. To reduce noise and redundancy, we used GSEA (Gene Set Enrichment Analysis)  to investigate hallmark gene sets (“Hallmark gene sets summarize and represent specific well-defined biological states or processes and display coherent expression”, as defined by GSEA).
Weighted gene co-expression network analysis
We conducted the weighted gene co-expression network analysis (WGCNA) using the R package WGCNA  to identify co-expressed genes. We used the normalized expression values from aforementioned analysis as input to construct a signed gene co-expression network. Then we looked at module preservation to find the most robust and generalized modules using a composite Z-summary statistic. For highly preserved modules, we calculated the module eigengene (the first principal component). Labels that denote the disease stages (normal, hyperplasia, adenoma/MIN, early carcinoma and late carcinoma) were permuted. We permuted the labels 10,000 times, and used an R package coin to test the statistical significance of Spearman’s correlation coefficient between tumor progression stages and module eigengenes. The results were exported to VisANT  for network visualization.
Time-course data analysis
To model the transitional changes in gene expression, we calculated a log2 ratio change (LRC) of each gene between stages i + 1 and i using results from the previous DEGs analysis: LRC = log2R(i+1) – log2R(i) (i = 1, 2, 3; stage 1 as hyperplasia). After normalizing LRCs to z-scores, we categorized the transitional changes into positive (+1), negative (–1), and constant (0) transition statuses by their deviance from the mean. Since there are three transitions and three transition statuses, altogether we have 33 = 27 possible groups of transition patterns.
To find gene expression trends along tumor progression, we selected the top 1,000 genes with highest PCA loadings, and divided the means of normalized expression at week 8, 10 and 12 by that of week 6 to get relative expression ratios. We then used k-means clustering to group genes based on the ratio. We plotted the number of clusters against the total within cluster SS (sums of squares) and selected the cluster value at the first inflection point, which is 5.
Quantitative real-time PCR (qPCR)
The qPCR was performed using the 2−ΔΔCt method. The RNA was converted to cDNA using SuperScript III First-Strand Synthesis Kit (Invitrogen), followed by reactions using ABI StepOnePlus Real-Time PCR System instrument (Applied Biosystems). Each 10μl reaction included 5μl Taqman Fast Master Mix, 2.5μl RNAse-free water, 2μl cDNA (10 nM) and 0.5μl Taqman primer set. Actb was used as the internal control gene to give genes tested a relative fold change using the 2−ΔΔCt method. The qPCR primers are by Taqman Gene Expression Assays, catalog number 4448892.
Differentially expressed genes are enriched in cancer-related pathways and E2F targets
Differentially expressed genes at four stages
We found genes in P53 signaling pathway with increased expression at early carcinoma stage (Fig. 2a), for example Rrm2, which functions in DNA repair and damage prevention. Rrm2 encodes a ribonucleotide reductase and was reported as an indicator of tamoxifen resistant in luminal patients as well as decreased survival in patients of all breast cancer subtypes . In our study, the expression of Rrm2 increased throughout the breast cancer development in PyMT mice (log2 fold change in hyperplasia, adenoma/MIN, early carcinoma, late carcinoma is 2.5, 4.0, 2.7, and 2.8, respectively) and may serve as potential biomarker for diagnosis as well as drug target.
In addition to terms related to cell cycle and proliferation, we also observed other significant enrichment such as PPAR (peroxisome proliferator-activated receptor) signaling (KEGG:03320) and drug metabolism- cytochrome P450 (KEGG:00982) (Fig. 2b). In late carcinoma, Pparg was among the 26 down-regulated PPAR pathway genes (Cd36, Sorbs1, Pparg, Angptl4, Lpl, Fabp4, Acox1, Acadm, Me1, Cpt2, Slc27a1, Scp2, Adipoq, Acsl1, Fabp3, Ppara, Plin1, Cyp27a1, Acox3, Pltp, Pck1, Ilk, Aqp7, Rxrg, Cpt1c, and Nr1h3); Cyp2e1 was among the 17 down-regulated drug metabolism-related cytochrome P450 genes (Mgst3, Ugt1a9, Adh1c, Fmo1, Gstz1, Aox1, Cyp2e1, Gsta3, Gstt2, Maob, Cyp2d6, Fmo2, Gstm1, Ugt1a5, Ugt1a3, Gstm2, and Aldh3b2).
Hub genes from co-expression network modules predict metastasis in human breast cancer datasets
WGCNA modules with high preservation
GO enrichment (FDR < 0.05)
KEGG pathways enrichment (FDR < 0.05)
Oxidative phosphorylation, TCA cycle
Extracellular matrix, cell adhesion
Focal adhesion, ECM-receptor interaction
Metabolism pathways, insulin pathway
Cell cycle, M phase, mitotic
DNA replication, cell cycle
The PyMT mice provide a very aggressive breast cancer model with a metastatic rate of over 90%. We hypothesized that genes whose expressions correlate with tumor progression may mark more aggressive tumors and thus can be used to predict metastasis risk. To test this hypothesis and to translate our model-based findings to the human disease, we assessed the impact of hub genes on distant metastasis-free survival (DMFS) using the Kaplan-Meier analysis by querying the Gene expression-based Outcome for Breast cancer Online (GOBO) . Based on their high or low expression levels, ten hub genes – H2AFX, KIF11, MCM6, MKI67, NCAPD2, PRC1, SMC2, TACC3, TOP2A, TPX2 – separated breast cancer patients into high- and low-risk groups, which showed statistically significant difference in DMFS (Additional file 4: Figure S3). Using the ten genes together as a panel also gave a significant prediction but showed no increase on predicting power. In the purple module, among top 15 genes that we selected as module hubs, 12 (PEX19, DHDH, ACSL1, ALDH6A1, PNPLA2, PPP2R5A, IVD, MME, ADIPOR2, ALAS1, CDO1, and CHPT1) are non-cell cycle genes and as a gene panel can separate patients with different risks by Kaplan-Meier analysis (p-value < 0.05), while the other three cell cycle-related genes (AKT2, INSR, and YWHAG) cannot. Hubs from the other two models showed no significant ability for separation.
Note that many hub genes as E2F targets are involved in tumorigenesis. As an E2F target, Dnmt1 (DNA methytrasferase 1) showed differential expression throughout tumor progression (expression fold changes in PyMT tumors vs. FVB controls during hyperplasia, adenoma/MIN, early carcinoma, and late carcinoma are 2.0, 2,5, 2.3, and 3.2 (FDR < 0.05), respectively). DNMT1 is crucial for the maintenance of DNA methylation, which causes gene silencing and is known to regulate many cancer-related genes including many tumor suppressors [35, 36]. These gene expression changes suggest a global methylation deregulation and a potential key role of DNA methylation in breast cancer. In addition, we have carried out an enhanced reduced representation bisulphite sequencing (ERRBS) of the same PyMT mouse model to examine DNA methylathion changes during its cancer progression process in another study to be published. Mki67 (Ki67) is a widely recognized marker for cell proliferation and is associated with breast cancer prognosis statistically . Top2a codes for topoisomerase II Alpha, a key enzyme in DNA replication that regulates gene expression and cell division and whose amplification is a predictor for anthracycline treatment in breast cancer . Another hub gene Tacc3 plays roles in cell cycle, immune system development, and microtubule cytoskeleton organization. An up-regulation of Tacc3 was observed in breast cancer and it was suggested that Tacc3 might be a deregulator of DNA damage response and a predictor of survival for breast cancer patients . Mcm6 and Ncapd2 are E2F targets as well, and little is known about their roles in breast tumor progression. Mcm6 is one of the mini-chromosome maintenance proteins and has been suggested as a prognostic marker in melanoma . In the black module, we found Mcm2-6 and Mcm9. The Mcm2-7 and Mcm8-9 form protein complexes respectively to function in DNA replication initiation and DNA recombination repair [41, 42]. An early study showed a strong association between the over expression of Mcm2-7 and the short survival in a breast cancer cohort  and suggested their roles in breast tumor progression. Here we added an additional layer of evidence. Ncapd2 (non-SMC condensing I complex subunit D2) is involved in mitotic sister chromosome segregation  without notable involvement in breast cancer. We think its alteration may represent elevated cell mitosis in tumors like that of Mki67.
The green module showed a slight negative correlation between eigengene and tumor progression compared to controls. 20 genes from the TCA cycle pathway (Acly, Aco2, Csl, Cs, Dlat, Dlst, Dld, Fh1, Idh2, Idh3a, Idh3b, Idh3g, Mdh1, Ogdh, Pdhb, Pdha1, Sdha, Sdhd, Suclg1, and Sucla2) were found in the green module. This enrichment indicates a significant involvement of energy metabolism in tumor progression. In the purple module, we found Eif4ebp1 with reduced expression in tumor. This gene is a eukaryotic translation initiation factor biding protein that functions in the insulin-signaling pathway. It inhibits the oncogene Eif4e, whose overexpression leads to human hepatocellular carcinoma development together with Ras activation .
Time-course analysis identifies expression patterns of tumor progression and malignant transition
As a compliment to the transition analysis, we explored the continuous monotonic pattern in tumor samples to study genes with possible tumor maintenance functions (Methods). We first isolated 1,000 genes with highest PCA loading and divided the mean of normalized expression at each time point by that of week 6 to obtain the relative expression ratio. Using k-means clustering, we then separated these genes into five groups with three major ones (Fig. 4c). Genes in red group mainly function in cell adhesion (GO: 0007155) and biological adhesion (GO: 0022610), maintained high expression as tumor progress. Genes in the blue group play roles in the immune system such as antigen processing and presentation, showed decreased expression especially at early carcinoma stage. Genes in the grey group, enriched in response to external stimulus (GO: 0009605), showed no obvious change as tumor progressed (Fig. 4d).
Novel genes involved in PyMT mouse model of breast cancer
Using R package RISmed, we mined PubMed for the functional relevance of the newly discovered genes in tumor development. We used RISmed to extract bibliographic data from PubMed, and examined if the official symbols of genes that we are interested in concurred with the key word “breast cancer” in article titles or abstracts. If so, such genes are considered to have functions related to breast cancer. Of the 4,391 DEGs, about half have no reported roles in breast cancer. Of the 79 E2F targets we identified in late carcinoma, 14 genes (H2afz, Spc24, Cenpm, Gins1, Pola1, Prim1, Nop56, Kif18b, Mms22l, Lyar, Tcf19, Cdca3, Trip13, Ncapd2) are not previously described in breast cancer. Among those five E2F targets correlated with malignant transition, Cdk1 and Plk1 were well studied in breast cancer, but much less so for Ccnb2, Aurkb, and Spag5 (7, 10, and 3 reports, respectively). Among the 79 E2F targets, Hells, Hmgb3, and Cit are of special interests because of their functions. HELLS interacts with DNMT3A and DNMT3B in the STRING protein-protein interaction database [48, 49], and is known to be involved in survival, de novo DNA methylation, and DNA methylation maintenance . It also cooperates with HDAC1, HDAC2, and DNMTs to silence transcription . This is consistent with our finding of increased Dnmt1 expression and global gene down-regulation. HMGB3 promotes cell proliferation in bladder cancer  and interacts with TOP2A and TOP2B, two targets for some anticancer agents. CIT, a serine/threonine-protein kinase, functions together with Kif14 (its expression increased by more than four folds in the TCGA breast cancer samples) in cell division. In summary, we propose some less studied genes as potential new players in breast cancer progression.
To see how the differences in expression are reflected on the gene product level, we queried The Human Protein Atlas  for RNA expression and protein staining scores in the metastatic breast cancer cell line MCF7. Of the ten hub genes tested by the Kaplan-Meier analysis, seven showed significant protein expression (Additional file 2: Table S6).
We also compared our results with other studies of PyMT breast cancer. In a previous study, Andrechek et al  compiled various PyMT datasets from the GEO database and identified 4,184 DEGs between PyMT mouse tumor and normal mammary gland. Comparing this set of genes with our 3,476 DEGs from the late carcinoma stage, we found an overlap of 2,069 genes with expression changes in the same direction. To see how our DEGs behave in other mouse models, we compared the same 3,476 DEGs from our study with 4,018 DEGs identified in the same study  between carcinoma in MMTV-Neu mice (a mouse model of breast cancer with pulmonary metastasis) and normal mammary gland. We found an overlap of 1,642 genes with expression change in the same direction.
To verify our RNA-sequencing results, we performed qPCR on the same samples used for sequencing. We confirmed increased expression of Top2a and Dnmt1 in all four stages (p-value < 0.05, Additional file 6: Figure S5A). To validate that our findings are reproducible, we also carried our qPCR on independent samples from PyMT and FVB mice not used for sequencing. We selected six genes showing increased expression in the PyMT mouse from all stages, and observed similar expression in the qPCR results (Additional file 6: Figure S5B).
A major challenge in breast cancer study is to understand the molecular mechanisms of malignant progression. Mouse models have been widely used to characterize the molecular events in breast cancer. Compared to early microarray-based gene expression studies of PyMT mouse carcinoma [58–60], we demonstrated for the first time transcriptional changes associated with distinct stages in breast cancer progression using RNA-sequencing and observed temporal global mRNA expression deregulation in the PyMT mouse model. Some previous studies used PyMT mouse model to identify signature or biomarkers of tumor virulence , or residual tumors and lung metastases . Another study  has indicated that highly expressed genes in PyMT mouse are enriched for gene that predict metastasis using the van’t Veer gene set . Many more studies compared various breast cancer mouse models as well as human samples to get mouse model classes and general molecular features of the disease [9, 57, 59]. But none focused on the changes during tumor progression.
We observed that most DEGs identified in the late carcinoma stage first appeared in the much earlier hyperplasia stage. This is consistent with previous human cancer studies, which reported substantial similarity of transcriptional profiles in various pathological stages  and suggested that expression profiles of early primary tumors may tell progressive potential [29, 30]. Similar to our findings, earlier investigations with microarrays also observed that more genes were down-regulated in both PyMT mouse  and human luminal B subtype . Many of the significantly enriched pathways are actually part of the big “pathway in cancer” (KEGG:05200), such as ECM-receptor interaction and focal adhesion. Although not significant enriched in our data, the PyMT oncogenesis causal pathway PI3K-AKT signaling is also included in the “pathway of cancer” so is linked to many of those enriched pathways. Metabolism deregulation are widely recognized in cancer, such that high calorie diet elevates cancer risk, especially on glucose metabolism . In both DEGs and WGCNA module genes, we found enrichment of many metabolic related terms and pathways, including the TCA cycle pathway.
With the co-expression network analysis, we identified genes that function together. Because co-expression networks derived from RNA-sequencing data exhibit higher correlation and hence higher connectivity than ones derived from microarrays , gene modules found in this study could better explain molecular mechanisms of the disease progression. One module with increased gene expression in tumor is enriched with DNA replication and cell cycle genes. Module hub genes are not only potential prognosis biomarkers of the human disease demonstrated by our Kaplan-Meier analysis with data from human patients but also potential risk factors underlying the molecular mechanisms of human breast cancer. Some hubs in this module are known cancer-related genes, however some have not been identified as oncogenes or cancer biomarkers so far, such as Mcm6 and Ncapd2. Although many hub genes can predict future metastasis with statistical significance, the high and low risk groups show only small difference on risk for metastasis. It is possible that individual genes have only limited predictive power and more genes need to be combined into gene panels for finer separation of different risk groups as shown in many gene signature studies [29, 30]. The module enriched in oxidative phosphorylation and the TCA cycle showed decreased expression in tumor, which is consistent with the facts that cancers switch from oxidative phosphorylation to glycolysis due to hypoxia and rely less on the TCA cycle . The reduced expression of TCA cycle genes in carcinoma stage has been reported before in the TAG mouse model . It is very interesting that we observed metabolic changes in the early stages; they can be potential biomarkers for diagnosis. Nonetheless, we may not entirely rule out the possibility that the changes were due to hypoxia induced by the aggressive development for this particular mouse model. Focal adhesion, ECM-receptor interaction, and insulin pathway in other modules were also discovered in previous cancer DEGs studies.
E2F factors, including E2F1 to E2F3A, are well-known transcription activators in cell cycle [65–67], apoptosis, and proliferation . Although we did not detect significant variation in their expression, possibly due to low expression levels, their target genes were enriched in DEGs as well as in co-expression modules. Recent studies have unveiled other functions of E2Fs beyond cell cycle . E2F1 expression is suggested to correlate with esophageal squamous carcinoma progression . Its knockdown reduced invasion potential but not proliferation . E2Fs have been associated with relapse-free survival time in Myc-induced tumors . They mediate tumor development and metastasis by regulating the expression of genes involved in angiogenesis , such as Vegfa, Cyr61, and Angpt2, and thus remodeling cell survival and ECM . Besides, among 200 E2F target genes, there are 34 cell cycle genes and 22 DNA replication genes. Because the overlap between E2F targets and proliferation genes is only marginal, our enrichment result on E2F targets reflects more than a mere effect on proliferation. Although the involvement of E2Fs and their targets are well reported, here we proposed some new targets whose functions in breast cancer have not been documented. More importantly, we studied the expression of those E2F target in the early stages of tumor progression to investigate their potentials for cancer early diagnosis. The important roles that E2Fs play in metastasis warrant further investigation of their novel target genes that were discovered in this study. Moreover, the differential expression of Hells and Dnmt1 suggests a crucial role of DNA methylation in tumor development.
The enrichment analysis of DEGs confirmed the importance of cell cycle process in the PyMT breast cancer tumorigenesis. Indeed, cell cycle and proliferation genes are strong predictors of metastasis . Nevertheless, we also found DEGs enriched in other pathways such as PPAR signaling throughout tumor progression. PPAR signaling commonly regulates fatty acids metabolism and energy homeostasis and is part of the extensive “pathways in cancer” as defined by KEGG. There is emerging evidence of new cellular functions such as cell differentiation and tumorigenesis for PPAR signaling [74–76]. PPARγ, encoded by PPARG, is involved in multiple types of cancer including colonic tumor and breast tumor [74–76]. In breast cancer, PPARγ promotes terminal differentiation of malignant breast epithelial cells, and its activation triggered by drugs like antidiabetic thiazolidinedione (TZD) is associated with reduced cell growth and less malignancy . Another down-regulated pathway in PyMT was drug metabolism-cytochrome P450. Cytochrome P450 (CYP) genes are responsible for phase I drug metabolism and are key enzymes contributing to tumorigenesis through metabolic activation of precaricinogenes . As drug metabolizing enzymes, they also have crucial impact on anticancer drug treatment. CYP2E1, a feature gene in this pathway, was reported to inhibit cell migration in breast cancer cell lines on ectopic expression . Taken together, some key genes from these non-cell proliferation pathways may also regulate tumorigenesis, and thus other genes in these pathways are worth further investigation.
Even though PyMT mouse mimic the human disease, expression levels of some oncogenes vary between the two species . The PyMT oncoprotein is not expressed in human and the lung-specific metastasis nature in PyMT mice may lead to a not generalized molecular signature of breast cancer metastasis. Despite all the limitations, PyMT mouse model is widely appreciated in breast cancer research.
We found remarkable expression similarity among hyperplasia, adenoma/MIN, early and late carcinoma samples, meaning genes altered in the late stages showed trace in the beginning of tumor progression. It is suggested that we may use early stage expression profile to help diagnosis and treatment. In addition, some E2F target genes identified by differential expression analysis, co-expression network analysis and time-course analysis may promote tumor progression and are new candidate for breast cancer that worth further investigation.
Differentially expressed gene
Epithelia mesenchymal transition
log2 ratio change
Polyoma middle T
We thank Drs. John Greally, Jan Vijg and Yousin Suh for their expert advice. We also thank the Einstein Center for Epigenomics, the Epigenomics Shared Facility, the Histotechnology and Comparative Pathology Facility (supported by the Albert Einstein Cancer Center Support Grant from the National Institutes of Health under award number P30CA013330), and the Computational Epigenomics Group for technical support.
This work was supported by an NIH Pathway to Independence Award from National Library of Medicine (R00LM009770-06) to Z.D.Z.
Availability of data
The PyMT RNA-sequencing data is available in the Gene Expression Omnibus (GEO) database as GSE76772.
ZDZ conceived and designed this study. YC carried out the analyses and also participated in the study design. YC and ZDZ wrote the manuscript. KB performed animal experiments and sample preparation. RNC, JRL, QZ, WZ, and CM provided useful input for the analyses and helped edit the manuscript. All authors read and approved the final manuscript.
The authors declare that they have no competing interests.
Ethics approval and consent to participate
This study of cancer in mice was approved by the Institutional Animal Care and Use Committee (IACUC) of Albert Einstein College of Medicine. All procedures involving mice were conducted in accordance with the National Institutes of Health guidelines concerning the use and care of experimental animals.
Concent for Pulication
Consent for publication is not applicable in this study, no individual person’s data was used.
Open AccessThis article is distributed under the terms of the Creative Commons Attribution 4.0 International License (http://creativecommons.org/licenses/by/4.0/), which permits unrestricted use, distribution, and reproduction in any medium, provided you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made. The Creative Commons Public Domain Dedication waiver (http://creativecommons.org/publicdomain/zero/1.0/) applies to the data made available in this article, unless otherwise stated.
- Weigelt B, Peterse JL, van’t Veer LJ. Breast cancer metastasis: Markers and models. Nat Rev Cancer. 2005;5(8):591–602.View ArticlePubMedGoogle Scholar
- Jemal A, Bray F, Center MM, Ferlay J, Ward E, Forman D. Global cancer statistics. CA Cancer J Clin. 2011;61(2):69–90.View ArticlePubMedGoogle Scholar
- Ferlay J, Soerjomataram I, Dikshit R, Eser S, Mathers C, Rebelo M, Parkin DM, Forman D, Bray F. Cancer incidence and mortality worldwide: sources, methods and major patterns in GLOBOCAN 2012. Int J Cancer. 2015;136(5):E359–86.View ArticlePubMedGoogle Scholar
- Ma XJ, Salunga R, Tuggle JT, Gaudet J, Enright E, McQuary P, Payette T, Pistone M, Stecker K, Zhang BM, et al. Gene expression profiles of human breast cancer progression. Proc Natl Acad Sci U S A. 2003;100(10):5974–9.View ArticlePubMedPubMed CentralGoogle Scholar
- Rappaport N, Nativ N, Stelzer G, Twik M, Guan-Golan Y, Stein TI, Bahir I, Belinky F, Morrey CP, Safran M, et al. MalaCards: an integrated compendium for diseases and their annotation. Database. 2013;2013:bat018.View ArticlePubMedPubMed CentralGoogle Scholar
- Lakhani SR. The transition from hyperplasia to invasive carcinoma of the breast. J Pathol. 1999;187(3):272–8.View ArticlePubMedGoogle Scholar
- Fantozzi A, Christofori G. Mouse models of breast cancer metastasis. Breast Cancer Res: BCR. 2006;8(4):212.
- Hollern DP, Andrechek ER. A genomic analysis of mouse models of breast cancer reveals molecular features of mouse models and relationships to human breast cancer. Breast cancer research: BCR. 2014;16(3):R59.
- Herschkowitz JI, Simin K, Weigman VJ, Mikaelian I, Usary J, Hu ZY, Rasmussen KE, Jones LP, Assefnia S, Chandrasekharan S, et al. Identification of conserved gene expression features between murine mammary carcinoma models and human breast tumors. Genome Biol. 2007;8(5):R76.
- Lin EY, Jones JG, Li P, Zhu UY, Whitney KD, Muller WJ, Pollard JW. Progression to malignancy in the polyoma middle T oncoprotein mouse breast cancer model provides a reliable model for human diseases. Am J Pathol. 2003;163(5):2113–26.View ArticlePubMedPubMed CentralGoogle Scholar
- Rodriguez-Viciana P, Collins C, Fried M. Polyoma and SV40 proteins differentially regulate PP2A to activate distinct cellular signaling pathways involved in growth control. Proc Natl Acad Sci U S A. 2006;103(51):19290–5.View ArticlePubMedPubMed CentralGoogle Scholar
- Fluck MM, Schaffhausen BS. Lessons in signaling and tumorigenesis from polyomavirus middle T antigen. Microbiol Mol Biol Rev. 2009;73(3):542–63. Table of Contents.View ArticlePubMedPubMed CentralGoogle Scholar
- Mitra K, Carvunis AR, Ramesh SK, Ideker T. Integrative approaches for finding modular structure in biological networks. Nat Rev Genet. 2013;14(10):719–32.View ArticlePubMedPubMed CentralGoogle Scholar
- Koboldt DC, Fulton RS, McLellan MD, Schmidt H, Kalicki-Veizer J, McMichael JF, Fulton LL, Dooling DJ, Ding L, Mardis ER, et al. Comprehensive molecular portraits of human breast tumours. Nature. 2012;490(7418):61–70.View ArticleGoogle Scholar
- Prevarskaya N, Skryma R, Shuba Y. Calcium in tumour metastasis: new roles for known actors. Nat Rev Cancer. 2011;11(8):609–18.View ArticlePubMedGoogle Scholar
- McLellan AS, Dubin RA, Jing Q, Broin PO, Moskowitz D, Suzuki M, Calder RB, Hargitai J, Golden A, Greally JM. The Wasp System: an open source environment for managing and analyzing genomic data. Genomics. 2012;100(6):345–51.View ArticlePubMedGoogle Scholar
- Andrews S. FastQC: a quality control tool for high throughput sequence data. 2010.Google Scholar
- Wu TD, Nacu S. Fast and SNP-tolerant detection of complex variants and splicing in short reads. Bioinformatics. 2010;26(7):873–81.View ArticlePubMedPubMed CentralGoogle Scholar
- Anders S, Pyl PT, Huber W. HTSeq — A Python framework to work with high-throughput sequencing data. 2014.Google Scholar
- Team RC. R: A Language and Environment for Statistical Computing. 2014.Google Scholar
- Anders S, Huber W. Differential expression analysis for sequence count data. Genome Biol. 2010;11(10):R106.View ArticlePubMedPubMed CentralGoogle Scholar
- Robinson MD, McCarthy DJ, Smyth GK. edgeR: a Bioconductor package for differential expression analysis of digital gene expression data. Bioinformatics. 2010;26(1):139–40.View ArticlePubMedGoogle Scholar
- Croft D, Mundo AF, Haw R, Milacic M, Weiser J, Wu G, Caudy M, Garapati P, Gillespie M, Kamdar MR, et al. The Reactome pathway knowledgebase. Nucleic Acids Res. 2014;42(Database issue):D472–7.View ArticlePubMedGoogle Scholar
- Kanehisa M, Goto S, Sato Y, Kawashima M, Furumichi M, Tanabe M. Data, information, knowledge and principle: back to metabolism in KEGG. Nucleic Acids Res. 2014;42(Database issue):D199–205.View ArticlePubMedGoogle Scholar
- Young MD, Wakefield MJ, Smyth GK, Oshlack A. Gene ontology analysis for RNA-seq: accounting for selection bias. Genome Biol. 2010;11(2):R14.View ArticlePubMedPubMed CentralGoogle Scholar
- Subramanian A, Tamayo P, Mootha VK, Mukherjee S, Ebert BL, Gillette MA, Paulovich A, Pomeroy SL, Golub TR, Lander ES, et al. Gene set enrichment analysis: a knowledge-based approach for interpreting genome-wide expression profiles. Proc Natl Acad Sci U S A. 2005;102(43):15545–50.View ArticlePubMedPubMed CentralGoogle Scholar
- Langfelder P, Horvath S. WGCNA: an R package for weighted correlation network analysis. Bmc Bioinformatics. 2008;9.
- Hu Z, Snitkin ES, DeLisi C. VisANT: an integrative framework for networks in systems biology. Brief Bioinform. 2008;9(4):317–25.View ArticlePubMedPubMed CentralGoogle Scholar
- Wang Y, Klijn JG, Zhang Y, Sieuwerts AM, Look MP, Yang F, Talantov D, Timmermans M, Meijer-van Gelder ME, Yu J, et al. Gene-expression profiles to predict distant metastasis of lymph-node-negative primary breast cancer. Lancet. 2005;365(9460):671–9.View ArticlePubMedGoogle Scholar
- van’t Veer LJ, Dai HY, van de Vijver MJ, He YDD, Hart AAM, Mao M, Peterse HL, van der Kooy K, Marton MJ, Witteveen AT, et al. Gene expression profiling predicts clinical outcome of breast cancer. Nature. 2002;415(6871):530–6.View ArticleGoogle Scholar
- Cam H, Dynlacht BD. Emerging roles for E2F: beyond the G1/S transition and DNA replication. Cancer Cell. 2003;3(4):311–6.View ArticlePubMedGoogle Scholar
- Hollern DP, Honeysett J, Cardiff RD, Andrechek ER. The E2F transcription factors regulate tumor development and metastasis in a mouse model of metastatic breast cancer. Mol Cell Biol. 2014;34(17):3229–43.View ArticlePubMedPubMed CentralGoogle Scholar
- Putluri N, Maity S, Kommagani R, Creighton CJ, Putluri V, Chen F, Nanda S, Bhowmik SK, Terunuma A, Dorsey T, et al. Pathway-Centric Integrative Analysis Identifies RRM2 as a Prognostic Marker in Breast Cancer Associated with Poor Survival and Tamoxifen Resistance. Neoplasia. 2014;16(5):390–402.View ArticlePubMedPubMed CentralGoogle Scholar
- Ringner M, Fredlund E, Hakkinen J, Borg A, Staaf J. GOBO: gene expression-based outcome for breast cancer online. PLoS One. 2011;6(3), e17911.View ArticlePubMedPubMed CentralGoogle Scholar
- Esteller M. Cancer epigenomics: DNA methylomes and histone-modification maps. Nat Rev Genet. 2007;8(4):286–98.View ArticlePubMedGoogle Scholar
- Robert MF, Morin S, Beaulieu N, Gauthier F, Chute IC, Barsalou A, MacLeod AR. DNMT1 is required to maintain CpG methylation and aberrant gene silencing in human cancer cells. Nat Genet. 2003;33(1):61–5.View ArticlePubMedGoogle Scholar
- Dowsett M, Nielsen TO, A’Hern R, Bartlett J, Coombes RC, Cuzick J, Ellis M, Henry NL, Hugh JC, Lively T, et al. Assessment of Ki67 in breast cancer: recommendations from the International Ki67 in Breast Cancer working group. J Natl Cancer Inst. 2011;103(22):1656–64.View ArticlePubMedPubMed CentralGoogle Scholar
- Wang J, Xu B, Yuan P, Zhang P, Li Q, Ma F, Fan Y. TOP2A amplification in breast cancer is a predictive marker of anthracycline-based neoadjuvant chemotherapy efficacy. Breast Cancer Res Treat. 2012;135(2):531–7.View ArticlePubMedGoogle Scholar
- Ha GH, Kim JL, Petersson A, Oh S, Denning MF, Patel T, Breuer EK. TACC3 deregulates the DNA damage response and confers sensitivity to radiation and PARP inhibition. Oncogene. 2015;34(13):1667–78.View ArticlePubMedGoogle Scholar
- Winnepenninckx V, Lazar V, Michiels S, Dessen P, Stas M, Alonso SR, Avril MF, Ortiz Romero PL, Robert T, Balacescu O, et al. Gene expression profiling of primary cutaneous melanoma and clinical outcome. J Natl Cancer Inst. 2006;98(7):472–82.View ArticlePubMedGoogle Scholar
- Eppig JT, Blake JA, Bult CJ, Kadin JA, Richardson JE, Mouse Genome Database G. The Mouse Genome Database (MGD): facilitating mouse as a model for human biology and disease. Nucleic Acids Res. 2015;43(Database issue):D726–36.View ArticlePubMedGoogle Scholar
- Forsburg SL. Eukaryotic MCM Proteins: Beyond Replication Initiation. Microbiol Mol Biol Rev. 2004;68(1):109–31.View ArticlePubMedPubMed CentralGoogle Scholar
- Kwok HF, Zhang SD, McCrudden CM, Yuen HF, Ting KP, Wen Q, Khoo US, Chan KY. Prognostic significance of minichromosome maintenance proteins in breast cancer. Am J Cancer Res. 2015;5(1):52–71.PubMedGoogle Scholar
- Watrin E, Legagneux V. Contribution of hCAP-D2, a non-SMC subunit of condensin I, to chromosome and chromosomal protein dynamics during mitosis. Mol Cell Biol. 2005;25(2):740–50.View ArticlePubMedPubMed CentralGoogle Scholar
- Wang C, Cigliano A, Jiang L, Li X, Fan B, Pilo MG, Liu Y, Gui B, Sini M, Smith JW, et al. 4EBP1/eIF4E and p70S6K/RPS6 axes play critical and distinct roles in hepatocarcinogenesis driven by AKT and N-Ras proto-oncogenes in mice. Hepatology. 2015;61(1):200–13.View ArticlePubMedGoogle Scholar
- Mill CP, Zordan MD, Rothenberg SM, Settleman J, Leary JF, Riese 2nd DJ. ErbB2 Is Necessary for ErbB4 Ligands to Stimulate Oncogenic Activities in Models of Human Breast Cancer. Genes Cancer. 2011;2(8):792–804.View ArticlePubMedPubMed CentralGoogle Scholar
- Tang CK, Concepcion XZ, Milan M, Gong X, Montgomery E, Lippman ME. Ribozyme-mediated down-regulation of ErbB-4 in estrogen receptor-positive breast cancer cells inhibits proliferation both in vitro and in vivo. Cancer Res. 1999;59(20):5315–22.PubMedGoogle Scholar
- Szklarczyk D, Franceschini A, Wyder S, Forslund K, Heller D, Huerta-Cepas J, Simonovic M, Roth A, Santos A, Tsafou KP, et al. STRING v10: protein-protein interaction networks, integrated over the tree of life. Nucleic Acids Res. 2015;43(Database issue):D447–52.View ArticlePubMedGoogle Scholar
- Zhu H, Geiman TM, Xi S, Jiang Q, Schmidtmann A, Chen T, Li E, Muegge K. Lsh is involved in de novo methylation of DNA. EMBO J. 2006;25(2):335–45.View ArticlePubMedPubMed CentralGoogle Scholar
- Myant K, Termanis A, Sundaram AY, Boe T, Li C, Merusi C, Burrage J, de Las Heras JI, Stancheva I. LSH and G9a/GLP complex are required for developmentally programmed DNA methylation. Genome Res. 2011;21(1):83–94.View ArticlePubMedPubMed CentralGoogle Scholar
- Myant K, Stancheva I. LSH cooperates with DNA methyltransferases to repress transcription. Mol Cell Biol. 2008;28(1):215–26.View ArticlePubMedGoogle Scholar
- Li M, Cai Y, Zhao H, Xu Z, Sun Q, Luo M, Gu L, Meng M, Han X, Sun H. Overexpression of HMGB3 protein promotes cell proliferation, migration and is associated with poor prognosis in urinary bladder cancer patients. Tumour Biol. 2015;36(6):4785–92.View ArticlePubMedGoogle Scholar
- Parker JS, Mullins M, Cheang MC, Leung S, Voduc D, Vickery T, Davies S, Fauron C, He X, Hu Z, et al. Supervised risk predictor of breast cancer based on intrinsic subtypes. J Clin Oncol. 2009;27(8):1160–7.View ArticlePubMedPubMed CentralGoogle Scholar
- Begley DA, Krupke DM, Neuhauser SB, Richardson JE, Bult CJ, Eppig JT, Sundberg JP. The Mouse Tumor Biology Database (MTB): A Central Electronic Resource for Locating and Integrating Mouse Tumor Pathology Data. Vet Pathol. 2012;49(1):218–23.View ArticlePubMedGoogle Scholar
- Smith CM, Finger JH, Hayamizu TF, McCright IJ, Xu J, Berghout J, Campbell J, Corbani LE, Forthofer KL, Frost PJ, et al. The mouse Gene Expression Database (GXD): 2014 update. Nucleic Acids Res. 2014;42(Database issue):D818–24.View ArticlePubMedGoogle Scholar
- Uhlen M, Fagerberg L, Hallstrom BM, Lindskog C, Oksvold P, Mardinoglu A, Sivertsson A, Kampf C, Sjostedt E, Asplund A, et al. Tissue-based map of the human proteome. Science. 2015;347(6220):1260419.View ArticlePubMedGoogle Scholar
- Hollern DP, Andrechek ER. A genomic analysis of mouse models of breast cancer reveals molecular features of mouse models and relationships to human breast cancer. Breast Cancer Res. 2014;16(3):R59.View ArticlePubMedPubMed CentralGoogle Scholar
- Franci C, Zhou J, Jiang Z, Modrusan Z, Good Z, Jackson E, Kouros-Mehr H. Biomarkers of Residual Disease, Disseminated Tumor Cells, and Metastases in the MMTV-PyMT Breast Cancer Model. PloS One. 2013;8(3):e58183.
- Pfefferle AD, Herschkowitz JI, Usary J, Harrell JC, Spike BT, Adams JR, Torres-Arzayus MI, Brown M, Egan SE, Wahl GM, et al. Transcriptomic classification of genetically engineered mouse models of breast cancer identifies human subtype counterparts. Genome Biol. 2013;14(11):R125.
- Zhu M, Yi M, Kim CH, Deng C, Li Y, Medina D, Stephens RM, Green JE. Integrated miRNA and mRNA expression profiling of mouse mammary tumor models identifies miRNA signatures associated with mammary tumor lineage. Genome Biol. 2011;12(8):R77.View ArticlePubMedPubMed CentralGoogle Scholar
- Qiu TH, Chandramouli GV, Hunter KW, Alkharouf NW, Green JE, Liu ET. Global expression profiling identifies signatures of tumor virulence in MMTV-PyMT-transgenic mice: correlation to human disease. Cancer Res. 2004;64(17):5973–81.View ArticlePubMedGoogle Scholar
- Dang CV. Links between metabolism and cancer. Genes Dev. 2012;26(9):877–90.View ArticlePubMedPubMed CentralGoogle Scholar
- Iancu OD, Kawane S, Bottomly D, Searles R, Hitzemann R, McWeeney S. Utilizing RNA-Seq data for de novo coexpression network inference. Bioinformatics. 2012;28(12):1592–7.View ArticlePubMedPubMed CentralGoogle Scholar
- Levine AJ, Puzio-Kuter AM. The Control of the Metabolic Switch in Cancers by Oncogenes and Tumor Suppressor Genes. 2010.Google Scholar
- Attwooll C, Lazzerini Denchi E, Helin K. The E2F family: specific functions and overlapping interests. EMBO J. 2004;23(24):4709–4716.
- Nevins JR. The Rb/E2F pathway and cancer. Hum Mol Genet. 2001;10(7):699–703.View ArticlePubMedGoogle Scholar
- Trimarchi JM, Lees JA. Sibling rivalry in the E2F family. Nat Rev Mol Cell Biol. 2002;3(1):11–20.View ArticlePubMedGoogle Scholar
- Hallstrom TC, Mori S, Nevins JR. An E2F1-dependent gene expression program that determines the balance between proliferation and cell death. Cancer Cell. 2008;13(1):11–22.View ArticlePubMedPubMed CentralGoogle Scholar
- Chen HZ, Tsai SY, Leone G. Emerging roles of E2Fs in cancer: an exit from cell cycle control. Nat Rev Cancer. 2009;9(11):785–97.View ArticlePubMedPubMed CentralGoogle Scholar
- Ebihara Y, Miyamoto M, Shichinohe T, Kawarada Y, Cho Y, Fukunaga A, Murakami S, Uehara H, Kaneko H, Hashimoto H, et al. Over-expression of E2F-1 in esophageal squamous cell carcinoma correlates with tumor progression. Dis Esophagus. 2004;17(2):150–4.View ArticlePubMedGoogle Scholar
- Alla V, Engelmann D, Niemetz A, Pahnke J, Schmidt A, Kunz M, Emmrich S, Steder M, Koczan D, Putzer BM. E2F1 in melanoma progression and metastasis. J Natl Cancer Inst. 2010;102(2):127–33.View ArticlePubMedGoogle Scholar
- Fujiwara K, Yuwanita I, Hollern DP, Andrechek ER. Prediction and genetic demonstration of a role for activator E2Fs in Myc-induced tumors. Cancer Res. 2011;71(5):1924–32.View ArticlePubMedGoogle Scholar
- Merdzhanova G, Gout S, Keramidas M, Edmond V, Coll JL, Brambilla C, Brambilla E, Gazzeri S, Eymin B. The transcription factor E2F1 and the SR protein SC35 control the ratio of pro-angiogenic versus antiangiogenic isoforms of vascular endothelial growth factor-A to inhibit neovascularization in vivo. Oncogene. 2010;29(39):5392–403.View ArticlePubMedGoogle Scholar
- Tachibana K, Yamasaki D, Ishimoto K, Doi T. The Role of PPARs in Cancer. PPAR Res. 2008;2008:102737.View ArticlePubMedPubMed CentralGoogle Scholar
- Sarraf P, Mueller E, Jones D, King FJ, DeAngelo DJ, Partridge JB, Holden SA, Chen LB, Singer S, Fletcher C, et al. Differentiation and reversal of malignant changes in colon cancer through PPARgamma. Nat Med. 1998;4(9):1046–52.View ArticlePubMedGoogle Scholar
- Mueller E, Sarraf P, Tontonoz P, Evans RM, Martin KJ, Zhang M, Fletcher C, Singer S, Spiegelman BM. Terminal differentiation of human breast cancer through PPAR gamma. Mol Cell. 1998;1(3):465–70.View ArticlePubMedGoogle Scholar
- Rodriguez-Antona C, Ingelman-Sundberg M. Cytochrome P450 pharmacogenetics and cancer. Oncogene. 2006;25(11):1679–91.View ArticlePubMedGoogle Scholar
- Leung T, Rajendran R, Singh S, Garva R, Krstic-Demonacos M, Demonacos C. Cytochrome P450 2E1 (CYP2E1) regulates the response to oxidative stress and migration of breast cancer cells. Breast Cancer Res. 2013;15(6):R107.View ArticlePubMedPubMed CentralGoogle Scholar
- Wagner KU. Models of breast cancer: quo vadis, animal modeling? Breast Cancer Res. 2004;6(1):31–8.View ArticlePubMedGoogle Scholar
- Wickham H. ggplot2: elegant graphics for data analysis. New York: Springer; 2009.View ArticleGoogle Scholar