Last updated: 2019-06-12

Checks: 6 0

Knit directory: Porello-heart-snRNAseq/

This reproducible R Markdown analysis was created with workflowr (version 1.3.0). The Checks tab describes the reproducibility checks that were applied when the results were created. The Past versions tab lists the development history.

R Markdown file: up-to-date

Great! Since the R Markdown file has been committed to the Git repository, you know the exact version of the code that produced these results.

Environment: empty

Great job! The global environment was empty. Objects defined in the global environment can affect the analysis in your R Markdown file in unknown ways. For reproduciblity it’s best to always run the code in an empty environment.

Seed: set.seed(20190603)

The command set.seed(20190603) was run prior to running the code in the R Markdown file. Setting a seed ensures that any results that rely on randomness, e.g. subsampling or permutations, are reproducible.

Session information: recorded

Great job! Recording the operating system, R version, and package versions is critical for reproducibility.

Cache: none

Nice! There were no cached chunks for this analysis, so you can be confident that you successfully produced the results during this run.

Repository version: 0b36828

Great! You are using Git for version control. Tracking code development and connecting the code version to the results is critical for reproducibility. The version displayed above was the version of the Git repository at the time these results were generated.

Note that you need to be careful to ensure that all relevant files for the analysis have been committed to Git prior to generating the results (you can use wflow_publish or wflow_git_commit). workflowr only checks the R Markdown file, but you know if there are other scripts or data files that it depends on. Below is the status of the Git repository when the results were generated: 


Ignored files:
    Ignored:    .Rhistory
    Ignored:    .Rproj.user/

Note that any generated files, e.g. HTML, png, CSS, etc., are not included in this status report because it is ok for generated content to have uncommitted changes.

These are the previous versions of the R Markdown and HTML files. If you’ve configured a remote Git repository (see ?wflow_git_remote), click on the hyperlinks in the table below to view them.

File Version Author Date Message
Rmd 0b36828 Belinda Phipson 2019-06-12 update QC
html 6f2bf8d Belinda Phipson 2019-06-09 Build site.
Rmd 225f001 Belinda Phipson 2019-06-09 Add MDSplots to QC
html 5839497 Belinda Phipson 2019-06-07 Build site.
Rmd 565d32b Belinda Phipson 2019-06-07 Adding cluster of fetal samples
Rmd 74d2e89 Belinda Phipson 2019-06-07 added new plots to QC report
html 2b103a6 Belinda Phipson 2019-06-06 Build site.
Rmd ce619de Belinda Phipson 2019-06-06 updated QC analysis and added fetal clustering ananalysis (incomplete)
Rmd 5a0d4e1 Belinda Phipson 2019-06-05 Added QC analysis and plots

Introduction

The first step of the analysis is to perform quality control of the cells to make sure that low quality cells are removed prior to further analysis.

Load libraries and functions

library(edgeR)
Loading required package: limma
library(RColorBrewer)
library(org.Hs.eg.db)
Loading required package: AnnotationDbi
Loading required package: stats4
Loading required package: BiocGenerics
Loading required package: parallel

Attaching package: 'BiocGenerics'
The following objects are masked from 'package:parallel':

    clusterApply, clusterApplyLB, clusterCall, clusterEvalQ,
    clusterExport, clusterMap, parApply, parCapply, parLapply,
    parLapplyLB, parRapply, parSapply, parSapplyLB
The following object is masked from 'package:limma':

    plotMA
The following objects are masked from 'package:stats':

    IQR, mad, sd, var, xtabs
The following objects are masked from 'package:base':

    anyDuplicated, append, as.data.frame, basename, cbind,
    colnames, dirname, do.call, duplicated, eval, evalq, Filter,
    Find, get, grep, grepl, intersect, is.unsorted, lapply, Map,
    mapply, match, mget, order, paste, pmax, pmax.int, pmin,
    pmin.int, Position, rank, rbind, Reduce, rownames, sapply,
    setdiff, sort, table, tapply, union, unique, unsplit, which,
    which.max, which.min
Loading required package: Biobase
Welcome to Bioconductor

    Vignettes contain introductory material; view with
    'browseVignettes()'. To cite Bioconductor, see
    'citation("Biobase")', and for packages 'citation("pkgname")'.
Loading required package: IRanges
Loading required package: S4Vectors

Attaching package: 'S4Vectors'
The following object is masked from 'package:base':

    expand.grid
library(limma)
library(Seurat)
Registered S3 method overwritten by 'R.oo':
  method        from       
  throw.default R.methodsS3
