Identify pure trisomy12 related protein expression signature (onlyTri12 vs none)

#only trisomy12 vs none
protTest <- protSub[,protSub$group != "both"]
testMat <- assay(protTest)
group <-  factor(protTest$group, levels = c("none","onlyTri12"))
fit <- proDA(testMat, design = ~  group)
resTab.tri12 <- test_diff(fit, contrast = "grouponlyTri12") %>%
  dplyr::rename(id = name, logFC = diff, t=t_statistic,
                P.Value = pval, adj.P.Val = adj_pval) %>% 
  mutate(name = rowData(protTest[id,])$hgnc_symbol) %>%
  select(name, id, logFC, t, P.Value, adj.P.Val, n_obs) %>%  
  arrange(P.Value) %>%
  as_tibble()
resTab.tri12.sig <- filter(resTab.tri12, adj.P.Val <0.1)

For proteins associated with trisomy12 (10%), test the impact of trisomy19 on them.

protTest <- protSub[resTab.tri12.sig$id,protSub$group != "none"]
testMat <- assay(protTest)
group <-  factor(protTest$group, levels = c("both","onlyTri12"))
fit <- proDA(testMat, design = ~ group)
resTab.tri19 <- test_diff(fit, contrast = "grouponlyTri12") %>%
  dplyr::rename(id = name, logFC = diff, t=t_statistic,
                P.Value = pval, adj.P.Val = adj_pval) %>% 
  mutate(name = rowData(protTest[id,])$hgnc_symbol) %>%
  select(name, id, logFC, t, P.Value, adj.P.Val,n_obs) %>%  
  arrange(P.Value) %>%
  as_tibble()

Compare the two results

resList <- list()
resList[["tri12_up"]] <- filter(resTab.tri12.sig,  logFC >0)$name
resList[["tri12_down"]] <- filter(resTab.tri12.sig,  logFC <0)$name
resList[["tri19_up"]] <- filter(resTab.tri19, P.Value <= 0.01, logFC >0)$name
resList[["tri19_down"]] <- filter(resTab.tri19, P.Value <= 0.01, logFC <0)$name

upset(fromList(resList))

Plot the expressions of proteins that show additional change in trisomy19 samples (raw P-value < 0.01)

proteinList <- intersect(filter(resTab.tri19, P.Value <=0.01)$name, resTab.tri12.sig$name)

protTab <- sumToTiday(protSub,"patID") %>% mutate(name = hgnc_symbol)

plotList <- lapply(proteinList, function(symbol) {
  plotTab <- filter(protTab, name == symbol) %>%
    mutate(group = factor(group, levels = c("none","onlyTri12","both")))
  chrLoc <- unique(filter(protTab, hgnc_symbol == symbol)$chromosome_name)
  ggplot(plotTab, aes(x=group, y = count)) + geom_boxplot(aes(fill = group)) + geom_point() +
    ggtitle(sprintf("%s (on chr%s)",symbol, chrLoc)) +
    ggrepel::geom_text_repel(aes(label = colID))
})
plot_grid(plotlist= plotList, ncol =3)

As we don’t have samples with only trisomy19 but not trisomy12, we don’t have a signature of “pure” trisomy19 effect on expression. So those are not “epistatic interactions”. But can be considered as additional effect of trisomy19 on the genes related to trisomy12.

Sample subsetting

Subset sample: select M-CLLs with annotation for trisomy12 and trisomy19

dds$IGHV <- patMeta[match(colnames(dds),patMeta$Patient.ID),]$IGHV.status
dds$trisomy12 <- patMeta[match(colnames(dds),patMeta$Patient.ID),]$trisomy12
dds$trisomy19 <- patMeta[match(colnames(dds),patMeta$Patient.ID),]$trisomy19
ddsSub <- dds[,dds$diag == "CLL" & dds$IGHV %in% "M" & !is.na(dds$trisomy12) & !is.na(dds$trisomy19)]

Summary of trisomy12 and trisomy19 status

