RNA-seq Data Analysis (Limma)

Mechanism of Doxorubicin-induced Cardiotoxicity

Organism Mus musculus

To compare expression profiles in the cardiomyocytes with wild type top2b and those with top2b deletion after in vivo treatment of mice with doxorubicin or drug vehicle ¹

Doxorubicin is widely used in modern cancer treatments, despite the advent of targeted therapy. However, a dose-dependent cardiotoxicity often limits its clinical use. The prevailing theory hypothesizes that doxorubicin-induced cardiotoxicity is the result of reactive oxygen species (ROS) generation due to redox-cycling of doxorubicin. Here is shown that cardiomyocyte-specific deletion of Topoisomerase II beta (Top2b) markedly reduced DNA double-strand breaks, apoptosis, and functional damages in doxorubicin-treated hearts. To investigate transcriptomic changes after doxorubicin treatment in wild type mouse and mouse with cardiac specific deletion of Top2b, researchers have examined the expression profiles in 4 groups of mice (3/group), ie. wildtype mice with or without doxorubicin treatment and mice with Top2b deletion in the cardiomyocytes with or without doxorubicin treatment. Mice were treated with doxorubicin (25mg/kg, i.p.) or PBS (drug vehicle) for 16 hr or 72 hr. The heart was removed and cardiomyocytes were isolated by using a Langendorff apparatus. After purification, total RNA was extracted from the cardiomyocytes, purified, and used for gene expression analysis. Compared with that in control cardiomyocytes or cardiomyocytes with Top2b deletion, doxorubicin caused a significant expression change in the genome of cardiomyocytes from the wildtype mice. Among the changes, multiple genes encoding mitochondrial structural protein and components of the respiratory chain complexes were down-regulated 72 hr after treatment while multiple genes in the p53 pathway were up-regulated 16 hr after treatment in the wildtype cardiomyocytes.

Overall design: Expression changes were examined in 2 groups of mice (wild type and conditional knockout of top2b in the cardiomyocytes) treated with doxorubicin or PBS for 16 or 72 hours

---

library(limma)
library(tidyverse)
library(dplyr)
library(magrittr)

Features  Samples 
   29532       12 

        treatment
genotype dox pbs
   top2b   3   3
   wt      3   3

---

Feature Inspection from Density-plot and Pre-processing

eset <- doxiru
exprs(eset) <- log(exprs(eset))
plotDensities(eset,  group = pData(eset)[,"genotype"], legend = "topright")

# Quantile normalize
exprs(eset) <- normalizeBetweenArrays(exprs(eset))
plotDensities(eset,  group = pData(eset)[,"genotype"], legend = "topright")
abline(v=0)

# Determining the genes with mean expression level greater than 0
keep <- rowMeans(exprs(eset)) > 0
sum(keep)

[1] 18361

eset <- eset[keep,]
plotDensities(eset, group = pData(eset)[,"genotype"], legend = "topright")

NA
NA

Boxplot of top2b

As expected, the expresssion of Top2b is much lower in the null mice (top2b) compared to wild type (wt).

Batch Effects & Checking Sources of Variations

Visualized the variation effect in PCA with Multidimensional Scaling. Both treatment and genotypic variates are checked.

plotMDS(eset, labels = pData(eset)[,"genotype"], gene.selection = "common")

plotMDS(eset, labels = pData(eset)[,"treatment"], gene.selection = "common")

Reassuringly, the samples cluster by their genotype and treatment. Interestingly, the Top2b null samples cluster more tightly compared to the wild type samples.

This supports the hypothesis that, top2b null mice are resistant to cardiotoxic effect of Doxorubicin.

---

Data Modelling

Factorial Design Using Group-mean Parameterization Model for Doxorubicine Study;

\(Y \:= \:\beta_1 X_1 \:+\: \beta_2 X_2 \:+\:\beta_3 X_3\:+\: \beta_4 X_4\:+\: \epsilon\)

• \(β_1\) - Mean expression level in top2b mice treated with dox
• \(β_2\) - Mean expression level in top2b treated with pbs
• \(β_3\) - Mean expression level in wt mice treated with dox
• \(β_4\) - Mean expression level in wt mice treated with pbs

