-
Notifications
You must be signed in to change notification settings - Fork 28
/
sCellGBM.R
80 lines (61 loc) · 4.19 KB
/
sCellGBM.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
library(CaSpER)
## "sCell_gbm_data.rda" contains the following objects:
## data: normalized gene expression matrix
## loh.name.mapping: data.frame for mapping loh files to expression files
## annotation
## loh
data("scell_gbm")
data <- scell_gbm$data
loh <- scell_gbm$loh
loh.name.mapping <- scell_gbm$loh.name.mapping
## generate annotation data.frame
annotation <- generateAnnotation(id_type="hgnc_symbol", genes=rownames(data), ishg19=T, centromere)
data <- data[match(annotation$Gene,rownames(data)), ]
## create CaSpER object
object <- CreateCasperObject(raw.data=data,loh.name.mapping=loh.name.mapping,
sequencing.type="single-cell",
cnv.scale=3, loh.scale=3,
annotation=annotation, method="iterative", loh=loh,
control.sample.ids="REF", cytoband=cytoband)
## runCaSpER
final.objects <- runCaSpER(object, removeCentromere=T, cytoband=cytoband, method="iterative")
## plot median filtered gene expression matrix
plotHeatmap(object, fileName="heatmap.png", cnv.scale= 3, cluster_cols = F, cluster_rows = T, show_rownames = T, only_soi = T)
## summarize large scale events
finalChrMat <- extractLargeScaleEvents (final.objects, thr=0.75)
#### VISUALIZATION
## plot large scale events using event summary matrix 1: amplification, -1:deletion, 0: neutral
plotLargeScaleEvent2 (finalChrMat, fileName="large.scale.events.summarized.pdf")
## plot BAF deviation for all samples together in one plot (can be used only with small sample size)
plotBAFAllSamples (loh = final.objects[[1]]@loh.median.filtered.data, fileName="LOHAllSamples.png")
## plot BAF signal in different scales for all samples
plotBAFOneSample (object, fileName="LohPlotsAllScales.pdf")
## plot large scale event summary for selected sample and chromosomes
plotSingleCellLargeScaleEventHeatmap(finalChrMat, sampleName="MGH31", chrs=c("5p", "14q"))
plotSingleCellLargeScaleEventHeatmap(finalChrMat, sampleName="MGH29", chrs=c("10p", "4p"))
plotSingleCellLargeScaleEventHeatmap(finalChrMat, sampleName="MGH30", chrs=c("6p", "7p"))
plotSingleCellLargeScaleEventHeatmap(finalChrMat, sampleName="MGH31", chrs=c("5p", "14q"))
plotSingleCellLargeScaleEventHeatmap(finalChrMat, sampleName="MGH28", chrs=c("10p", "14q"))
plotSingleCellLargeScaleEventHeatmap(finalChrMat, sampleName="MGH28", chrs=c("5q", "19q"))
plotSingleCellLargeScaleEventHeatmap(finalChrMat, sampleName="MGH28", chrs=c("8q", "20p"))
plotSingleCellLargeScaleEventHeatmap(finalChrMat, sampleName="MGH26", chrs=c("5q", "7p"))
plotSingleCellLargeScaleEventHeatmap(finalChrMat, sampleName="MGH26", chrs=c("1q", "22q"))
## calculate significant mutual exclusive and co-occurent events
results <- extractMUAndCooccurence (finalChrMat, loh, loh.name.mapping)
## visualize mutual exclusive and co-occurent events
plotMUAndCooccurence (results)
## visualize CNV pyhlogenetic tree
## first please download pyhlip package from http://evolution.genetics.washington.edu/phylip/getme-new1.html
#plotSCellCNVTree (finalChrMat, sampleName="MGH31", path="C:\\Users\\aharmanci\\Downloads\\phylip-3.695\\phylip-3.695\\exe", fileName="CNVTree.pdf")
all.summary <- getDiffExprGenes(final.objects, sampleName="MGH31", chrs=c("5q", "14q"), event.type=c(1, -1))
#enrichment <- generateEnrichmentSummary (results=all.summary)
all.summary <- getDiffExprGenes(final.objects, sampleName="MGH31", chrs=c("1p", "13q"), event.type=c(1, -1))
all.summary <- getDiffExprGenes(final.objects, sampleName="MGH31", chrs=c("5q", "13q"), event.type=c(1, -1))
all.summary <- getDiffExprGenes(final.objects, sampleName="MGH28", chrs=c("8q", "20p"), event.type=c(1, -1))
all.summary <- getDiffExprGenes(final.objects, sampleName="MGH30", chrs=c("7p", "6p"), event.type=c(1, -1))
all.summary <- getDiffExprGenes(final.objects, sampleName="MGH26", chrs=c("1q", "22q"), event.type=c(1, -1))
genes <- as.character(all.summary$ID[all.summary$adj.P.Val<0.05])
all.summary <- getDiffExprGenes(final.objects, sampleName="MGH28", chrs=c("5q", "19q"), event.type=c(1, -1))
genes <- as.character(all.summary$ID[all.summary$adj.P.Val<0.05])
all.summary <- getDiffExprGenes(final.objects, sampleName="MGH29", chrs=c("10p", "4p"), event.type=c(-1, -1))
genes <- as.character(all.summary$ID[all.summary$adj.P.Val<0.05])