plot_EM_fit: Visualize EM fitting for each cell.

In rujinwang/SCOPE: A normalization and copy number estimation method for single-cell DNA sequencing

Visualize EM fitting for each cell.

Description

A pdf file containing EM fitting results and plots is generated.

Usage

plot_EM_fit(Y_qc, gc_qc, norm_index, T, ploidyInt, beta0,
                minCountQC = 20, filename)

Arguments

Y_qc	read depth matrix across all cells after quality control
gc_qc	vector of GC content for each bin after quality control
norm_index	indices of normal/diploid cells
T	a vector of integers indicating number of CNV groups. Use BIC to select optimal number of CNV groups. If T = 1, assume all reads are from normal regions so that EM algorithm is not implemented. Otherwise, we assume there is always a CNV group of heterozygous deletion and a group of null region. The rest groups are representative of different duplication states.
ploidyInt	a vector of initialized ploidy return from initialize_ploidy
beta0	a vector of initialized bin-specific biases returned from CODEX2 without latent factors
minCountQC	the minimum read coverage required for EM fitting. Defalut is 20
filename	the name of output pdf file

Value

pdf file with EM fitting results and two plots: log likelihood, and BIC versus the number of CNV groups.

Author(s)

Rujin Wang [e-mail]

Examples

Gini <- get_gini(Y_sim)
# first-pass CODEX2 run with no latent factors
normObj.sim <- normalize_codex2_ns_noK(Y_qc = Y_sim,
                                        gc_qc = ref_sim$gc,
                                        norm_index = which(Gini<=0.12))
Yhat.noK.sim <- normObj.sim$Yhat
beta.hat.noK.sim <- normObj.sim$beta.hat
fGC.hat.noK.sim <- normObj.sim$fGC.hat
N.sim <- normObj.sim$N

# Ploidy initialization
ploidy.sim <- initialize_ploidy(Y = Y_sim, 
                                Yhat = Yhat.noK.sim, 
                                ref = ref_sim)
ploidy.sim

plot_EM_fit(Y_qc = Y_sim, gc_qc = ref_sim$gc,
        norm_index = which(Gini<=0.12), T = 1:7,
        ploidyInt = ploidy.sim,
        beta0 = beta.hat.noK.sim,
        filename = 'plot_EM_fit_demo.pdf')

rujinwang/SCOPE documentation built on Jan. 1, 2023, 5:40 a.m.