Commit ad25df78 authored by tekath's avatar tekath
Browse files

One to one mapping function. Plot prop bar overhaul.

parent 510a37ab

DESCRIPTION

+1 −1
Original line number Diff line number Diff line
@@ -10,6 +10,6 @@ License: GPL (>= 3) 
Encoding: UTF-8

LazyData: true

Depends: R (>= 3.4.0)

Imports: Matrix, tximport, sparseDRIMSeq, stageR, rtracklayer, formattable, DT, Gviz, assertthat

Imports: Matrix, tximport, sparseDRIMSeq, stageR, rtracklayer, formattable, DT, Gviz, assertthat, scales, ggplot2, reshape2, BiocParallel

RoxygenNote: 7.1.0

Roxygen: list(markdown = TRUE)

NAMESPACE

+2 −0
Original line number Diff line number Diff line
@@ -7,7 +7,9 @@ export(get_by_partition) 
export(import_counts)

export(import_gtf)

export(move_columns_to_front)

export(one_to_one_mapping)

export(plot_dtu_table)

export(plot_prop_bar)

export(posthoc_and_stager)

export(rm_version)

export(run_drimseq)

R/dtu_analysis.R

+24 −18
Original line number Diff line number Diff line
@@ -268,7 +268,9 @@ run_drimseq <- function(counts, tx2gene, pd, id_col=NULL, cond_col, cond_levels= 
    tx2gene <- rapply(tx2gene, as.character, classes="factor", how="replace")

    tx2gene <- tx2gene[match(rownames(counts), tx2gene[[1]]),]

    assertthat::assert_that(nrow(tx2gene)==nrow(counts))

    message("Using tx2gene columns:\n\t",colnames(tx2gene)[[1]]," ---> 'feature_id'\n\t",colnames(tx2gene)[[2]]," ---> 'gene_id'")

    colnames(tx2gene)[c(1,2)] <- c("feature_id", "gene_id")

    colnames(tx2gene) <- make.names(colnames(tx2gene), unique = T)



    if(is.null(cond_levels)){

        if(length(unique(pd[[cond_col]]))==2){
@@ -279,7 +281,7 @@ run_drimseq <- function(counts, tx2gene, pd, id_col=NULL, cond_col, cond_levels= 
        }

    }

    assertthat::assert_that(length(cond_levels)==2)

    message("Comparing ", cond_levels[1], " vs ", cond_levels[2])

    message("\nComparing ", cond_levels[1], " vs ", cond_levels[2])



    if(is.null(id_col)){

        samp <- data.frame("sample_id"=row.names(pd), "condition"=pd[[cond_col]], stringsAsFactors = F)
@@ -294,7 +296,7 @@ run_drimseq <- function(counts, tx2gene, pd, id_col=NULL, cond_col, cond_levels= 
        samp <- samp[!is.na(samp$condition),]

        counts <- counts[ , !(colnames(counts) %in% exclude)]

    }

    message("Proceed with cells/samples: ",paste0(capture.output(table(samp$condition)), collapse = "\n"))

    message("\nProceed with cells/samples: ",paste0(capture.output(table(samp$condition)), collapse = "\n"))



    #TODO: Reevaluate!

    message("\nFiltering...\n")
@@ -350,14 +352,14 @@ run_drimseq <- function(counts, tx2gene, pd, id_col=NULL, cond_col, cond_levels= 
    tictoc::tic("Metadata")

    #carry over metadata

    if(carry_over_metadata&ncol(tx2gene)>2){

        temp <- tx2gene[match(rownames(exp_in_tx), tx2gene$feature_id),]

        temp <- temp[,!apply(temp,2,function(x) any(is.na(x)))]

        if(nrow(temp)==nrow(exp_in_tx)&ncol(temp)>2){

            exp_in_tx <- cbind(exp_in_tx, temp[,-c(1,2)])

            add_to_gene_columns <- check_unique_by_partition(temp[,-c(1,2)], drim@counts@partitioning)

        tx2gene <- tx2gene[match(rownames(exp_in_tx), tx2gene$feature_id),]

        tx2gene <- tx2gene[,!apply(tx2gene,2,function(x) any(is.na(x)))]

        if(nrow(tx2gene)==nrow(exp_in_tx)&ncol(tx2gene)>2){

            exp_in_tx <- cbind(exp_in_tx, tx2gene[,-c(1,2)], stringsAsFactors = F)

            add_to_gene_columns <- check_unique_by_partition(tx2gene[,-c(1,2)], drim@counts@partitioning)

            if(!is.null(add_to_gene_columns)){

                temp <- get_by_partition(df = temp, columns =add_to_gene_columns, partitioning = drim@counts@partitioning, FUN=unique, BPPARAM = BPPARAM)

                exp_in_gn <- cbind(exp_in_gn, temp[match(rownames(exp_in_gn), temp[[1]]),-c(1)])

                tx2gene <- get_by_partition(df = tx2gene, columns =add_to_gene_columns, partitioning = drim@counts@partitioning, FUN=unique, BPPARAM = BPPARAM)

                exp_in_gn <- cbind(exp_in_gn, tx2gene[match(rownames(exp_in_gn), tx2gene[[1]]),-c(1)], stringsAsFactors = F)

            }

        }

    }
@@ -378,6 +380,10 @@ run_drimseq <- function(counts, tx2gene, pd, id_col=NULL, cond_col, cond_levels= 


    group <- factor(samp$condition, levels = cond_levels, ordered = T)

    tictoc::toc(log = T)



    exp_in_gn <- rapply(exp_in_gn, as.character, classes="factor", how="replace")

    exp_in_tx <- rapply(exp_in_tx, as.character, classes="factor", how="replace")



    return_obj <- list("meta_table_gene"=exp_in_gn, "meta_table_tx"=exp_in_tx,

                       "drim"=drim_test, "design_full"=design_full, "group"=group,

                       "used_filtering_options"=list("DRIM"=filter_opt_list), "tictoc"=tictoc::tic.log(format = T))
@@ -402,8 +408,8 @@ run_drimseq <- function(counts, tx2gene, pd, id_col=NULL, cond_col, cond_levels= 
#' @param ofdr Overall false discovery rate (OFDR) threshold.

#' @param posthoc_filt Specify the minimal standard deviation of a transcripts porportion level that should be kept when performing posthoc filtering. To disbale poshoc filtering 0 or `FALSE` can be provided.

#' @return An extended `dturtle` object. Additional slots include:

#' - `DE_gene`: A character vector of all genes where the first stageR step was significant. Basically the significant genes, that showed signs of DTU.

#' - `DE_tx` : A named character vector of transcripts where the second stageR step was significant. Basically the significant transcripts of the significant genes.

#' - `sig_gene`: A character vector of all genes where the first stageR step was significant. Basically the significant genes, that showed signs of DTU.

#' - `sig_tx` : A named character vector of transcripts where the second stageR step was significant. Basically the significant transcripts of the significant genes.

#' - `FDR_table` : A data frame of the stage-wise adjusted p-values for all genes/transcripts. Might contain NA-values, as transcript level p-values are not avaible when the gene level test was not significant.

#'

#' @family DTUrtle
@@ -438,18 +444,18 @@ posthoc_and_stager <- function(dturtle, ofdr=0.05, posthoc=0.1){ 
    stageRObj <- stageR::stageRTx(pScreen = pscreen, pConfirmation = pconfirm, pScreenAdjusted = F, tx2gene = tx2gene)

    stageRObj <- stageR::stageWiseAdjustment(stageRObj, method = "dtu", alpha = ofdr)

    fdr_table <- stageR::getAdjustedPValues(stageRObj, order = F, onlySignificantGenes = F)

    DE_gene <- unique(as.character(fdr_table$geneID[fdr_table$gene<ofdr]))

    sig_gene <- unique(as.character(fdr_table$geneID[fdr_table$gene<ofdr]))



    if(length(DE_gene)>0){

    if(length(sig_gene)>0){

        temp <- fdr_table[fdr_table$gene<ofdr&fdr_table$transcript<ofdr,]

        DE_tx <- setNames(as.character(temp$txID), temp$geneID)

        message("Found ",length(DE_gene)," significant genes with ",length(DE_tx)," significant transcripts (OFDR: ",ofdr,")")

        sig_tx <- setNames(as.character(temp$txID), temp$geneID)

        message("Found ",length(sig_gene)," significant genes with ",length(sig_tx)," significant transcripts (OFDR: ",ofdr,")")

    }else{

        DE_gene <- NULL

        DE_tx <- NULL

        sig_gene <- NULL

        sig_tx <- NULL

        message("No gene passed the screening test. If applicable try to adjust the OFDR level.")

    }

    return_obj <- append(list("DE_gene" = DE_gene, "DE_tx" = DE_tx,

    return_obj <- append(list("sig_gene" = sig_gene, "sig_tx" = sig_tx,

                                       "FDR_table" = fdr_table), dturtle)

    return_obj$used_filtering_options$posthoc_stager <- list("ofdr" = ofdr,

                                       "posthoc"=ifelse(posthoc==F, 0, posthoc))

R/utility_functions.R

+43 −9
Original line number Diff line number Diff line
@@ -241,11 +241,11 @@ check_unique_by_partition <- function(df, partitioning, columns=NULL){ 
    assertthat::assert_that(is.list(partitioning))

    if(!is.null(columns)){

        assertthat::assert_that(all(columns %in% colnames(df)))

        df <- df[,columns]

        df <- df[,columns, drop=F]

    }

    cols <- colnames(df)

    for(part in partitioning){

            dat <- df[part,cols]

        dat <- df[part, cols, drop=F]

        cols <- cols[apply(dat, 2, function(x){all(x == x[1])})]

        if(length(cols)==0){

            return(NULL)
@@ -283,3 +283,37 @@ get_by_partition <- function(df, partitioning, FUN, columns=NULL, simplify=T, dr 
    return(BiocParallel::bpaggregate(df, by=list(rep(names(partitioning), lengths(partitioning))), FUN=FUN,  simplify=simplify, drop=T, BPPARAM = BPPARAM))

}





#' Ensure one-to-one mapping

#'

#' Ensure one-to-one mapping of the two specified columns in the data frame.

#'

#' First checks if every unique value in `name` corresponds with a unique value in `id`.

#' If not, changes the disagreeing values in `name` by extending the label with the `ext` character and a number.

#'

#' @param df A data frame, containing the two provided columns `name` and `id`

#' @param name A column of `df`. If no one-to-one mapping exists, the values of this column will be changed (by extending with `ext`)!

#' @param id A column of df. This is the preferred place for identifier columns, as this column will not be touched in case of a disagreement.

#' @param ext The extension character.

#'

#' @return A data frame, where one to one mapping for the two columns is ensured.

#' @export

one_to_one_mapping <- function(df, name, id, ext="_"){

    assertthat::assert_that(is.data.frame(df))

    assertthat::assert_that(all(c(name,id) %in% colnames(df)))

    assertthat::assert_that(is(ext, "character"))



    not_correct <- lapply(split(df[[id]], df[[name]]), unique)

    not_correct <- not_correct[lengths(not_correct)!=1]



    if(length(not_correct)>0){

        lapply(not_correct, function(x)

            lapply(seq(from = 2, along.with = x[-1]), function(i)

                df[[name]][df[[id]]==x[[i]]] <<- paste0(df[[name]][df[[id]]==x[[i]]],ext,i)))

        message("Changed ", sum(lengths(not_correct))-length(not_correct), " names.")

        return(df)

    }else{

        message("No changes needed - already one to one mapping.")

        return(df)

    }

}

R/visualization.R

+159 −107
Original line number Diff line number Diff line
@@ -20,22 +20,22 @@ 
#' @examples

create_dtu_table <- function(dturtle, add_gene_metadata = list("pct_gene_expr"="exp_in"), add_tx_metadata = list("max_pct_tx_expr"=c("exp_in", max))){

  assertthat::assert_that(is(dturtle,"dturtle"), msg = "The provided dturtle object is not of class 'dturtle'.")

    assertthat::assert_that(!is.null(dturtle$DE_gene), msg = "The provided dturtle object does not contain all the needed information. Have you run 'posthoc_and_stager()'?")

    assertthat::assert_that(!is.null(dturtle$DE_tx), msg = "The provided dturtle object does not contain all the needed information. Have you run 'posthoc_and_stager()'?")

  assertthat::assert_that(!is.null(dturtle$sig_gene), msg = "The provided dturtle object does not contain all the needed information. Have you run 'posthoc_and_stager()'?")

  assertthat::assert_that(!is.null(dturtle$sig_tx), msg = "The provided dturtle object does not contain all the needed information. Have you run 'posthoc_and_stager()'?")

  assertthat::assert_that(!is.null(dturtle$FDR_table), msg = "The provided dturtle object does not contain all the needed information. Have you run 'posthoc_and_stager()'?")

  assertthat::assert_that(!is.null(dturtle$group), msg = "The provided dturtle object does not contain all the needed information. Have you run 'posthoc_and_stager()'?")

  assertthat::assert_that(!is.null(dturtle$drim), msg = "The provided dturtle object does not contain all the needed information. Have you run 'posthoc_and_stager()'?")

    assertthat::assert_that(length(dturtle$DE_gene)>0, msg = "The provided dturtle object does not contain any significant gene. Maybe try to rerun the pipeline with more relaxes thresholds.")

    assertthat::assert_that(is.null(add_gene_metadata)|(is(add_gene_metadata, "list")&&all(lengths(add_tx_metadata)>0)), msg = "The add_gene_metadata object must be a list of non-empty elements or NULL.")

    assertthat::assert_that(is.null(add_tx_metadata)|(is(add_tx_metadata, "list")&&all(lengths(add_tx_metadata)>0)), msg = "The add_tx_metadata object must be a list of non-empty elements or or NULL.")

  assertthat::assert_that(length(dturtle$sig_gene)>0, msg = "The provided dturtle object does not contain any significant gene. Maybe try to rerun the pipeline with more relaxes thresholds.")

  assertthat::assert_that(is.null(add_gene_metadata)||(is(add_gene_metadata, "list")&&all(lengths(add_tx_metadata)>0)), msg = "The add_gene_metadata object must be a list of non-empty elements or NULL.")

  assertthat::assert_that(is.null(add_tx_metadata)||(is(add_tx_metadata, "list")&&all(lengths(add_tx_metadata)>0)), msg = "The add_tx_metadata object must be a list of non-empty elements or or NULL.")



  max_delta_col <- paste0("max(",levels(dturtle$group)[1], "-",levels(dturtle$group)[2],")")

    dtu_table <- data.frame("geneID" = dturtle$DE_gene, stringsAsFactors = F)

  dtu_table <- data.frame("geneID" = dturtle$sig_gene, stringsAsFactors = F)



  dtu_table$gene_qval <- sapply(dtu_table$geneID, FUN = function(x) min(dturtle$FDR_table$gene[dturtle$FDR_table$geneID == x]))

  dtu_table$min_tx_qval <- sapply(dtu_table$geneID, FUN = function(x) min(dturtle$FDR_table$transcript[dturtle$FDR_table$geneID == x]))

  dtu_table$n_tx <- sapply(dtu_table$geneID, FUN = function(x) length(dturtle$FDR_table$geneID[dturtle$FDR_table$geneID == x]))

    dtu_table$n_sig_tx <- sapply(dtu_table$geneID, FUN = function(x) length(dturtle$DE_tx[names(dturtle$DE_tx) == x]))

  dtu_table$n_sig_tx <- sapply(dtu_table$geneID, FUN = function(x) length(dturtle$sig_tx[names(dturtle$sig_tx) == x]))

  dtu_table[[max_delta_col]] <- as.numeric(mapply(dtu_table$geneID, FUN = getmax, MoreArgs = list(dturtle = dturtle)))

  # dtu_table$pct_expr_gene <- dturtle$meta_table_gene$pct_exp[match(dtu_table$geneID, dturtle$pct_exp_gene$gene)]

  # dtu_table$max_pct_expr_tx <- sapply(dtu_table$geneID, FUN = function(x) max(dturtle$pct_exp_tx$pct_exp[dturtle$pct_exp_tx$gene == x]))
@@ -340,27 +340,75 @@ plot_dtu_table <- function(dtu, txdf, title, folder, cores=1, image_folder="imag 






#TODO: group parameter? different shape (44?) as current one does not scale

plotProp_bar <- function(dtu, txdf=txdf, gene_id, path){

#' Title

#'

#' @param dturtle Result object of `posthoc_and_stager()`. Must be of class `dturtle`.

#' @param genes Character vector of genes to plot. If `NULL`, defaults to all found significant genes (`sig_genes`).

#' @param meta_gene_id

#' @param group_colors

#' @param dash_color

#' @param savepath

#'

#' @return

#' @family DTUrtle visualization

#' @export

#'

#' @examples

plot_prop_bar <- function(dturtle, genes=NULL, meta_gene_id=NULL, group_colors=NULL, dash_color="red", savepath=NULL){

  assertthat::assert_that(is(dturtle,"dturtle"), msg = "The provided dturtle object is not of class 'dturtle'.")

  assertthat::assert_that(is.null(genes)||(is(genes,"character")&&length(genes)>0), msg = "The genes object must be a non-empty character vector or NULL.")

  assertthat::assert_that(is.null(meta_gene_id)||(is(meta_gene_id, "character")&&meta_gene_id %in% colnames(dturtle$meta_table_gene)), msg = "The provided meta_gene_id column could not be found or is of wrong format.")

  assertthat::assert_that(is.null(group_colors)||(is(group_colors, "list")&&!is.null(names(group_colors))), msg = "The provided group colors must be a named list or NULL.")

  assertthat::assert_that(is.null(dash_color)||is(dash_color,"character"), msg = "The provided dash_color must be of type character or NULL.")

  assertthat::assert_that(is.null(savepath)||is(savepath,"character"), msg = "The provided savepath must be of type character or NULL.")

  assertthat::assert_that(!is.null(dturtle$sig_gene), msg = "The provided dturtle object does not contain all the needed information. Have you run 'posthoc_and_stager()'?")

  assertthat::assert_that(length(dturtle$sig_gene)>0, msg = "The provided dturtle object does not contain any significant gene. Maybe try to rerun the pipeline with more relaxes thresholds.")

  assertthat::assert_that(!is.null(dturtle$meta_table_gene), msg = "The provided dturtle object does not contain all the needed information. Have you run 'posthoc_and_stager()'?")

  assertthat::assert_that(!is.null(dturtle$drim), msg = "The provided dturtle object does not contain all the needed information. Have you run 'posthoc_and_stager()'?")

  assertthat::assert_that(!is.null(dturtle$sig_tx), msg = "The provided dturtle object does not contain all the needed information. Have you run 'posthoc_and_stager()'?")

  assertthat::assert_that(!is.null(dturtle$group), msg = "The provided dturtle object does not contain all the needed information. Have you run 'posthoc_and_stager()'?")



  if(is.null(genes)){

    genes <- dturtle$sig_gene

  }



  valid_genes <- genes[genes %in% dturtle$drim@results_gene$gene_id]



  if(length(genes)!=length(valid_genes)){

    message("Removed ", length(genes)-length(valid_genes), " genes, which where not present in the drimseq analysis.")

  }



  if(!is.null(meta_gene_id)){

    gene_ids <- dturtle$meta_table_gene[match(valid_genes, dturtle$meta_table_gene$gene),meta_gene_id]

    names(gene_ids) <- valid_genes

  }



  plot_list <- lapply(valid_genes, function(gene){

    counts_gene <- as.matrix(dturtle$drim@counts[[gene]])

    group <- dturtle$group



    gene_id <- as.character(gene_id)

    group

    counts_gene <- dtu$drim@counts[[gene_id]]

    sum1 <- colSums(counts_gene[,which(group==levels(group)[1])])

    sum2 <- colSums(counts_gene[,which(group==levels(group)[2])])

    mean_1 <- mean(sum1, na.rm = T)

    mean_2 <- mean(sum2, na.rm = T)

    sd_1 <- sd(sum1, na.rm = T)

    sd_2 <- sd(sum2, na.rm = T)



    gene_name <- txdf$external_gene_name[txdf$GENEID==gene_id][1]

    mean_1 <- mean(colSums(counts_gene[,which(group==levels(group)[1])]), na.rm = TRUE)

    mean_2 <- mean(colSums(counts_gene[,which(group==levels(group)[2])]), na.rm = TRUE)

    sd_1 <- sd(colSums(counts_gene[,which(group==levels(group)[1])]), na.rm = TRUE)

    sd_2 <- sd(colSums(counts_gene[,which(group==levels(group)[2])]), na.rm = TRUE)

    main <- paste0(gene_name," (",gene_id,")",

    if(!is.null(meta_gene_id)){

      gene_id <- gene_ids[[gene]]

      main <- paste0(gene," (",gene_id,")")

    }else{

      main <- gene

    }

    main <- paste0(main,

                   "\n\tMean expression (CV) ", levels(group)[1], " = ", round(mean_1), " (" ,round((sd_1/mean_1)*100, digits=1), "%)",

                   "\n\tMean expression (CV) ", levels(group)[2], " = ", round(mean_2), " (" ,round((sd_2/mean_2)*100, digits=1), "%)")



    fit_full <- dtu$drim@fit_full[[gene_id]]

    md <- dtu$drim@samples

    fit_full <- dturtle$drim@fit_full[[gene]]

    md <- dturtle$drim@samples



    proportions <- prop.table(counts_gene, 2)

    proportions[proportions == "NaN"] <- NA

    proportions <- sweep(counts_gene, 2, colSums(counts_gene), "/")

    proportions[is.nan(proportions)] <- NA

    prop_samp <- data.frame(feature_id = rownames(proportions), proportions, stringsAsFactors = F, check.names = F)

    prop_fit <- data.frame(feature_id = rownames(fit_full), fit_full, stringsAsFactors = F, check.names = F)
@@ -372,11 +420,6 @@ plotProp_bar <- function(dtu, txdf=txdf, gene_id, path){ 
    o <- order(group, -prop_samp[feature_levels[1],-1])

    sample_levels <- colnames(counts_gene)[o]



    #txname labels with transcript support level

    labels <- txdf$external_transcript_name[match(feature_levels, txdf$TXNAME)]

    tsl <- gsub("tsl", "", txdf$transcript_tsl[match(feature_levels, txdf$TXNAME)])

    x_label <- paste0(labels, rep(" (", length(labels)), tsl, rep(")", length(labels)))



    prop_samp <- reshape2::melt(prop_samp, id.vars = "feature_id", variable.name = "sample_id",

                                value.name = "proportion", factorsAsStrings = FALSE)

    prop_samp$feature_id <- factor(prop_samp$feature_id, levels = feature_levels)
@@ -400,26 +443,35 @@ plotProp_bar <- function(dtu, txdf=txdf, gene_id, path){ 
    }



    #colours

    if(is.null(group_colors)){

      group_colors <- scales::hue_pal()(nlevels(group))

      names(group_colors) <- levels(group)

    text_colour <- ifelse(feature_levels %in% dtu$final_q_tx$txID, "red", "dimgrey")

    }

    text_colour <- ifelse(feature_levels %in% dturtle$sig_tx, "red", "dimgrey")



    #barplot

    library(ggplot2)

    ggp <- ggplot() +

        geom_bar(data = prop_samp,

                 aes_string(x = "feature_id", y = "proportion", group = "sample_id", fill = "group"),

                 stat = "identity", position = position_dodge(width = 0.9)) +

        theme_bw() +

        theme(axis.text.x = element_text(angle = 25, vjust = 1, hjust=1, colour = text_colour), axis.text = element_text(size = 12),

              axis.title = element_text(size = 12, face = "bold"), plot.title = element_text(size = 12),

              legend.position = "right", legend.title = element_text(size = 12), legend.text = element_text(size = 12)) +

        scale_fill_manual(name = "Groups", values = group_colors, breaks = names(group_colors)) +

        ggtitle(main) + xlab("Features") + ylab("Proportions") +

        geom_point(data = prop_fit, aes_string(x = "feature_id", y = "proportion", group = "sample_id"),

                   position = position_dodge(width = 0.9), size = 3, shape = 46, alpha = 0.2, colour = "red") +

        scale_x_discrete(labels = x_label)

    ggplot2::ggsave(paste0(path,make.names(gene_name),"_bar.png"),ggp, dpi=250, scale=1.5)

    ggp <- ggplot2::ggplot() +

      ggplot2::geom_bar(data = prop_samp, ggplot2::aes_string(x = "feature_id", y = "proportion", group = "sample_id", fill = "group"),

               stat = "identity", position = ggplot2::position_dodge(width = 0.9)) + ggplot2::theme_bw() +

      ggplot2::theme(axis.text.x = ggplot2::element_text(angle = 25, vjust = 1, hjust=1, colour = text_colour), axis.text = ggplot2::element_text(size = 12),

                     axis.title = ggplot2::element_text(size = 12, face = "bold"), plot.title = ggplot2::element_text(size = 12),

                     legend.position = "right", legend.title = ggplot2::element_text(size = 12), legend.text = ggplot2::element_text(size = 12)) +

      ggplot2::scale_fill_manual(name = "Groups", values = group_colors, breaks = names(group_colors)) +

      ggplot2::labs(title = main, x = "Features", y = "Proportions")



    if(!is.null(dash_color)){

      ggp <- ggp + ggplot2::geom_errorbar(data = prop_fit, ggplot2::aes_string(x = "feature_id", ymin = "proportion", ymax = "proportion", group = "sample_id"),

                                 position = ggplot2::position_dodge(width = 0.9), size=0.5, linetype = "dashed", inherit.aes = F, width = 1, colour = dash_color)

    }

    if(!is.null(savepath)){

      if(!dir.exists(savepath)){

        dir.create(savepath)

      }

      ggplot2::ggsave(filename = paste0(savepath,make.names(gene),"_bar.png"), plot = ggp, ...)

    }

    return(setNames(list(ggp),gene))

  })

  return(plot_list)

}



plot_heat_per_gene <- function(gene_id, txdf, all_counts, mut, mut_genes, path){

CRA Git | Maintained and supported by SUSTech CRA and CCSE