Are those three miRNA in our RNAseq dataset?

mirList <- c("MIR15","MIR16","MIR34A")
sapply(mirList, function(n) rowData(dds)$symbol[grepl(n,rowData(dds)$symbol)])

Loading required package: DESeq2

$MIR15
 [1] "MIR153-1" "MIR15B"   "MIR150"   "MIR152"   "MIR153-2" "MIR154"  
 [7] "MIR1537"  "MIR1538"  "MIR155HG" "MIR151A"  "MIR1587"  "MIR151B" 

$MIR16
[1] "MIR16-2" "MIR16-1"

$MIR34A
[1] "MIR34AHG"

miRNAs selected for test

mirSele <- c("MIR15B","MIR16-1","MIR16-2","MIR34AHG")
mirSele

[1] "MIR15B"   "MIR16-1"  "MIR16-2"  "MIR34AHG"

Expression level of those 4 miRNAs

ddsMIR <- dds[rowData(dds)$symbol %in% mirSele,dds$PatID %in% colnames(protCLL)]
plotTab <- counts(ddsMIR,normalized = TRUE) %>% data.frame() %>% rownames_to_column("id") %>%
  gather(key = "patID", value = "expr", -id) %>%
  mutate(symbol = rowData(dds[id,])$symbol,
         IGHV = protCLL[,patID]$IGHV.status,
         trisomy12 = protCLL[,patID]$trisomy12)

ggplot(plotTab, aes(x=symbol, y = expr)) + geom_boxplot() + geom_point(aes(col = IGHV, shape = trisomy12)) +
  scale_y_log10() + ylab("Normalized expression")

Warning: Transformation introduced infinite values in continuous y-axis

Warning: Transformation introduced infinite values in continuous y-axis

Warning: Removed 158 rows containing non-finite values (stat_boxplot).

Past versions of unnamed-chunk-4-1.png

The expressions are pretty low

Do the expressions of miRNA correlate with IGHV or trisomy12?

testRes <- group_by(plotTab, symbol) %>% nest() %>%
  mutate(m = map(data, ~wilcox.test(expr~IGHV, .))) %>%
  mutate(res = map(m, broom::tidy)) %>%
  unnest(res) %>%
  select(symbol, p.value)
testRes

# A tibble: 4 x 2
# Groups:   symbol [4]
  symbol    p.value
  <chr>       <dbl>
1 MIR16-2    0.540 
2 MIR15B   NaN     
3 MIR16-1    0.568 
4 MIR34AHG   0.0636

testRes <- group_by(plotTab, symbol) %>% nest() %>%
  mutate(m = map(data, ~wilcox.test(expr~trisomy12, .))) %>%
  mutate(res = map(m, broom::tidy)) %>%
  unnest(res) %>%
  select(symbol, p.value)
testRes

# A tibble: 4 x 2
# Groups:   symbol [4]
  symbol   p.value
  <chr>      <dbl>
1 MIR16-2    0.137
2 MIR15B   NaN    
3 MIR16-1    0.127
4 MIR34AHG   0.244

Correlated miRNA expression with protein expression

Subset dataset and vst

#remove gene that are not expressed 
ddsMIR <- ddsMIR[rowSums(counts(ddsMIR, normalized = TRUE)) > 0,]
protSub <- protCLL[,colnames(ddsMIR)]
rnaMat <- counts(ddsMIR,normalized = TRUE)
protMat <- assay(protSub)

corTab <- lapply(rownames(rnaMat), function(i) {
  lapply(rownames(protMat), function(j) {
    rna <- rnaMat[i,]
    prot <- protMat[j,]
    res <- cor.test(rna, prot, method ="spearman")
    tibble(geneID = i, proteinID = j, P.Value = res$p.value, coef = res$estimate)
  }) %>% bind_rows()
}) %>% bind_rows() %>% mutate(geneSymbol = rowData(dds[geneID,])$symbol,
                                 proteinSymbol = rowData(protCLL[proteinID,])$hgnc_symbol) %>%
  mutate(adj.P.Val = p.adjust(P.Value, method = "BH")) %>%
  arrange(P.Value)

P-value histogram

hist(corTab$P.Value)

Overall, there are not many significant associations

Table of significant associations (raw P value < 0.01, no assocation passed 10% FDR)

corTab.sig <- filter(corTab, P.Value < 0.01) %>% arrange(P.Value)
corTab.sig  %>% select(geneSymbol, proteinSymbol, P.Value, adj.P.Val) %>%
  arrange(P.Value) %>%
  mutate_if(is.numeric, formatC, digits=2, format= "e") %>%
  DT::datatable()

Number of significant associations per miRNA (10% FDR)

sumTab <- corTab %>% filter(P.Value < 0.01) %>%
  mutate(direction = ifelse(coef >0, "pos","neg")) %>%
  group_by(geneSymbol, direction) %>% summarise(n=length(geneID)) %>%
  arrange(desc(n))

Bar plot

ggplot(sumTab, aes(x=geneSymbol, y = n, fill = direction)) + geom_bar(stat = "identity", position = position_dodge2()) +
  theme(axis.text.x = element_text(angle =90, hjust=1, vjust=0.5))

