inst/doc/SCATE.R

In SCATE: SCATE: Single-cell ATAC-seq Signal Extraction and Enhancement

## ----setup, include=FALSE-----------------------------------------------------
knitr::opts_chunk$set(echo = TRUE)

## -----------------------------------------------------------------------------
# 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="SCATEData"),"/extdata/"),full.names = TRUE,pattern='.bam$')
bamlist <- bamlist[grepl('.bam$', bamlist)]
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[seq(1,3), seq(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=156,ncores=1,verbose=TRUE)
# check the 10000-10005th row of the matrix
res[seq(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

## -----------------------------------------------------------------------------
userclusterres <- SCATE(satac,genome="hg19",cluster=usercellcluster,clunum=156,ncores=1,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

## -----------------------------------------------------------------------------
extractres <- extractfeature(res,region,mode='overlap',folder='destination folder')

## -----------------------------------------------------------------------------
peakres <- peakcall(res)
# check the result for the first cluster
head(peakres[[1]])

## -----------------------------------------------------------------------------
write.table(peakres[[1]],file='your file.bed',sep='\t',quote=FALSE,col.names = FALSE,row.names = FALSE)

## -----------------------------------------------------------------------------
piperes <- SCATEpipeline(bamlist,genome="hg19",cellclunum=2,CREclunum=156,perplexity=5,ncores=1)
# 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']]

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

## -----------------------------------------------------------------------------
sessionInfo()