-
Notifications
You must be signed in to change notification settings - Fork 0
/
Copy pathmy-diff-map.R
90 lines (70 loc) · 2.93 KB
/
my-diff-map.R
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
# Michelle Li
# diffusion mapping using destiny package
# last updated: July 20, 2018
# Got destiny package working with Nestorowa data and paul data
# Added interactive 3D plot
# Sorted out how to choose dimensions when converting to ExpressionSet
# required packages
library(destiny)
library(Biobase)
library(rgl)
#### PREPARING DATA ----------------------------------------------
## NOTES ON DATA:
## Nestorowa data: 1645 cells x 4290 genes (need transpose for diff mapping)
## Paul data: 8716 genes x 2730 cells
## Grover data: 46175 genes x 135 cells
## Rockne data: 12674 genes x 101 cells
## read in data and convert to matrix (Nestorowa, Grover, non-.RData file)
# rawdata <- read.delim("~/GitHub/cell-diff-reu-2018/data/coordinates_gene_counts_flow_cytometry.txt", row.names=1)
rawdata <- read.csv("~/GitHub/cell-diff-reu-2018/data/grover_expression.txt", row.names=1, sep="")
cleandata <- na.omit(rawdata) # remove rows with NA
# data <- cleandata[,14:ncol(cleandata)] # choose relevant flow cytometry data
data <- cleandata
data <- as.matrix(data)
## load .RData file
# load("~/GitHub/cell-diff-reu-2018/data/Paul_Cell_MARSseq_GSE72857.RData")
dataname <- "Grover Data" # set name
## create factor for grouping
rownames <- row.names(X)
groups <- as.character(rownames)
groups[grepl("HSPC",groups)] <- "HSPC" # replace all labels with hspc
groups[grepl("LT",groups)] <- "LT.HSC" # replace with lt.hsc
groups[grepl("Prog",groups)] <- "PROG" # same for prog
groupsf <- factor(groups)
# groupsf <- factor(cluster.id) # paul data
nclusters <- nlevels(groupsf) # number of groups to cluster
colors <- rainbow(nclusters) # choose colors
## convert to expression set
## EXPRESSION SET DIMENSIONS: cells = columns ; genes = rows
# data <- t(data) # may need to transpose to have correct rows/columns
dataES <- ExpressionSet(assayData=data) # you MUST include 'assayData=' to be able to choose your dimensions properly
#### DIFFUSION MAP -------------------------------------------------
## diffusion map (un-normalized data)
dm <- DiffusionMap(dataES)
#### PLOTTING ------------------------------------------------------
# par(mar = c(5.1,4.1,4.1,8.1), xpd = "TRUE") # use if you need extra space in margins for legend
## plot 2d
plot(x=eigenvectors(dm)[,1], y=eigenvectors(dm)[,2],
xlab="DC1", ylab="DC2",
col = colors[groupsf],
pch=20,
main = paste("Diffusion Map (destiny) of", dataname))
legend("bottomleft", inset=c(0,0),
legend=levels(groupsf),
col=colors,
pch=20,
ncol=1, cex=0.4)
## plot 3d
plot3d(x=eigenvectors(dm)[,1], y=eigenvectors(dm)[,2], z=eigenvectors(dm)[,3],
xlab="DC1", ylab="DC2", zlab="DC3",
col = colors[groupsf],
pch=20,
main = paste("Diffusion Map (destiny) of", dataname))
legend3d("bottomright", inset=c(0.1,0.1),
legend = levels(groupsf),
pch = 20,
col = colors,
ncol=3, cex=1)
## sigmas
sigmas <- find_sigmas(dataES)
plot(sigmas)