annoTab <- tibble(patID = colnames(ddsSub),
                  tri12 = ifelse(ddsSub$trisomy12 == 1, "tri12","wt"),
                  tri19 = ifelse(ddsSub$trisomy19 ==1, "tri19","wt")) %>%
  mutate(group = ifelse(tri12 == "tri12", ifelse(tri19 == "tri19","both","onlyTri12"),"none"))
table(annoTab$tri12, annoTab$tri19)

       
        tri19 wt
  tri12     5  7
  wt        0 97

There are too many samples without either trisomy19 and trisomy12. The sample imbalance will cause bias in hypothesis testing. In addition, the other recurrent alterations in wildtype group may also complicate the problem.

To ensure the compatibility and reduce noise, in the below analysis, I will:

1. Check the distribution of several recurrent mutations that impact gene expression (based on Almut’s paper), namely: del13q, SF3B1, del17p, TP53, del11q, BRAF, gain8q, del8p, NOTCH1, MED12, ATM.

2. Get mutations that occurred in non-wildtype samples (with trisomy12, trisomy19 or both).

3. Remove wild type samples that contain the recurrent mutations that do not occur in non-wildtype samples.

4. Block for the mutations that occur in both wild type and non-wildtype samples.

#Get the distribution of recurrent mutations in RNAseq samples
patSum <- select(patMeta, Patient.ID, del13q, SF3B1, del17p, TP53, del11q, BRAF, gain8q, del8p, NOTCH1, MED12, ATM) %>%
  filter(Patient.ID %in% annoTab$patID) %>%
  gather(key = "gene",value = "status",-Patient.ID) %>%
  mutate(status = as.integer(status)) %>% filter(status %in% 1)

Which mutations occur in non-wildtype samples?

geneKeep <- unique(filter(patSum, Patient.ID %in% filter(annoTab, group != "none")$patID)$gene)
geneKeep

[1] "del13q" "SF3B1"  "del17p" "TP53"  

Remove wildtype samples with recurrent mutations except for those four: del13q, SF3B1, del17p and TP53. Those genes will be blocked when doing hypothesis test

patSum.other <- filter(patSum, !gene %in% geneKeep)
annoTab <- filter(annoTab, !(group == "none" & patID %in% patSum.other$Patient.ID))

Sample summary after filtering

table(annoTab$tri12, annoTab$tri19)

       
        tri19 wt
  tri12     5  7
  wt        0 82

Filtering transcripts

ddsSub <- ddsSub[,annoTab$patID]
ddsSub$group <- annoTab$group
geneTab <- patMeta[match(ddsSub$PatID, patMeta$Patient.ID),c("Patient.ID", "del13q","SF3B1","del17p","TP53")] %>%
  data.frame() %>% column_to_rownames("Patient.ID")
colData(ddsSub) <- cbind(colData(ddsSub),geneTab)

#filter out none protein coding genes and gene on sex chromosome
ddsSub<-ddsSub[rowData(ddsSub)$biotype %in% "protein_coding",]
ddsSub <- ddsSub[! rowData(ddsSub)$symbol %in% c("",NA),]

##vst
ddsSub.vst <- varianceStabilizingTransformation(ddsSub)
dim(ddsSub)

[1] 20074    94

Table of associations with 10% FDR

corRes.sig <- resTab.rna %>% filter(adj.P.Val <0.1)
corRes.sig %>% select(id, name, logFC, P.Value, adj.P.Val) %>%
  mutate_if(is.numeric, formatC, digits=2, format="e") %>% DT::datatable()

Volcano plot (1% FDR)

plotTab <- resTab.rna %>%
  mutate(onChr19 = ifelse(chromosome == "19","yes","no"))
plotVolcano(plotTab, fdrCut =0.01, x_lab="log2FoldChange", 
            plotTitle = "trisomy19", ifLabel = TRUE, colLabel = "onChr19")

Warning: Removed 50 rows containing missing values (geom_point).

Plot top 9 most differentially expressed proteins