Correlation plot of top 9 most correlated miRNA protein pairs

plotList <- lapply(seq(1,9), function(i) {
  idProt <- corTab.sig[i,]$proteinID
  idGene <- corTab.sig[i,]$geneID
  plotTab <- tibble(pro = protMat[idProt,], gene = rnaMat[idGene,])
  ggplot(plotTab, aes(x=gene, y=pro)) + geom_point() + geom_smooth(method = "lm") +
    ylab(corTab.sig[i,]$proteinSymbol) + xlab(corTab.sig[i,]$geneSymbol)
})

cowplot::plot_grid(plotlist = plotList, ncol =3)

Get miRNA profile from RNAseq data

Not that the mRNAseq is optimized to capture protein coding genes, not miRNA.

ddsMIR <- dds[grepl("MIR",rowData(dds)$symbol),dds$PatID %in% colnames(protCLL)]
#remove gene that are not expressed 
ddsMIR <- ddsMIR[rowSums(counts(ddsMIR, normalized = TRUE)) > 10,]
ddsMIR.vst <- varianceStabilizingTransformation(ddsMIR)
sds <- genefilter::rowSds(assay(ddsMIR.vst))
ddsMIR.vst <- ddsMIR.vst[sds > genefilter::shorth(sds),]
dim(ddsMIR.vst)

[1] 228  46

protSub <- protCLL[,colnames(ddsMIR.vst)]

Here I use all the miRNAs that were detected by RNAseq and show variations of expression, if we have a more focused list, we may increase the power by reducing number of multilpe tests.

Calculate correlations between protein abundance and miRNA expression

In this test, I blocked for ighv and trisomy12, which means the test will identify correlations that are not confounded by IGHV and trisomy12

rnaMat <- assay(ddsMIR.vst)
proMat <- assays(protSub)[["count"]]
ighv <- protSub$IGHV.status
trisomy12 <- protSub$trisomy12

corTab <- lapply(rownames(rnaMat), function(i) {
    rna <- rnaMat[i,]

    designMat <- model.matrix(~1+ighv + trisomy12 + rna) #block for IGHV and trisomy12
    fit <- lmFit(proMat,  designMat)
    fit2 <- eBayes(fit)
    corRes <- topTable(fit2, number ="all", adjust.method = "BH", coef = "rna") %>% rownames_to_column("proteinID") %>%
      mutate(geneID = i)
    corRes
}) %>% bind_rows() %>% mutate(geneSymbol = rowData(dds[geneID,])$symbol,
                                 proteinSymbol = rowData(protCLL[proteinID,])$hgnc_symbol) %>%
  mutate(adj.P.All = p.adjust(P.Value, method = "BH"))

P-value histogram

hist(corTab$P.Value)

Overall, there are not many significant associations

Table of significant associations (10% FDR)

corTab.sig <- filter(corTab, adj.P.Val < 0.1) %>% arrange(P.Value)
corTab.sig  %>% select(geneSymbol, proteinSymbol, P.Value, adj.P.Val) %>%
  arrange(P.Value) %>%
  mutate_if(is.numeric, formatC, digits=2, format= "e") %>%
  DT::datatable()

Number of significant associations per miRNA (10% FDR)

sumTab <- corTab %>% filter(adj.P.Val < 0.1) %>%
  mutate(direction = ifelse(logFC >0, "pos","neg")) %>%
  group_by(geneSymbol, direction) %>% summarise(n=length(geneID)) %>%
  arrange(desc(n))

Bar plot

ggplot(sumTab, aes(x=geneSymbol, y = n, fill = direction)) + geom_bar(stat = "identity", position = position_dodge2()) +
  theme(axis.text.x = element_text(angle =90, hjust=1, vjust=0.5))

Correlation plot of top 9 most correlated miRNA protein pairs

plotList <- lapply(seq(1,9), function(i) {
  idProt <- corTab.sig[i,]$proteinID
  idGene <- corTab.sig[i,]$geneID
  plotTab <- tibble(pro = proMat[idProt,], gene = rnaMat[idGene,])
  ggplot(plotTab, aes(x=gene, y=pro)) + geom_point() + geom_smooth(method = "lm") +
    ylab(corTab.sig[i,]$proteinSymbol) + xlab(corTab.sig[i,]$geneSymbol)
})

cowplot::plot_grid(plotlist = plotList, ncol =3)

Session information

sessionInfo()

R version 3.6.0 (2019-04-26)
Platform: x86_64-apple-darwin15.6.0 (64-bit)
Running under: macOS Mojave 10.14.6

Matrix products: default
BLAS:   /Library/Frameworks/R.framework/Versions/3.6/Resources/lib/libRblas.0.dylib
LAPACK: /Library/Frameworks/R.framework/Versions/3.6/Resources/lib/libRlapack.dylib

locale:
[1] en_US.UTF-8/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8

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

