Subset samples and define 11q deleted region.

Only genes on chr11 will be considered in this analysis

protCLL$del11q <- patMeta[match(colnames(protCLL),patMeta$Patient.ID),]$del11q
patTab <- colData(protCLL) %>% data.frame() %>% rownames_to_column("patID") %>%
  filter(!is.na(del11q))
dds <- dds[,dds$PatID %in% patTab$patID]
allRnaTab <- filter(allRnaTab, patID %in% patTab$patID, ChromID %in% "chr11" )
allProtTab <- filter(allProtTab, patID %in% patTab$patID, ChromID %in% "chr11")

I will define the region of 11q22.3 to 11q23.2 as the commonly deleted region, as most del11q samples show deletion in this region. Let me know if this is appropriate. This region will be refered to as 11q in the analysis

start <- 102.9
end <- 114.5
allRnaTab <- mutate(allRnaTab, ChromID = ifelse(start_position > start & end_position < end, "11q", "other"))
allProtTab <- mutate(allProtTab, ChromID = ifelse(start_position > start & end_position < end, "11q", "other"))

Compare gene dosage effect on RNA and protein level

#remove genes never expressed
noExpTab <- group_by(allRnaTab, id) %>% summarise(sumExpr = sum(expr)) %>%
  filter(sumExpr == 0)

rnaExprTab <- allRnaTab %>% filter(!id%in% noExpTab$id) %>%
  mutate(del11q = patMeta[match(patID, patMeta$Patient.ID),]$del11q) %>%
  filter(!is.na(del11q)) %>% mutate(cnv = ifelse(del11q %in% 1, "del11q","wt"))

meanExprChr11 <- rnaExprTab %>% filter(ChromID %in% "11q") %>%
  group_by(id, symbol, cnv) %>% summarise(meanExpr = mean(expr, na.rm=TRUE)) %>%
  ungroup()

ggplot(meanExprChr11, aes(x=meanExpr, fill = cnv)) + geom_histogram(position = "identity", alpha=0.5) 

meanExprChr11 <- rnaExprTab %>% filter(ChromID %in% "11q") %>%
  group_by(id) %>% mutate(med=mean(expr),sd = sd(expr)) %>% mutate(expr = (expr-med)) %>%
  group_by(id, symbol, cnv) %>% summarise(meanExpr = mean(expr, na.rm=TRUE)) %>%
  ungroup()
ggplot(meanExprChr11, aes(x=meanExpr, fill = cnv)) + geom_histogram(position = "identity", alpha=0.5) 

plotTab <- rnaExprTab %>%
  group_by(id, cnv, ChromID) %>% summarise(meanExpr = mean(expr, na.rm=TRUE)) %>%
  ungroup() %>%
  spread(key = cnv, value = meanExpr) %>%
  mutate(diff = del11q-wt)

ggplot(plotTab, aes(x=diff, fill = ChromID, y = ..density..)) + geom_histogram(position = "identity", alpha=0.5) +
  xlab("Mean expression difference")

protExprTab <- allProtTab %>% 
  mutate(del11q = patMeta[match(patID, patMeta$Patient.ID),]$del11q) %>%
  filter(!is.na(del11q)) %>% mutate(cnv = ifelse(del11q %in% 1, "del11q","wt"))

meanExprChr11 <- protExprTab %>% filter(ChromID %in% "11q") %>%
  group_by(id, symbol, cnv) %>% summarise(meanExpr = mean(expr, na.rm=TRUE)) %>%
  ungroup()

ggplot(meanExprChr11, aes(x=meanExpr, fill = cnv)) + geom_histogram(position = "identity", alpha=0.5) 

meanExprChr11 <- protExprTab %>% filter(ChromID %in% "11q") %>%
  group_by(id) %>% mutate(med=mean(expr),sd = sd(expr)) %>% mutate(expr = (expr-med)) %>%
  group_by(id, symbol, cnv) %>% summarise(meanExpr = mean(expr, na.rm=TRUE)) %>%
  ungroup()
ggplot(meanExprChr11, aes(x=meanExpr, fill = cnv)) + geom_histogram(position = "identity", alpha=0.5) 