exprTab <- assay(ddsSub.vst[corRes.sig$id[1:9],colnames(ddsTest)]) %>% data.frame() %>%
  rownames_to_column("id") %>% 
  mutate(name = rowData(ddsSub[id])$symbol) %>%
  gather(key = "patID", value = "expression", -id, -name) %>%
  mutate(trisomy19 = ddsSub[,patID]$trisomy19)

plotTab <- filter(exprTab, name %in% corRes.sig$name[1:9]) 
ggplot(plotTab, aes(x=trisomy19, y =  expression)) + geom_boxplot(aes(fill = trisomy19)) + geom_point() +
  facet_wrap(~name, scale = "free") +
  ggrepel::geom_text_repel(aes(label = patID))

Heatmap of differentially expressed protein (1% FDR)

proList <- filter(corRes.sig, adj.P.Val < 0.01) %>% arrange(desc(logFC)) %>%
  filter(!is.na(name)) %>% distinct(name, .keep_all = TRUE) %>% pull(id)
plotMat <- assay(ddsSub.vst)[proList,colnames(ddsTest)]
rownames(plotMat) <- rowData(ddsSub.vst[proList,])$symbol

colAnno <- colData(ddsSub.vst)[,c("trisomy19","trisomy12","del13q","SF3B1","TP53", "del17p")] %>%
  data.frame()

rowAnno <- filter(corRes.sig, adj.P.Val < 0.01) %>% 
  mutate(onChr19 = ifelse(chromosome == "19","yes","no")) %>%
  select(name, onChr19) %>% data.frame() %>% column_to_rownames("name")

plotMat <- jyluMisc::mscale(plotMat, censor = 6)
pheatmap(plotMat, scale = "none", annotation_col = colAnno,  annotation_row = rowAnno, cluster_rows = FALSE,
         clustering_method = "ward.D2",
         color = colorRampPalette(c("navy","white","firebrick"))(100),
         breaks = seq(-6,6, length.out = 101))

Trisomy19 samples form its own cluster.

Enrichment analysis (10% FDR)

gmts = list(H= "../data/gmts/h.all.v6.2.symbols.gmt",
            KEGG = "../data/gmts/c2.cp.kegg.v6.2.symbols.gmt")
inputTab <- resTab.rna %>% filter(P.Value <0.05) %>%
  distinct(name, .keep_all = TRUE) %>%
  select(name, stat) %>% data.frame() %>% column_to_rownames("name")
enRes <- list()
enRes[["HALLMARK"]] <- runGSEA(inputTab, gmts$H, "page")
enRes[["KEGG"]] <- runGSEA(inputTab, gmts$KEGG, "page")
p <- plotEnrichmentBar(enRes, pCut =0.1, ifFDR= TRUE)
#pdf("tri12Enrich.pdf", height = 15, width = 6)
plot(p)

#dev.off()

Heatmap plots for enriched pathways (genes with raw P <0.05)

HALLMARK_TNFA_SIGNALING_VIA_NFKB

plotSetHeatmap <- function(geneSigTab, setDir, setName, ddsObj, colAnno, scale = TRUE, annoChrom = NULL) {
  geneList <- loadGSC(setDir)[["gsc"]][[setName]]
  sigGene <- filter(geneSigTab, name %in% geneList) %>%
    arrange(desc(logFC))
  colAnno <- colAnno[order(colAnno[,1]),,drop = FALSE]
  plotMat <- assay(ddsObj)[sigGene[["id"]],rownames(colAnno)]
  if (scale) {
    #calculate z-score and sensor
    plotMat <- t(scale(t(plotMat)))
    plotMat[plotMat >= 6] <- 6
    plotMat[plotMat <= -6] <- -6
  }
  
  if (!is.null(annoChrom)) {
    rowAnno <- mutate(geneSigTab, onChr = ifelse(chromosome %in% annoChrom, "yes","no")) %>%
      select(id, onChr) %>% data.frame() %>% column_to_rownames("id")
    colnames(rowAnno) <- paste0("onChr",annoChrom)
    
  } else rowAnno <- NULL
  
  pheatmap(plotMat, color = colorRampPalette(c("navy","white","firebrick"))(100),
         cluster_cols = FALSE, cluster_rows = FALSE,
         annotation_col = colAnno, labels_row = sigGene$name, annotation_row = rowAnno,
         show_colnames = TRUE, fontsize_row = 8, breaks = seq(-6,6, length.out = 101), treeheight_row = 0,
         border_color = NA, main = setName)

}
colAnnoSub <- colAnno[colnames(ddsTest),c("trisomy19"),drop=FALSE]
resTab.rna.sig <- filter(resTab.rna, P.Value < 0.05)
plotSetHeatmap(resTab.rna.sig, gmts$H, "HALLMARK_TNFA_SIGNALING_VIA_NFKB", ddsSub.vst[,colnames(ddsTest)], 
               colAnno = colAnnoSub, annoChrom = "19")

