ROSMAP snRNA-seq: trajectory inference

Example pipeline to perform trajectory inference

rm(list = ls())

print(getwd())

source("./scripts/R_functions/MiscPreprocessing.R")
source("./scripts/R_functions/LineageFunctions.R")
require(Matrix)

## Loading required package: Matrix

## Warning: package 'Matrix' was built under R version 3.6.3

require(Seurat)

## Loading required package: Seurat

## Warning: package 'Seurat' was built under R version 3.6.3

Read in Seurat Object file.

obj = readRDS("Output/SeuratObject.RDS")
#urlstr <- "https://www.dropbox.com/s/w0mzab5i3dwxwlp/SeuratObject.RDS?raw=1"
#obj = readRDS(url(urlstr))

Prepare output directory.

outdir = 'Trajectory'
if(!file.exists(outdir)) dir.create(outdir, recursive = TRUE)
setwd(outdir)

As an example, we will perform trajectory inference for AD cells in cluster 5.

clusterid = 5
obj2 = obj[, Idents(obj) == clusterid & obj@meta.data$pathologic.diagnosis.of.AD == 'YES']
Dat = GetAssayData(obj2)
meta = obj2@meta.data
GeneNames = rownames(Dat)
Names = colnames(Dat)
cNames = meta$Barcode
l = length(Names)

#deleting columns not in the covariate list
temp = rep(T, l)
for (i in 1:l){
  if (!(Names[i] %in% cNames)){
    temp[i] = F
  }
}

In = which(temp)
#print(temp)
Dat = Dat[,In]

#deleting extra rows in covariate list
Names = Names[In]
l = length(cNames)
temp = rep(T,l)
for (i in 1:l){
  if (!(cNames[i] %in% Names)){
    temp[i] = F
  }
}
In = which(temp)
meta = meta[In, ]

In_genes = rownames(obj2)[rowSums(obj2@assays$SCT@counts>0) > 0.2*ncol(obj2)]
DatNorm = Dat
DatNorm2 = DatNorm[match(In_genes, rownames(DatNorm)), ]
GeneNamesAD = GeneNames[match(In_genes, GeneNames)]

DatNorm4 = DatNorm2
Dat4 = meta

temp = DatNorm4
temp2 = Dat4

#converting ENSG to gene symbols
gene_short_name = In_genes

Run monocle pipeline.

rownames(temp) = NULL
colnames(temp) = NULL
MonRun = RunMonocleTobit(temp, temp2, C_by = 'Pseudotime', gene_short_name = gene_short_name)

## Loading required package: Biobase

## 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:Matrix':
## 
##     which

## 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

## Welcome to Bioconductor
## 
##     Vignettes contain introductory material; view with
##     'browseVignettes()'. To cite Bioconductor, see
##     'citation("Biobase")', and for packages 'citation("pkgname")'.

## Loading required package: ggplot2

## Warning: package 'ggplot2' was built under R version 3.6.3

## Loading required package: VGAM

## Loading required package: stats4

## Loading required package: splines

## Loading required package: DDRTree

## Loading required package: irlba

## Warning in if (cds@expressionFamily@vfamily %in% c("negbinomial",
## "negbinomial.size")) {: the condition has length > 1 and only the first element
## will be used

## Warning in if (cds@expressionFamily@vfamily == "binomialff") {: the condition
## has length > 1 and only the first element will be used

## Warning in if (cds@expressionFamily@vfamily == "Tobit") {: the condition has
## length > 1 and only the first element will be used

## Warning in if (cds@expressionFamily@vfamily == "uninormal") {: the condition has
## length > 1 and only the first element will be used

## Warning: `select_()` is deprecated as of dplyr 0.7.0.
## Please use `select()` instead.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_warnings()` to see where this warning was generated.

Save outputs.

saveRDS(MonRun, file = paste0("cluster", clusterid, "-monocle.RDS"))
saveRDS(DatNorm4, file = paste0("cluster", clusterid, "-input_mat.RDS"))