Run AlteredPRQ algorithm to detect protein complex ratio complex changes

source ("../code/AlteredPQR.R")
quant_data_all = assays(protCLL)[["QRILC"]]
cols_with_reference_data = seq(ncol(protCLL))[protCLL$trisomy12 %in% 0]
RepresentativePairs = Altered_PQR(modif_z_score_threshold = 3.0, fraction_of_samples_threshold = 0.3)

[1] "Running"
[1] "..."
[1] "..."
[1] "Top 0.1, 1 and 5% upper and lower z-score values are: 8.26742846535531 3.93569392733038 2.19885719576535 and -6.27632295669363 -3.42621202641534 -1.92712652454253."
[1] "Top 1% of the absolute values for the modified z-scores is 4.55241553106446."

Re-format output

protRes.pqr <- lapply(RepresentativePairs, function(x) x) %>% bind_cols() %>%
  separate(Protein_pair, into = c("idA","idB"),"-") %>%
  mutate(protA = rowData(protCLL[idA,])$hgnc_symbol,
         protB = rowData(protCLL[idB,])$hgnc_symbol,
         chrA = rowData(protCLL[idA,])$chromosome_name,
         chrB = rowData(protCLL[idB,])$chromosome_name) %>%
  mutate(withChr12 = ifelse(chrA %in% "12" | chrB%in% "12", "yes", "no"))%>% mutate(idx = seq(nrow(.))) %>%
  mutate(pair=map2_chr(idA, idB, ~paste0(sort(c(.x,.y)), collapse = "-")))

Run the same algorithm on RNA expression data

dds$trisomy12 <- patMeta[match(dds$PatID,patMeta$Patient.ID),]$trisomy12
rowData(dds)$uniprotID <- rownames(protCLL)[match(rownames(dds), rowData(protCLL)$ensembl_gene_id)]
ddsSub <- dds[!is.na(rowData(dds)$uniprotID), dds$PatID %in% colnames(protCLL)]
rownames(ddsSub) <- rowData(ddsSub)$uniprotID
ddsSub.vst <- varianceStabilizingTransformation(ddsSub)

source ("../code/AlteredPQR.R")
quant_data_all = assay(ddsSub.vst)
cols_with_reference_data = seq(ncol(ddsSub.vst))[ddsSub.vst$trisomy12 %in% 0]
RepresentativePairs = Altered_PQR(modif_z_score_threshold = 3.0, fraction_of_samples_threshold = 0.3)

[1] "Running"
[1] "..."
[1] "..."
[1] "Top 0.1, 1 and 5% upper and lower z-score values are: 9.02665546130853 3.95530891659894 2.44967137434484 and -6.48456606482675 -3.77874942568752 -2.3012824895991."
[1] "Top 1% of the absolute values for the modified z-scores is 4.61537110973003."

rnaRes.pqr <- lapply(RepresentativePairs, function(x) x) %>% bind_cols() %>%
  separate(Protein_pair, into = c("idA","idB"),"-") %>%
  dplyr::rename(rnaChange = Change, rnaScore= Score) %>%
  mutate(pair=map2_chr(idA, idB, ~paste0(sort(c(.x,.y)), collapse = "-"))) %>%
  select(pair, rnaScore, rnaChange)

Combine protein and RNA result

comRes.pqr <- left_join(protRes.pqr, rnaRes.pqr, by ="pair") %>%
  mutate(explainedByRNA = ifelse(is.na(rnaChange), "no", 
                                 ifelse(Change == rnaChange, "yes", "no")))

Exploring results

Whether protein pairs involve genes on Chr12 have higher score?

ggplot(comRes.pqr, aes(x=Score, fill = withChr12, col = withChr12)) + geom_histogram(alpha=0.5, position = "identity") + theme_bw()

Not really, but there are some chr12 pairs have extremely high score.

How many of those changes can be explained at RNA level

table(comRes.pqr$explainedByRNA)

 no yes 
183   2 

Only 2

List of detected pairs

comRes.pqr %>% select(protA, protB, Score, Change, chrA, chrB, withChr12, explainedByRNA) %>%
  mutate(Score = format(Score, digits = 1)) %>%
  DT::datatable() 

