mosaicsFitHMM {mosaics} | R Documentation

## Fit MOSAiCS-HMM model

### Description

Fit MOSAiCS-HMM model.

### Usage

    mosaicsFitHMM( object, ... )
    ## S4 method for signature 'MosaicsFit'
    mosaicsFitHMM( object, signalModel="2S", binsize=NA,
    	init="mosaics", init.FDR=0.05,
    	init.maxgap=200, init.minsize=50, init.thres=10, init.piMat=as.matrix(NA),
    	max.iter=100, eps=1e-20, parallel=FALSE, nCore=8 )

### Arguments

`object` | Object of class `MosaicsFit`, a fitted MOSAiCS model obtained using function `mosaicsFit`.
`signalModel` | Signal model. Possible values are "1S" (one-signal-component model) and "2S" (two-signal-component model). Default is "2S".
`binsize` | Size of each bin. Value should be positive integer. If `binsize=NA`, `mosaicsFitHMM` function calcuates the value from data. Default is `NA`.
`init` | Approach to initialize MOSAiCS-HMM. Possible values are `"mosaics"` (use MOSAiCS peak calling results for initialization) or `"specify"` (explicitly specify transition matrix). Default is `"mosaics"`.
`init.FDR` | Parameter for the MOSAiCS-HMM initialization. False discovery rate. Default is 0.05. Related only if `init="mosaics"`.
`init.maxgap` | Parameter for the MOSAiCS-HMM initialization. Initial nearby peaks are merged if the distance (in bp) between them is less than `init.maxgap`. Default is 200. Related only if `init="mosaics"`.
`init.minsize` | Parameter for the MOSAiCS-HMM initialization. An initial peak is removed if its width is narrower than `init.minsize`. Default is 50. Related only if `init="mosaics"`.
`init.thres` | Parameter for the MOSAiCS-HMM initialization. A bin within initial peak is removed if its ChIP tag counts are less than `init.thres`. Default is 10. Related only if `init="mosaics"`.
`init.piMat` | Initial value for transition matrix. The first rows/columns correspond to the non-binding state while the second rows/columns correspond to the binding state. Related only if `init="specify"`. If `init="specify"` but `init.piMat` is not specified, `mosaicsFitHMM()` uses its default for the MOSAiCS-HMM initialization.
`max.iter` | Number of iterations for fitting MOSAiCS-HMM. Default is 100.
`eps` | Criterion to stop iterations for fitting MOSAiCS-HMM. Default is 1e-20.
`parallel` | Utilize multiple CPUs for parallel computing using `"parallel"` package? Possible values are `TRUE` (utilize multiple CPUs) or `FALSE` (do not utilize multiple CPUs). Default is `FALSE` (do not utilize multiple CPUs).
`nCore` | Number of CPUs when parallel computing is utilized.
`...` | Other parameters to be passed through to generic `mosaicsFitHMM`.

### Details

`mosaicsFitHMM` and `mosaicsPeakHMM` are developed to identify broad peaks such as histone modifications, using Hidden Markov Model (HMM) approach, as proposed in Chung et al. (2014). If you are interested in identifying narrow peaks such as transcription factor binding sites, please use `mosaicsPeak` instead of `mosaicsFitHMM` and `mosaicsPeakHMM`.

When peaks are called, proper signal model needs to be specified. The optimal choice for the number of signal components depends on the characteristics of ChIP-seq data. In order to support users in the choice of optimal signal model, Bayesian Information Criterion (BIC) values and Goodness of Fit (GOF) plot are provided for the fitted MOSAiCS model. BIC values and GOF plot can be obtained by applying `show` and `plot` methods, respectively, to the `MosaicsFit` class object, which is a fitted MOSAiCS model.

`init.FDR`, `init.maxgap`, `init.minsize`, and `init.thres` are the parameters for MOSAiCS-HMM initialization when MOSAiCS peak calling results are used for initialization (`init="mosaics"`). If user specifies transition matrix (`init="specify"`), only `init.piMat` is used for initialization. If you use a bin size shorter than the average fragment length of the experiment, we recommend to set `init.maxgap` to the average fragment length and `init.minsize` to the bin size. If you set the bin size to the average fragment length or if bin size is larger than the average fragment length, set `init.maxgap` to the average fragment length and `init.minsize` to a value smaller than the average fragment length. See the vignette for further details.