plotTab <- protExprTab %>%
  group_by(id, cnv, ChromID) %>% summarise(meanExpr = mean(expr, na.rm=TRUE)) %>%
  ungroup() %>%
  spread(key = cnv, value = meanExpr) %>%
  mutate(diff = del11q-wt)

ggplot(plotTab, aes(x=diff, fill = ChromID, y = ..density..)) + geom_histogram(position = "identity", alpha=0.5) 

varRna <- filter(rnaExprTab, ChromID == "11q") %>%
  group_by(id,del11q) %>% summarise(sd = sd(expr)) %>%
  mutate(set = "rna")
varProt <- filter(protExprTab, ChromID == "11q") %>%
  group_by(id,del11q) %>% summarise(sd = sd(expr)) %>%
  mutate(set = "protein")
plotTab <- bind_rows(varRna, varProt) %>% ungroup()

ggplot(plotTab, aes(x=sd, fill = set, y = ..density..)) + 
  geom_histogram(position = "identity", alpha=0.5, bins = 100)+
  facet_wrap(~del11q) + xlab("Expression variance")

expRna <- filter(rnaExprTab, ChromID == "11q") %>%
  group_by(id) %>% mutate(meanVal = mean(expr)) %>%
  mutate(expr = expr-meanVal) %>%
  group_by(id,del11q) %>%
  summarise(meanExpr = mean(expr)) %>%
  mutate(set = "rna")
expProt <- filter(protExprTab, ChromID == "11q") %>%
  group_by(id) %>% mutate(meanVal = mean(expr)) %>%
  mutate(expr = expr-meanVal) %>%
  group_by(id,del11q) %>%
  summarise(meanExpr = mean(expr)) %>%
  mutate(set = "protein")
plotTab <- bind_rows(expRna, expProt)

ggplot(plotTab, aes(x=meanExpr, fill = set, y = ..density..)) + 
  geom_histogram(position = "identity", alpha=0.5, bins = 100) +
  facet_wrap(~del11q) +
  xlab("Mean expression")

protDiffTab <- protExprTab %>%
  group_by(id, cnv, ChromID) %>% summarise(meanExpr = mean(expr, na.rm=TRUE)) %>%
  ungroup() %>%
  spread(key = cnv, value = meanExpr) %>%
  mutate(diffProt = log(del11q/wt)) %>%
  mutate(chr = ifelse(ChromID == "11q","11q","Other")) %>%
  select(id, diffProt, chr)

rnaDiffTab <- rnaExprTab %>%
  group_by(id, cnv, ChromID) %>% summarise(meanExpr = mean(expr, na.rm=TRUE)) %>%
  ungroup() %>%
  spread(key = cnv, value = meanExpr) %>%
  mutate(diffRna = log(del11q/wt)) %>% select(id, diffRna)

compareTab <- left_join(protDiffTab, rnaDiffTab, by = "id") %>%
  filter(!is.na(diffProt),!is.na(diffRna)) %>%
  gather(key = dataset, value=diff,-id,-chr) %>%
  mutate(group =paste0(chr,"_",dataset)) %>%
  filter(chr == "11q")
  

ggplot(compareTab, aes(x=diff, fill = group, y = ..density..)) + 
  geom_histogram(position = "identity", alpha=0.5, bins = 100) +
  geom_vline(xintercept = 0, color = "blue", linetype ="dashed") +
  xlab("Expression different (del11q - wt)")

Analysis of the buffering effect

gene11q <- intersect(filter(allProtTab, ChromID == "11q")$id, filter(allRnaTab, ChromID =="11q")$id)
#subset samples
overSample <- intersect(dds$PatID, colnames(protCLL))
protSub <- protCLL[rowData(protCLL)$ensembl_gene_id %in% gene11q ,overSample]
rownames(protSub) <- rowData(protSub)$ensembl_gene_id
designMat <- data.frame(row.names = colnames(protSub),
                        del11q = protSub$del11q)
exprMat <- assays(protSub)[["count"]]

#testing