Registered S3 methods overwritten by 'ggplot2':
  method         from 
  [.quosures     rlang
  c.quosures     rlang
  print.quosures rlang
library(monocle)
Loading required package: Matrix

Attaching package: 'Matrix'
The following object is masked from 'package:S4Vectors':

    expand
Loading required package: ggplot2
Loading required package: VGAM
Loading required package: splines
Loading required package: DDRTree
Loading required package: irlba
library(cowplot)

Attaching package: 'cowplot'
The following object is masked from 'package:ggplot2':

    ggsave
library(DelayedArray)
Loading required package: matrixStats

Attaching package: 'matrixStats'
The following objects are masked from 'package:Biobase':

    anyMissing, rowMedians
Loading required package: BiocParallel

Attaching package: 'DelayedArray'
The following objects are masked from 'package:matrixStats':

    colMaxs, colMins, colRanges, rowMaxs, rowMins, rowRanges
The following objects are masked from 'package:base':

    aperm, apply, rowsum
library(scran)
Loading required package: SingleCellExperiment
Loading required package: SummarizedExperiment
Loading required package: GenomicRanges
Loading required package: GenomeInfoDb

Attaching package: 'SingleCellExperiment'
The following object is masked from 'package:edgeR':

    cpm
library(NMF)
Loading required package: pkgmaker
Loading required package: registry

Attaching package: 'pkgmaker'
The following object is masked from 'package:S4Vectors':

    new2
The following object is masked from 'package:base':

    isFALSE
Loading required package: rngtools
Loading required package: cluster
NMF - BioConductor layer [OK] | Shared memory capabilities [NO: synchronicity] | Cores 31/32
  To enable shared memory capabilities, try: install.extras('
NMF
')

Attaching package: 'NMF'
The following object is masked from 'package:DelayedArray':

    seed
The following object is masked from 'package:BiocParallel':

    register
The following object is masked from 'package:S4Vectors':

    nrun
library(workflowr)
This is workflowr version 1.3.0
Run ?workflowr for help getting started
source("/group/bioi1/belinda/SingleCell/Normalisation/normCounts.R")
source("/misc/card2-single_cell_nuclei_rnaseq/Porello-heart-snRNAseq/code/findModes.R")
source("/misc/card2-single_cell_nuclei_rnaseq/Porello-heart-snRNAseq/code/ggplotColors.R")

Read in data

I have stored the sample-level information in a targets file under the data directory. The cellranger output files are read in using the Seurat function Read10X.

targets <- read.delim("/misc/card2-single_cell_nuclei_rnaseq/Porello-heart-snRNAseq/data/targets.txt",header=TRUE, stringsAsFactors = FALSE)
targets$FileName2 <- paste(targets$FileName,"/",sep="")
targets$Group_ID2 <- gsub("LV_","",targets$Group_ID)
group <- c("Fetal_1","Fetal_2","Fetal_3",
           "Young_1","Young_2","Young_3",
           "Adult_1","Adult_2","Adult_3", 
           "Diseased_1","Diseased_2",
           "Diseased_3","Diseased_4")
m <- match(group, targets$Group_ID2)
targets <- targets[m,]


f1 <- Read10X(data.dir = targets$FileName2[targets$Group_ID2=="Fetal_1"])
f2 <- Read10X(data.dir = targets$FileName2[targets$Group_ID2=="Fetal_2"])
f3 <- Read10X(data.dir = targets$FileName2[targets$Group_ID2=="Fetal_3"])

y1 <- Read10X(data.dir = targets$FileName2[targets$Group_ID2=="Young_1"])
y2 <- Read10X(data.dir = targets$FileName2[targets$Group_ID2=="Young_2"])
y3 <- Read10X(data.dir = targets$FileName2[targets$Group_ID2=="Young_3"])

a1 <- Read10X(data.dir = targets$FileName2[targets$Group_ID2=="Adult_1"])
a2 <- Read10X(data.dir = targets$FileName2[targets$Group_ID2=="Adult_2"])
a3 <- Read10X(data.dir = targets$FileName2[targets$Group_ID2=="Adult_3"])

d1 <- Read10X(data.dir = targets$FileName2[targets$Group_ID2=="Diseased_1"])
d2 <- Read10X(data.dir = targets$FileName2[targets$Group_ID2=="Diseased_2"])
d3 <- Read10X(data.dir = targets$FileName2[targets$Group_ID2=="Diseased_3"])
d4 <- Read10X(data.dir = targets$FileName2[targets$Group_ID2=="Diseased_4"])

Basic QC

I will QC each sample individually. I examined the distribution of library sizes per cell, as well as the proportions of zeroes in each cell and the number of genes detected per cell. I also looked at the ribosomal and mitochondrial content per cell. The ribosomal and mitochondrial content is very low compared to single cell protocols and mostly the number of genes detected is at least 500 per cell. Considering that cellranger has already performed cell filtering using their default setting, I’m keeping all cells for now and will proceed with clustering.

Get gene annotation

I’m using the org.Hs.eg.db annotation package to extract gene symbols, entrez ID, ensembl ID and the genename. This is useful for identifying ribosomal and mitochondrial related genes for QC purposes, as well as for downstream analysis.

columns(org.Hs.eg.db)
 [1] "ACCNUM"       "ALIAS"        "ENSEMBL"      "ENSEMBLPROT" 
 [5] "ENSEMBLTRANS" "ENTREZID"     "ENZYME"       "EVIDENCE"    
 [9] "EVIDENCEALL"  "GENENAME"     "GO"           "GOALL"       
[13] "IPI"          "MAP"          "OMIM"         "ONTOLOGY"    
[17] "ONTOLOGYALL"  "PATH"         "PFAM"         "PMID"        
[21] "PROSITE"      "REFSEQ"       "SYMBOL"       "UCSCKG"      
[25] "UNIGENE"      "UNIPROT"     
ann <- select(org.Hs.eg.db,keys=rownames(f1),columns=c("SYMBOL","ENTREZID","ENSEMBL","GENENAME"),keytype = "SYMBOL")
'select()' returned 1:many mapping between keys and columns
m <- match(rownames(f1),ann$SYMBOL)
ann <- ann[m,]
table(ann$SYMBOL==rownames(f1))

 TRUE 
33939 
mito <- grep("mitochondrial",ann$GENENAME)
length(mito)
[1] 226
ribo <- grep("ribosomal",ann$GENENAME)
length(ribo)
[1] 198

Fetal 1

# Fetal 1
par(mfrow=c(1,2))
libsize <- colSums(f1)
pz <- colMeans(f1==0)
numgene <- colSums(f1!=0)

find_modes(libsize)
 [1]   748.4111  9962.9469 24040.7100 27368.1812 29287.8762 36070.7984
 [7] 40166.1477 44773.4156 46437.1512 53604.0124 57827.3414 62818.5483
find_modes(numgene)
[1]  444.6676 3607.0781 5604.3900 6048.2371 7379.7784 8563.3707 8914.7496
find_modes(pz)
[1] 0.7373302 0.7476835 0.7825576 0.8217909 0.8348687 0.8937188 0.9868980
smoothScatter(libsize,pz,xlab="Library size",ylab="Proportion zeroes per cell",main="Fetal 1")
smoothScatter(libsize,numgene,xlab="Library size",ylab="Detected genes",main="Fetal 1")

Version Author Date
2b103a6 Belinda Phipson 2019-06-06 
par(mfrow=c(1,3))
plot(density(libsize),main="Distribution of library sizes")
#abline(v=3000,col=2)
plot(density(pz), main="Distribution of proportion of zeroes per cell")
#abline(v=0.95,col=2)
plot(density(numgene), main="Distribution of detected genes per cell")
abline(v=500,col=4, lty=3)
legend("topright",lty=2,col=4,legend="500 genes")

Version Author Date
2b103a6 Belinda Phipson 2019-06-06 
mycol <- rep(1,length(libsize))
mycol[pz>0.95] <- 2
par(mfrow=c(1,2))
propmito <- colSums(f1[mito,])/libsize
propribo <- colSums(f1[ribo,])/libsize
plot(libsize,propribo,xlab="Library size",ylab="Proportion ribosomal counts",main="Fetal 1: Ribosomal content",col=mycol)
legend("topright",fill=c(1,2),legend=c("pz <= 0.95","pz > 0.95"))
plot(libsize,propmito,xlab="Library size",ylab="Proportion mitochondrial counts",main="Fetal 1: Mitochondrial content",col=mycol)
legend("topright",fill=c(1,2),legend=c("pz <= 0.95","pz > 0.95"))

Version Author Date
2b103a6 Belinda Phipson 2019-06-06

Fetal 2

# Fetal 2
par(mfrow=c(1,2))
libsize <- colSums(f2)
pz <- colMeans(f2==0)
numgene <- colSums(f2!=0)

find_modes(libsize)
 [1]   914.8093  8566.1059 19217.9110 28069.4110 31970.0720 34070.4280
 [7] 35870.7331 38721.2161 41421.6737 44722.2331 49222.9958 52673.5805
[13] 58224.5212 64825.6398 65275.7161 65725.7924 66475.9195 67376.0720
[19] 67676.1229 68126.1992 69176.3771 69926.5042 70526.6059 74427.2670
find_modes(numgene)
[1]  468.3299 3223.0165 4543.7567 6826.7504 7449.3850 7656.9299 8109.7551
[8] 9260.6858
find_modes(pz)
[1] 0.7271373 0.7610491 0.7743914 0.7805066 0.7988523 0.8661199 0.9050350
[8] 0.9862008
smoothScatter(libsize,pz,xlab="Library size",ylab="Proportion zeroes per cell",main="Fetal 2")
smoothScatter(libsize,numgene,xlab="Library size",ylab="Detected genes",main="Fetal 2")

Version Author Date
2b103a6 Belinda Phipson 2019-06-06 
par(mfrow=c(1,3))
plot(density(libsize),main="Distribution of library sizes")
#abline(v=3000,col=2)
plot(density(pz), main="Distribution of proportion of zeroes per cell")
#abline(v=0.95,col=2)
plot(density(numgene), main="Distribution of detected genes per cell")
abline(v=500,col=4, lty=3)
legend("topright",lty=2,col=4,legend="500 genes")

Version Author Date
2b103a6 Belinda Phipson 2019-06-06 
mycol <- rep(1,length(libsize))
mycol[pz>0.95] <- 2

par(mfrow=c(1,2))
propmito <- colSums(f2[mito,])/libsize
propribo <- colSums(f2[ribo,])/libsize
plot(libsize,propribo,xlab="Library size",ylab="Proportion ribosomal counts",main="Fetal 2: Ribosomal content",col=mycol)
legend("topright",fill=c(1,2),legend=c("pz <= 0.95","pz > 0.95"))
plot(libsize,propmito,xlab="Library size",ylab="Proportion mitochondrial counts",main="Fetal 2: Mitochondrial content",col=mycol)
legend("topright",fill=c(1,2),legend=c("pz <= 0.95","pz > 0.95"))

Version Author Date
2b103a6 Belinda Phipson 2019-06-06

Fetal 3

# Fetal 3
par(mfrow=c(1,2))
libsize <- colSums(f3)
pz <- colMeans(f3==0)
numgene <- colSums(f3!=0)

find_modes(libsize)
 [1]   855.0211 13126.2327 35437.5266 39527.9305 47336.8833 49568.0127
 [7] 52170.9970 55889.5460 60351.8047 72623.0164 76341.5654 84150.5182
[13] 91029.8338
find_modes(numgene)
[1]   464.8103  3765.3342  6609.1818  8103.7587  8560.4350  9100.1433
[7]  9328.4814 10138.0439
find_modes(pz)
[1] 0.7012863 0.7251398 0.7318677 0.7477700 0.7612258 0.8052629 0.8890558
[8] 0.9863045
smoothScatter(libsize,pz,xlab="Library size",ylab="Proportion zeroes per cell",main="Fetal 3")
smoothScatter(libsize,numgene,xlab="Library size",ylab="Detected genes",main="Fetal 3")

Version Author Date
2b103a6 Belinda Phipson 2019-06-06 
par(mfrow=c(1,3))
plot(density(libsize),main="Distribution of library sizes")
#abline(v=3000,col=2)
plot(density(pz), main="Distribution of proportion of zeroes per cell")
#abline(v=0.95,col=2)
plot(density(numgene), main="Distribution of detected genes per cell")
abline(v=500,col=4, lty=3)
legend("topright",lty=2,col=4,legend="500 genes")

Version Author Date
2b103a6 Belinda Phipson 2019-06-06 
mycol <- rep(1,length(libsize))
mycol[pz>0.95] <- 2

par(mfrow=c(1,2))
propmito <- colSums(f3[mito,])/libsize
propribo <- colSums(f3[ribo,])/libsize
plot(libsize,propribo,xlab="Library size",ylab="Proportion ribosomal counts",main="Fetal 3: Ribosomal content",col=mycol)
legend("topright",fill=c(1,2),legend=c("pz <= 0.95","pz > 0.95"))
plot(libsize,propmito,xlab="Library size",ylab="Proportion mitochondrial counts",main="Fetal 3: Mitochondrial content",col=mycol)
legend("topright",fill=c(1,2),legend=c("pz <= 0.95","pz > 0.95"))

Version Author Date
2b103a6 Belinda Phipson 2019-06-06

Young 1

# Young 1
par(mfrow=c(1,2))
libsize <- colSums(y1)
pz <- colMeans(y1==0)
numgene <- colSums(y1!=0)

find_modes(libsize)
 [1]  1128.536  8307.830 22405.353 25016.005 34283.821 38069.267 39896.723
 [8] 43943.234 48250.811 54385.844 55299.572 56604.898 57649.159 63523.127
find_modes(numgene)
[1]  511.8448 1043.3141 3302.0583 5494.3689 6391.2232 6739.9999 8068.6730
find_modes(pz)
[1] 0.7622596 0.8014084 0.8116850 0.8381105 0.9027061 0.9692591 0.9849187
smoothScatter(libsize,pz,xlab="Library size",ylab="Proportion zeroes per cell",main="Young 1")
smoothScatter(libsize,numgene,xlab="Library size",ylab="Detected genes",main="Young 1")

Version Author Date
2b103a6 Belinda Phipson 2019-06-06 
par(mfrow=c(1,3))
plot(density(libsize),main="Distribution of library sizes")
#abline(v=3000,col=2)
plot(density(pz), main="Distribution of proportion of zeroes per cell")
#abline(v=0.95,col=2)
plot(density(numgene), main="Distribution of detected genes per cell")
abline(v=500,col=4, lty=3)
legend("topright",lty=2,col=4,legend="500 genes")

Version Author Date
2b103a6 Belinda Phipson 2019-06-06 
mycol <- rep(1,length(libsize))
mycol[pz>0.95] <- 2

par(mfrow=c(1,2))
propmito <- colSums(y1[mito,])/libsize
propribo <- colSums(y1[ribo,])/libsize
plot(libsize,propribo,xlab="Library size",ylab="Proportion ribosomal counts",main="Young 1: Ribosomal content",col=mycol)
legend("topright",fill=c(1,2),legend=c("pz <= 0.95","pz > 0.95"))
plot(libsize,propmito,xlab="Library size",ylab="Proportion mitochondrial counts",main="Young 1: Mitochondrial content",col=mycol)
legend("topright",fill=c(1,2),legend=c("pz <= 0.95","pz > 0.95"))

Version Author Date
2b103a6 Belinda Phipson 2019-06-06

Young 2

# Young 2
par(mfrow=c(1,2))
libsize <- colSums(y2)
pz <- colMeans(y2==0)
numgene <- colSums(y2!=0)

find_modes(libsize)
 [1]   6020.806  54822.374  69591.269  76333.591  83396.976  97523.745
 [7] 112934.766 118392.836 121282.403 122887.718 123850.907 124814.095
[13] 125777.284 126740.473 128345.788 129308.977 130272.165 131235.354
[19] 132519.606 133161.732 134445.984 135730.235 136693.424 137977.676
[25] 138940.865 139904.054 140546.180 141188.305 143435.746 144077.872
[31] 145041.061 156278.264
find_modes(numgene)
[1]  3049.700  9748.472 12208.177
find_modes(pz)
[1] 0.6402906 0.7127649 0.9101417
smoothScatter(libsize,pz,xlab="Library size",ylab="Proportion zeroes per cell",main="Young 2")
smoothScatter(libsize,numgene,xlab="Library size",ylab="Detected genes",main="Young 2")

Version Author Date
2b103a6 Belinda Phipson 2019-06-06 
par(mfrow=c(1,3))
plot(density(libsize),main="Distribution of library sizes")
#abline(v=3000,col=2)
plot(density(pz), main="Distribution of proportion of zeroes per cell")
#abline(v=0.95,col=2)
plot(density(numgene), main="Distribution of detected genes per cell")
abline(v=500,col=4, lty=3)
legend("topright",lty=2,col=4,legend="500 genes")

Version Author Date
2b103a6 Belinda Phipson 2019-06-06 
mycol <- rep(1,length(libsize))
mycol[pz>0.95] <- 2

par(mfrow=c(1,2))
propmito <- colSums(y2[mito,])/libsize
propribo <- colSums(y2[ribo,])/libsize
plot(libsize,propribo,xlab="Library size",ylab="Proportion ribosomal counts",main="Young 2: Ribosomal content",col=mycol)
legend("topright",fill=c(1,2),legend=c("pz <= 0.95","pz > 0.95"))
plot(libsize,propmito,xlab="Library size",ylab="Proportion mitochondrial counts",main="Young 2: Mitochondrial content",col=mycol)
legend("topright",fill=c(1,2),legend=c("pz <= 0.95","pz > 0.95"))

Version Author Date
2b103a6 Belinda Phipson 2019-06-06

Young 3

# Young 3
par(mfrow=c(1,2))
libsize <- colSums(y3)
pz <- colMeans(y3==0)
numgene <- colSums(y3!=0)

find_modes(libsize)
 [1]  1092.059  4834.037 16416.351 19623.761 24791.255 27820.476 30671.507
 [8] 33522.538 34948.054 38333.653 45817.610 48668.641 51519.673 54192.514
[15] 57756.303 61498.282 64705.692 67913.102 70051.375 72189.649 75753.438
[22] 81990.068 88048.509
find_modes(numgene)
[1]  703.5671 2840.3481
find_modes(pz)
[1] 0.9163102 0.9792697
smoothScatter(libsize,pz,xlab="Library size",ylab="Proportion zeroes per cell",main="Young 3")
smoothScatter(libsize,numgene,xlab="Library size",ylab="Detected genes",main="Young 3")

Version Author Date
2b103a6 Belinda Phipson 2019-06-06 
par(mfrow=c(1,3))
plot(density(libsize),main="Distribution of library sizes")
#abline(v=3000,col=2)
plot(density(pz), main="Distribution of proportion of zeroes per cell")
#abline(v=0.95,col=2)
plot(density(numgene), main="Distribution of detected genes per cell")
abline(v=500,col=4, lty=3)
legend("topright",lty=2,col=4,legend="500 genes")

Version Author Date
2b103a6 Belinda Phipson 2019-06-06 
mycol <- rep(1,length(libsize))
mycol[pz>0.95] <- 2

par(mfrow=c(1,2))
propmito <- colSums(y3[mito,])/libsize
propribo <- colSums(y3[ribo,])/libsize
plot(libsize,propribo,xlab="Library size",ylab="Proportion ribosomal counts",main="Young 3: Ribosomal content",col=mycol)
legend("topright",fill=c(1,2),legend=c("pz <= 0.95","pz > 0.95"))
plot(libsize,propmito,xlab="Library size",ylab="Proportion mitochondrial counts",main="Young 3: Mitochondrial content",col=mycol)
legend("topright",fill=c(1,2),legend=c("pz <= 0.95","pz > 0.95"))

Version Author Date
2b103a6 Belinda Phipson 2019-06-06

Adult 1

# Adult 1
par(mfrow=c(1,2))
libsize <- colSums(a1)
pz <- colMeans(a1==0)
numgene <- colSums(a1!=0)

find_modes(libsize)
 [1]   7571.146  52089.152  57964.625  64518.037  74461.145  87341.989
 [7]  87793.949  89375.807  90279.726  91183.644  94573.340 107002.225
find_modes(numgene)
[1]   520.9095  3410.7457  7681.5657  9190.4182 12080.2545
find_modes(pz)
[1] 0.6440598 0.7292077 0.7736655 0.8995036 0.9846516
smoothScatter(libsize,pz,xlab="Library size",ylab="Proportion zeroes per cell",main="Adult 1")
smoothScatter(libsize,numgene,xlab="Library size",ylab="Detected genes",main="Adult 1")

Version Author Date
2b103a6 Belinda Phipson 2019-06-06 
par(mfrow=c(1,3))
plot(density(libsize),main="Distribution of library sizes")
#abline(v=3000,col=2)
plot(density(pz), main="Distribution of proportion of zeroes per cell")
#abline(v=0.95,col=2)
plot(density(numgene), main="Distribution of detected genes per cell")
abline(v=500,col=4, lty=3)
legend("topright",lty=2,col=4,legend="500 genes")

Version Author Date
2b103a6 Belinda Phipson 2019-06-06 
mycol <- rep(1,length(libsize))
mycol[pz>0.95] <- 2

par(mfrow=c(1,2))
propmito <- colSums(a1[mito,])/libsize
propribo <- colSums(a1[ribo,])/libsize
plot(libsize,propribo,xlab="Library size",ylab="Proportion ribosomal counts",main="Adult 1: Ribosomal content",col=mycol)
legend("topright",fill=c(1,2),legend=c("pz <= 0.95","pz > 0.95"))
plot(libsize,propmito,xlab="Library size",ylab="Proportion mitochondrial counts",main="Adult 1: Mitochondrial content",col=mycol)
legend("topright",fill=c(1,2),legend=c("pz <= 0.95","pz > 0.95"))

Version Author Date
2b103a6 Belinda Phipson 2019-06-06

Adult 2

# Adult 2
par(mfrow=c(1,2))
libsize <- colSums(a2)
pz <- colMeans(a2==0)
numgene <- colSums(a2!=0)

find_modes(libsize)
 [1]   2188.160   7238.031  32992.375  39557.207  51171.911  60261.680
 [7]  67331.499  72886.358  78441.216  88540.959  94600.804 107730.470
[13] 113285.328 115810.264 121870.109 126919.981 130454.891 139039.672
[19] 150654.376 152674.324 154694.273 156714.221 158734.170 160754.118
[25] 162269.080 164794.015 166813.964 168328.925 184993.501 186508.462
[31] 188528.410 190548.359 192568.307 195598.230 196608.204 199638.127
[37] 200648.102 203678.024 212262.805 216302.703 218322.651 220342.600
[43] 222362.548 224382.497 225897.458 226907.432 228927.381 251651.802
find_modes(numgene)
[1]   889.0499  2908.6962  7241.1634 11052.4314 15254.5988
find_modes(pz)
[1] 0.5505289 0.6743442 0.7866418 0.9142963 0.9738045
smoothScatter(libsize,pz,xlab="Library size",ylab="Proportion zeroes per cell",main="Adult 2")
smoothScatter(libsize,numgene,xlab="Library size",ylab="Detected genes",main="Adult 2")

Version Author Date
2b103a6 Belinda Phipson 2019-06-06 
par(mfrow=c(1,3))
plot(density(libsize),main="Distribution of library sizes")
#abline(v=3000,col=2)
plot(density(pz), main="Distribution of proportion of zeroes per cell")
#abline(v=0.95,col=2)
plot(density(numgene), main="Distribution of detected genes per cell")
abline(v=500,col=4, lty=3)
legend("topright",lty=2,col=4,legend="500 genes")

Version Author Date
2b103a6 Belinda Phipson 2019-06-06 
mycol <- rep(1,length(libsize))
mycol[pz>0.95] <- 2

par(mfrow=c(1,2))
propmito <- colSums(a2[mito,])/libsize
propribo <- colSums(a2[ribo,])/libsize
plot(libsize,propribo,xlab="Library size",ylab="Proportion ribosomal counts",main="Adult 2: Ribosomal content",col=mycol)
legend("topright",fill=c(1,2),legend=c("pz <= 0.95","pz > 0.95"))
plot(libsize,propmito,xlab="Library size",ylab="Proportion mitochondrial counts",main="Adult 2: Mitochondrial content",col=mycol)
legend("topright",fill=c(1,2),legend=c("pz <= 0.95","pz > 0.95"))

Version Author Date
2b103a6 Belinda Phipson 2019-06-06

Adult 3

# Adult 3
par(mfrow=c(1,2))
libsize <- colSums(a3)
pz <- colMeans(a3==0)
numgene <- colSums(a3!=0)

find_modes(libsize)
 [1]   2212.058   9412.938  30168.417  39487.204  50076.734  54736.127
 [7]  62360.589  73373.700  86504.718  89893.367  96247.085 109378.103
[13] 122509.120 137334.462 155972.035 165714.403 177151.095 179692.583
[19] 180539.745 181810.489 183081.232 184775.557 186469.882 188587.788
[25] 189858.532 191976.438 193247.181 195365.087 196635.831 198330.156
[31] 209343.267
find_modes(numgene)
[1]   897.0151  3625.1939  6353.3728  7540.9330 10397.4967 13607.1189
find_modes(pz)
[1] 0.5990713 0.6936416 0.7778092 0.8128002 0.8931850 0.9735698
smoothScatter(libsize,pz,xlab="Library size",ylab="Proportion zeroes per cell",main="Adult 3")
smoothScatter(libsize,numgene,xlab="Library size",ylab="Detected genes",main="Adult 3")

Version Author Date
2b103a6 Belinda Phipson 2019-06-06 
par(mfrow=c(1,3))
plot(density(libsize),main="Distribution of library sizes")
#abline(v=3000,col=2)
plot(density(pz), main="Distribution of proportion of zeroes per cell")
#abline(v=0.95,col=2)
plot(density(numgene), main="Distribution of detected genes per cell")
abline(v=500,col=4, lty=3)
legend("topright",lty=2,col=4,legend="500 genes")

Version Author Date
2b103a6 Belinda Phipson 2019-06-06 
mycol <- rep(1,length(libsize))
mycol[pz>0.95] <- 2

par(mfrow=c(1,2))
propmito <- colSums(a3[mito,])/libsize
propribo <- colSums(a3[ribo,])/libsize
plot(libsize,propribo,xlab="Library size",ylab="Proportion ribosomal counts",main="Adult 3: Ribosomal content",col=mycol)
legend("topright",fill=c(1,2),legend=c("pz <= 0.95","pz > 0.95"))
plot(libsize,propmito,xlab="Library size",ylab="Proportion mitochondrial counts",main="Adult 3: Mitochondrial content",col=mycol)
legend("topright",fill=c(1,2),legend=c("pz <= 0.95","pz > 0.95"))

Version Author Date
2b103a6 Belinda Phipson 2019-06-06

Disease 1

# Disease 1
par(mfrow=c(1,2))
libsize <- colSums(d1)
pz <- colMeans(d1==0)
numgene <- colSums(d1!=0)

find_modes(libsize)
 [1]   741.9607  1973.8215  8133.1256  9892.9268 11652.7280 12708.6087
 [7] 13940.4695 15524.2905 16228.2110 17284.0917 19043.8929 20275.7537
[13] 21507.6145 22387.5151 24499.2765 25907.1175 27490.9385 28546.8192
[19] 29954.6602 31890.4415 32946.3222 34354.1631 35234.0637 38225.7257
[25] 39281.6064 41041.4076 42273.2684 43681.1094 45088.9503 46848.7515
[31] 48960.5129 50192.3737 51776.1948 53008.0556 55295.7971 56703.6381
[37] 58463.4392 58991.3796 60751.1808 64270.7831 66558.5246 67438.4252
[43] 70430.0872 75533.5106 77997.2323 80988.8943 82924.6756 88732.0194
find_modes(numgene)
[1]  512.4921 1325.6812 4563.6522 5228.9887 5790.8284 6160.4598 6751.8700
[8] 7121.5014
find_modes(pz)
[1] 0.7901676 0.8010587 0.8184843 0.8293754 0.8459298 0.8655337 0.9609393
[8] 0.9848996
smoothScatter(libsize,pz,xlab="Library size",ylab="Proportion zeroes per cell",main="Disease 1")
smoothScatter(libsize,numgene,xlab="Library size",ylab="Detected genes",main="Disease 1")

Version Author Date
2b103a6 Belinda Phipson 2019-06-06 
par(mfrow=c(1,3))
plot(density(libsize),main="Distribution of library sizes")
#abline(v=3000,col=2)
plot(density(pz), main="Distribution of proportion of zeroes per cell")
#abline(v=0.95,col=2)
plot(density(numgene), main="Distribution of detected genes per cell")
abline(v=500,col=4, lty=3)
legend("topright",lty=2,col=4,legend="500 genes")

Version Author Date
2b103a6 Belinda Phipson 2019-06-06 
mycol <- rep(1,length(libsize))
mycol[pz>0.95] <- 2

par(mfrow=c(1,2))
propmito <- colSums(d1[mito,])/libsize
propribo <- colSums(d1[ribo,])/libsize
plot(libsize,propribo,xlab="Library size",ylab="Proportion ribosomal counts",main="Disease 1: Ribosomal content",col=mycol)
legend("topright",fill=c(1,2),legend=c("pz <= 0.95","pz > 0.95"))
plot(libsize,propmito,xlab="Library size",ylab="Proportion mitochondrial counts",main="Disease 1: Mitochondrial content",col=mycol)
legend("topright",fill=c(1,2),legend=c("pz <= 0.95","pz > 0.95"))

Version Author Date
2b103a6 Belinda Phipson 2019-06-06

Disease 2

# Disease 2
par(mfrow=c(1,2))
libsize <- colSums(d2)
pz <- colMeans(d2==0)
numgene <- colSums(d2!=0)

find_modes(libsize)
 [1]   1098.901   3315.519  17881.866  20098.484  21681.783  24215.060
 [7]  25798.359  27381.658  28964.956  32131.553  35298.150  37198.109
[13]  39414.727  41314.685  42897.984  46064.581  50181.157  51447.796
[19]  53031.094  56831.011  59047.629  63797.525  65064.164  66647.462
[25]  68864.080  69814.059  72980.656  76147.254  79313.851  80897.149
[31]  82797.108  84697.066  86913.684  88813.642  91980.239  93563.538
[37]  97680.114  99896.732 101796.690 103063.329 106229.926 108129.885
[43] 110663.162 113196.440 115729.718 118896.315 121429.593 125229.509
[49] 126496.148 127129.468 128396.106 129662.745 132196.023 140429.176
[55] 146129.050 155312.182 159745.418
find_modes(numgene)
[1]   535.0048  2029.4649  6149.9047  6854.4359  8391.5948 10377.0918
find_modes(pz)
[1] 0.6942429 0.7527448 0.7980366 0.8187953 0.9402026 0.9842363
smoothScatter(libsize,pz,xlab="Library size",ylab="Proportion zeroes per cell",main="Disease 2")
smoothScatter(libsize,numgene,xlab="Library size",ylab="Detected genes",main="Disease 2")

Version Author Date
2b103a6 Belinda Phipson 2019-06-06 
par(mfrow=c(1,3))
plot(density(libsize),main="Distribution of library sizes")
#abline(v=3000,col=2)
plot(density(pz), main="Distribution of proportion of zeroes per cell")
#abline(v=0.95,col=2)
plot(density(numgene), main="Distribution of detected genes per cell")
abline(v=500,col=4, lty=3)
legend("topright",lty=2,col=4,legend="500 genes")

Version Author Date
2b103a6 Belinda Phipson 2019-06-06 
mycol <- rep(1,length(libsize))
mycol[pz>0.95] <- 2

par(mfrow=c(1,2))
propmito <- colSums(d2[mito,])/libsize
propribo <- colSums(d2[ribo,])/libsize
plot(libsize,propribo,xlab="Library size",ylab="Proportion ribosomal counts",main="Disease 2: Ribosomal content",col=mycol)
legend("topright",fill=c(1,2),legend=c("pz <= 0.95","pz > 0.95"))
plot(libsize,propmito,xlab="Library size",ylab="Proportion mitochondrial counts",main="Disease 2: Mitochondrial content",col=mycol)
legend("topright",fill=c(1,2),legend=c("pz <= 0.95","pz > 0.95"))

Version Author Date
2b103a6 Belinda Phipson 2019-06-06

Disease 3

# Disease 3
par(mfrow=c(1,2))
libsize <- colSums(d3)
pz <- colMeans(d3==0)
numgene <- colSums(d3!=0)

find_modes(libsize)
[1]   4073.714  15918.699  73959.127  87936.209  94332.501 110204.781
find_modes(numgene)
[1] 2348.607 4571.828
find_modes(pz)
[1] 0.8652928 0.9307992
smoothScatter(libsize,pz,xlab="Library size",ylab="Proportion zeroes per cell",main="Disease 3")
smoothScatter(libsize,numgene,xlab="Library size",ylab="Detected genes",main="Disease 3")

Version Author Date
2b103a6 Belinda Phipson 2019-06-06 
par(mfrow=c(1,3))
plot(density(libsize),main="Distribution of library sizes")
#abline(v=3000,col=2)
plot(density(pz), main="Distribution of proportion of zeroes per cell")
#abline(v=0.95,col=2)
plot(density(numgene), main="Distribution of detected genes per cell")
abline(v=500,col=4, lty=3)
legend("topright",lty=2,col=4,legend="500 genes")

Version Author Date
2b103a6 Belinda Phipson 2019-06-06 
mycol <- rep(1,length(libsize))
mycol[pz>0.95] <- 2

par(mfrow=c(1,2))
propmito <- colSums(d3[mito,])/libsize
propribo <- colSums(d3[ribo,])/libsize
plot(libsize,propribo,xlab="Library size",ylab="Proportion ribosomal counts",main="Disease 3: Ribosomal content",col=mycol)
legend("topright",fill=c(1,2),legend=c("pz <= 0.95","pz > 0.95"))
plot(libsize,propmito,xlab="Library size",ylab="Proportion mitochondrial counts",main="Disease 3: Mitochondrial content",col=mycol)
legend("topright",fill=c(1,2),legend=c("pz <= 0.95","pz > 0.95"))

Version Author Date
2b103a6 Belinda Phipson 2019-06-06

Disease 4

# Disease 4
par(mfrow=c(1,2))
libsize <- colSums(d4)
pz <- colMeans(d4==0)
numgene <- colSums(d4!=0)

find_modes(libsize)
 [1]   5500.593  19798.398  23264.533  28030.468  30196.802  32796.403
 [7]  36262.538  40595.207  49693.810  53593.212  56626.079  60958.748
[13]  63991.616  69190.818  75256.553  81755.556  90420.893  96919.895
[19] 100386.030 108618.100 119016.504 125082.240 134614.110 137213.711
[25] 138513.511 140246.579 143279.447 153244.584 160610.120 161476.654
[31] 162343.187 164076.255 165809.322 167542.389 169275.457 173174.858
[37] 174907.926 178374.060 179240.594 182273.462 187039.397 188772.464
[43] 190505.532 192238.599 193971.666 195271.467 197004.534 198737.602
[49] 199604.135 200470.669 201337.203 202203.736 203936.804 205669.871
[55] 206536.405 207402.938 208269.472 210869.073 218234.609
find_modes(numgene)
[1]  2249.813  4643.692  5338.689  9972.003 11413.478 12262.919
find_modes(pz)
[1] 0.6386777 0.6637061 0.7061787 0.8426975 0.8631754 0.9337101
smoothScatter(libsize,pz,xlab="Library size",ylab="Proportion zeroes per cell",main="Disease 4")
smoothScatter(libsize,numgene,xlab="Library size",ylab="Detected genes",main="Disease 4")

Version Author Date
2b103a6 Belinda Phipson 2019-06-06 
par(mfrow=c(1,3))
plot(density(libsize),main="Distribution of library sizes")
#abline(v=3000,col=2)
plot(density(pz), main="Distribution of proportion of zeroes per cell")
#abline(v=0.95,col=2)
plot(density(numgene), main="Distribution of detected genes per cell")
abline(v=500,col=4, lty=3)
legend("topright",lty=2,col=4,legend="500 genes")

Version Author Date
2b103a6 Belinda Phipson 2019-06-06 
mycol <- rep(1,length(libsize))
mycol[pz>0.95] <- 2

par(mfrow=c(1,2))
propmito <- colSums(d4[mito,])/libsize
propribo <- colSums(d4[ribo,])/libsize
plot(libsize,propribo,xlab="Library size",ylab="Proportion ribosomal counts",main="Disease 4: Ribosomal content",col=mycol)
legend("topright",fill=c(1,2),legend=c("pz <= 0.95","pz > 0.95"))
plot(libsize,propmito,xlab="Library size",ylab="Proportion mitochondrial counts",main="Disease 4: Mitochondrial content",col=mycol)
legend("topright",fill=c(1,2),legend=c("pz <= 0.95","pz > 0.95"))

Version Author Date
2b103a6 Belinda Phipson 2019-06-06

Comparing between samples

Here I am looking at the library sizes, number of detected genes, first modes of both library size distributions and number of detected genes per cell, as well as numbers of cells for each sample. This will help highlight if any of the samples look vastly different to the others in the dataset.

Calculating statistics

libmode <- rep(NA,13)
numgmode <- rep(NA,13)
numcells <- rep(NA,13)
names(libmode) <- names(numgmode) <- names(numcells) <- c("Fetal_1","Fetal_2","Fetal_3",
                                       "Young_1","Young_2","Young_3",
                                       "Adult_1","Adult_2","Adult_3",
                                       "Diseased_1","Diseased_2",
                                       "Diseased_3","Diseased_4")

mylibs <- mynumg <- vector("list", 13) 

mylibs[[1]] <- colSums(f1)
mynumg[[1]] <- colSums(f1!=0)
numcells[1] <- ncol(f1)

libmode[1] <- find_modes(mylibs[[1]])[1]
numgmode[1] <- find_modes(mynumg[[1]])[1]

mylibs[[2]] <- colSums(f2)
mynumg[[2]] <- colSums(f2!=0)
numcells[2] <- ncol(f2)

libmode[2] <- find_modes(mylibs[[2]])[1]
numgmode[2] <- find_modes(mynumg[[2]])[1]

mylibs[[3]] <- colSums(f3)
mynumg[[3]] <- colSums(f3!=0)
libmode[3] <- find_modes(mylibs[[3]])[1]
numgmode[3] <- find_modes(mynumg[[3]])[1]
numcells[3] <- ncol(f3)

mylibs[[4]] <- colSums(y1)
mynumg[[4]] <- colSums(y1!=0)
libmode[4] <- find_modes(mylibs[[4]])[1]
numgmode[4] <- find_modes(mynumg[[4]])[1]
numcells[4] <- ncol(y1)

mylibs[[5]] <- colSums(y2)
mynumg[[5]] <- colSums(y2!=0)
libmode[5] <- find_modes(mylibs[[5]])[1]
numgmode[5] <- find_modes(mynumg[[5]])[1]
numcells[5] <- ncol(y2)

mylibs[[6]] <- colSums(y3)
mynumg[[6]] <- colSums(y3!=0)
libmode[6] <- find_modes(mylibs[[6]])[1]
numgmode[6] <- find_modes(mynumg[[6]])[1]
numcells[6] <- ncol(y3)

mylibs[[7]] <- colSums(a1)
mynumg[[7]] <- colSums(a1!=0)
libmode[7] <- find_modes(mylibs[[7]])[1]
numgmode[7] <- find_modes(mynumg[[7]])[1]
numcells[7] <- ncol(a1)

mylibs[[8]] <- colSums(a2)
mynumg[[8]] <- colSums(a2!=0)
libmode[8] <- find_modes(mylibs[[8]])[1]
numgmode[8] <- find_modes(mynumg[[8]])[1]
numcells[8] <- ncol(a2)

mylibs[[9]] <- colSums(a3)
mynumg[[9]] <- colSums(a3!=0)
libmode[9] <- find_modes(mylibs[[9]])[1]
numgmode[9] <- find_modes(mynumg[[9]])[1]
numcells[9] <- ncol(a3)

mylibs[[10]] <- colSums(d1)
mynumg[[10]] <- colSums(d1!=0)
libmode[10] <- find_modes(mylibs[[10]])[1]
numgmode[10] <- find_modes(mynumg[[10]])[1]
numcells[10] <- ncol(d1)

mylibs[[11]] <- colSums(d2)
mynumg[[11]] <- colSums(d2!=0)
libmode[11] <- find_modes(mylibs[[11]])[1]
numgmode[11] <- find_modes(mynumg[[11]])[1]
numcells[11] <- ncol(d2)

mylibs[[12]] <- colSums(d3)
mynumg[[12]] <- colSums(d3!=0)
libmode[12] <- find_modes(mylibs[[12]])[1]
numgmode[12] <- find_modes(mynumg[[12]])[1]
numcells[12] <- ncol(d3)

mylibs[[13]] <- colSums(d4)
mynumg[[13]] <- colSums(d4!=0)
libmode[13] <- find_modes(mylibs[[13]])[1]
numgmode[13] <- find_modes(mynumg[[13]])[1]
numcells[13] <- ncol(d4)

Number of cells per sample

par(mar=c(6,4,2,2))
mycols <- rep(ggplotColors(4),c(3,3,3,4))
barplot(numcells,col=mycols,las=2,main="Number of cells per sample")
abline(h=1000,lty=2,lwd=2)

Version Author Date
5839497 Belinda Phipson 2019-06-07
2b103a6 Belinda Phipson 2019-06-06

Library size distributions (sequencing depth)

plot(density(mylibs[[1]]),lwd=2,col=mycols[1],ylim=c(0,0.0003),main="Library size distributions")
for(i in 2:13) lines(density(mylibs[[i]]),lwd=2,col=mycols[i])
legend("topright",legend=c("Fetal","Young","Adult","Diseased"),fill=ggplotColors(4))

Version Author Date
5839497 Belinda Phipson 2019-06-07
2b103a6 Belinda Phipson 2019-06-06 
#abline(v=1000,lty=2,lwd=2)
par(mfrow=c(1,1))
par(mar=c(6,4,2,2))
boxplot(mylibs,col=mycols,names=names(libmode),las=2, main="Library size distributions",ylim=c(0,150000))

Version Author Date
2b103a6 Belinda Phipson 2019-06-06 
boxplot(mylibs,col=ggplotColors(2)[factor(targets$Batch)],names=names(libmode),las=2, main="Library size distributions by batch",ylim=c(0,150000))
legend("topleft",fill=ggplotColors(2),legend=levels(factor(targets$Batch)))

Number of genes detected

plot(density(mynumg[[1]]),lwd=2,col=mycols[1],ylim=c(0,0.0007),main="Number of genes detected distributions")
for(i in 2:13) lines(density(mynumg[[i]]),lwd=2,col=mycols[i])
legend("topright",legend=c("Fetal","Young","Adult","Diseased"),fill=ggplotColors(4))

Version Author Date
5839497 Belinda Phipson 2019-06-07 
par(mfrow=c(1,1))
par(mar=c(6,4,2,2))
boxplot(mynumg,col=mycols,names=names(libmode),las=2, main="Number of genes detected")

Version Author Date
5839497 Belinda Phipson 2019-06-07 
boxplot(mynumg,col=ggplotColors(2)[factor(targets$Batch)],names=names(libmode),las=2, main="Number of genes detected by batch")
legend("topleft",fill=ggplotColors(2),legend=levels(factor(targets$Batch)))

Version Author Date
5839497 Belinda Phipson 2019-06-07

Examining the first modes from density ditributions

barplot(libmode,col=mycols,las=2)
title("First mode: library size distributions")
abline(h=1000,lty=2,lwd=2)

barplot(numgmode,col=mycols,las=2)
title("First mode: number of genes detected distributions")
abline(h=500,lty=2,lwd=2)

MDS plot of all samples

To get a high-level idea of the overall sources of variability in the dataset, I have summed the counts over all cells within a sample to obtain a “pseudobulk” sample and made MDS plots using functions in edgeR.

pseudobulk <- matrix(NA,ncol=13,nrow=nrow(f1))
colnames(pseudobulk) <- group
rownames(pseudobulk) <- rownames(f1)

pseudobulk[,1] <- rowSums(f1)
pseudobulk[,2] <- rowSums(f2)
pseudobulk[,3] <- rowSums(f3)
pseudobulk[,4] <- rowSums(y1)
pseudobulk[,5] <- rowSums(y2)
pseudobulk[,6] <- rowSums(y3)
pseudobulk[,7] <- rowSums(a1)
pseudobulk[,8] <- rowSums(a2)
pseudobulk[,9] <- rowSums(a3)
pseudobulk[,10] <- rowSums(d1)
pseudobulk[,11] <- rowSums(d2)
pseudobulk[,12] <- rowSums(d3)
pseudobulk[,13] <- rowSums(d4)
y <- DGEList(pseudobulk)
keep <- rowSums(y$counts)>10
y.keep <- y[keep,]
#pdf(file="/group/bioi1/belinda/SingleCell/Cardiac/MDSplots.pdf",
#    width = 12,height = 5)
par(mfrow=c(1,2))
par(mar=c(4,4,2,2))
plotMDS(y.keep,col=ggplotColors(4)[factor(targets$Group)])
legend("bottom",fill=ggplotColors(4),legend=levels(factor(targets$Group)),bty="n")

plotMDS(y.keep,dim=c(3,4),col=ggplotColors(2)[factor(targets$Sex)])
legend("bottomright",fill=ggplotColors(2),legend=levels(factor(targets$Sex)),bty="n")

#dev.off()

sessionInfo()
R version 3.6.0 (2019-04-26)
Platform: x86_64-pc-linux-gnu (64-bit)
Running under: CentOS release 6.7 (Final)

Matrix products: default
BLAS:   /usr/local/installed/R/3.6.0/lib64/R/lib/libRblas.so
LAPACK: /usr/local/installed/R/3.6.0/lib64/R/lib/libRlapack.so

locale:
 [1] LC_CTYPE=en_US.UTF-8       LC_NUMERIC=C              
 [3] LC_TIME=en_US.UTF-8        LC_COLLATE=en_US.UTF-8    
 [5] LC_MONETARY=en_US.UTF-8    LC_MESSAGES=en_US.UTF-8   
 [7] LC_PAPER=en_US.UTF-8       LC_NAME=C                 
 [9] LC_ADDRESS=C               LC_TELEPHONE=C            
[11] LC_MEASUREMENT=en_US.UTF-8 LC_IDENTIFICATION=C       

attached base packages:
 [1] splines   parallel  stats4    stats     graphics  grDevices utils    
 [8] datasets  methods   base     

other attached packages:
 [1] workflowr_1.3.0             NMF_0.21.0                 
 [3] bigmemory_4.5.33            cluster_2.0.9              
 [5] rngtools_1.3.1.1            pkgmaker_0.27              
 [7] registry_0.5-1              scran_1.12.0               
 [9] SingleCellExperiment_1.6.0  SummarizedExperiment_1.14.0
[11] GenomicRanges_1.36.0        GenomeInfoDb_1.20.0        
[13] DelayedArray_0.10.0         BiocParallel_1.18.0        
[15] matrixStats_0.54.0          cowplot_0.9.4              
[17] monocle_2.12.0              DDRTree_0.1.5              
[19] irlba_2.3.3                 VGAM_1.1-1                 
[21] ggplot2_3.1.1               Matrix_1.2-17              
[23] Seurat_3.0.1                org.Hs.eg.db_3.8.2         
[25] AnnotationDbi_1.46.0        IRanges_2.18.0             
[27] S4Vectors_0.22.0            Biobase_2.44.0             
[29] BiocGenerics_0.30.0         RColorBrewer_1.1-2         
[31] edgeR_3.26.3                limma_3.40.2               

loaded via a namespace (and not attached):
  [1] backports_1.1.4          plyr_1.8.4              
  [3] igraph_1.2.4.1           lazyeval_0.2.2          
  [5] densityClust_0.3         listenv_0.7.0           
  [7] gridBase_0.4-7           scater_1.12.2           
  [9] fastICA_1.2-1            digest_0.6.18           
 [11] foreach_1.4.4            htmltools_0.3.6         
 [13] viridis_0.5.1            gdata_2.18.0            
 [15] magrittr_1.5             memoise_1.1.0           
 [17] doParallel_1.0.14        ROCR_1.0-7              
 [19] globals_0.12.4           R.utils_2.8.0           
 [21] docopt_0.6.1             colorspace_1.4-1        
 [23] blob_1.1.1               ggrepel_0.8.1           
 [25] xfun_0.6                 dplyr_0.8.1             
 [27] bigmemory.sri_0.1.3      sparsesvd_0.1-4         
 [29] crayon_1.3.4             RCurl_1.95-4.12         
 [31] jsonlite_1.6             iterators_1.0.10        
 [33] survival_2.44-1.1        zoo_1.8-5               
 [35] ape_5.3                  glue_1.3.1              
 [37] gtable_0.3.0             zlibbioc_1.30.0         
 [39] XVector_0.24.0           BiocSingular_1.0.0      
 [41] future.apply_1.2.0       scales_1.0.0            
 [43] pheatmap_1.0.12          DBI_1.0.0               
 [45] bibtex_0.4.2             Rcpp_1.0.1              
 [47] metap_1.1                xtable_1.8-4            
 [49] viridisLite_0.3.0        dqrng_0.2.1             
 [51] reticulate_1.12          bit_1.1-14              
 [53] rsvd_1.0.0               SDMTools_1.1-221.1      
 [55] tsne_0.1-3               htmlwidgets_1.3         
 [57] httr_1.4.0               FNN_1.1.3               
 [59] gplots_3.0.1.1           ica_1.0-2               
 [61] pkgconfig_2.0.2          R.methodsS3_1.7.1       
 [63] dynamicTreeCut_1.63-1    locfit_1.5-9.1          
 [65] tidyselect_0.2.5         rlang_0.3.4             
 [67] reshape2_1.4.3           munsell_0.5.0           
 [69] tools_3.6.0              RSQLite_2.1.1           
 [71] ggridges_0.5.1           evaluate_0.13           
 [73] stringr_1.4.0            yaml_2.2.0              
 [75] npsurv_0.4-0             knitr_1.22              
 [77] bit64_0.9-7              fs_1.3.1                
 [79] fitdistrplus_1.0-14      caTools_1.17.1.2        
 [81] purrr_0.3.2              RANN_2.6.1              
 [83] pbapply_1.4-0            future_1.12.0           
 [85] nlme_3.1-139             whisker_0.3-2           
 [87] slam_0.1-45              R.oo_1.22.0             
 [89] compiler_3.6.0           beeswarm_0.2.3          
 [91] plotly_4.9.0             png_0.1-7               
 [93] lsei_1.2-0               statmod_1.4.30          
 [95] tibble_2.1.1             stringi_1.4.3           
 [97] lattice_0.20-38          HSMMSingleCell_1.4.0    
 [99] pillar_1.3.1             combinat_0.0-8          
[101] Rdpack_0.11-0            lmtest_0.9-37           
[103] BiocNeighbors_1.2.0      data.table_1.11.6       
[105] bitops_1.0-6             gbRd_0.4-11             
[107] R6_2.4.0                 KernSmooth_2.23-15      
[109] gridExtra_2.3            vipor_0.4.5             
[111] codetools_0.2-16         MASS_7.3-51.4           
[113] gtools_3.8.1             assertthat_0.2.1        
[115] rprojroot_1.3-2          withr_2.1.2             
[117] qlcMatrix_0.9.7          sctransform_0.2.0       
[119] GenomeInfoDbData_1.2.1   grid_3.6.0              
[121] tidyr_0.8.3              DelayedMatrixStats_1.6.0
[123] rmarkdown_1.12.6         Rtsne_0.15              
[125] git2r_0.25.2             ggbeeswarm_0.6.0