R/base_lineStrandData.R

In carter-allen/dpeak: dPeak (Deconvolution of Peaks in ChIP-seq Analysis)

# plot fitting results

.lineStrandData <- function( stackedFragment, fragSet, adjustment, plot_title, xlim,
    PET, extension=1, smoothing=FALSE ) {

    xvar <- stackedFragment[,1]

    # forward strand

    if ( PET == FALSE ) {
        # SET

        indF <- which( fragSet[,3] == "F" )

        if ( length(indF) > 0 ) {
            Fstart <- fragSet[ indF, 1 ]
            Fend <- pmin( fragSet[ indF, 1 ] + extension - 1, max(xvar) )

            xvarq <- c( min(Fstart), max(Fend) )
            xvar_org <- c(seq(from = xvarq[1], to = xvarq[2]))
            yvar_org <- .ff_stack( Fstart, Fend, xvarq[1], xvarq[2] )

            yvarF <- rep( 0, length(xvar) )
            yvarF[ match( xvar_org, xvar ) ] <- yvar_org
        } else {
            yvarF <- c()
        }
    } else {
        # PET

        Fstart <- fragSet[,1]
        Fend <- pmin( fragSet[,1] + extension - 1, max(xvar) )

        xvarq <- c( min(Fstart), max(Fend) )
        xvar_org <- c(seq(from = xvarq[1], to = xvarq[2]))
        yvar_org <- .ff_stack( Fstart, Fend, xvarq[1], xvarq[2] )

        yvarF <- rep( 0, length(xvar) )
        yvarF[ match( xvar_org, xvar ) ] <- yvar_org
    }

    # reverse strand

    if ( PET == FALSE ) {
        # SET

        if ( nrow(fragSet) - length(indF) > 0 ) {
            if ( length(indF) > 0 ) {
                Rstart <- pmax( fragSet[ -indF, 2 ] - extension + 1, min(xvar) )
                Rend <- fragSet[ -indF, 2 ]
            } else {
                Rstart <- pmax( fragSet[ , 2 ] - extension + 1, min(xvar) )
                Rend <- fragSet[ , 2 ]
            }

            xvarq <- c( min(Rstart), max(Rend) )
            xvar_org <- c(seq(from = xvarq[1], to = xvarq[2]))
            yvar_org <- .ff_stack( Rstart, Rend, xvarq[1], xvarq[2] )

            yvarR <- rep( 0, length(xvar) )
            yvarR[ match( xvar_org, xvar ) ] <- yvar_org
        } else {
            yvarR <- c()
        }
    } else {
        # PET

        Rstart <- pmax( fragSet[,2] - extension + 1, min(xvar) )
        Rend <- fragSet[,2]

        xvarq <- c( min(Rstart), max(Rend) )
        xvar_org <- c(seq(from = xvarq[1], to = xvarq[2]))
        yvar_org <- .ff_stack( Rstart, Rend, xvarq[1], xvarq[2] )

        yvarR <- rep( 0, length(xvar) )
        yvarR[ match( xvar_org, xvar ) ] <- yvar_org
    }

    # plot

    if ( smoothing ) {

        if ( length(yvarF) > 0 ) {
            smoothF <- smooth.spline( xvar, yvarF * adjustment )
        }
        if ( length(yvarR) > 0 ) {
            smoothR <- smooth.spline( xvar, yvarR * adjustment )
        }

        if ( length(yvarF) > 0 & length(yvarR) > 0 ) {
            lines( smoothF, col="black", lty=2 )
            lines( smoothR, col="red", lty=2 )
        } else if ( length(yvarF) > 0 ) {
            lines( smoothF, col="black", lty=2 )
        } else if ( length(yvarR) > 0 ) {
            lines( smoothR, col="red", lty=2 )
        }
    } else {
        if ( length(yvarF) > 0 & length(yvarR) > 0 ) {
            lines( xvar, yvarF * adjustment, col="black", lty=2 )
            lines( xvar, yvarR * adjustment, col="red", lty=2 )
        } else if ( length(yvarF) > 0 ) {
            lines( xvar, yvarF * adjustment, col="black", lty=2 )
        } else if ( length(yvarR) > 0 ) {
            lines( xvar, yvarR * adjustment, col="red", lty=2 )
        }
    }
}