fit <- proDA(exprMat, design = ~ ., col_data = designMat )
  
resTab.prot <- test_diff(fit, contrast = "del11q1") %>%
  select(name, pval, diff, t_statistic, adj_pval) %>%
  dplyr::rename(id = name, logFC.prot = diff, stat.prot = t_statistic, pval.prot = pval, padj.prot = adj_pval)

#subset samples
overSample <- intersect(dds$PatID, colnames(protCLL))
ddsSub <- dds[gene11q,overSample]
ddsSub$IGHV <- protSub[,ddsSub$PatID]$IGHV.status
ddsSub$del11q <- protSub[,ddsSub$PatID]$del11q

design(ddsSub) <- ~ del11q
#testing
deRes <- DESeq(ddsSub)
  
resTab.rna <- results(deRes, name = "del11q_1_vs_0", tidy = TRUE) %>%
  select(row, log2FoldChange, stat, pvalue, padj) %>%
  dplyr::rename(id = row, logFC.rna = log2FoldChange, stat.rna = stat, pval.rna = pvalue, padj.rna = padj) %>%
  mutate(logFC.rna = logFC.rna*log(2)) #change to log fold change, same as proteomic

comTab <- left_join(resTab.prot, resTab.rna, by = "id") %>%
  mutate(symbol = rowData(dds[id,])$symbol)

bufferTab <- comTab %>% filter(stat.rna < 0, stat.prot<0) %>%
  ungroup() %>%
  mutate(stat.prot.sqrt = sqrt(abs(stat.prot)),
         stat.prot.center = stat.prot.sqrt - mean(stat.prot.sqrt)) %>%
  mutate(score = stat.prot.center*stat.rna) %>%
  mutate(ifBuffer = case_when(
    padj.prot < 0.1 & padj.rna < 0.1 ~ "noBuffer",
    padj.prot > 0.1 & padj.rna < 0.1 ~ "Buffered",
    padj.prot < 0.1 & padj.rna > 0.1 ~ "Enhanced",
    TRUE ~ "undetermined"
  )) %>%
  arrange(desc(score))

ggplot(bufferTab, aes(x=ifBuffer,y=score, fill = ifBuffer)) + geom_boxplot() + geom_point()

table(bufferTab$ifBuffer)

    Buffered     noBuffer undetermined 
           3            8            1 

ggplot(bufferTab, aes(x=logFC.rna, y=logFC.prot)) + geom_point(aes(col = ifBuffer)) +
  xlab("RNA log fold change") + ylab("Protein log fold change") +
  geom_smooth(method = "lm") 

select(bufferTab, symbol, ifBuffer, score, padj.prot, padj.rna, logFC.prot,logFC.rna) %>%
  mutate_if(is.numeric, formatC, digits=2) %>%
  DT::datatable()

geneList <- filter(bufferTab, ifBuffer == "Buffered")$id
pList <- lapply(geneList, function(i) {
  tabProt <- allProtTab %>% filter(id == i) %>%
    select(id, patID, symbol,expr) %>% dplyr::rename(protExpr = expr)
  tabRna <- allRnaTab %>% filter(id == i) %>%
    select(id, patID, expr) %>% dplyr::rename(rnaExpr = expr)
  plotTab <- left_join(tabProt, tabRna, by = c("id","patID")) %>% 
    filter(!is.na(protExpr), !is.na(rnaExpr)) %>%
    mutate(del11q = patMeta[match(patID, patMeta$Patient.ID),]$del11q)
  p <- ggplot(plotTab, aes(x=rnaExpr, y = protExpr)) +
    geom_point(aes(col=del11q)) + geom_smooth(method="lm") + ggtitle(unique(plotTab$symbol)) +
    theme(legend.position = "bottom")
  ggMarginal(p, type = "histogram", groupFill = TRUE)
  })
cowplot::plot_grid(plotlist = pList, ncol=3)

