vignettes/MACPET.R

In IoannisVardaxis/MACPET: Model based analysis for paired-end data

## ----style,eval=TRUE,echo=FALSE,results='hide'-----------------------------
BiocStyle::latex

## ----eval=TRUE,echo=TRUE---------------------------------------------------
#Create a temporary test folder, or anywhere you want:
SA_AnalysisDir=file.path(tempdir(),"MACPETtest")
dir.create(SA_AnalysisDir)#where you will save the results.

## --------------------------------------------------------------------------
library(MACPET)

## --------------------------------------------------------------------------
load(system.file("extdata", "MACPET_pselfData.rda", package = "MACPET"))
class(MACPET_pselfData) #example name
MACPET_pselfData #print method

## --------------------------------------------------------------------------
metadata(MACPET_pselfData)

## --------------------------------------------------------------------------
seqinfo(MACPET_pselfData)

## --------------------------------------------------------------------------
load(system.file("extdata", "MACPET_psfitData.rda", package = "MACPET"))
class(MACPET_psfitData) #example name
MACPET_psfitData #print method

## --------------------------------------------------------------------------
head(metadata(MACPET_psfitData)$Peaks.Info)

## --------------------------------------------------------------------------
load(system.file("extdata", "MACPET_pinterData.rda", package = "MACPET"))
class(MACPET_pinterData) #example name
MACPET_pinterData #print method

## ----eval=TRUE-------------------------------------------------------------
metadata(MACPET_pinterData)

## --------------------------------------------------------------------------
load(system.file("extdata", "MACPET_pintraData.rda", package = "MACPET"))
class(MACPET_pintraData)#example name
MACPET_pintraData#print method

## --------------------------------------------------------------------------
metadata(MACPET_pintraData)

## --------------------------------------------------------------------------
load(system.file("extdata", "MACPET_GenomeMapData.rda", package = "MACPET"))
class(MACPET_GenomeMapData) #example name
MACPET_GenomeMapData #print method

## --------------------------------------------------------------------------
metadata(MACPET_GenomeMapData)

## --------------------------------------------------------------------------
class(MACPET_pselfData)
summary(MACPET_pselfData)

## --------------------------------------------------------------------------
class(MACPET_psfitData)
summary(MACPET_psfitData)

## --------------------------------------------------------------------------
class(MACPET_pintraData)
requireNamespace("ggplot2")
requireNamespace("reshape2")
summary(MACPET_pintraData,heatmap=TRUE)

## --------------------------------------------------------------------------
class(MACPET_pinterData)
requireNamespace("ggplot2")
requireNamespace("reshape2")
summary(MACPET_pinterData,heatmap=TRUE)

## --------------------------------------------------------------------------
class(MACPET_GenomeMapData)
summary(MACPET_GenomeMapData)

## --------------------------------------------------------------------------
requireNamespace("ggplot2")
class(MACPET_pselfData)
# PET counts plot
plot(MACPET_pselfData)

## --------------------------------------------------------------------------
class(MACPET_psfitData)
#binding site couts:
plot(MACPET_psfitData,kind="PeakCounts")

## --------------------------------------------------------------------------
# region example with binding sites:
plot(MACPET_psfitData,kind="PeakPETs",RegIndex=1)

## --------------------------------------------------------------------------
class(MACPET_pintraData)
#plot counts:
plot(MACPET_pintraData)

## --------------------------------------------------------------------------
class(MACPET_pinterData)
requireNamespace("igraph")
#network plot:
plot(MACPET_pinterData)

## --------------------------------------------------------------------------
class(MACPET_GenomeMapData)
requireNamespace("igraph")
#network plot:
plot(MACPET_GenomeMapData,Type='network-circle')

## ----eval=TRUE,echo=TRUE---------------------------------------------------
class(MACPET_psfitData)#PSFit class
exportPeaks(object=MACPET_psfitData,file.out="Peaks",threshold=1e-5,savedir=SA_AnalysisDir)

## ----eval=TRUE,echo=TRUE---------------------------------------------------
class(MACPET_psfitData)#PSFit class
object=PeaksToGRanges(object=MACPET_psfitData,threshold=1e-5)
object

## ----eval=TRUE,echo=TRUE---------------------------------------------------
class(MACPET_psfitData)#PSFit class
TagsToGInteractions(object=MACPET_psfitData,threshold=1e-5)


## ----eval=TRUE,echo=TRUE---------------------------------------------------
class(MACPET_psfitData)#PSFit class
PeaksToNarrowPeak(object=MACPET_psfitData,threshold=1e-5,
                  file.out="MACPET_peaks.narrowPeak",savedir=SA_AnalysisDir)

## ----eval=TRUE,echo=TRUE---------------------------------------------------
 #--remove information and convert to GInteractions:
object=MACPET_pselfData
#--remove information and convert to GInteractions:
S4Vectors::metadata(object)=list(NULL)
class(object)='GInteractions'
#----input parameters
S2_BlackList=TRUE
SA_prefix="MACPET"
S2_AnalysisDir=SA_AnalysisDir