HALLMARK_MTORC1_SIGNALING

plotSetHeatmap(resTab.rna.sig, gmts$H, "HALLMARK_MTORC1_SIGNALING", ddsSub.vst[,colnames(ddsTest)], colAnno = colAnnoSub, annoChrom = "19")

HALLMARK_ALLOGRAFT_REJECTION

plotSetHeatmap(resTab.rna.sig, gmts$H, "HALLMARK_ALLOGRAFT_REJECTION", ddsSub.vst[,colnames(ddsTest)], colAnno = colAnnoSub, annoChrom = "19")

HALLMARK_INTERFERON_GAMMA_RESPONSE

plotSetHeatmap(resTab.rna.sig, gmts$H, "HALLMARK_INTERFERON_GAMMA_RESPONSE", ddsSub.vst[,colnames(ddsTest)], colAnno = colAnnoSub, annoChrom = 19)

HALLMARK_APOPTOSIS

plotSetHeatmap(resTab.rna.sig, gmts$H, "HALLMARK_APOPTOSIS", ddsSub.vst[,colnames(ddsTest)], 
               colAnno = colAnnoSub, annoChrom = "19")

Compare the results with proteomic analysis

Upset plot (I use cut-off of raw p value < 0.01 here)

compareList <- list(Protein_up = filter(resTab, P.Value < 0.01, logFC >0)$name,
                    Protein_down = filter(resTab, P.Value < 0.01, logFC <0)$name,
                    RNA_up = filter(resTab.rna, P.Value < 0.01, logFC >0)$name,
                    RNA_down = filter(resTab.rna, P.Value < 0.01, logFC <0)$name)
upset(fromList(compareList))

intersect(compareList$Protein_up, compareList$RNA_up)

 [1] "GNAO1"   "IGF2BP3" "RRAS"    "MEF2C"   "CYB5A"   "SYNE2"   "RANBP3" 
 [8] "TMEM205" "FCER2"   "DPP9"    "GPX4"   

intersect(compareList$Protein_down, compareList$RNA_down)

[1] "BMF"     "SAMSN1"  "GPM6A"   "ALDH3A2" "BASP1"   "IRF5"    "EZR"    

Identify pure trisomy12 related protein expression signature (onlyTri12 vs none)

#only trisomy12 vs none
ddsTest <- ddsSub[,ddsSub$group != "both"]
ddsTest$group <-  factor(ddsTest$group, levels = c("none","onlyTri12"))

##filter out low count genes
minrs <- 100
rs  <- rowSums(counts(ddsTest, normalized = TRUE))
ddsTest<-ddsTest[ rs >= minrs, ]

design(ddsTest) <- ~ del13q + SF3B1 + del17p + TP53 + group
ddsTest <- DESeq(ddsTest)

using pre-existing size factors

estimating dispersions

gene-wise dispersion estimates

mean-dispersion relationship

final dispersion estimates

fitting model and testing

-- replacing outliers and refitting for 482 genes
-- DESeq argument 'minReplicatesForReplace' = 7 
-- original counts are preserved in counts(dds)

estimating dispersions

fitting model and testing

