PAGODA2 walkthrough

This is a quick walkthrough covering key PAGODA2 usage patterns to date.

We will analyze an 8k PBMC dataset from 10x as an example (starting with unfiltered read count matrices).

Try executing this chunk by clicking the Run button within the chunk or by placing your cursor inside it and pressing Cmd+Shift+Enter.

First we’ll load the library

library(pagoda2)

Now let’s load the 10x matrices:

cd <- read.10x.matrices(list(PBMC8K='pbmc8k/raw_gene_bc_matrices/GRCh38'))[[1]]

reading 1 dataset(s) . done

str(cd)

Formal class 'dgCMatrix' [package "Matrix"] with 6 slots
  ..@ i       : int [1:17171166] 1614 21625 110 153 231 484 548 550 633 673 ...
  ..@ p       : int [1:737281] 0 2 2 2 116 120 120 125 126 126 ...
  ..@ Dim     : int [1:2] 33694 737280
  ..@ Dimnames:List of 2
  .. ..$ : chr [1:33694] "RP11-34P13.3" "FAM138A" "OR4F5" "RP11-34P13.7" ...
  .. ..$ : chr [1:737280] "PBMC8K_AAACCTGAGAAACCAT-1" "PBMC8K_AAACCTGAGAAACCGC-1" "PBMC8K_AAACCTGAGAAACCTA-1" "PBMC8K_AAACCTGAGAAACGAG-1" ...
  ..@ x       : num [1:17171166] 1 1 1 1 1 2 1 1 1 1 ...
  ..@ factors : list()

Look at the summary counts across genes and across cells:

par(mfrow=c(1,2), mar = c(3.5,3.5,2.0,0.5), mgp = c(2,0.65,0), cex = 1.0)
hist(log10(colSums(cd)+1),main='molecules per cell',col='wheat',xlab='log10[ molecules per cell]')
hist(log10(rowSums(cd)+1),main='molecules per gene',col='wheat',xlab='log10[ molecules per gene]')

This is unfiltered data, so most of the barcodes are empty/background. Let’s require at least 500 molecules per cell, and use our quick umi-gene relationship function to filter the cells.

counts <- gene.vs.molecule.cell.filter(cd,min.cell.size=500)

str(counts)

Formal class 'dgCMatrix' [package "Matrix"] with 6 slots
  ..@ i       : int [1:11965872] 33 45 59 153 335 375 408 440 484 495 ...
  ..@ p       : int [1:8787] 0 871 1677 2993 3891 5417 6912 8165 9598 11230 ...
  ..@ Dim     : int [1:2] 33694 8786
  ..@ Dimnames:List of 2
  .. ..$ : chr [1:33694] "RP11-34P13.3" "FAM138A" "OR4F5" "RP11-34P13.7" ...
  .. ..$ : chr [1:8786] "PBMC8K_AAACCTGAGCATCATC-1" "PBMC8K_AAACCTGAGCTAACTC-1" "PBMC8K_AAACCTGAGCTAGTGG-1" "PBMC8K_AAACCTGCACATTAGC-1" ...
  ..@ x       : num [1:11965872] 1 1 1 5 1 1 1 1 18 1 ...
  ..@ factors : list()

We can also looka at the number of molecules per gene, and omit low-expressed genes to save computational time:

hist(log10(rowSums(counts)+1),main='Molecules per gene',xlab='molecules (log10)',col='wheat')
abline(v=1,lty=2,col=2)

counts <- counts[rowSums(counts)>=10,]
str(counts)

Formal class 'dgCMatrix' [package "Matrix"] with 6 slots
  ..@ i       : int [1:11947057] 11 17 30 84 174 192 207 222 247 257 ...
  ..@ p       : int [1:8787] 0 869 1675 2991 3889 5412 6907 8158 9590 11219 ...
  ..@ Dim     : int [1:2] 15556 8786
  ..@ Dimnames:List of 2
  .. ..$ : chr [1:15556] "RP11-34P13.7" "FO538757.2" "AP006222.2" "RP4-669L17.10" ...
  .. ..$ : chr [1:8786] "PBMC8K_AAACCTGAGCATCATC-1" "PBMC8K_AAACCTGAGCTAACTC-1" "PBMC8K_AAACCTGAGCTAGTGG-1" "PBMC8K_AAACCTGCACATTAGC-1" ...
  ..@ x       : num [1:11947057] 1 1 1 5 1 1 1 1 18 1 ...
  ..@ factors : list()