Contrasts for Doxorubicine Study

• Response of wild type mice to dox treatment: \(β_3 − β_4\) = 0
• Response of Top2b null mice to dox treatment: \(β_1 − β_2\) = 0
• Differences between Top2b null and wild type mice in response to dox treatment: \((β_1 − β_2) - (β_3 − β_4) = 0\)

group <- with(pData(eset), paste(genotype, treatment, sep = "."))
group <- factor(group)


design <- model.matrix(~0 + group)
colnames(design) <- levels(group)


colSums(design)

top2b.dox top2b.pbs    wt.dox    wt.pbs 
        3         3         3         3 

This resulted in 4 coefficients that each model 3 samples.

To test for the effect of doxorubicin on the hearts of wild type and Top2b null mice (and any interaction between treatment and genotype), it’s needed to contrast the coefficients from the design matrix.

cm <- makeContrasts(dox_wt = wt.dox - wt.pbs,
                    dox_top2b = top2b.dox - top2b.pbs,
                    interaction = (top2b.dox - top2b.pbs) - (wt.dox - wt.pbs),
                    levels = design)

# View the contrasts matrix
cm

           Contrasts
Levels      dox_wt dox_top2b interaction
  top2b.dox      0         1           1
  top2b.pbs      0        -1          -1
  wt.dox         1         0          -1
  wt.pbs        -1         0           1

Model Fitting

fit <- lmFit(eset, design)
fit2 <- contrasts.fit(fit, contrasts = cm)
fit2 <- eBayes(fit2)

results <- decideTests(fit2)
summary(results)

       dox_wt dox_top2b interaction
Down     3180         0        1254
NotSig  12265     18361       15621
Up       2916         0        1486

# Create a Venn diagram
vennDiagram(results)

As expected, the doxorubucin only had an effect on the wild type mice.

Contrast Specified P-value Histogram

stats_dox_wt <- topTable(fit2, coef = "dox_wt", number = nrow(fit2),
                         sort.by = "none")

stats_dox_top2b <- topTable(fit2, coef = "dox_top2b", number = nrow(fit2),
                            sort.by = "none")

stats_interaction <- topTable(fit2, coef = "interaction", number = nrow(fit2),
                              sort.by = "none")


hist(stats_dox_wt[,"P.Value"])

hist(stats_dox_top2b[,"P.Value"])

hist(stats_interaction[,"P.Value"])

The contrasts dox_wt and interaction were enriched for low p-values, and the p-values for the contrast dox_top2b were uniformly distributed.

Contrast Specified Volcano-Plot

# Extract the gene symbols
gene_symbols <- fit2$genes[,"symbol"]

# Create a volcano plot for the contrast dox_wt
volcanoplot(fit2, coef = "dox_wt", highlight = 5, names = gene_symbols)

# Create a volcano plot for the contrast dox_top2b
volcanoplot(fit2, coef = "dox_top2b", highlight = 5, names = gene_symbols)

# Create a volcano plot for the contrast interaction
volcanoplot(fit2, coef = "interaction", highlight = 5, names =gene_symbols)

The difference in the x- and y-axis ranges for the dox_top2b contrast comapred to the other two are noted.

---

Pathway Enrichment Analysis

To better understand the effect of the differentially expressed genes in the doxorubicin study, enrichment of known biological pathways curated in the KEGG database are tested.

entrez <- fit2$genes[,"entrez"]

enrich_dox_wt <- kegga(fit2, coef = "dox_wt", geneid = entrez, species = "Mm")
topKEGG(enrich_dox_wt)

enrich_interaction <- kegga(fit2, coef = "interaction", geneid = entrez, species = "Mm")
topKEGG(enrich_interaction)

One of the top hits for both contrasts was a pathway for cardiomyopathy, so the genes in this pathway would be worth investigating further.

R version 3.6.1 (2019-07-05)
Platform: x86_64-w64-mingw32/x64 (64-bit)
Running under: Windows 10 x64 (build 18363)

Matrix products: default

locale:
[1] LC_COLLATE=English_United States.1252  LC_CTYPE=English_United States.1252   
[3] LC_MONETARY=English_United States.1252 LC_NUMERIC=C                          
[5] LC_TIME=English_United States.1252    