geneList <- filter(bufferTab, ifBuffer == "noBuffer")$id
pList <- lapply(geneList, function(i) {
  tabProt <- allProtTab %>% filter(id == i) %>%
    select(id, patID, symbol,expr) %>% dplyr::rename(protExpr = expr)
  tabRna <- allRnaTab %>% filter(id == i) %>%
    select(id, patID, expr) %>% dplyr::rename(rnaExpr = expr)
  plotTab <- left_join(tabProt, tabRna, by = c("id","patID")) %>% 
    filter(!is.na(protExpr), !is.na(rnaExpr)) %>%
    mutate(del11q = patMeta[match(patID, patMeta$Patient.ID),]$del11q)
  p <- ggplot(plotTab, aes(x=rnaExpr, y = protExpr)) +
    geom_point(aes(col=del11q)) + geom_smooth(method="lm") +
    ggtitle(unique(plotTab$symbol)) +
    theme(legend.position = "bottom")
  ggMarginal(p, type = "histogram", groupFill = TRUE)
  })
cowplot::plot_grid(plotlist = pList, ncol=3)

Whether buffering is related to protein complexes

int_pairs = read_delim("../data/proteins_in_complexes", delim = "\t") %>%
  mutate(symbolA = rowData(protCLL)[match(ProtA, rownames(protCLL)),]$hgnc_symbol,
         symbolB = rowData(protCLL)[match(ProtB, rownames(protCLL)),]$hgnc_symbol) %>%
  filter(!is.na(symbolA),!is.na(symbolB))

bufferTab <- mutate(bufferTab, inComplex = ifelse(symbol %in% c(int_pairs$symbolA, int_pairs$symbolB), TRUE, FALSE))

ggplot(bufferTab, aes(x=inComplex, y=score)) + geom_boxplot() + geom_point()

t.test(score~inComplex, bufferTab, var.equal= TRUE)

    Two Sample t-test

data:  score by inComplex
t = -0.19566, df = 10, p-value = 0.8488
alternative hypothesis: true difference in means is not equal to 0
95 percent confidence interval:
 -6.761654  5.669999
sample estimates:
mean in group FALSE  mean in group TRUE 
          -1.899900           -1.354073 

Plot buffering effect on genomic coordinates

