Are thechnical variables confounded with major genomic variabes?

# A tibble: 3 x 4
  genomic     technical      pval  p.adj
  <chr>       <chr>         <dbl>  <dbl>
1 trisomy12   proteinConc 0.00185 0.0778
2 IGHV.status processDate 0.0277  0.581 
3 del13q      operator    0.0494  0.581 

Trisomy12 has some impact on total protein concentration, but not very strong

Trisomy12 versus protein concentration

plotTab <- tibble(gene = geneTab$trisomy12, proteinConc = techTab$proteinConc)
ggplot(plotTab, aes(x=gene, y = proteinConc)) + geom_boxplot() + geom_point()

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

There’s a small trend that trisomy12 samples tend to have higher protein concentration.

Association between technical variables and priciple components of protein expression

plotMat <- assays(protCLL)[["QRILC"]]
pcRes <- prcomp(t(plotMat), center =TRUE, scale. = FALSE)$x

testRes <- lapply(colnames(pcRes), function(pc) {
  lapply(colnames(techTab), function(tech) {
    pcVar <- pcRes[,pc]
    techVar <- techTab[[tech]]
    res <- summary(aov(pcVar ~ techVar))
    p <- res[[1]][1,4]
    tibble(component = pc, technical = tech, pval = p)
  }) %>% bind_rows()
}) %>% bind_rows() %>% mutate(p.adj = p.adjust(pval, method = "BH")) %>%
  arrange(pval)
filter(testRes, p.adj < 0.1)

# A tibble: 1 x 4
  component technical        pval  p.adj
  <chr>     <chr>           <dbl>  <dbl>
1 PC26      proteinConc 0.0000741 0.0218

There are some principle components correlated with technical variables. But the PCs are not top PCs, suggesting the know technical factor do not have impact on the major trends of the dataset.

Associations between technical variables and individual protein expressions

Association test

techTab <- colData(protCLL)[,c("operator", "viability","batch","processDate","proteinConc","freeThawCycle")] %>%
  data.frame() %>%rownames_to_column("patID") %>% as_tibble() %>% mutate(processDate = as.character(processDate)) %>%
  mutate_if(is.character, as.factor)%>% mutate_at(vars(-patID), as.numeric)
testTab <- assays(protCLL)[["QRILC"]] %>% data.frame() %>% 
  rownames_to_column("id") %>% mutate(name = rowData(protCLL)[id,]$hgnc_symbol) %>%
  gather(key = "patID", value = "expr", -id, -name) %>%
  left_join(techTab, by ="patID") %>% gather(key = technical, value = value, -id, -name, -patID, -expr)

Warning: Column `patID` joining character vector and factor, coercing into
character vector

testRes <- filter(testTab, !is.na(value)) %>% 
  group_by(name, technical) %>% nest() %>%
  mutate(m = map(data, ~lm(expr~value,.))) %>%
  mutate(res = map(m, broom::tidy)) %>%
  unnest(res) %>% filter(term=="value") %>%
  ungroup() %>%
  mutate(p.adj = p.adjust(p.value, method = "BH"))

sumTab <- filter(testRes, p.adj < 0.1) %>% group_by(technical) %>% summarise(n=length(name)) 
ggplot(sumTab, aes(x=technical, y = n)) + geom_bar(stat = "identity") + coord_flip() +
  xlab("") + ylab("Number of significantly associated proteins")

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

There are more than 750 proteins whose abundance may be significantly affected by overall protein concentration.

Assocation P value histogram for each technical factor

ggplot(testRes, aes(x=p.value)) + geom_histogram() + facet_wrap(~technical) + xlim(0,1)

`stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

Warning: Removed 12 rows containing missing values (geom_bar).

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

Also based on the p-value histogram, only overall protein concentration may have potential impact on protein abundance detection.

List of proteins that show significant correlation with total protein concentration (10% FDR)

filter(testRes, technical == "proteinConc") %>% select(name, p.value, p.adj, estimate) %>%
  arrange(p.value) %>% filter(p.adj <= 0.1) %>%
  mutate_if(is.numeric, formatC, digits=2, format="e") %>%
  DT::datatable()

Will the correlation with RNA expression improve if we adjust for total protein concentration?

proteinConc <- techTab[match(colnames(proMat.qrilc), techTab$patID),]$proteinConc
corTab <- lapply(geneOverlap, function(n) {
  rna <- rnaMat[n,]
  pro.q <- proMat.qrilc[n,]
  p.single <- anova(lm(rna ~ pro.q))$`Pr(>F)`[1]
  p.multi <- car::Anova(lm(rna ~ pro.q + proteinConc))$`Pr(>F)`[1]
  tibble(name = n, corrected = c("no","yes"),
         p = c(p.single, p.multi))
}) %>% bind_rows() %>% mutate(p.adj = p.adjust(p, method = "BH")) %>% arrange(p)

Number of significant correlations VS FDR cut-off

plotTab <- lapply(seq(0,0.1, length.out = 100), function(fdr) {
  filTab <- dplyr::filter(corTab, p.adj < fdr) %>%
    group_by(corrected) %>% summarise(n = length(name)) %>% mutate(fdr = fdr)
}) %>% bind_rows()

ggplot(plotTab, aes(x=fdr, y = n, col = corrected))+ geom_line() +
  ylab("Number of significant correlations") +
  xlab("FDR cut-off")

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

Seems to improve the correlation a little, but not much. We can include this factor in association test later.

List of proteins that show significant correlation with viability (10% FDR)

filter(testRes, technical == "viability") %>% select(name, p.value, p.adj, estimate) %>%
  arrange(p.value) %>% filter(p.adj <= 0.1) %>%
  mutate_if(is.numeric, formatC, digits=2, format="e") %>%
  DT::datatable()

Hierarchical clustering

plotMat <- assays(protCLL)[["QRILC"]]

colAnno <- colData(protCLL)[,c("gender","IGHV.status","trisomy12")] %>%
  data.frame()

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

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

Some separation can be observed for trisomy12

PCA

pcOut <- prcomp(t(plotMat), center =TRUE, scale. = FALSE)
pcRes <- pcOut$x
eigs <- pcOut$sdev^2
varExp <- structure(eigs/sum(eigs),names = colnames(pcRes))

plotTab <- pcRes[,1:2] %>% data.frame() %>% cbind(colAnno[rownames(.),]) %>%
  rownames_to_column("patID") %>% as_tibble()

ggplot(plotTab, aes(x=PC1, y=PC2, col = IGHV.status, shape = trisomy12)) + geom_point(size=4) +
  xlab(sprintf("PC1 (%1.2f%%)",varExp[["PC1"]]*100)) +
  ylab(sprintf("PC2 (%1.2f%%)",varExp[["PC2"]]*100))

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

No strong separation can be observed.

Correlation PCs with trisomy12 and IGHV status

corTab <- lapply(colnames(pcRes),  function(pc) {
  ighvCor <- t.test(pcRes[,pc] ~ colAnno$IGHV.status)
  tri12Cor <- t.test(pcRes[,pc] ~ colAnno$trisomy12)
  tibble(PC = pc, 
         feature=c("IGHV", "trisomy12"),
         p = c(ighvCor$p.value, tri12Cor$p.value))
}) %>% bind_rows() %>% mutate(p.adj = p.adjust(p, method = "BH")) %>%
  filter(p <= 0.05) %>% arrange(p)
corTab

# A tibble: 4 x 4
  PC    feature          p        p.adj
  <chr> <chr>        <dbl>        <dbl>
1 PC3   trisomy12 1.72e-10 0.0000000168
2 PC4   IGHV      6.50e- 6 0.000319    
3 PC2   trisomy12 1.83e- 2 0.598       
4 PC1   IGHV      4.53e- 2 0.998       

PCA plot using PC3 and PC4

plotTab <- pcRes[,1:10] %>% data.frame() %>% cbind(colAnno[rownames(.),]) %>%
  rownames_to_column("patID") %>% as_tibble()

ggplot(plotTab, aes(x=PC3, y=PC4, col = IGHV.status, shape = trisomy12, label = patID)) + geom_point() + ggrepel::geom_text_repel() +
  xlab(sprintf("PC3 (%1.2f%%)",varExp[["PC3"]]*100)) +
  ylab(sprintf("PC4 (%1.2f%%)",varExp[["PC4"]]*100))

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

Using PC3 and PC4 better seperates IGHV and trisomy12.

Annotate PC1

proMat <- assays(protCLL)[["QRILC"]]
pc <- pcRes[,1][colnames(proMat)]
designMat <- model.matrix(~1+pc)
fit <- lmFit(proMat,  designMat)
fit2 <- eBayes(fit)
corRes.pc1 <- topTable(fit2, number ="all", adjust.method = "BH", coef = "pc") %>% rownames_to_column("id") %>%
    mutate(symbol = rowData(protCLL[id,])$hgnc_symbol)

hist(corRes.pc1$P.Value,breaks=100)

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

colAnno <- tibble(patID = colnames(proMat), PC1 = pcRes[colnames(proMat),1]) %>%
  arrange(PC1) %>% data.frame() %>% column_to_rownames("patID")

plotMat <- proMat[resTab.sig$id[1:100],rownames(colAnno)]
plotMat <- jyluMisc::mscale(plotMat, censor = 6)

pheatmap(plotMat, scale = "none", annotation_col = colAnno, clustering_method = "ward.D2",
         cluster_cols = FALSE,
         labels_row = resTab.sig$symbol[1:100], color = colorRampPalette(c("navy","white","firebrick"))(100),
         breaks = seq(-6,6, length.out = 101))

gmts = list(H= "../data/gmts/h.all.v6.2.symbols.gmt",
            C6 = "../data/gmts/c6.all.v6.2.symbols.gmt",
            KEGG = "../data/gmts/c2.cp.kegg.v6.2.symbols.gmt")

res <- runCamera(proMat, designMat, gmts$H, id = rowData(protCLL[rownames(proMat),])$hgnc_symbol)
res$enrichPlot

res <- runCamera(proMat, designMat, gmts$C6, id = rowData(protCLL[rownames(proMat),])$hgnc_symbol)
res$enrichPlot

Annotate PC2

proMat <- assays(protCLL)[["QRILC"]]
pc <- pcRes[,2][colnames(proMat)]
designMat <- model.matrix(~1+pc)
fit <- lmFit(proMat,  designMat)
fit2 <- eBayes(fit)
corRes.pc2 <- topTable(fit2, number ="all", adjust.method = "BH", coef = "pc") %>% rownames_to_column("id") %>%
    mutate(symbol = rowData(protCLL[id,])$hgnc_symbol)

hist(corRes.pc2$P.Value,breaks=100)

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

colAnno <- tibble(patID = colnames(proMat), PC2 = pcRes[colnames(proMat),2]) %>%
  arrange(PC2) %>% data.frame() %>% column_to_rownames("patID")

plotMat <- proMat[resTab.sig$id[1:100],rownames(colAnno)]
plotMat <- jyluMisc::mscale(plotMat, censor = 6)

pheatmap(plotMat, scale = "none", annotation_col = colAnno, clustering_method = "ward.D2",
         cluster_cols = FALSE,
         labels_row = resTab.sig$symbol[1:100], color = colorRampPalette(c("navy","white","firebrick"))(100),
         breaks = seq(-6,6, length.out = 101))

res <- runCamera(proMat, designMat, gmts$H, id = rowData(protCLL[rownames(proMat),])$hgnc_symbol)
res$enrichPlot

Association with lymphocyte count (potential contamination?)

load("~/CLLproject_jlu/ShinyApps/sampleTimeline//timeline.RData")
plotTab <- pcRes[,1:9] %>% data.frame()  %>%
  rownames_to_column("patID") %>% as_tibble() %>%
  mutate(sampleID = protCLL[,patID]$sampleID) %>%
  mutate(lymCount = sampleTab[match(sampleID, sampleTab$sampleID),]$leukCount) %>%
  gather(key = "pc", value = "val",-patID,-sampleID,-lymCount)

ggplot(plotTab, aes(x=val, y=lymCount)) + geom_point() + geom_smooth(method = "lm") +
  facet_wrap(~pc, ncol =3, scale = "free")

`geom_smooth()` using formula 'y ~ x'

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

corRes <- plotTab %>% group_by(pc) %>% nest() %>%
  mutate(m= map(data, ~cor.test(~ val+ lymCount,.))) %>%
  mutate(res = map(m, broom::tidy)) %>%
  unnest(res) %>% select(pc, estimate, p.value) %>%
  arrange(p.value)
corRes

# A tibble: 9 x 3
# Groups:   pc [9]
  pc    estimate p.value
  <chr>    <dbl>   <dbl>
1 PC1    -0.253   0.0789
2 PC3     0.228   0.115 
3 PC6     0.219   0.130 
4 PC8     0.206   0.156 
5 PC2    -0.194   0.182 
6 PC7    -0.183   0.209 
7 PC5     0.110   0.451 
8 PC9     0.0840  0.566 
9 PC4     0.0284  0.847 

Association wit %Lymphocyte PB

library(DBI)

con <- dbConnect(RPostgreSQL::PostgreSQL(),
                 dbname = "tumorbank",
                 host = "huber-vm01.embl.de",
                 user = "admin",
                 password = "bloodcancertumorbank")

PBtab <- tbl(con, "patient") %>% 
  left_join(tbl(con, "sample"), by = c(patid = "smppatidref")) %>% 
  left_join(tbl(con, "analysis"), by = c(smpid = "anlsmpidref")) %>%
  collect() %>%
  select(patpatientid, smpsampleid, smpleukocytes, smpsampledate, smppblymphocytes)

dbDisconnect(con)

[1] TRUE

plotTab <- pcRes[,1:9] %>% data.frame()  %>%
  rownames_to_column("patID") %>% as_tibble() %>%
  mutate(sampleID = protCLL[,patID]$sampleID) %>%
  mutate(perPB = PBtab[match(sampleID, PBtab$smpsampleid),]$smppblymphocytes) %>%
  gather(key = "pc", value = "val",-patID,-sampleID,-perPB)

ggplot(plotTab, aes(x=val, y=perPB)) + geom_point() + geom_smooth(method = "lm") +
  facet_wrap(~pc, ncol =3, scale = "free")

`geom_smooth()` using formula 'y ~ x'

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

corRes <- plotTab %>% group_by(pc) %>% nest() %>%
  mutate(m= map(data, ~cor.test(~ val+ perPB,.))) %>%
  mutate(res = map(m, broom::tidy)) %>%
  unnest(res) %>% select(pc, estimate, p.value) %>%
  arrange(p.value)
corRes

# A tibble: 9 x 3
# Groups:   pc [9]
  pc    estimate p.value
  <chr>    <dbl>   <dbl>
1 PC2   -0.388   0.00586
2 PC1   -0.360   0.0110 
3 PC3    0.171   0.240  
4 PC8    0.165   0.256  
5 PC4    0.0617  0.674  
6 PC9   -0.0415  0.777  
7 PC7    0.0235  0.872  
8 PC5   -0.00258 0.986  
9 PC6   -0.00207 0.989  

Will the correlation with transcriptomics increase if PC1&2 are regressed out?

pc1 <- pcRes[colnames(rnaMat),1]
pc2 <- pcRes[colnames(rnaMat),2]

corTab <- lapply(geneOverlap, function(n) {
  rna <- rnaMat[n,]
  pro.q <- proMat.qrilc[n,]
  p.single <- anova(lm(rna ~ pro.q))$`Pr(>F)`[1]
  p.multi <- car::Anova(lm(rna ~ pro.q + pc1 + pc2))$`Pr(>F)`[1]
  tibble(name = n, corrected = c("no","yes"),
         p = c(p.single, p.multi))
}) %>% bind_rows() %>% mutate(p.adj = p.adjust(p, method = "BH")) %>% arrange(p)

Number of significant correlations VS FDR cut-off

plotTab <- lapply(seq(0,0.1, length.out = 100), function(fdr) {
  filTab <- dplyr::filter(corTab, p.adj < fdr) %>%
    group_by(corrected) %>% summarise(n = length(name)) %>% mutate(fdr = fdr)
}) %>% bind_rows()

ggplot(plotTab, aes(x=fdr, y = n, col = corrected))+ geom_line() +
  ylab("Number of significant correlations") +
  xlab("FDR cut-off")

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

Seems to improve the correlation a little, but not much. We can include this factor in association test later.

Hierarchical clustering using reconstructed matrix

plotMat <- protMat.new

colAnno <- colData(protCLL)[,c("gender","IGHV.status","trisomy12")] %>%
  data.frame()

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

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

Better separation for trisomy12 can be observed

PCA using reconstructed matrix

pcOut <- prcomp(t(plotMat), center =TRUE, scale. = FALSE)
pcRes.new <- pcOut$x
eigs <- pcOut$sdev^2
varExp <- structure(eigs/sum(eigs),names = colnames(pcRes.new))

plotTab <- pcRes.new[,1:2] %>% data.frame() %>% cbind(colAnno[rownames(.),]) %>%
  rownames_to_column("patID") %>% as_tibble()

ggplot(plotTab, aes(x=PC1, y=PC2, col = IGHV.status, shape = trisomy12)) + geom_point(size=4) +
  xlab(sprintf("PC1 (%1.2f%%)",varExp[["PC1"]]*100)) +
  ylab(sprintf("PC2 (%1.2f%%)",varExp[["PC2"]]*100))

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