## ----style-knitr, eval=TRUE, echo=FALSE, results="asis"-------------------- BiocStyle::latex() ## ----data_preparation1----------------------------------------------------- library(BUSseq) library(SingleCellExperiment) #Input data is should be a SingleCellExperiment object or a list CountData <- assay(BUSseqfits_example, "counts") batch_ind <- unlist(colData(BUSseqfits_example)) # Construct a SingleCellExperiment object with colData as batch indicators sce_input <- SingleCellExperiment(assays = list(counts = CountData), colData = DataFrame(Batch_ind = factor(batch_ind))) # Or, construct a list with each element represents the count data matrix # of a batch num_cell_batch <- table(batch_ind) list_input <- list(Batch1 = CountData[,1:num_cell_batch[1]], Batch2 = CountData[,1:num_cell_batch[2] + num_cell_batch[1]]) # Cell numbers within each batch print(num_cell_batch) #Peek at the read counts print(CountData[1:5,1:5]) ## ----BUSgibbs-------------------------------------------------------------- # Conduct MCMC sampling and posterior inference for BUSseq model BUSseqfits_res <- BUSseq_MCMC(ObservedData = sce_input, seed = 1234, n.cores = 2, n.celltypes = 4, n.iterations = 500) ## ----BUSseq_output--------------------------------------------------------- # The output is a SingleCellExperiment object class(BUSseqfits_res) # Peek at the output BUSseqfits_res ## ----Output_extract-------------------------------------------------------- # Extract the imputed read counts Imputed_count <- assay(BUSseqfits_res, "imputed_data") ## ----Celltypes------------------------------------------------------------- celltyes_est <- celltypes(BUSseqfits_res) ## ----BatchEffects---------------------------------------------------------- location_batch_effects_est <- location_batch_effects(BUSseqfits_res) head(location_batch_effects_est) overdispersion_est <- overdispersions(BUSseqfits_res) head(overdispersion_est) ## ----CellEffects----------------------------------------------------------- cell_effects_est <- cell_effect_values(BUSseqfits_res) head(cell_effects_est) ## ----CelltypeEffects------------------------------------------------------- celltype_mean_expression_est <- celltype_mean_expression(BUSseqfits_example) head(celltype_mean_expression_est) celltype_effects_est <- celltype_effects(BUSseqfits_res) head(celltype_effects_est) ## ----IG-------------------------------------------------------------------- #obtain the intrinsic gene indicators intrinsic_gene_indicators <- intrinsic_genes_BUSseq(BUSseqfits_res) print(intrinsic_gene_indicators) #The estimated FDR, the first 240 genes are known as intrinsic #genes in the simulation setting. index_intri <- which(unlist(intrinsic_gene_indicators) == "Yes") false_discovery_ind <- !(index_intri %in% 1:240) fdr_est <- sum(false_discovery_ind)/length(index_intri) print(fdr_est) ## ----adjusted-------------------------------------------------------------- # Obtain the corrected read count data BUSseqfits_res <- corrected_read_counts(BUSseqfits_res) # An new assay "corrected_data" is incorporated BUSseqfits_res ## ----visualize1------------------------------------------------------------ #generate the heatmap of raw read count data heatmap_data_BUSseq(BUSseqfits_res, data_type = "Raw", project_name = "BUSseq_raw_allgenes", image_dir = "./heatmap") #display only the first 100 genes in the heatmap heatmap_data_BUSseq(BUSseqfits_res, data_type = "Raw", gene_set = 1:100, project_name = "BUSseq_raw_100genes", image_dir = "./heatmap") ## ----visualize2------------------------------------------------------------ #generate the heatmap of imputed read count data heatmap_data_BUSseq(BUSseqfits_res, data_type = "Imputed", project_name = "BUSseq_imputed_allgenes", image_dir = "./heatmap") #generate the heatmap of corrected read count data heatmap_data_BUSseq(BUSseqfits_res, data_type = "Corrected", project_name = "BUSseq_corrected_allgenes", image_dir = "./heatmap") ## ----visualize3------------------------------------------------------------ #Only show the identified intrinsic genes heatmap_data_BUSseq(BUSseqfits_res, data_type = "Corrected", gene_set = index_intri, project_name = "BUSseq_corrected_intrinsic_genes", image_dir = "./heatmap") ## ----RSession-------------------------------------------------------------- sessionInfo()