attached base packages:
[1] parallel  stats4    stats     graphics  grDevices utils     datasets  methods   base     

other attached packages:
 [1] limma_3.40.6                plotly_4.9.2.1              annotables_0.1.91          
 [4] forcats_0.5.0               stringr_1.4.0               dplyr_0.8.3                
 [7] purrr_0.3.3                 readr_1.3.1                 tidyr_1.0.0                
[10] tibble_2.1.3                ggplot2_3.3.2               tidyverse_1.3.0            
[13] magrittr_1.5                SummarizedExperiment_1.14.1 DelayedArray_0.10.0        
[16] BiocParallel_1.18.1         matrixStats_0.57.0          Biobase_2.44.0             
[19] GenomicRanges_1.36.1        GenomeInfoDb_1.20.0         IRanges_2.18.3             
[22] S4Vectors_0.22.1            BiocGenerics_0.30.0        

loaded via a namespace (and not attached):
 [1] colorspace_1.4-1       ellipsis_0.3.1         rsconnect_0.8.16       htmlTable_2.1.0       
 [5] XVector_0.24.0         base64enc_0.1-3        fs_1.3.1               rstudioapi_0.11       
 [9] remotes_2.2.0          bit64_4.0.5            AnnotationDbi_1.46.1   fansi_0.4.1           
[13] lubridate_1.7.4        xml2_1.2.2             splines_3.6.1          geneplotter_1.62.0    
[17] knitr_1.30             Formula_1.2-4          jsonlite_1.7.1         broom_0.5.3           
[21] annotate_1.62.0        GO.db_3.8.2            cluster_2.1.0          dbplyr_1.4.4          
[25] png_0.1-7              BiocManager_1.30.10    compiler_3.6.1         httr_1.4.2            
[29] backports_1.1.5        assertthat_0.2.1       Matrix_1.2-18          lazyeval_0.2.2        
[33] cli_2.1.0              org.Mm.eg.db_3.8.2     htmltools_0.5.0        tools_3.6.1           
[37] gtable_0.3.0           glue_1.3.1             GenomeInfoDbData_1.2.1 Rcpp_1.0.3            
[41] cellranger_1.1.0       vctrs_0.3.4            nlme_3.1-143           xfun_0.18             
[45] rvest_0.3.6            lifecycle_0.2.0        XML_3.99-0.3           zlibbioc_1.30.0       
[49] scales_1.1.1           hms_0.5.3              RColorBrewer_1.1-2     yaml_2.2.1            
[53] curl_4.3               memoise_1.1.0          gridExtra_2.3          rpart_4.1-15          
[57] latticeExtra_0.6-29    stringi_1.4.4          RSQLite_2.2.1          genefilter_1.66.0     
[61] checkmate_2.0.0        rlang_0.4.8            pkgconfig_2.0.3        bitops_1.0-6          
[65] evaluate_0.14          lattice_0.20-38        htmlwidgets_1.5.2      bit_4.0.4             
[69] tidyselect_0.2.5       R6_2.4.1               generics_0.0.2         Hmisc_4.4-1           
[73] DBI_1.1.0              pillar_1.4.6           haven_2.3.1            foreign_0.8-74        
[77] withr_2.3.0            survival_3.1-8         RCurl_1.98-1.2         nnet_7.3-12           
[81] modelr_0.1.8           crayon_1.3.4           rmarkdown_2.5          jpeg_0.1-8.1          
[85] locfit_1.5-9.4         grid_3.6.1             readxl_1.3.1           data.table_1.13.0     
[89] blob_1.2.1             reprex_0.3.0           digest_0.6.26          xtable_1.8-4          
[93] munsell_0.5.0          viridisLite_0.3.0     

---

1. “Zhang S, Liu X, Bawa-Khalfe T, Lu LS et al. Identification of the molecular basis of doxorubicin-induced cardiotoxicity. Nat Med 2012 Nov;18(11):1639-42. PMID: 23104132”↩

 

---

A work by Md. Tabassum Hossain Emon

emon.biotech.10th@gmail.com