ConvertToPSelf(object=object,
               S2_BlackList=S2_BlackList,
               SA_prefix=SA_prefix,
               S2_AnalysisDir=S2_AnalysisDir)
#load object:
rm(MACPET_pselfData)#old object
load(file.path(S2_AnalysisDir,"MACPET_pselfData"))
class(MACPET_pselfData)


## ----eval=TRUE,echo=TRUE---------------------------------------------------
class(MACPET_GenomeMapData)#GenomeMap class
GetSignInteractions(object=MACPET_GenomeMapData,
                     threshold = NULL,
                     ReturnedAs='GInteractions')

## ----eval=TRUE,echo=TRUE---------------------------------------------------
class(MACPET_GenomeMapData)#GenomeMap class
GetShortestPath(object=MACPET_GenomeMapData,
                     threshold = NULL,
                     ChrFrom="chr1",
                     ChrTo="chr1",
                     SummitFrom=10000,
                     SummitTo=1000000)

## ----echo=TRUE,eval=TRUE---------------------------------------------------
AnalysisStatistics(x.self=MACPET_psfitData,
                   x.intra=MACPET_pintraData,
                   x.inter=MACPET_pinterData,
                   file.out='AnalysisStats',
                   savedir=SA_AnalysisDir,
                   threshold=1e-5)


## ----echo=TRUE,eval=TRUE---------------------------------------------------
requireNamespace('ggplot2')

#Create a temporary forder, or anywhere you want:
S1_AnalysisDir=SA_AnalysisDir

#directories of the BAM files:
BAM_file_1=system.file('extdata', 'SampleChIAPETDataRead_1.bam', package = 'MACPET')
BAM_file_2=system.file('extdata', 'SampleChIAPETDataRead_2.bam', package = 'MACPET')
SA_prefix="MACPET"

#convert to paired-end BAM:
ConvertToPE_BAM(S1_AnalysisDir=S1_AnalysisDir,
             SA_prefix=SA_prefix,
             S1_BAMStream=2000000,S1_image=TRUE,
             S1_genome="hg19",BAM_file_1=BAM_file_1,
             BAM_file_2=BAM_file_2)

#test if the resulted BAM is paired-end:
PairedBAM=file.path(S1_AnalysisDir,paste(SA_prefix,"_Paired_end.bam",sep=""))
Rsamtools::testPairedEndBam(file = PairedBAM, index = PairedBAM)

bamfile = Rsamtools::BamFile(file = PairedBAM,asMates = TRUE)
GenomicAlignments::readGAlignmentPairs(file = bamfile,use.names = FALSE,
                                with.which_label = FALSE,
                                strandMode = 1)


## ----echo=TRUE,eval=TRUE---------------------------------------------------

#give directory of the BAM file:
S2_PairedEndBAMpath=system.file('extdata', 'SampleChIAPETData.bam', package = 'MACPET')

#give prefix name:
SA_prefix="MACPET"

#parallel backhead can be created using the BiocParallel package
#parallel backhead can be created using the BiocParallel package
#requireNamespace('BiocParallel')
#snow <- BiocParallel::SnowParam(workers = 4, type = 'SOCK', progressbar=FALSE)
#BiocParallel::register(snow, default=TRUE)

# packages for plotting:
requireNamespace('ggplot2')

#-run for the whole binding site analysis:
MACPETUlt(SA_AnalysisDir=SA_AnalysisDir,
       SA_stages=c(2:4),
       SA_prefix=SA_prefix,
       S2_PairedEndBAMpath=S2_PairedEndBAMpath,
       S2_image=TRUE,
       S2_BlackList=TRUE,
       S3_image=TRUE,
       S4_image=TRUE,
       S4_FDR_peak=1)# the data is small so use all the peaks found.

#load results:
SelfObject=paste(SA_prefix,"_pselfData",sep="")
load(file.path(SA_AnalysisDir,"S2_results",SelfObject))
SelfObject=get(SelfObject)
class(SelfObject) # see methods for this class

IntraObject=paste(SA_prefix,"_pintraData",sep="")
load(file.path(SA_AnalysisDir,"S2_results",IntraObject))
IntraObject=get(IntraObject)
class(IntraObject) # see methods for this class

InterObject=paste(SA_prefix,"_pinterData",sep="")
load(file.path(SA_AnalysisDir,"S2_results",InterObject))
InterObject=get(InterObject)
class(InterObject) # see methods for this class

SelfFitObject=paste(SA_prefix,"_psfitData",sep="")
load(file.path(SA_AnalysisDir,"S3_results",SelfFitObject))
SelfFitObject=get(SelfFitObject)
class(SelfFitObject) # see methods for this class

GenomeMapObject=paste(SA_prefix,"_GenomeMapData",sep="")
load(file.path(SA_AnalysisDir,"S4_results",GenomeMapObject))
GenomeMapObject=get(GenomeMapObject)
class(GenomeMapObject) # see methods for this class

#-----delete test directory:
unlink(SA_AnalysisDir,recursive=TRUE)