Now we have a clean/lean count matrix and are ready to start analysis. First we’ll create pagoda2 object that will maintain all of the results. It will also provide handles for running all operations on the data. Use of log scale is recommended, as it makes some statistics smoother.

r <- Pagoda2$new(counts,log.scale=TRUE)

Error in setCountMatrix(x, min.cells.per.gene = min.cells.per.gene, trim = trim,  : 
  duplicate gene names are not allowed - please reduce

Yes, we need to make gene names unique:

rownames(counts) <- make.unique(rownames(counts))
r <- Pagoda2$new(counts,log.scale=TRUE, n.cores=20)

8786 cells, 13004 genes; normalizing ... using plain model winsorizing ... log scale ... done.

Note, that above we’ve also specified the default number of processing cores that will be used to speed up some of the analysis steps below. Please match it to your system.

Next, we’ll adjust the variance, to normalize the extent to which genes with (very) different expression magnitudes will contribute to the downstream anlaysis:

r$adjustVariance(plot=T,gam.k=10)

calculating variance fit ... using gam 

Loading required package: mgcv
Loading required package: nlme
This is mgcv 1.8-18. For overview type 'help("mgcv-package")'.

538 overdispersed genes ... 538 persisting ... done.

There are many alternative ways of proceeding with the downstream analysis. Below we’ll use the simplest, default scenario, where we first reduce the dataset dimensions by running PCA, and then move into k-nearest neighbor graph space for clustering and visualization calculations. First, the PCA reduction:

r$calculatePcaReduction(nPcs=50,n.odgenes=3e3)

Loading required package: irlba

The next few steps will make kNN graph, find clusters and generate a quick largeVis embedding to visualize the subpopulations:

r$makeKnnGraph(k=40,type='PCA',center=T,distance='cosine');

Loading required package: igraph

Attaching package: ‘igraph’

The following objects are masked from ‘package:stats’:

    decompose, spectrum

The following object is masked from ‘package:base’:

    union

creating space of type angular done
adding data ... done
building index ... done
querying ... done

r$getKnnClusters(method=infomap.community,type='PCA')
M <- 30; r$getEmbedding(type='PCA',M=M,perplexity=30,gamma=1/M,alpha=1)

Estimating embeddings.

0%   10   20   30   40   50   60   70   80   90   100%
|----|----|----|----|----|----|----|----|----|----|
**************************************************|

Now we can visualize the embedding using the determined clusters:

r$plotEmbedding(type='PCA',show.legend=F,mark.clusters=T,min.group.size=50,shuffle.colors=F,mark.cluster.cex=1,alpha=0.1,main='clusters (largeVis)')

Alternatively, we can generate tSNE embedding (takes longer)

r$getEmbedding(type='PCA',embeddingType='tSNE',perplexity=50,verbose=F)

Loading required package: Rtsne

calculating distance ... pearson ...running tSNE using 20 cores:

r$plotEmbedding(type='PCA',embeddingType='tSNE',show.legend=F,mark.clusters=T,min.group.size=1,shuffle.colors=F,mark.cluster.cex=1,alpha=0.1,main='clusters (tSNE)')

We can use the same plotEmbedding() function to show all kinds of other values. For instance, let’s look at depth, or an expresson pattern of a gene:

str(r$depth)

 Named num [1:8786] 2376 1660 4517 2779 4656 ...
 - attr(*, "names")= chr [1:8786] "PBMC8K_AAACCTGAGCATCATC-1" "PBMC8K_AAACCTGAGCTAACTC-1" "PBMC8K_AAACCTGAGCTAGTGG-1" "PBMC8K_AAACCTGCACATTAGC-1" ...

par(mfrow=c(1,2))
r$plotEmbedding(type='PCA',embeddingType='tSNE',colors=r$depth,shuffle.colors=F,mark.cluster.cex=1,alpha=0.1,main='depth')

treating colors as a gradient with zlim: 1440 8903.75 

gene <-"LYZ"
r$plotEmbedding(type='PCA',embeddingType='tSNE',colors=r$counts[,gene],shuffle.colors=F,mark.cluster.cex=1,alpha=0.1,main=gene)

treating colors as a gradient with zlim: 0 2.879075 

We can generate multiple potential clusterings, with different names. Here we’ll use multilevel clustering:

r$getKnnClusters(method=multilevel.community,type='PCA',name='multilevel')
str(r$clusters)

List of 1
 $ PCA:List of 2
  ..$ community : Factor w/ 36 levels "1","2","3","4",..: 27 2 5 3 6 2 14 3 28 3 ...
  .. ..- attr(*, "names")= chr [1:8786] "PBMC8K_AAACCTGAGCATCATC-1" "PBMC8K_AAACCTGAGCTAACTC-1" "PBMC8K_AAACCTGAGCTAGTGG-1" "PBMC8K_AAACCTGCACATTAGC-1" ...
  ..$ multilevel: Factor w/ 16 levels "1","2","3","4",..: 6 9 4 1 9 9 15 1 9 1 ...
  .. ..- attr(*, "names")= chr [1:8786] "PBMC8K_AAACCTGAGCATCATC-1" "PBMC8K_AAACCTGAGCTAACTC-1" "PBMC8K_AAACCTGAGCTAGTGG-1" "PBMC8K_AAACCTGCACATTAGC-1" ...

Compare with infomap:

par(mfrow=c(1,2))
r$plotEmbedding(type='PCA',embeddingType='tSNE',groups=r$clusters$PCA$community,show.legend=F,mark.clusters=T,min.group.size=1,shuffle.colors=F,mark.cluster.cex=1,alpha=0.1,main='infomap clusters (tSNE)')

using provided groups as a factor

r$plotEmbedding(type='PCA',embeddingType='tSNE',clusterType='multilevel',show.legend=F,mark.clusters=T,min.group.size=1,shuffle.colors=F,mark.cluster.cex=1,alpha=0.1,main='multlevel clusters (tSNE)')

Run differential expression on the infomap clusters:

r$getDifferentialGenes(type='PCA',verbose=T,clusterType='community')

running differential expression with  36  clusters ... adjusting p-values ... done.

Visualize top genes:

names(r$diffgenes)

[1] "PCA"

de <- r$diffgenes$PCA[[1]][['4']];
r$plotGeneHeatmap(genes=rownames(de)[1:15],groups=r$clusters$PCA[[1]])

Spot-check a gene

gene <-"IGHM"
r$plotEmbedding(type='PCA',embeddingType='tSNE',colors=r$counts[,gene],shuffle.colors=F,mark.cluster.cex=1,alpha=0.1,main=gene)

treating colors as a gradient with zlim: 0 1.078653 

Pathway overdispersion analysis (a la PAGODA1)

First, build GO->gene environment:

suppressMessages(library(org.Hs.eg.db))
# translate gene names to ids
ids <- unlist(lapply(mget(colnames(r$counts),org.Hs.egALIAS2EG,ifnotfound=NA),function(x) x[1]))
# reverse map
rids <- names(ids); names(rids) <- ids;
# list all the ids per GO category
go.env <- list2env(eapply(org.Hs.egGO2ALLEGS,function(x) as.character(na.omit(rids[x]))))

Now run overdispersion anlaysis

r$testPathwayOverdispersion(go.env,verbose=T,correlation.distance.threshold=0.95,recalculate.pca=F,top.aspects=15)

We’ll use hierarchical differential expression results instead:

r$getHierarchicalDiffExpressionAspects(type='PCA',clusterName='community',z.threshold=3)

using community  clustering for PCA space

1 function calls resulted in an error

We’ll make an app with that, ordering the “differential expression aspects” explicitly (otherwise if row clustering is omitted they’ll be clustered by similarity)

app <- p2.make.pagoda1.app(r,inner.clustering=TRUE,embeddingType='tSNE',clusterType='multilevel',min.group.size=50,row.clustering=list(order=rev(1:nrow(r$misc$pathwayOD$xv))))

Loading required package: fastcluster

Attaching package: ‘fastcluster’

The following object is masked from ‘package:stats’:

    hclust

Show app:

show.app(app,'pbmc',browse=T)

Alternatively, we can create a more advanced p2 app, which can be saved as a standalone binary file and viewed from the browser either by accessing that binary file on the local drive, or over regular HTTP (in other words, without requiring a server or an R session).

First, compile gene sets and metadata that we want to make available in the offline web app:

library(GO.db)
termDescriptions <- Term(GOTERM[names(go.env)]); # saves a good minute or so compared to individual lookups
sn <- function(x) { names(x) <- x; x}
geneSets <- lapply(sn(names(go.env)),function(x) {
  list(properties=list(locked=T,genesetname=x,shortdescription=as.character(termDescriptions[x])),genes=c(go.env[[x]]))
})

# Make metadata
additionalMetadata <- list(multilevel = p2.metadata.from.factor(r$clusters$PCA$multilevel, displayname = 'Multilevel', s = 0.8))

Now generate the web app and the binary output.

# Generate the web application
# dendrogramCellGroups specifies the clustering that will be used to generate the main dendrogram
p2app <- make.p2.app(r, dendrogramCellGroups = r$clusters$PCA$community, additionalMetadata = additionalMetadata, geneSets = geneSets);
# Optional showing of app
#show.app(p2w, name='newPagoda')
# Save serialised web object, RDS app and session image
p2app$serializeToStaticFast(binary.filename = 'pbmc8k.bin',verbose=T)

Writing... matsparse
        p array size: 13005 [First entry value: 0]
        i array size: 11901901 [First entry value: 3]
        x array size: 11901901 [First entry value: 0.86584]
    Sparse matrix header information
        dim1=8786
        dim2=13004
        pStartOffset=32
        iStartOffset=52052
        xStartOffset=47659656
        dimnames1StartOffset=95267260
        dimnames2StartOffset=95513270
        dimnames2EndOffset=95630240
Writing... mataspect
        p array size: 34 [First entry value: 0]
        i array size: 52132 [First entry value: 0]
        x array size: 52132 [First entry value: -19.18]
    Sparse matrix header information
        dim1=8786
        dim2=33
        pStartOffset=32
        iStartOffset=168
        xStartOffset=208696
        dimnames1StartOffset=417224
        dimnames2StartOffset=663234
        dimnames2EndOffset=663923
Writing... sparseMatrixTransp
        p array size: 8787 [First entry value: 0]
        i array size: 11901901 [First entry value: 8]
        x array size: 11901901 [First entry value: 0.351273]
    Sparse matrix header information
        dim1=13004
        dim2=8786
        pStartOffset=32
        iStartOffset=35180
        xStartOffset=47642784
        dimnames1StartOffset=95250388
        dimnames2StartOffset=95367358
        dimnames2EndOffset=95613368
Making File from payload...
    File format information
        Index entry size is 140 bytes
        File header size is 48 bytes
    Preparing header...
    Total index size is: 1960 bytes
    Constructing index...
        Writing entry 0 of size 1 blocks (or 32768 bytes)
        Writing entry 1 of size 25 blocks (or 819200 bytes)
        Writing entry 2 of size 8 blocks (or 262144 bytes)
        Writing entry 3 of size 30 blocks (or 983040 bytes)
        Writing entry 4 of size 1 blocks (or 32768 bytes)
        Writing entry 5 of size 1 blocks (or 32768 bytes)
        Writing entry 6 of size 1 blocks (or 32768 bytes)
        Writing entry 7 of size 67 blocks (or 2195456 bytes)
        Writing entry 8 of size 536 blocks (or 17563648 bytes)
        Writing entry 9 of size 14 blocks (or 458752 bytes)
        Writing entry 10 of size 14 blocks (or 458752 bytes)
        Writing entry 11 of size 2919 blocks (or 95649792 bytes)
        Writing entry 12 of size 21 blocks (or 688128 bytes)
        Writing entry 13 of size 2918 blocks (or 95617024 bytes)
Free up entry payloads
NULL

Now the pbmc8k.bin file can be viewed locally using http://pklab.med.harvard.edu/nikolas/pagoda2/frontend/current/pagodaLocal/index.html (which will ask you to select a file), or by placing the .bin file into some HTTP-accessible location and using a URL like this, where ‘fileURL’ argument would specify the location at which the .bin file can be accessed: http://pklab.med.harvard.edu/nikolas/pagoda2/frontend/current/pagodaURL/index.html?fileURL=http://pklab.med.harvard.edu/nikolas/pagoda2/staticDemo/pagodaPublicDemo.bin

The index.html file (and associated scripts) are avaialble in the “inst” folder of the pagoda2 package, in case you want to set up your own copy.

Here’s a view of the binary file that was produced by this tutorial: http://pklab.med.harvard.edu/peterk/p2/r/index.html?fileURL=http://pklab.med.harvard.edu/peterk/p2/pbmc8k.bin