R Notebook

This is an R Markdown Notebook. When you execute code within the notebook, the results appear beneath the code.

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

library(dplyr)
library(Seurat)
data_dir <- 'HIV'
HIV.data <- Read10X(data_dir)

library(dplyr)
library(Seurat)
data_dir <- 'Control'
control.data <- Read10X(data_dir)
control <- CreateSeuratObject(counts = control.data, project = "Control1", min.cells = 3, min.features = 200)

Feature names cannot have underscores ('_'), replacing with dashes ('-')

control[["percent.mt"]] <- PercentageFeatureSet(control, pattern = "^MT-")
control <- subset(control, subset = nFeature_RNA > 200 & nFeature_RNA < 2500 & percent.mt < 5)
control

An object of class Seurat 
17297 features across 5511 samples within 1 assay 
Active assay: RNA (17297 features)

library(dplyr)
library(Seurat)
data.dir <- 'NF'
NF.data <- Read10X(data.dir)
NF <- CreateSeuratObject(counts = NF.data, project = "NF1", min.cells = 3, min.features = 200)

Feature names cannot have underscores ('_'), replacing with dashes ('-')

NF[["percent.mt"]] <- PercentageFeatureSet(NF, pattern = "^MT-")
NF <- subset(NF, subset = nFeature_RNA > 200 & nFeature_RNA < 2500 & percent.mt < 5)
NF

An object of class Seurat 
17943 features across 6141 samples within 1 assay 
Active assay: RNA (17943 features)

control@meta.data$treatment <- "CTRL"
HIV@meta.data$treatment <- "HIV"
NF@meta.data$treatment <- "HIV/NF499"

tonsil <- merge(control, y = c(HIV, NF), add.cell.ids = c("control", "HIV", "NF"), project = "tonsil")
tonsil

An object of class Seurat 
19191 features across 15373 samples within 1 assay 
Active assay: RNA (19191 features)

head(colnames(tonsil))

[1] "control_AAACCTGAGTACCGGA" "control_AAACCTGAGTGTCCAT" "control_AAACCTGAGTTACCCA" "control_AAACCTGAGTTGCAGG"
[5] "control_AAACCTGCACGGACAA" "control_AAACCTGCACTCTGTC"

tonsil <- SplitObject (tonsil, split.by = "treatment")
tonsil <- tonsil[c("CTRL", "HIV", "HIV/NF499")]

tonsil <- lapply(X = tonsil, FUN = function(x) {
  x <- NormalizeData(x)
  x <- FindVariableFeatures(x, selection.method = "vst", nfeatures = 2000)
})