resTab.tri12 <- results(ddsTest, contrast = c("group","onlyTri12","none"), tidy = TRUE) %>%
  dplyr::rename(id = row, logFC = log2FoldChange, t=stat,
                P.Value = pvalue, adj.P.Val = padj) %>% 
  mutate(name = rowData(dds[id,])$symbol) %>%
  select(name, id, logFC, t, P.Value, adj.P.Val) %>%  
  arrange(P.Value) %>%
  as_tibble()
resTab.tri12.sig <- filter(resTab.tri12, adj.P.Val < 0.1)

For proteins associated with trisomy12 (10%), test the impact of trisomy19 on them.

#trisomy19+trisomy12 signature
ddsTest <- ddsSub[,ddsSub$group != "none"]
ddsTest$group <-  factor(ddsTest$group, levels = c("onlyTri12","both"))

##filter out low count genes
minrs <- 100
rs  <- rowSums(counts(ddsTest, normalized = TRUE))
ddsTest<-ddsTest[ rs >= minrs, ]

design(ddsTest) <- ~ del13q + SF3B1 + del17p + TP53 + group
ddsTest <- DESeq(ddsTest)

using pre-existing size factors

estimating dispersions

gene-wise dispersion estimates

mean-dispersion relationship

final dispersion estimates

fitting model and testing

resTab.tri19 <- results(ddsTest, contrast = c("group","both","onlyTri12"), tidy = TRUE) %>%
  dplyr::rename(id = row, logFC = log2FoldChange, t=stat,
                P.Value = pvalue, adj.P.Val = padj) %>% 
  mutate(name = rowData(dds[id,])$symbol) %>%
  select(name, id, logFC, t, P.Value, adj.P.Val) %>%  
  arrange(P.Value) %>%
  as_tibble()
resTab.tri19.sig <- filter(resTab.tri19, adj.P.Val < 0.1)

Compare the two results

resList <- list()
resList[["tri12_up"]] <- filter(resTab.tri12.sig,  logFC >0)$name
resList[["tri12_down"]] <- filter(resTab.tri12.sig,  logFC <0)$name
resList[["tri19_up"]] <- filter(resTab.tri19.sig, logFC >0)$name
resList[["tri19_down"]] <- filter(resTab.tri19.sig, logFC <0)$name
upset(fromList(resList))

compareTab <- full_join(resTab.tri12.sig, resTab.tri19.sig, by =c("name","id")) %>%
  select(name, id, logFC.x, P.Value.x, adj.P.Val.x, logFC.y, P.Value.y, adj.P.Val.y) %>%
  dplyr::rename(logFC.tri12 = logFC.x, logFC.tri19 = logFC.y,
         P.Value.tri12 = P.Value.x, P.Value.tri19 = P.Value.y,
         adj.P.Val.tri12 = adj.P.Val.x, adj.P.Val.tri19 = adj.P.Val.y) %>%
  mutate(chr = rowData(ddsSub[id,])$chromosome)

showList <- filter(compareTab, logFC.tri12 > 0 & logFC.tri19 >0) %>%
  mutate_if(is.numeric, formatC, digits=2)
DT::datatable(showList)

plotList <- lapply(seq(nrow(showList)), function(i) {
  id <- showList[i,]$id
  plotTab <- tibble(patID = ddsSub.vst$PatID,
                    count = assay(ddsSub.vst)[id,],
                    group = factor(ddsSub.vst$group,levels = c("none","onlyTri12","both")))
  
  chrLoc <- showList[i,]$chr
  ggplot(plotTab, aes(x=group, y = count)) + geom_boxplot(aes(fill = group)) + geom_point() +
    ggtitle(sprintf("%s (on chr%s)",showList[i,]$name, chrLoc)) +
    ggrepel::geom_text_repel(data = filter(plotTab, group %in% c("onlyTri12","both")),aes(label = patID))+
    theme(legend.position ="none")
})
plot_grid(plotlist= plotList, ncol =2)

showList <- filter(compareTab, logFC.tri12 < 0 & logFC.tri19 <0) %>%
  mutate_if(is.numeric, formatC, digits=2)