Visualization using network plot

Only pairs involving genes on chr12 are used for now Build network

comRes.filt <- filter(comRes.pqr, Score > 0, withChr12 == "yes")

#get node list
allNodes <- union(comRes.filt$protA, comRes.filt$protB) 

nodeList <- data.frame(id = seq(length(allNodes))-1, name = allNodes, stringsAsFactors = FALSE) %>%
  mutate(onChr12 = ifelse(rowData(protCLL[match(name, rowData(protCLL)$hgnc_symbol),])$chromosome_name %in% "12", 
                          "chr12","otherChr")) %>% 
  mutate(group = corRes.sig[match(name, corRes.sig$name),]$direction) %>%
  mutate(group = ifelse(is.na(group),"ns",group))

#get edge list
edgeList <- select(comRes.filt, protA, protB, Change, explainedByRNA) %>%
  dplyr::rename(Source = protA, Target = protB) %>% 
  mutate(Source = nodeList[match(Source,nodeList$name),]$id,
         Target = nodeList[match(Target, nodeList$name),]$id) %>%
  data.frame(stringsAsFactors = FALSE)

net <- graph_from_data_frame(vertices = nodeList, d=edgeList, directed = FALSE)

Visualize using ggraph

tidyNet <- as_tbl_graph(net)
ggraph(tidyNet) + geom_edge_link(aes(color = Change, edge_linetype = explainedByRNA), width=1) + 
  geom_node_point(aes(color =group, shape = onChr12), size=4) + 
  geom_node_text(aes(label = name), repel = TRUE) +
  scale_color_manual(values = c(Up = "pink",Down = "lightblue", ns="grey"))+
  theme_graph() 

In this plot, proteins on Chr12 are shown as solid circles and others are shown as triangles. The proteins up-regulated in trisomy12 samples are colored by red, down-regulated proteins are colored by cyan and the proteins with no significant changes are colored by grey. The color of the edges indication the whether the ratio of two proteins in pairs are increased or decreased in the reference group. But I don’t think this value has real biological meaning. It depends which protein is placed as the first one when calculating the ratio.

Inspecting some interesting pairs

plotPair <- function(comRes, protList, protCLL, gene) { 
  pairList <- filter(comRes, protA %in% protList | protB %in% protList)
  plotList <- lapply(seq(nrow(pairList)), function(i) {
    idA <- pairList[i,]$idA
    idB <- pairList[i,]$idB
    protA <- pairList[i,]$protA
    protB  <- pairList[i,]$protB
    idPair <- c(idA, idB)
    protPair <- c(protA, protB)
    ord <- order(protPair)
    idPair <- idPair[ord]
    protPair <- protPair[ord]
    plotTab <- assays(protCLL)[["count"]][idPair,] %>%
      t() %>% data.frame()
    colnames(plotTab) <- protPair
    plotTab$logRatio <- log2(plotTab[,1]) - log2(plotTab[,2])
    plotTab <- rownames_to_column(plotTab,"patID") %>%
      mutate(status = factor(protCLL[,patID][[gene]])) %>%
      filter(!is.na(logRatio))
    histP <- ggplot(plotTab, aes(x=logRatio, fill = status, col = status)) +
      geom_histogram(position = "identity", alpha=0.5) +
      ggtitle(sprintf("Stoichiometry: %s ~ %s",protPair[1], protPair[2]))
    corP <- ggplot(plotTab, aes_string(x=protPair[1], y=protPair[2], col="status")) + 
      geom_point() + geom_smooth(formula = y~x, method = "lm") +
      scale_color_discrete(name = gene)
    plot_grid(histP, corP)
  })
  return(plotList)
}