plotFoldGenome <- function(bufferTab, allBand, allProtTab, chr, region = c(-Inf,Inf),
                           ifTrend = FALSE, maxVal =1, minVal=-1) {
  
  #table for cyto band
  bandTab <- filter(allBand, ChromID == chr, chromStart >= region[1], chromEnd <= region[2]) %>%
    mutate(chromMid = chromMid)
  
  #table for fold change
  protCoordTab <- allProtTab %>% distinct(symbol, start_position, end_position, mid_position)
  foldTab <- bufferTab %>% select(symbol, logFC.prot, logFC.rna, score, ifBuffer) %>%
    gather(key = "set", value = "logFC", -symbol, -score,-ifBuffer) %>%
    left_join(protCoordTab)
  bufferLine <- filter(bufferTab, ifBuffer %in% c("Buffered","noBuffer")) %>%
    left_join(protCoordTab) %>%
    distinct(symbol, mid_position, logFC.prot, logFC.rna, ifBuffer) %>%
    mutate(minY = ifelse(logFC.prot > logFC.rna, logFC.rna, logFC.prot),
           maxY= ifelse(logFC.prot > logFC.rna, logFC.prot, logFC.rna))
  
  xMax <- max(bandTab$chromEnd, na.rm = T)
  xMin <- min(bandTab$chromStart, na.rm=T)
  #main plot for Protein
  gPro <- ggplot() + 
    geom_rect(data=bandTab, mapping=aes(xmin=chromStart, xmax=chromEnd, ymin=minVal, ymax=maxVal, 
                                        fill=Colour, label = band), alpha=0.1) +
    geom_text(data=bandTab, mapping=aes(label=band, x=chromMid), y=maxVal, hjust =1, angle = 90, size=2.5) +
    geom_rect(data = foldTab, 
            mapping=aes(xmin=start_position,
                        xmax=end_position, ymin=logFC, ymax=logFC+0.1,
                        fill = set)) +
    geom_segment(data = bufferLine, aes(x=mid_position, xend = mid_position, 
                                        y=minY + 0.1, yend = maxY, col = ifBuffer),
                 linetype = "dashed") +
    scale_x_continuous(expand=c(0,0),limits = c(xMin,xMax)) +
    xlab("Genomic position [Mb]") + 
    ylab("Log Fold Change") + 
    scale_fill_manual(values = c(even = "white",odd = "grey50",
                                logFC.rna = "orange", logFC.prot = "darkblue")) +
    scale_color_manual(values = c(logFC.rna = "orange",logFC.prot = "darkblue",
                                  Buffered = "red",noBuffer = "green")) +
    ggtitle(paste0("Log fold change comparison","(",chr,")")) +
    ggrepel::geom_text_repel(data = bufferLine,
                                 aes(x=mid_position, y=logFC.prot, label = symbol, col = ifBuffer)) +
    theme(plot.title = element_text(face = "bold", size = 10, hjust = 0.3),
        legend.position = "none",
        panel.background = element_blank(),
        panel.grid.major = element_line(colour="grey90", size=0.1))
  
    if (ifTrend) {
      gPro <- gPro + stat_smooth(data =foldTab, geom="line",
                mapping = aes(y=logFC, x= mid_position,
                              color = set), 
                formula = y ~ x, method = "loess", se=FALSE, span=0.2,
                size =0.5, alpha=0.5)
    }
    
   
    #for legend
    ## if the patient has CNV data
    lgTab <- tibble(x= seq(6),y=seq(6),
                    Dataset = c(rep("logFC.rna",3), rep("logFC.prot",3)),
                    ifBuffer = c(rep("Buffered",3), rep("noBuffer",3)))
  
    lg <- ggplot(lgTab, aes(x=x,y=y)) +
      geom_point(aes(fill = Dataset, color = ifBuffer), shape =22,size=3) + 
      scale_fill_manual(values = c(logFC.rna = "orange", logFC.prot = "darkblue")) +
      scale_color_manual(values = c(logFC.rna = "orange",logFC.prot = "darkblue",
                                  Buffered = "red",noBuffer = "green")) + 
      theme(legend.position = "bottom")
    
    lg <- get_legend(lg)
    
    return(list(plot = gPro, legend = lg))
}

band11q <- filter(allBand, ChromID == "chr11", band %in% c("q22.3","q23.1","q23.2")) %>%
  mutate(ChromID = "11q")
g <- plotFoldGenome(bufferTab, band11q, allProtTab, "11q", region = c(-Inf,Inf),
                           ifTrend = FALSE, maxVal =0.5, minVal=-2)
pg <- plot_grid(g$plot, g$legend, ncol = 1, rel_heights = c(1,0.1))
pg

ggsave(filename = "../public/del11q_buffer_plot.pdf", plot = pg, device = "pdf", height = 10, width = 18)

PDF version: del11q_buffer_plot.pdf

In this plot, the y axis in the log fold change of either protein (blue) or RNA (orange) expression (del11q vs WT). If there’s a “Buffering” effect, the protein and rna is connected by a red dotted line. The noBuffer effect will be joined by a green dotted line.

Summary:

1. Similar as trisomy12, the gene dosage effect of del11q is visible in both RNA expression and protein expression, as compared to WT samples. And it seems there’s less buffering based on the ratio of buffered proteins, compared to trisomy12 and trisomy19. This maybe because it’s more difficult to buffer a deletion than a gain.

2. Same as trisomy12 and trisomy19, there’s no significant association between the buffering effect and whether the protein is in complex or not.

