Compare CVs between individuals, remove PCs

Objective

We filted the first PC from the data and computed per gene CVs based on the filtered data. We observed that the adjusted CV profiles differed significantly between the three individuals.

Also, interestingly, pairwise comparisons revealed that the adjusted CV profile of individual 19239 is is significantly different from the other two individuals (19239 vs. 19098: p < 1e-9; 19239 vs. 19101: p < .007; 19098 vs. 19101: p = .0474. Note that individual 19239 is the one outlier in the PCA plots. Coincidently, PC2 separates 19239 from the other two individuals. However, we observed little correlations between PC2 gene loadings and data-wide CV (before/after filter the first PC), suggesting that factors other than coefficient of variation affects PC2.

Set up

library("data.table")
library("dplyr")
library("limma")
library("edgeR")
library("ggplot2")
library("grid")
theme_set(theme_bw(base_size = 12))
source("functions.R")

Prepare data

Input annotation of only QC-filtered single cells. Remove NA19098.r2

anno_qc_filter <- read.table("../data/annotation-filter.txt", header = TRUE,
                             stringsAsFactors = FALSE)

Import endogeneous gene molecule counts that are QC-filtered, CPM-normalized, ERCC-normalized, and also processed to remove unwanted variation from batch effet. ERCC genes are removed from this file.

molecules_ENSG <- read.table("../data/molecules-final.txt", header = TRUE, stringsAsFactors = FALSE)

Input moleclule counts before log2 CPM transformation. This file is used to compute percent zero-count cells per sample.

molecules_sparse <- read.table("../data/molecules-filter.txt", header = TRUE, stringsAsFactors = FALSE)

molecules_sparse <- molecules_sparse[grep("ENSG", rownames(molecules_sparse)), ]
stopifnot( all.equal(rownames(molecules_ENSG), rownames(molecules_sparse)) )

SVD/PCA

library(matrixStats)
centered_ENSG <- molecules_ENSG - rowMeans(molecules_ENSG)
svd_all <- svd( centered_ENSG )
head( (svd_all$d^2)/(sum(svd_all$d^2) ) )

[1] 0.137852451 0.018414108 0.009950908 0.005890242 0.003358104 0.003204841

# reproduce SVD results with PCA
# pca_results <- prcomp( t(centered_ENSG) )
# head((pca_results$sdev^2)/(sum(pca_results$sdev^2)) )

par(mfrow = c(2,2))
plot(x = rowMeans(molecules_ENSG), y = svd_all$u[ ,1],
     xlab = "mean molecule counts", ylab = "first eigen-cell", cex = .4)
plot(x = rowMeans(molecules_ENSG), y = svd_all$u[ ,2],
     xlab = "mean molecule counts", ylab = "second eigen-cell", cex = .4)
plot(x = rowMeans(molecules_ENSG), y = svd_all$u[ ,3],
     xlab = "mean molecule counts", ylab = "third eigen-cell", cex = .4)
plot(x = rowMeans(molecules_ENSG), y = svd_all$u[ ,4],
     xlab = "mean molecule counts", ylab = "third eigen-cell", cex = .4)

cor(rowMeans(molecules_ENSG), y = svd_all$u[ ,1], method ="spearman")

[1] -0.820101

1. First PC ~ mean gene expression levels/percent genes detected
2. Second PC ?
3. Third PC ~ differences between the three individuals…

library(gridExtra)
p1 <- ggplot(data.frame(pc1 = svd_all$v[,1],
                        pc2 = svd_all$v[,2],
                        individual = anno_qc_filter$individual),
              aes(x = pc1, y = pc2, col = individual)) +
        geom_point()

p2 <- ggplot(data.frame(pc2 = svd_all$v[,2],
                        pc3 = svd_all$v[,3],
                        individual = anno_qc_filter$individual),
            aes(x = pc3, y = pc2, col = individual) ) +
        geom_point() 

grid.arrange(p1, p2, ncol = 2)

Remove the first PC

filtered_data <- with(svd_all, u %*% diag( c(0, d[-1]) ) %*% t(v))

Compute CV of the filtered data

cv_filtered <- lapply(1:3, function(ii_individual) {
    individuals <- unique(anno_qc_filter$individual)
    counts <- filtered_data[ , anno_qc_filter$individual == individuals [ii_individual]]
    means <- apply(counts, 1, mean)
    sds <- apply(counts, 1, sd)
    cv <- sds/means
    return(cv)
})
names(cv_filtered) <- unique(anno_qc_filter$individual)
cv_filtered <- do.call(cbind, cv_filtered)
rownames(cv_filtered) <- rownames(molecules_ENSG)
friedman.test(cv_filtered)

    Friedman rank sum test

data:  cv_filtered
Friedman chi-squared = 14.327, df = 2, p-value = 0.0007744

friedman.test(cv_filtered[ , c(1,2)])

    Friedman rank sum test

data:  cv_filtered[, c(1, 2)]
Friedman chi-squared = 4.2544, df = 1, p-value = 0.03915

friedman.test(cv_filtered[ , c(1,3)])

    Friedman rank sum test

data:  cv_filtered[, c(1, 3)]
Friedman chi-squared = 15.758, df = 1, p-value = 7.2e-05

friedman.test(cv_filtered[ , c(2,3)])

    Friedman rank sum test

data:  cv_filtered[, c(2, 3)]
Friedman chi-squared = 0.045816, df = 1, p-value = 0.8305

Cell-cycle genes

Import cell-cycle genes.

cellcycle_genes <- read.table("../data/cellcyclegenes.txt", sep = "\t",
                              header = TRUE)
dim(cellcycle_genes)

[1] 151   5

ii_cellcycle <- rownames(cv_filtered) %in% unlist(cellcycle_genes)
friedman.test(cv_filtered[which(ii_cellcycle), ])

    Friedman rank sum test

data:  cv_filtered[which(ii_cellcycle), ]
Friedman chi-squared = 15.471, df = 2, p-value = 0.0004369

Non-cell-cycle genes

friedman.test(cv_filtered[which(!ii_cellcycle), ] )

    Friedman rank sum test

data:  cv_filtered[which(!ii_cellcycle), ]
Friedman chi-squared = 9.4449, df = 2, p-value = 0.008893

Pluripotent genes

pluripotent_genes <- read.table("../data/pluripotency-genes.txt", sep = "\t",
                                header = TRUE)$To

Friedman test.

ii_pluripotent <- (rownames(cv_filtered) %in% unlist(pluripotent_genes))
pca_all_pluripotent <- cv_filtered[which(ii_pluripotent), ]
friedman.test(pca_all_pluripotent)

    Friedman rank sum test

data:  pca_all_pluripotent
Friedman chi-squared = 9.1429, df = 2, p-value = 0.01034

Find driving genes of PC2

# compute data-wide cv
data_cv <- apply(molecules_ENSG, 1, sd)/apply(molecules_ENSG, 1, mean)
par(mfrow = c(2,2))
plot(x = data_cv, y = svd_all$u[ ,1])
plot(x = data_cv, y = svd_all$u[ ,2])
plot(x = data_cv, y = svd_all$u[ ,3])
plot(x = data_cv, y = svd_all$u[ ,4])

# compute data-wide cv based on the filtered data
filtered_data_cv <- apply(filtered_data, 1, sd)/apply(filtered_data, 1, mean)
par(mfrow = c(2,2))
plot(x = filtered_data_cv, y = svd_all$u[ ,1])
plot(x = filtered_data_cv, y = svd_all$u[ ,2])
plot(x = filtered_data_cv, y = svd_all$u[ ,3])
plot(x = filtered_data_cv, y = svd_all$u[ ,4])

pc2_genes <- svd_all$u[,2]

# find genes outside of 3 mean +/- sd
mean_pc2 <- mean(pc2_genes)
sd_pc2 <- sd(pc2_genes)
ii_extreme <- which( (pc2_genes > mean_pc2 + 3*sd_pc2)|(pc2_genes < mean_pc2 - 3*sd_pc2) )
length(ii_extreme)

[1] 107

library(broman)
par(mfrow = c(1,1))
xy <- cbind(c(1:length(pc2_genes)), pc2_genes)
plot(xy, cex = .4, pch = 16,
     main = paste("PC2 gene loadings \n", 
                  length(ii_extreme), "genes"), xlab = "genes", ylab = "PC2 loadings")
points(xy[ii_extreme, ], pch = 1, cex = .7, col = brocolors("crayon")["Red Orange"])

Annotation

library("biomaRt")
ensembl <- useMart(host = "grch37.ensembl.org",
                   biomart = "ENSEMBL_MART_ENSEMBL",
                   dataset = "hsapiens_gene_ensembl")
pc2_gene_info <- getBM(attributes = c("ensembl_gene_id", "chromosome_name",
                                       "external_gene_name", "transcript_count",
                                       "description"),
                        filters = "ensembl_gene_id",
                        values = rownames(molecules_ENSG[ii_extreme, ]),
                        mart = ensembl)
kable(pc2_gene_info[ , c(1, 2, 3, 5)])

if (file.exists("rda/svd-filtered/pc2-cv.rda")) {
  load("rda/svd-filtered/pc2-cv.rda")  
} else {
  library(Humanzee)
  go_pc2 <- GOtest(my_ensembl_gene_universe = rownames(molecules_ENSG),
                my_ensembl_gene_test = rownames(molecules_ENSG)[ii_extreme],
                pval_cutoff = 1, ontology=c("BP","CC","MF") )
  save(go_pc2, file = "rda/svd-filtered/pc2-cv.rda")
}

Extract terms

if (file.exists("rda/svd-filtered/go-pc2-terms.rda")) {
  load("rda/svd-filtered/go-pc2-terms.rda")  
} else {
  # Biological process
  goterms_bp <- summary(go_pc2$GO$BP, pvalue = 1)
  goterms_bp <- data.frame(ID = goterms_bp[[1]],
                           Pvalue = goterms_bp[[2]],
                           Terms = goterms_bp[[7]])
  goterms_bp <- goterms_bp[order(goterms_bp$Pvalue), ]
  
  # Cellular component
  goterms_cc <- summary(go_pc2$GO$CC, pvalue = 1)
  goterms_cc <- data.frame(ID = goterms_cc[[1]],
                           Pvalue = goterms_cc[[2]],
                           Terms = goterms_cc[[7]])
  goterms_cc <- goterms_cc[order(goterms_cc$Pvalue), ]

  # Molecular function
  goterms_mf <- summary(go_pc2$GO$MF, pvalue = 1)
  goterms_mf <- data.frame(ID = goterms_mf[[1]],
                           Pvalue = goterms_mf[[2]],
                           Terms = goterms_mf[[7]])
  goterms_mf <- goterms_mf[order(goterms_mf$Pvalue), ]

  save(goterms_bp, goterms_cc, goterms_mf, 
       file = "rda/svd-filtered/go-pc2-terms.rda")
}

head(goterms_bp, 20)

           ID       Pvalue
1  GO:0010043 4.233878e-05
2  GO:0071294 6.411339e-05
3  GO:1902305 1.795128e-04
4  GO:1990267 2.651255e-04
5  GO:0048669 3.380583e-04
6  GO:1903273 3.380583e-04
7  GO:1903275 3.380583e-04
8  GO:1903276 3.380583e-04
9  GO:1903278 3.380583e-04
10 GO:1903286 3.380583e-04
11 GO:1903288 3.380583e-04
12 GO:0045926 3.525005e-04
13 GO:0050794 3.846393e-04
14 GO:0050789 4.175013e-04
15 GO:0002028 6.899663e-04
16 GO:0006814 1.072938e-03
17 GO:0007379 1.111173e-03
18 GO:2000649 1.376235e-03
19 GO:0035725 1.428822e-03
20 GO:0003008 1.477490e-03
                                                         Terms
1                                         response to zinc ion
2                                cellular response to zinc ion
3             regulation of sodium ion transmembrane transport
4                    response to transition metal nanoparticle
5                    collateral sprouting in absence of injury
6                              regulation of sodium ion export
7                     positive regulation of sodium ion export
8                    regulation of sodium ion export from cell
9           positive regulation of sodium ion export from cell
10                          regulation of potassium ion import
11                 positive regulation of potassium ion import
12                               negative regulation of growth
13                              regulation of cellular process
14                            regulation of biological process
15                          regulation of sodium ion transport
16                                        sodium ion transport
17                                       segment specification
18 regulation of sodium ion transmembrane transporter activity
19                          sodium ion transmembrane transport
20                                              system process

head(goterms_cc, 20)

           ID      Pvalue                                            Terms
1  GO:0005886 0.001486595                                  plasma membrane
2  GO:0005890 0.001518698       sodium:potassium-exchanging ATPase complex
3  GO:0016323 0.002148167                      basolateral plasma membrane
4  GO:0071944 0.002224033                                   cell periphery
5  GO:0044459 0.003264812                             plasma membrane part
6  GO:0090533 0.003572430               cation-transporting ATPase complex
7  GO:0031226 0.003674196           intrinsic component of plasma membrane
8  GO:0016324 0.004219984                           apical plasma membrane
9  GO:0045177 0.005759799                              apical part of cell
10 GO:0098533 0.007536376 ATPase dependent transmembrane transport complex
11 GO:0005887 0.007724379            integral component of plasma membrane
12 GO:0098590 0.008426372                           plasma membrane region
13 GO:0032591 0.010249974                         dendritic spine membrane
14 GO:0036477 0.017893454                      somatodendritic compartment
15 GO:0005584 0.020395947                           collagen type I trimer
16 GO:0005967 0.020395947     mitochondrial pyruvate dehydrogenase complex
17 GO:0008043 0.020395947                   intracellular ferritin complex
18 GO:0033162 0.020395947                              melanosome membrane
19 GO:0045009 0.020395947                                        chitosome
20 GO:0045025 0.020395947                        mitochondrial degradosome

head(goterms_mf, 20)

           ID       Pvalue
1  GO:0004551 0.0006593199
2  GO:0005391 0.0016255779
3  GO:0008556 0.0016255779
4  GO:0001671 0.0068617248
5  GO:0035035 0.0080535610
6  GO:0004081 0.0106095053
7  GO:0004654 0.0106095053
8  GO:0005275 0.0106095053
9  GO:0008504 0.0106095053
10 GO:0008796 0.0106095053
11 GO:0008803 0.0106095053
12 GO:0015222 0.0106095053
13 GO:0035240 0.0106095053
14 GO:0035643 0.0106095053
15 GO:0072544 0.0106095053
16 GO:0060590 0.0136606917
17 GO:0015662 0.0152630549
18 GO:0031072 0.0194406746
19 GO:0003835 0.0211075854
20 GO:0004704 0.0211075854
                                                                                   Terms
1                                                      nucleotide diphosphatase activity
2                                            sodium:potassium-exchanging ATPase activity
3                                                 potassium-transporting ATPase activity
4                                                              ATPase activator activity
5                                                      histone acetyltransferase binding
6                             bis(5'-nucleosyl)-tetraphosphatase (asymmetrical) activity
7                                     polyribonucleotide nucleotidyltransferase activity
8                                               amine transmembrane transporter activity
9                                           monoamine transmembrane transporter activity
10                                           bis(5'-nucleosyl)-tetraphosphatase activity
11                             bis(5'-nucleosyl)-tetraphosphatase (symmetrical) activity
12                                          serotonin transmembrane transporter activity
13                                                                      dopamine binding
14                                                              L-DOPA receptor activity
15                                                                        L-DOPA binding
16                                                             ATPase regulator activity
17 ATPase activity, coupled to transmembrane movement of ions, phosphorylative mechanism
18                                                            heat shock protein binding
19                                 beta-galactoside alpha-2,6-sialyltransferase activity
20                                                    NF-kappaB-inducing kinase activity

Find driving genes of PC3

pc3_genes <- svd_all$u[ ,3]

# find genes outside of 3 mean +/- sd
mean_pc3 <- mean(pc3_genes)
sd_pc3 <- sd(pc3_genes)
ii_extreme <- which( (pc3_genes > mean_pc2 + 3*sd_pc2)|(pc3_genes < mean_pc2 - 3*sd_pc2) )
length(ii_extreme)

[1] 90

library(broman)
par(mfrow = c(1,1))
xy <- cbind(c(1:length(pc3_genes)), pc3_genes)
plot(xy, cex = .4, pch = 16,
     main = paste("PC3 gene loadings \n", 
                  length(ii_extreme), "genes"), xlab = "genes", ylab = "PC3 loadings")
points(xy[ii_extreme, ], pch = 1, cex = .7, col = brocolors("crayon")["Shamrock"])

Annotation

library("biomaRt")
ensembl <- useMart(host = "grch37.ensembl.org",
                   biomart = "ENSEMBL_MART_ENSEMBL",
                   dataset = "hsapiens_gene_ensembl")
pc3_gene_info <- getBM(attributes = c("ensembl_gene_id", "chromosome_name",
                                       "external_gene_name", "transcript_count",
                                       "description"),
                        filters = "ensembl_gene_id",
                        values = rownames(molecules_ENSG[ii_extreme, ]),
                        mart = ensembl)
kable(pc3_gene_info[ , c(1, 2, 3, 5)])

if (file.exists("rda/svd-filtered/pc3-cv.rda")) {
  load("rda/svd-filtered/pc3-cv.rda")  
} else {
  library(Humanzee)
  go_pc3 <- GOtest(my_ensembl_gene_universe = rownames(molecules_ENSG),
                my_ensembl_gene_test = rownames(molecules_ENSG)[ii_extreme],
                pval_cutoff = 1, ontology=c("BP","CC","MF") )
  save(go_pc3, file = "rda/svd-filtered/pc2-cv.rda")
}

Extract terms

if (file.exists("rda/svd-filtered/go-pc3-terms.rda")) {
  load("rda/svd-filtered/go-pc3-terms.rda")  
} else {
  # Biological process
  goterms_bp <- summary(go_pc3$GO$BP, pvalue = 1)
  goterms_bp <- data.frame(ID = goterms_bp[[1]],
                           Pvalue = goterms_bp[[2]],
                           Terms = goterms_bp[[7]])
  goterms_bp <- goterms_bp[order(goterms_bp$Pvalue), ]
  
  # Cellular component
  goterms_cc <- summary(go_pc3$GO$CC, pvalue = 1)
  goterms_cc <- data.frame(ID = goterms_cc[[1]],
                           Pvalue = goterms_cc[[2]],
                           Terms = goterms_cc[[7]])
  goterms_cc <- goterms_cc[order(goterms_cc$Pvalue), ]

  # Molecular function
  goterms_mf <- summary(go_pc3$GO$MF, pvalue = 1)
  goterms_mf <- data.frame(ID = goterms_mf[[1]],
                           Pvalue = goterms_mf[[2]],
                           Terms = goterms_mf[[7]])
  goterms_mf <- goterms_mf[order(goterms_mf$Pvalue), ]

  save(goterms_bp, goterms_cc, goterms_mf, 
       file = "rda/svd-filtered/go-pc3-terms.rda")
}

head(goterms_bp, 20)

           ID       Pvalue
1  GO:0071294 3.878761e-07
2  GO:0071276 1.148025e-06
3  GO:0071241 1.117054e-05
4  GO:0010043 1.994167e-05
5  GO:0071248 7.431671e-05
6  GO:0046686 9.915044e-05
7  GO:1901685 4.151958e-04
8  GO:1901687 4.151958e-04
9  GO:0031589 4.247346e-04
10 GO:0036149 7.602917e-04
11 GO:0006936 9.161115e-04
12 GO:0002544 1.133887e-03
13 GO:0045989 1.133887e-03
14 GO:1990267 1.286417e-03
15 GO:0006941 1.474185e-03
16 GO:0007263 1.578329e-03
17 GO:0010038 1.624272e-03
18 GO:0072507 2.141512e-03
19 GO:0048469 2.281605e-03
20 GO:0033555 2.281980e-03
                                                Terms
1                       cellular response to zinc ion
2                    cellular response to cadmium ion
3            cellular response to inorganic substance
4                                response to zinc ion
5                      cellular response to metal ion
6                             response to cadmium ion
7            glutathione derivative metabolic process
8         glutathione derivative biosynthetic process
9                             cell-substrate adhesion
10         phosphatidylinositol acyl-chain remodeling
11                                 muscle contraction
12                      chronic inflammatory response
13 positive regulation of striated muscle contraction
14          response to transition metal nanoparticle
15                        striated muscle contraction
16          nitric oxide mediated signal transduction
17                              response to metal ion
18              divalent inorganic cation homeostasis
19                                    cell maturation
20        multicellular organismal response to stress

head(goterms_cc, 20)

           ID      Pvalue                              Terms
1  GO:0005589 0.008785692            collagen type VI trimer
2  GO:0030485 0.008785692    smooth muscle contractile fiber
3  GO:0044228 0.008785692                  host cell surface
4  GO:0098647 0.008785692           collagen beaded filament
5  GO:0031012 0.010730464               extracellular matrix
6  GO:0005577 0.017495106                 fibrinogen complex
7  GO:0097443 0.017495106                   sorting endosome
8  GO:0005578 0.021674083 proteinaceous extracellular matrix
9  GO:0048471 0.022341863    perinuclear region of cytoplasm
10 GO:0000798 0.026128898            nuclear cohesin complex
11 GO:0034991 0.026128898    nuclear meiotic cohesin complex
12 GO:0097513 0.026128898                 myosin II filament
13 GO:0005882 0.031905152              intermediate filament
14 GO:0018995 0.034687715                               host
15 GO:0030893 0.034687715            meiotic cohesin complex
16 GO:0032982 0.034687715                    myosin filament
17 GO:0033643 0.034687715                     host cell part
18 GO:0043245 0.034687715              extraorganismal space
19 GO:0043657 0.034687715                          host cell
20 GO:0044215 0.034687715                     other organism

head(goterms_mf, 20)

           ID       Pvalue
1  GO:0016672 7.222821e-05
2  GO:0045174 7.222821e-05
3  GO:0050610 7.222821e-05
4  GO:0004364 1.619106e-04
5  GO:0019838 1.920784e-04
6  GO:0030613 2.154802e-04
7  GO:0030614 2.154802e-04
8  GO:0016765 4.334233e-04
9  GO:0015038 7.103117e-04
10 GO:0001948 7.841233e-04
11 GO:0015037 1.059556e-03
12 GO:0050840 1.358660e-03
13 GO:0005178 1.680954e-03
14 GO:0050839 1.913923e-03
15 GO:0001786 2.500895e-03
16 GO:0004620 2.731199e-03
17 GO:0016667 3.509092e-03
18 GO:0004623 3.778650e-03
19 GO:0005520 4.509316e-03
20 GO:0001968 5.299753e-03
                                                                                                  Terms
1  oxidoreductase activity, acting on a sulfur group of donors, quinone or similar compound as acceptor
2                                                        glutathione dehydrogenase (ascorbate) activity
3                                                                     methylarsonate reductase activity
4                                                                      glutathione transferase activity
5                                                                                 growth factor binding
6                                    oxidoreductase activity, acting on phosphorus or arsenic in donors
7             oxidoreductase activity, acting on phosphorus or arsenic in donors, disulfide as acceptor
8                           transferase activity, transferring alkyl or aryl (other than methyl) groups
9                                                         glutathione disulfide oxidoreductase activity
10                                                                                 glycoprotein binding
11                                                            peptide disulfide oxidoreductase activity
12                                                                         extracellular matrix binding
13                                                                                     integrin binding
14                                                                       cell adhesion molecule binding
15                                                                           phosphatidylserine binding
16                                                                               phospholipase activity
17                                          oxidoreductase activity, acting on a sulfur group of donors
18                                                                            phospholipase A2 activity
19                                                                   insulin-like growth factor binding
20                                                                                  fibronectin binding

Session information

sessionInfo()

R version 3.2.1 (2015-06-18)
Platform: x86_64-apple-darwin13.4.0 (64-bit)
Running under: OS X 10.10.5 (Yosemite)

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    grid      stats     graphics  grDevices utils    
 [8] datasets  methods   base     

other attached packages:
 [1] org.Hs.eg.db_3.1.2   GOstats_2.34.0       graph_1.46.0        
 [4] Category_2.34.2      Matrix_1.2-3         GO.db_3.1.2         
 [7] RSQLite_1.0.0        DBI_0.3.1            AnnotationDbi_1.30.1
[10] GenomeInfoDb_1.4.3   IRanges_2.2.9        S4Vectors_0.6.6     
[13] Biobase_2.28.0       BiocGenerics_0.14.0  Humanzee_0.1.0      
[16] biomaRt_2.24.1       broman_0.59-5        gridExtra_2.0.0     
[19] matrixStats_0.15.0   ggplot2_1.0.1        edgeR_3.10.5        
[22] limma_3.24.15        dplyr_0.4.3          data.table_1.9.6    
[25] knitr_1.11          

loaded via a namespace (and not attached):
 [1] genefilter_1.50.0      reshape2_1.4.1         splines_3.2.1         
 [4] lattice_0.20-33        colorspace_1.2-6       htmltools_0.2.6       
 [7] yaml_2.1.13            survival_2.38-3        chron_2.3-47          
[10] XML_3.98-1.3           RBGL_1.44.0            plyr_1.8.3            
[13] stringr_1.0.0          munsell_0.4.2          gtable_0.1.2          
[16] evaluate_0.8           labeling_0.3           highr_0.5.1           
[19] GSEABase_1.30.2        proto_0.3-10           Rcpp_0.12.2           
[22] xtable_1.8-0           scales_0.3.0           formatR_1.2.1         
[25] annotate_1.46.1        digest_0.6.8           stringi_1.0-1         
[28] tools_3.2.1            bitops_1.0-6           magrittr_1.5          
[31] RCurl_1.95-4.7         MASS_7.3-45            assertthat_0.1        
[34] rmarkdown_0.8.1        AnnotationForge_1.10.1 R6_2.1.1