plotPair.rna <- function(comRes, protList, ddsSub.vst, gene) { 
  pairList <- filter(comRes, protA %in% protList | protB %in% protList)
  plotList <- lapply(seq(nrow(pairList)), function(i) {
    idA <- pairList[i,]$idA
    idB <- pairList[i,]$idB
    protA <- pairList[i,]$protA
    protB  <- pairList[i,]$protB
    idPair <- c(idA, idB)
    protPair <- c(protA, protB)
    ord <- order(protPair)
    idPair <- idPair[ord]
    protPair <- protPair[ord]
    plotTab <- assay(ddsSub.vst)[idPair,] %>%
      t() %>% data.frame()
    colnames(plotTab) <- protPair
    plotTab$logRatio <- log2(plotTab[,1]) - log2(plotTab[,2])
    plotTab <- rownames_to_column(plotTab,"patID") %>%
      mutate(status = factor(ddsSub.vst[,patID][[gene]])) %>%
      filter(!is.na(logRatio))
    histP <- ggplot(plotTab, aes(x=logRatio, fill = status, col = status)) +
      geom_histogram(position = "identity", alpha=0.5) +
      ggtitle(sprintf("Stoichiometry: %s ~ %s",protPair[1], protPair[2]))
    corP <- ggplot(plotTab, aes_string(x=protPair[1], y=protPair[2], col="status")) + 
      geom_point() + geom_smooth(formula = y~x, method = "lm") +
      scale_color_discrete(name = gene)
    plot_grid(histP, corP)
  })
  return(plotList)
}

protList <- c("PTPN6","PTPN11")
plotPair(comRes.pqr, protList, protCLL, "trisomy12")

[[1]]

[[2]]

[[3]]

[[4]]

[[5]]

protList <- c("STAT2")
plotPair(comRes.pqr, protList, protCLL, "trisomy12")

[[1]]

[[2]]

protList <- c("PTPN6","PTPN11")
plotPair.rna(comRes.pqr, protList, ddsSub.vst, "trisomy12")

[[1]]

[[2]]

[[3]]

[[4]]

[[5]]

protList <- c("STAT2")
plotPair.rna(comRes.pqr, protList, ddsSub.vst, "trisomy12")

[[1]]

[[2]]

Differential correlation detection using DGCA package

quant_data_all = assays(protCLL)[["QRILC"]]
quant_data_all <- quant_data_all[order(rownames(quant_data_all)),]

tri12 <- protCLL$trisomy12
designMat <- model.matrix(~tri12+0 )
colnames(designMat) <- c("WT","Tri12")
ddcor_res = ddcorAll(inputMat = quant_data_all, design = designMat,
  compare = c("WT", "Tri12"),
  adjust = "BH", heatmapPlot = FALSE, nPerm = 0, nPairs = "all")

Reformat output

comTab <- int_pairs %>%
  mutate(pair = map2_chr(ProtA, ProtB, ~paste0(sort(c(.x, .y)),collapse = "-"))) %>%
  separate(pair, c("Gene1","Gene2"), "-",remove = FALSE) %>%
  select(Gene1, Gene2) %>% mutate(inComplex= TRUE)

allRes <- left_join(ddcor_res, comTab, by = c("Gene1","Gene2")) %>%
  mutate(inComplex = ifelse(is.na(inComplex),FALSE,TRUE))

Distribution of p-values for protein in complexes and not in complexes

ggplot(allRes, aes(x=pValDiff, fill = inComplex)) + geom_histogram() + facet_wrap(~inComplex, scale="free") +
  xlim(0,1)

Not much difference.

Differential correlation detection on RNA level

quant_data_all = assay(ddsSub.vst)
quant_data_all <- quant_data_all[order(rownames(quant_data_all)),]

tri12 <- ddsSub.vst$trisomy12
designMat <- model.matrix(~tri12+0 )
colnames(designMat) <- c("WT","Tri12")
ddcor_res = ddcorAll(inputMat = quant_data_all, design = designMat,
  compare = c("WT", "Tri12"),
  adjust = "BH", heatmapPlot = FALSE, nPerm = 0, nPairs = "all")

rnaRes.cor <- ddcor_res %>%
  select(Gene1, Gene2,  pValDiff, pValDiff_adj, Classes) %>%
  dplyr::rename(p.rna = pValDiff, padj.rna = pValDiff_adj, Classes.rna = Classes)

Exploring the results