Session information

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] cowplot_0.9.4               ggplot2_3.3.0              
[11] ggExtra_0.9                 proDA_1.1.2                
[13] jyluMisc_0.1.5              piano_2.0.2                
[15] DESeq2_1.24.0               SummarizedExperiment_1.14.0
[17] DelayedArray_0.10.0         BiocParallel_1.18.0        
[19] matrixStats_0.54.0          Biobase_2.44.0             
[21] GenomicRanges_1.36.0        GenomeInfoDb_1.20.0        
[23] IRanges_2.18.1              S4Vectors_0.22.0           
[25] 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] fastmatch_1.1-0        drc_3.0-1              workflowr_1.6.0       
  [7] igraph_1.2.4.1         shinydashboard_0.7.1   splines_3.6.0         
 [10] crosstalk_1.0.0        TH.data_1.0-10         digest_0.6.19         
 [13] htmltools_0.4.0        gdata_2.18.0           magrittr_1.5          
 [16] checkmate_2.0.0        memoise_1.1.0          cluster_2.1.0         
 [19] openxlsx_4.1.0.1       limma_3.40.2           annotate_1.62.0       
 [22] modelr_0.1.5           sandwich_2.5-1         colorspace_1.4-1      
 [25] ggrepel_0.8.1          rvest_0.3.5            blob_1.1.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] survival_2.44-1.1      zoo_1.8-6              glue_1.3.2            
 [37] survminer_0.4.4        gtable_0.3.0           zlibbioc_1.30.0       
 [40] XVector_0.24.0         car_3.0-3              abind_1.4-5           
 [43] scales_1.1.0           mvtnorm_1.0-11         DBI_1.0.0             
 [46] relations_0.6-8        miniUI_0.1.1.1         Rcpp_1.0.1            
 [49] plotrix_3.7-6          cmprsk_2.2-8           xtable_1.8-4          
 [52] htmlTable_1.13.1       foreign_0.8-71         bit_1.1-14            
 [55] km.ci_0.5-2            Formula_1.2-3          DT_0.7                
 [58] httr_1.4.1             htmlwidgets_1.3        fgsea_1.10.0          
 [61] gplots_3.0.1.1         RColorBrewer_1.1-2     acepack_1.4.1         
 [64] ellipsis_0.2.0         farver_2.0.3           pkgconfig_2.0.2       
 [67] XML_3.98-1.20          dbplyr_1.4.2           nnet_7.3-12           
 [70] locfit_1.5-9.1         labeling_0.3           tidyselect_1.0.0      
 [73] rlang_0.4.5            later_0.8.0            AnnotationDbi_1.46.0  
 [76] munsell_0.5.0          cellranger_1.1.0       tools_3.6.0           
 [79] visNetwork_2.0.7       cli_1.1.0              generics_0.0.2        
 [82] RSQLite_2.1.1          broom_0.5.2            evaluate_0.14         
 [85] yaml_2.2.0             knitr_1.23             bit64_0.9-7           
 [88] fs_1.4.0               zip_2.0.2              survMisc_0.5.5        
 [91] caTools_1.17.1.2       nlme_3.1-140           mime_0.7              
 [94] slam_0.1-45            xml2_1.2.2             compiler_3.6.0        
 [97] rstudioapi_0.10        curl_3.3               ggsignif_0.5.0        
[100] marray_1.62.0          reprex_0.3.0           geneplotter_1.62.0    
[103] stringi_1.4.3          lattice_0.20-38        Matrix_1.2-17         
[106] KMsurv_0.1-5           shinyjs_1.0            vctrs_0.2.4           
[109] pillar_1.4.3           lifecycle_0.2.0        data.table_1.12.2     
[112] bitops_1.0-6           httpuv_1.5.1           extraDistr_1.8.11     
[115] R6_2.4.0               latticeExtra_0.6-28    promises_1.0.1        
[118] KernSmooth_2.23-15     gridExtra_2.3          rio_0.5.16            
[121] codetools_0.2-16       MASS_7.3-51.4          gtools_3.8.1          
[124] exactRankTests_0.8-30  assertthat_0.2.1       rprojroot_1.3-2       
[127] withr_2.1.2            multcomp_1.4-10        GenomeInfoDbData_1.2.1
[130] mgcv_1.8-28            hms_0.5.2              grid_3.6.0            
[133] rpart_4.1-15           rmarkdown_1.13         carData_3.0-2         
[136] ggpubr_0.2.1           git2r_0.26.1           maxstat_0.7-25        
[139] sets_1.0-18            lubridate_1.7.4        shiny_1.3.2           
[142] base64enc_0.1-3