other attached packages:
 [1] DESeq2_1.24.0               forcats_0.4.0              
 [3] stringr_1.4.0               dplyr_0.8.3                
 [5] purrr_0.3.3                 readr_1.3.1                
 [7] tidyr_1.0.0                 tibble_2.1.3               
 [9] tidyverse_1.3.0             SummarizedExperiment_1.14.0
[11] DelayedArray_0.10.0         BiocParallel_1.18.0        
[13] matrixStats_0.54.0          Biobase_2.44.0             
[15] GenomicRanges_1.36.0        GenomeInfoDb_1.20.0        
[17] IRanges_2.18.1              S4Vectors_0.22.0           
[19] BiocGenerics_0.30.0         jyluMisc_0.1.5             
[21] pheatmap_1.0.12             limma_3.40.2               
[23] piano_2.0.2                 cowplot_0.9.4              
[25] ggplot2_3.2.1              

loaded via a namespace (and not attached):
  [1] utf8_1.1.4             shinydashboard_0.7.1   tidyselect_0.2.5      
  [4] RSQLite_2.1.1          AnnotationDbi_1.46.0   htmlwidgets_1.3       
  [7] grid_3.6.0             maxstat_0.7-25         munsell_0.5.0         
 [10] codetools_0.2-16       DT_0.7                 withr_2.1.2           
 [13] colorspace_1.4-1       knitr_1.23             rstudioapi_0.10       
 [16] ggsignif_0.5.0         labeling_0.3           git2r_0.26.1          
 [19] slam_0.1-45            GenomeInfoDbData_1.2.1 KMsurv_0.1-5          
 [22] bit64_0.9-7            rprojroot_1.3-2        vctrs_0.2.0           
 [25] generics_0.0.2         TH.data_1.0-10         xfun_0.8              
 [28] sets_1.0-18            R6_2.4.0               locfit_1.5-9.1        
 [31] bitops_1.0-6           fgsea_1.10.0           assertthat_0.2.1      
 [34] promises_1.0.1         scales_1.0.0           multcomp_1.4-10       
 [37] nnet_7.3-12            gtable_0.3.0           sandwich_2.5-1        
 [40] workflowr_1.6.0        rlang_0.4.1            zeallot_0.1.0         
 [43] genefilter_1.66.0      cmprsk_2.2-8           splines_3.6.0         
 [46] lazyeval_0.2.2         acepack_1.4.1          broom_0.5.2           
 [49] checkmate_1.9.3        yaml_2.2.0             abind_1.4-5           
 [52] modelr_0.1.5           crosstalk_1.0.0        backports_1.1.4       
 [55] httpuv_1.5.1           Hmisc_4.2-0            tools_3.6.0           
 [58] relations_0.6-8        ellipsis_0.2.0         gplots_3.0.1.1        
 [61] RColorBrewer_1.1-2     Rcpp_1.0.1             base64enc_0.1-3       
 [64] visNetwork_2.0.7       zlibbioc_1.30.0        RCurl_1.95-4.12       
 [67] ggpubr_0.2.1           rpart_4.1-15           zoo_1.8-6             
 [70] haven_2.2.0            cluster_2.1.0          exactRankTests_0.8-30 
 [73] fs_1.3.1               magrittr_1.5           data.table_1.12.2     
 [76] openxlsx_4.1.0.1       reprex_0.3.0           survminer_0.4.4       
 [79] mvtnorm_1.0-11         whisker_0.3-2          hms_0.5.2             
 [82] shinyjs_1.0            mime_0.7               evaluate_0.14         
 [85] xtable_1.8-4           XML_3.98-1.20          rio_0.5.16            
 [88] readxl_1.3.1           gridExtra_2.3          compiler_3.6.0        
 [91] KernSmooth_2.23-15     crayon_1.3.4           htmltools_0.3.6       
 [94] later_0.8.0            Formula_1.2-3          geneplotter_1.62.0    
 [97] lubridate_1.7.4        DBI_1.0.0              dbplyr_1.4.2          
[100] MASS_7.3-51.4          Matrix_1.2-17          car_3.0-3             
[103] cli_1.1.0              marray_1.62.0          gdata_2.18.0          
[106] igraph_1.2.4.1         pkgconfig_2.0.2        km.ci_0.5-2           
[109] foreign_0.8-71         xml2_1.2.2             annotate_1.62.0       
[112] XVector_0.24.0         drc_3.0-1              rvest_0.3.5           
[115] digest_0.6.19          rmarkdown_1.13         cellranger_1.1.0      
[118] fastmatch_1.1-0        survMisc_0.5.5         htmlTable_1.13.1      
[121] curl_3.3               shiny_1.3.2            gtools_3.8.1          
[124] lifecycle_0.1.0        nlme_3.1-140           jsonlite_1.6          
[127] carData_3.0-2          fansi_0.4.0            pillar_1.4.2          
[130] lattice_0.20-38        httr_1.4.1             plotrix_3.7-6         
[133] survival_2.44-1.1      glue_1.3.1             zip_2.0.2             
[136] bit_1.1-14             stringi_1.4.3          blob_1.1.1            
[139] latticeExtra_0.6-28    caTools_1.17.1.2       memoise_1.1.0