Only for proteins involved in known complexes

comRes.cor <- filter(allRes, inComplex) %>%
   mutate(protA = rowData(protCLL[Gene1,])$hgnc_symbol,
         protB = rowData(protCLL[Gene2,])$hgnc_symbol,
         chrA = rowData(protCLL[Gene1,])$chromosome_name,
         chrB = rowData(protCLL[Gene2,])$chromosome_name) %>%
  mutate(withChr12 = ifelse(chrA %in% "12" | chrB%in% "12", "yes", "no")) %>% mutate(idx = seq(nrow(.))) %>%
  mutate(p=pValDiff,padj = pValDiff_adj)

Add test results from RNA

comRes.cor <- left_join(comRes.cor, rnaRes.cor, by =c("Gene1","Gene2")) %>%
  mutate(explainedByRNA = ifelse(is.na(p.rna),"no",
                                 ifelse(padj.rna < 0.25 & Classes == Classes.rna,"yes","no")))

P values for pair involving or not involving complex

ggplot(comRes.cor, aes(x=p, fill = withChr12, col = withChr12)) + 
  geom_histogram(alpha=0.5, position = "identity", col = NA) + facet_wrap(~withChr12, scale = "free") +
  xlim(0,1)

List of significant pairs (25% FDR) As this test is very stringent, I use the looser FDR cut-off here.

comRes.sig <- filter(comRes.cor, withChr12 %in% c("yes")) %>%
  mutate(padj = p.adjust(p, method = "BH"),
         ifSig = padj <0.25) %>%
  filter(ifSig)
comRes.sig %>% select(protA, protB, p, padj, chrA, chrB, withChr12, Classes, explainedByRNA) %>%
  mutate_if(is.numeric, formatC, digits=2, format="e") %>%
  DT::datatable() 

comRes.filt <- filter(comRes.sig)
#comRes.filt <- comRes
#get node list
allNodes <- union(comRes.filt$protA, comRes.filt$protB) 

nodeList <- data.frame(id = seq(length(allNodes))-1, name = allNodes, stringsAsFactors = FALSE) %>%
  mutate(onChr12 = ifelse(rowData(protCLL[match(name, rowData(protCLL)$hgnc_symbol),])$chromosome_name %in% "12", "chr12","otherChr")) %>% 
  mutate(group = corRes.sig[match(name, corRes.sig$name),]$direction) %>%
  mutate(group = ifelse(is.na(group),"ns",group)) %>%
  mutate(nnSize = ifelse(onChr12 == "chr12",100, 1))

#get edge list
edgeList <- select(comRes.filt, protA, protB, p, Classes) %>%
  dplyr::rename(Source = protA, Target = protB) %>% 
  mutate(Source = nodeList[match(Source,nodeList$name),]$id,
         Target = nodeList[match(Target, nodeList$name),]$id,
         Classes = as.character(Classes)) %>%
  data.frame(stringsAsFactors = FALSE)

net <- graph_from_data_frame(vertices = nodeList, d=edgeList, directed = FALSE)

tidyNet <- as_tbl_graph(net)
ggraph(tidyNet) + geom_edge_link(aes(color = Classes), width=1) + 
  geom_node_point(aes(color =group, shape = onChr12), size=4) + 
  geom_node_text(aes(label = name), repel = TRUE) +
  scale_color_manual(values = c(Up = "pink",Down = "lightblue", ns="grey"))+
  theme_graph() 

sessionInfo()

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

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] forcats_0.4.0               stringr_1.4.0              
 [3] dplyr_0.8.5                 purrr_0.3.3                
 [5] readr_1.3.1                 tidyr_1.0.0                
 [7] tibble_3.0.0                tidyverse_1.3.0            
 [9] ggraph_1.0.2                igraph_1.2.4.1             
[11] cowplot_0.9.4               ggplot2_3.3.0              
[13] DESeq2_1.24.0               proDA_1.1.2                
[15] DGCA_1.0.2                  tidygraph_1.1.2            
[17] SummarizedExperiment_1.14.0 DelayedArray_0.10.0        
[19] BiocParallel_1.18.0         matrixStats_0.54.0         
[21] Biobase_2.44.0              GenomicRanges_1.36.0       
[23] GenomeInfoDb_1.20.0         IRanges_2.18.1             
[25] S4Vectors_0.22.0            BiocGenerics_0.30.0        

