-
Notifications
You must be signed in to change notification settings - Fork 1
/
Copy pathSCATE.R
102 lines (81 loc) · 4.39 KB
/
SCATE.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
90
91
92
93
94
95
96
97
98
99
100
101
## -----------------------------------------------------------------------------
# load in SCATE
options(warn=-1)
suppressMessages(library(SCATE))
set.seed(12345)
# set up locations to bam files
bamlist <- list.files(paste0(system.file(package="SCATE"),"/extdata/example"),full.names = T)
head(bamlist)
## -----------------------------------------------------------------------------
satac <- satacprocess(input=bamlist,type='bam',libsizefilter=1000)
# Number of elements in satac
length(satac)
# Content in the first element as an example
satac[[1]]
## -----------------------------------------------------------------------------
names(satac) <- sub('.*/','',names(satac))
## -----------------------------------------------------------------------------
clusterres <- cellcluster(satac,genome="hg19",clunum=2,perplexity=5)
## -----------------------------------------------------------------------------
# tSNE results
tsne <- clusterres[[1]]
head(tsne)
# clustering results
cluster <- clusterres[[2]]
cluster
# cell 'GSM1596840.bam' belongs to cluster 1, and cell 'SRR1779746.bam' belongs to cluster 2.
# aggregated signal for CRE cluster
aggsig <- clusterres[[3]]
aggsig[1:3,1:3]
## -----------------------------------------------------------------------------
library(ggplot2)
plotdata <- data.frame(tSNE1=tsne[,1],tSNE2=tsne[,2],Cluster=as.factor(cluster))
ggplot(plotdata,aes(x=tSNE1,y=tSNE2,col=Cluster)) + geom_point()
## -----------------------------------------------------------------------------
res <- SCATE(satac,genome="hg19",cluster=cluster,clusterid=NULL,clunum=5000,ncores=10,verbose=TRUE)
# check the 10000-10005th row of the matrix
res[10000:10005,]
## -----------------------------------------------------------------------------
# use similar ways to construct the cluster
usercellcluster <- rep(1:2,each=9)
names(usercellcluster) <- names(satac)
# check the contents of the cluster
usercellcluster
## ----eval=FALSE---------------------------------------------------------------
# userclusterres <- SCATE(satac,genome="hg19",cluster=usercellcluster,clunum=5000,ncores=10,verbose=TRUE)
## -----------------------------------------------------------------------------
region <- data.frame(chr=c('chr5','chr5'),start=c(50000,50700),end=c(50300,51000))
region
## -----------------------------------------------------------------------------
extractres <- extractfeature(res,region,mode='overlap')
extractres
## ----eval=FALSE---------------------------------------------------------------
# extractres <- extractfeature(res,region,mode='overlap',folder='destination folder')
## -----------------------------------------------------------------------------
peakres <- peakcall(res)
# check the result for the first cluster
head(peakres[[1]])
## ----eval=FALSE---------------------------------------------------------------
# write.table(peakres[[1]],file='your file.bed',sep='\t',quote=F,col.names = F,row.names = F)
## ----eval=FALSE---------------------------------------------------------------
# piperes <- SCATEpipeline(bamlist,genome="hg19",cellclunum=2,CREclunum=5000,perplexity=5,ncores=10)
# # get the cell cluster results, same as calling 'cellcluster' function.
# cluster <- piperes[['cellcluster']]
# # get the SCATE outputs, same as calling 'SCATE' function.
# SCATEres <- piperes[['SCATE']]
# # get the peak calling results, same as calling 'peakcall' function.
# peakres <- piperes[['peak']]
## ----eval=FALSE---------------------------------------------------------------
# extractres <- extractfeature(piperes[['SCATE']],region,mode='overlap',folder='destination folder')
## ----eval=FALSE---------------------------------------------------------------
# makedatabase(datapath,savepath,bamfile=bamfile,cre=cre,genome='hg19')
## ----eval=FALSE---------------------------------------------------------------
# makedatabase(datapath=NULL,savepath,bamfile=bamfile,cre=cre,genomerange=genomerange)
## ----eval=FALSE---------------------------------------------------------------
# piperes <- SCATEpipeline(bamlist,datapath='path to new database')
## ----eval=FALSE---------------------------------------------------------------
# clusterres <- cellcluster(satac,datapath='path to new database',clunum=2,perplexity=5)
# cluster <- clusterres[[2]]
# res <- SCATE(satac,datapath='path to new database',cluster=cluster,clusterid=NULL)
## ----eval=FALSE---------------------------------------------------------------
# sessionInfo()