Parallel computing can be utilized for faster computing if `parallel=TRUE` and `parallel` package is loaded. `nCore` determines number of CPUs used for parallel computing.

### Value

Construct `MosaicsHMM` class object.

### Author(s)

Dongjun Chung, Pei Fen Kuan, Rene Welch, Sunduz Keles

### References

Kuan, PF, D Chung, G Pan, JA Thomson, R Stewart, and S Keles (2011), "A Statistical Framework for the Analysis of ChIP-Seq Data", Journal of the American Statistical Association, Vol. 106, pp. 891-903.

Chung, D, Zhang Q, and Keles S (2014), "MOSAiCS-HMM: A model-based approach for detecting regions of histone modifications from ChIP-seq data", Datta S and Nettleton D (eds.), Statistical Analysis of Next Generation Sequencing Data, Springer.

### See Also

`mosaicsFit`, `mosaicsPeakHMM`, `MosaicsFit`, `MosaicsHMM`.

### Examples

    ## Not run:
    library(mosaicsExample)

    constructBins( infile=system.file( file.path("extdata","wgEncodeBroadHistoneGm12878H3k4me3StdAlnRep1_chr22_sorted.bam"), package="mosaicsExample"),
        fileFormat="bam", outfileLoc="./",
        byChr=FALSE, useChrfile=FALSE, chrfile=NULL, excludeChr=NULL,
        PET=FALSE, fragLen=200, binSize=200, capping=0 )
    constructBins( infile=system.file( file.path("extdata","wgEncodeBroadHistoneGm12878ControlStdAlnRep1_chr22_sorted.bam"), package="mosaicsExample"),
        fileFormat="bam", outfileLoc="./",
        byChr=FALSE, useChrfile=FALSE, chrfile=NULL, excludeChr=NULL,
        PET=FALSE, fragLen=200, binSize=200, capping=0 )

    binHM <- readBins( type=c("chip","input"),
        fileName=c( "./wgEncodeBroadHistoneGm12878H3k4me3StdAlnRep1_chr22_sorted.bam_fragL200_bin200.txt",
        "./wgEncodeBroadHistoneGm12878ControlStdAlnRep1_chr22_sorted.bam_fragL200_bin200.txt" ) )
    binHM
    plot(binHM)
    plot( binHM, plotType="input" )

    fitHM <- mosaicsFit( binHM, analysisType="IO", bgEst="rMOM" )
    fitHM
    plot(fitHM)

    hmmHM <- mosaicsFitHMM( fitHM, signalModel = "2S",
      init="mosaics", init.FDR = 0.05, parallel=TRUE, nCore=8 )
    hmmHM
    plot(hmmHM)

    peakHM <- mosaicsPeakHMM( hmmHM, FDR = 0.05, decoding="posterior",
      thres=10, parallel=TRUE, nCore=8 )

    peakHM <- extractReads( peakHM,
      chipFile=system.file( file.path("extdata","wgEncodeBroadHistoneGm12878H3k4me3StdAlnRep1_chr22_sorted.bam"), package="mosaicsExample"),
      chipFileFormat="bam", chipPET=FALSE, chipFragLen=200,
      controlFile=system.file( file.path("extdata","wgEncodeBroadHistoneGm12878ControlStdAlnRep1_chr22_sorted.bam"), package="mosaicsExample"),
      controlFileFormat="bam", controlPET=FALSE, controlFragLen=200, parallel=TRUE, nCore=8 )
    peakHM

    peakHM <- findSummit( peakHM, parallel=TRUE, nCore=8 )
    head(print(peakHM))
    plot( peakHM, filename="./peakplot_HM.pdf" )

    peakHM <- adjustBoundary( peakHM, parallel=TRUE, nCore=8 )
    peakHM
    head(print(peakHM))

    peakHM <- filterPeak( peakHM, parallel=TRUE, nCore=8 )
    peakHM
    head(print(peakHM))

    export( peakHM, type = "txt", filename = "./peakHM.txt" )
    export( peakHM, type = "bed", filename = "./peakHM.bed" )
    export( peakHM, type = "gff", filename = "./peakHM.gff" )
    export( peakHM, type = "narrowPeak", filename = "./peakHM.narrowPeak" )
    export( peakHM, type = "broadPeak", filename = "./peakHM.broadPeak" )

    ## End(Not run)

[Package mosaics version 2.18.0 Index]