Performing log-normalization
0%   10   20   30   40   50   60   70   80   90   100%
[----|----|----|----|----|----|----|----|----|----|
**************************************************|
Calculating gene variances
0%   10   20   30   40   50   60   70   80   90   100%
[----|----|----|----|----|----|----|----|----|----|
**************************************************|
Calculating feature variances of standardized and clipped values
0%   10   20   30   40   50   60   70   80   90   100%
[----|----|----|----|----|----|----|----|----|----|
**************************************************|
Performing log-normalization
0%   10   20   30   40   50   60   70   80   90   100%
[----|----|----|----|----|----|----|----|----|----|
**************************************************|
Calculating gene variances
0%   10   20   30   40   50   60   70   80   90   100%
[----|----|----|----|----|----|----|----|----|----|
**************************************************|
Calculating feature variances of standardized and clipped values
0%   10   20   30   40   50   60   70   80   90   100%
[----|----|----|----|----|----|----|----|----|----|
**************************************************|
Performing log-normalization
0%   10   20   30   40   50   60   70   80   90   100%
[----|----|----|----|----|----|----|----|----|----|
**************************************************|
Calculating gene variances
0%   10   20   30   40   50   60   70   80   90   100%
[----|----|----|----|----|----|----|----|----|----|
**************************************************|
Calculating feature variances of standardized and clipped values
0%   10   20   30   40   50   60   70   80   90   100%
[----|----|----|----|----|----|----|----|----|----|
**************************************************|

tonsil.anchors <- FindIntegrationAnchors(object.list = tonsil, dims = 1:30)

Computing 2000 integration features
Scaling features for provided objects

  |                                                  | 0 % ~calculating  
  |+++++++++++++++++                                 | 33% ~02s          
  |++++++++++++++++++++++++++++++++++                | 67% ~01s          
  |++++++++++++++++++++++++++++++++++++++++++++++++++| 100% elapsed=03s  

Finding all pairwise anchors

  |                                                  | 0 % ~calculating  

Running CCA
Merging objects
Finding neighborhoods
Finding anchors
    Found 13579 anchors
Filtering anchors
    Retained 3294 anchors
Extracting within-dataset neighbors

  |+++++++++++++++++                                 | 33% ~03m 26s      

Running CCA
Warning in doTryCatch(return(expr), name, parentenv, handler) :
  restarting interrupted promise evaluation
Warning in doTryCatch(return(expr), name, parentenv, handler) :
  restarting interrupted promise evaluation
Merging objects
Finding neighborhoods
Finding anchors
    Found 17084 anchors
Filtering anchors
    Retained 3332 anchors
Extracting within-dataset neighbors

  |++++++++++++++++++++++++++++++++++                | 67% ~02m 15s      

Running CCA

tonsil.integrated <- IntegrateData(anchorset = tonsil.anchors, dims = 1:30)

Merging dataset 2 into 1
Extracting anchors for merged samples
Finding integration vectors
Finding integration vector weights
0%   10   20   30   40   50   60   70   80   90   100%
[----|----|----|----|----|----|----|----|----|----|
**************************************************|
Integrating data
Merging dataset 3 into 1 2
Extracting anchors for merged samples
Finding integration vectors
Finding integration vector weights
0%   10   20   30   40   50   60   70   80   90   100%
[----|----|----|----|----|----|----|----|----|----|
**************************************************|
Integrating data

library(ggplot2)
library(cowplot)

********************************************************
Note: As of version 1.0.0, cowplot does not change the
  default ggplot2 theme anymore. To recover the previous
  behavior, execute:
  theme_set(theme_cowplot())
********************************************************

DefaultAssay(tonsil.integrated) <- "integrated"
tonsil.integrated <- ScaleData(tonsil.integrated, verbose = FALSE)
tonsil.integrated <- RunPCA(tonsil.integrated, npcs = 30, verbose = FALSE)
tonsil.integrated <- RunUMAP(tonsil.integrated, reduction = "pca", dims = 1:20)

The default method for RunUMAP has changed from calling Python UMAP via reticulate to the R-native UWOT using the cosine metric
To use Python UMAP via reticulate, set umap.method to 'umap-learn' and metric to 'correlation'
This message will be shown once per session23:38:59 UMAP embedding parameters a = 0.9922 b = 1.112
23:38:59 Read 15373 rows and found 20 numeric columns
23:38:59 Using Annoy for neighbor search, n_neighbors = 30
23:38:59 Building Annoy index with metric = cosine, n_trees = 50
0%   10   20   30   40   50   60   70   80   90   100%
[----|----|----|----|----|----|----|----|----|----|
**************************************************|
23:39:03 Writing NN index file to temp file /var/folders/v8/pd7f0mt15fj37kdvjqks4lcm0000gn/T//Rtmpm9qfS0/file121bc4273dba6
23:39:03 Searching Annoy index using 1 thread, search_k = 3000
23:39:10 Annoy recall = 100%
23:39:11 Commencing smooth kNN distance calibration using 1 thread
23:39:12 Initializing from normalized Laplacian + noise
23:39:13 Commencing optimization for 200 epochs, with 712006 positive edges
0%   10   20   30   40   50   60   70   80   90   100%
[----|----|----|----|----|----|----|----|----|----|
**************************************************|
23:39:28 Optimization finished

tonsil.integrated <- FindNeighbors(tonsil.integrated, reduction = "pca", dims = 1:20)

Computing nearest neighbor graph
Computing SNN

tonsil.integrated <- FindClusters(tonsil.integrated, resolution = 0.5)

Modularity Optimizer version 1.3.0 by Ludo Waltman and Nees Jan van Eck

Number of nodes: 15373
Number of edges: 636454

Running Louvain algorithm...

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

Maximum modularity in 10 random starts: 0.8832
Number of communities: 14
Elapsed time: 4 seconds

p1 <- DimPlot(tonsil.integrated, reduction = "umap", group.by = "treatment")
p2 <- DimPlot(tonsil.integrated, reduction = "umap", label = TRUE)

Using `as.character()` on a quosure is deprecated as of rlang 0.3.0.
Please use `as_label()` or `as_name()` instead.
This warning is displayed once per session.

plot_grid(p1, p2)

DimPlot(tonsil.integrated, reduction = "umap", split.by = "treatment")

DefaultAssay(tonsil.integrated) <- "RNA"
markers0 <- FindConservedMarkers(tonsil.integrated, ident.1 = 0, grouping.var = "treatment", verbose = FALSE)

head(markers0, n = 30)

markers1 <- FindConservedMarkers(tonsil.integrated, ident.1 = 1, grouping.var = "treatment", verbose = FALSE)

head(markers1, n = 30)

markers2 <- FindConservedMarkers(tonsil.integrated, ident.1 = 2, grouping.var = "treatment", verbose = FALSE)

head(markers2, n = 30)

markers3 <- FindConservedMarkers(tonsil.integrated, ident.1 = 3, grouping.var = "treatment", verbose = FALSE)

head(markers3, n = 30)

markers4 <- FindConservedMarkers(tonsil.integrated, ident.1 = 4, grouping.var = "treatment", verbose = FALSE)

head(markers4, n = 30)

markers5 <- FindConservedMarkers(tonsil.integrated, ident.1 = 5, grouping.var = "treatment", verbose = FALSE)

head(markers5, n = 30)

markers6 <- FindConservedMarkers(tonsil.integrated, ident.1 = 6, grouping.var = "treatment", verbose = FALSE)

head(markers6, n = 30)

markers7 <- FindConservedMarkers(tonsil.integrated, ident.1 = 7, grouping.var = "treatment", verbose = FALSE)

head(markers7, n = 30)

markers8 <- FindConservedMarkers(tonsil.integrated, ident.1 = 8, grouping.var = "treatment", verbose = FALSE)

head(markers8, n = 30)

markers9 <- FindConservedMarkers(tonsil.integrated, ident.1 = 9, grouping.var = "treatment", verbose = FALSE)

head(markers9, n = 30)

markers10 <- FindConservedMarkers(tonsil.integrated, ident.1 = 10, grouping.var = "treatment", verbose = FALSE)

head(markers10, n = 20)

markers11 <- FindConservedMarkers(tonsil.integrated, ident.1 = 11, grouping.var = "treatment", verbose = FALSE)

head(markers11, n = 20)

markers12 <- FindConservedMarkers(tonsil.integrated, ident.1 = 12, grouping.var = "treatment", verbose = FALSE)

head(markers12, n = 20)

markers13 <- FindConservedMarkers(tonsil.integrated, ident.1 = 13, grouping.var = "treatment", verbose = FALSE)

head(markers13, n = 20)

Add a new chunk by clicking the Insert Chunk button on the toolbar or by pressing Cmd+Option+I.

When you save the notebook, an HTML file containing the code and output will be saved alongside it (click the Preview button or press Cmd+Shift+K to preview the HTML file).

The preview shows you a rendered HTML copy of the contents of the editor. Consequently, unlike Knit, Preview does not run any R code chunks. Instead, the output of the chunk when it was last run in the editor is displayed.