loaded via a namespace (and not attached):
  [1] readxl_1.3.1           backports_1.1.4        Hmisc_4.2-0           
  [4] workflowr_1.6.0        plyr_1.8.4             splines_3.6.0         
  [7] crosstalk_1.0.0        robust_0.4-18.1        digest_0.6.19         
 [10] foreach_1.4.4          htmltools_0.4.0        viridis_0.5.1         
 [13] GO.db_3.8.2            magrittr_1.5           checkmate_2.0.0       
 [16] memoise_1.1.0          fit.models_0.5-14      cluster_2.1.0         
 [19] doParallel_1.0.14      fastcluster_1.1.25     annotate_1.62.0       
 [22] modelr_0.1.5           colorspace_1.4-1       rvest_0.3.5           
 [25] blob_1.1.1             rrcov_1.4-9            ggrepel_0.8.1         
 [28] haven_2.2.0            xfun_0.8               crayon_1.3.4          
 [31] RCurl_1.95-4.12        jsonlite_1.6           genefilter_1.66.0     
 [34] impute_1.58.0          survival_2.44-1.1      iterators_1.0.10      
 [37] glue_1.3.2             polyclip_1.10-0        gtable_0.3.0          
 [40] zlibbioc_1.30.0        XVector_0.24.0         DEoptimR_1.0-8        
 [43] scales_1.1.0           mvtnorm_1.0-11         DBI_1.0.0             
 [46] Rcpp_1.0.1             viridisLite_0.3.0      xtable_1.8-4          
 [49] htmlTable_1.13.1       foreign_0.8-71         bit_1.1-14            
 [52] preprocessCore_1.46.0  Formula_1.2-3          DT_0.7                
 [55] htmlwidgets_1.3        httr_1.4.1             RColorBrewer_1.1-2    
 [58] acepack_1.4.1          ellipsis_0.2.0         pkgconfig_2.0.2       
 [61] XML_3.98-1.20          farver_2.0.3           nnet_7.3-12           
 [64] dbplyr_1.4.2           locfit_1.5-9.1         dynamicTreeCut_1.63-1 
 [67] labeling_0.3           tidyselect_1.0.0       rlang_0.4.5           
 [70] later_0.8.0            AnnotationDbi_1.46.0   cellranger_1.1.0      
 [73] munsell_0.5.0          tools_3.6.0            cli_1.1.0             
 [76] generics_0.0.2         RSQLite_2.1.1          broom_0.5.2           
 [79] evaluate_0.14          yaml_2.2.0             knitr_1.23            
 [82] bit64_0.9-7            fs_1.4.0               robustbase_0.93-5     
 [85] nlme_3.1-140           mime_0.7               xml2_1.2.2            
 [88] compiler_3.6.0         rstudioapi_0.10        reprex_0.3.0          
 [91] tweenr_1.0.1           geneplotter_1.62.0     pcaPP_1.9-73          
 [94] stringi_1.4.3          lattice_0.20-38        Matrix_1.2-17         
 [97] vctrs_0.2.4            pillar_1.4.3           lifecycle_0.2.0       
[100] data.table_1.12.2      bitops_1.0-6           httpuv_1.5.1          
[103] R6_2.4.0               latticeExtra_0.6-28    promises_1.0.1        
[106] gridExtra_2.3          codetools_0.2-16       MASS_7.3-51.4         
[109] assertthat_0.2.1       rprojroot_1.3-2        withr_2.1.2           
[112] GenomeInfoDbData_1.2.1 mgcv_1.8-28            hms_0.5.2             
[115] grid_3.6.0             rpart_4.1-15           rmarkdown_1.13        
[118] git2r_0.26.1           ggforce_0.2.2          shiny_1.3.2           
[121] lubridate_1.7.4        WGCNA_1.68             base64enc_0.1-3