DT::datatable(showList)

plotList <- lapply(seq(nrow(showList)), function(i) {
  id <- showList[i,]$id
  plotTab <- tibble(patID = ddsSub.vst$PatID,
                    count = assay(ddsSub.vst)[id,],
                    group = factor(ddsSub.vst$group,levels = c("none","onlyTri12","both")))
  
  chrLoc <- showList[i,]$chr
  ggplot(plotTab, aes(x=group, y = count)) + geom_boxplot(aes(fill = group)) + geom_point() +
    ggtitle(sprintf("%s (on chr%s)",showList[i,]$name, chrLoc)) +
    ggrepel::geom_text_repel(data = filter(plotTab, group %in% c("onlyTri12","both")),aes(label = patID))+
    theme(legend.position ="none")
})
plot_grid(plotlist= plotList, ncol =2)

showList <- filter(compareTab, logFC.tri12 > 0 & logFC.tri19 <0) %>%
  arrange(desc(abs(logFC.tri19-logFC.tri12))) %>%
  mutate_if(is.numeric, formatC, digits=2)
DT::datatable(showList)

plotList <- lapply(seq(9), function(i) {
  id <- showList[i,]$id
  plotTab <- tibble(patID = ddsSub.vst$PatID,
                    count = assay(ddsSub.vst)[id,],
                    group = factor(ddsSub.vst$group,levels = c("none","onlyTri12","both")))
  
  chrLoc <- showList[i,]$chr
  ggplot(plotTab, aes(x=group, y = count)) + geom_boxplot(aes(fill = group)) + geom_point() +
    ggtitle(sprintf("%s (on chr%s)",showList[i,]$name, chrLoc)) +    
    ggrepel::geom_text_repel(data = filter(plotTab, group %in% c("onlyTri12","both")),aes(label = patID))+
    theme(legend.position ="none")
})
plot_grid(plotlist= plotList, ncol =3)

showList <- filter(compareTab, logFC.tri12 < 0 & logFC.tri19 > 0) %>%
  arrange(desc(abs(logFC.tri19-logFC.tri12))) %>%
  mutate_if(is.numeric, formatC, digits=2)
DT::datatable(showList)

plotList <- lapply(seq(9), function(i) {
  id <- showList[i,]$id
  plotTab <- tibble(patID = ddsSub.vst$PatID,
                    count = assay(ddsSub.vst)[id,],
                    group = factor(ddsSub.vst$group,levels = c("none","onlyTri12","both")))
  
  chrLoc <- showList[i,]$chr
  ggplot(plotTab, aes(x=group, y = count)) + geom_boxplot(aes(fill = group)) + geom_point() +
    ggtitle(sprintf("%s (on chr%s)",showList[i,]$name, chrLoc)) +
    ggrepel::geom_text_repel(data = filter(plotTab, group %in% c("onlyTri12","both")),aes(label = patID))+
    theme(legend.position ="none")
})
plot_grid(plotlist= plotList, ncol =3)

Table of associations with raw P-value < 0.05

corRes.sig <- resTab.lm %>% filter(P.Value <0.05)
corRes.sig %>% select(id, Drug, conc, logFC, P.Value, adj.P.Val) %>%
  mutate_if(is.numeric, formatC, digits=2) %>% DT::datatable()

Plot of significant assocations (raw P <0.05)

plotTab <- filter(viabTab, id %in% corRes.sig$id, trisomy12 %in% 1) %>%
  mutate(id = factor(id, levels = corRes.sig$id))
ggplot(plotTab, aes(x=trisomy19, y =  viab)) + geom_boxplot(aes(fill = trisomy19)) + geom_point() +
  facet_wrap(~id, scale = "free", ncol=3) + ggrepel::geom_text_repel(aes(label = patientID))

Although the individual P-Values are not very significant, but there’s a general trend that additional presence trisomy19 lead to drug resistance. But most of those associations happen at low concentrations and the overall effect size is not large. This maybe due to the baseline survival difference of samples with trisomy19.