Raw counts

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

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

There is no strong different, this is because the baseline expression values among genes are much larger than the relative expression difference between trisomy12 and wt samples.

Expression values centered by mean

meanExprChr12 <- rnaExprTab %>% filter(ChromID %in% "chr12") %>%
  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(meanExprChr12, aes(x=meanExpr, fill = cnv)) + geom_histogram(position = "identity", alpha=0.5) 

Now the difference is clearly visible.

Mean expression difference between trisomy12 and wt samples for the genes on chr12 and other chromosomes

plotTab <- rnaExprTab %>%
  group_by(id, cnv, ChromID) %>% summarise(meanExpr = mean(expr, na.rm=TRUE)) %>%
  ungroup() %>%
  spread(key = cnv, value = meanExpr) %>%
  mutate(diff = trisomy12-wt) %>%
  mutate(chr = ifelse(ChromID == "chr12","Chr12","Other"))

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

It’s also visible that most genes on chr12 tend to have higher expression in trisomy12 samples. While other genes follow normal distribution centered at 0.

Raw counts

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

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

Similar as in the RNA expression, the difference is not very visible when raw count is used.

Expression values centered by mean

meanExprChr12 <- protExprTab %>% filter(ChromID %in% "chr12") %>%
  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(meanExprChr12, aes(x=meanExpr, fill = cnv)) + geom_histogram(position = "identity", alpha=0.5) 

Now difference is clearly visible.

Mean protein expression difference between trisomy12 and wt samples for the genes on chr12 and other chromosomes

plotTab <- protExprTab %>%
  group_by(id, cnv, ChromID) %>% summarise(meanExpr = mean(expr, na.rm=TRUE)) %>%
  ungroup() %>%
  spread(key = cnv, value = meanExpr) %>%
  mutate(diff = trisomy12-wt) %>%
  mutate(chr = ifelse(ChromID == "chr12","Chr12","Other"))

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

The gene dosage effect is also clearly visible in the proteomic dataset

Here I use two ways to quantify the buffering effect:

1. A buffering score, which is based on the difference of log fold change between protein and rna dataset and the t-statistics of the differentially expressed RNAs. The purpose is to give the gene that show significant and strong RNA change, but little protein change a high buffering score. While the genes that do not show strong RNA expression change will have a score close to zero. And the genes that show both strong protein and RNA expression change a more negative score.

2. A categorical variable, “ifBuffer”, based on the the significance of differential expression. The genes that show both significant protein and RNA up-regulation are in the “noBuffer” group, while the genes that show significant RNA-up-regulation but no significant protein expression change are in the “Buffered” group. The “Enhanced” group contains the genes that do not show significant changes in RNA level but with significant changes in protein level. The buffering score can not differentiate this group and will categorize it as undetermined. But the genes in this group, although pretty rare, may also be potentially interesting. Other genes are in the “undetermined” group.

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

The buffering score and the categorical variable are related. Perhaps the buffering score can estimate more subtle effect, like the degree of buffering.

table(bufferTab$ifBuffer)

    Buffered     Enhanced     noBuffer undetermined 
          50           10          115           12 

Compare fold change in RNA expression and protein expression

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") 

Top 9 most buffered based on buffering score

geneList <- bufferTab$id[1:9]
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(trisomy12 = patMeta[match(patID, patMeta$Patient.ID),]$trisomy12)
  p <- ggplot(plotTab, aes(x=rnaExpr, y = protExpr)) +
    geom_point(aes(col=trisomy12)) + 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)

The 9 least buffered based on buffering score

geneList <- tail(bufferTab$id, n=9)
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(trisomy12 = patMeta[match(patID, patMeta$Patient.ID),]$trisomy12)
  p <- ggplot(plotTab, aes(x=rnaExpr, y = protExpr)) +
    geom_point(aes(col=trisomy12)) + 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 == "Enhanced")$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(trisomy12 = patMeta[match(patID, patMeta$Patient.ID),]$trisomy12)
  p <- ggplot(plotTab, aes(x=rnaExpr, y = protExpr)) +
    geom_point(aes(col=trisomy12)) + 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)