#' Import the quantification results for DGE analysis of many RNA-seq quantifiers, including `alevin` and `bustools` for single-cell data.
#' Most likely the first step in your DTUrtle DGE analysis.
#'
#' It is necessary to specify a `tx2gene` data frame as a parameter.
#' This data frame must be a a two-column data frame linking transcript id (column 1) to gene id/name (column 2).
#' Please see [import_gtf()], [move_columns_to_front()] and [one_to_one_mapping()] to help with tx2gene creation.
#' See also [combine_to_matrix()], when output is a list of single-cell runs.
#'
#' @param files Vector of files to be imported. Optionally can be named to keep the samples names.
#' @param type Type of the quantification data. All tools supported by \code{\link[tximport:tximport]{tximport}} can be selected, additionally to the newly implemented `bustools` support for single-cell data. If you have single-cell data, the use of `alevin` or `bustools` is proposed.
#' @return - For bulk data: A list containing a count matrix, a matrix of average effective transcript lengths and a flag how counts where inferred from abundance estimates.
#' - For single-cell data: A list of count matrices per sample. Should be combined and optionally added to a Seurat object with [combine_to_matrix()].
#' @family DTUrtle DGE
#' @export
#' @seealso Please see [import_gtf()], [move_columns_to_front()] and [one_to_one_mapping()] to help with tx2gene creation. See also [combine_to_matrix()], when output is a list of single-cell runs.
message("Reading in ",length(files)," ",type," runs.")
args=list(...)
if(methods::hasArg("countsFromAbundance")){
if(args$countsFromAbundance!="no"){
warning("DESeq2 expects 'un-normalized counts', as a correction for library size is performed by DESeq2 itself. Changing the default 'countsFromAbundance' value is therefore strongly discouraged.")
}
}
if(methods::hasArg("txOut")){
if(args$txOut==T){
warning("Unless you exactly know what you are doing, it is not recommended to set txOut to TRUE\nDownstream analysis may fail!")
}
}else{
assertthat::assert_that(methods::hasArg("tx2gene"),msg="Please provide a 'tx2gene' data frame for gene-level summarisation.")
args$txOut<-F
}
if(type=="alevin"||type=="bustools"){
return_obj=list()
if(type=="alevin"){
assertthat::assert_that(all(basename(files)=="quants_mat.gz"),msg="Expecting 'files' to point to 'quants_mat.gz' file in a directory 'alevin'\n also containing 'quants_mat_rows.txt' and 'quant_mat_cols.txt'.\n Please re-run alevin preserving output structure.")
#' Perform differential gene expression analysis with DESeq2
#'
#' Offers functionality to perform a DGE analysis for bulk and single-cell data with DESeq2, automatically applying recommended models and parameter settings.
#' It is strongly advised to provide 'raw' count data, as imported with [import_dge_counts()].
#' Installation of package 'apeglm' is recommended for LFC-shrinkage, for single-cell data the package 'glmGamPoi' is additionally recommended.
#' For questions about DESeq2, LFC-shrinkage or s-values, please refer to the excellent [DESeq2 vignette](http://bioconductor.org/packages/release/bioc/vignettes/DESeq2/inst/doc/DESeq2.html).
#'
#' @param counts Can be either:
#' 1. (sparse) matrix with gene counts, where the rows correspond to genes. One column per sample/cell with the count data for the specified genes must be present (The names of these columns must match with the identifiers in `id_col`).
#' 2. An results object from tximport (preferred if [import_dge_counts()] was used with bulk data).
#' 3. Seurat object with a gene level assay as `active.assay` (most likely result object from [combine_to_matrix()])
#' @param lfc_threshold Specify a log2 fold change threshold (on log2 scale) to test against. 0 implicates no threshold.
#' @param sig_threshold Specify a significance threshold for the results adjusted p-values or s-values. 1 implicates no threshold.
#' @param dge_calling_strategy Should be either 'bulk' or 'sc'. Specify the type of the provided data (bulk or single-cell RNA-seq) so that appropriate parameters can be applied.
#' @param deseq_opts Manually specify parameters for the \code{\link[DESeq2:DESeq]{DESeq}} function. Will overwrite recommended parameters if necessary.
#' @param lfc_shrink_opts Manually specify parameters for the \code{\link[DESeq2:lfcShrink]{lfcShrink}} function. Will overwrite recommended parameters if necessary.
#' @param return_dds Should the DESeqDataSet object be returned?
#' @inheritParams run_drimseq
#'
#' @return A list with the analysis results and parameters:
#' - `results_all`: Data frame of the DGE test results for all analyzed genes.
#' - `results_sig`: Data frame of the significant DEG test results, according to the specified parameters (`sig_threshold`, `lfc_threshold`).
#' - `dds`: The DESeqDataSet of the analysis, if `return_dds=TRUE`.
#' - `drim`: Results of the DRIMSeq statistical computations (`dmTest()`).
#' - `sval_threshold / adjp_threshold`: The given significance threshold used to either filter s-values or adjusted p-values.
#' - `comparison`: A string representation of the performed comparison.
#' - `condition1`: The first condition of the performed comparison.
#' - `condition2`: The second condition of the performed comparison.
#' - `sample_table`: The (filtered) sample table (`pd`) - including the `condition` column used for comparison.
#' - `deseq_opts`: A list of used \code{\link[DESeq2:DESeq]{DESeq}} parameters.
#' - `lfc_shrink_opts`: A list of used \code{\link[DESeq2:lfcShrink]{lfcShrink}} parameters.
#' @family DTUrtle DGE
#' @export
#' @seealso [import_dge_counts()] for correct import of gene-level counts. [combine_to_matrix()] to summarize scRNA counts to one matrix.
assertthat::assert_that(requireNamespace("Seurat",quietly=T),msg="The package Seurat is needed if a Seurat object is provided.")
assertthat::assert_that(utils::packageVersion("Seurat")>="3.0.0",msg="At least Version 3 of Seurat is needed. Currently only Seurat 3 objects are supported.")
assertthat::assert_that(counts@version>="3.0.0",msg="The provided 'counts' is not a Seurat 3 object. Currently only Seurat 3 objects are supported.")
counts<-Seurat::GetAssayData(counts)
}
assertthat::assert_that(methods::is(counts,"matrix")|methods::is(counts,"sparseMatrix")|is.list(counts),msg="'counts' must be a (sparse) Matrix or a tximport result object.")
assertthat::assert_that(is.null(id_col)||(is.character(id_col)&&length(id_col)==1&&id_col%in%colnames(pd)),msg="id_col should be a single column name of pd or NULL.")
assertthat::assert_that((is.character(cond_col)&&length(cond_col)==1&&cond_col%in%colnames(pd)),msg=paste0("Could not find",cond_col," in column names of pd."))
assertthat::assert_that(is.null(cond_levels)||length(cond_levels)==2,msg="'cond_levels' should be of length two or NULL.")
assertthat::assert_that(is.numeric(lfc_threshold),msg="'lfc_threshold' must be numeric. Use a theshold of 0 to disable subsetting.")
assertthat::assert_that(is.numeric(sig_threshold),msg="'sig_threshold' must be numeric. Use a threshold of 1 to disable subsetting.")
assertthat::assert_that(dge_calling_strategy%in%c("bulk","sc"),msg="Please select a valid dge_calling_strategy ('bulk' or 'sc').")
assertthat::assert_that(is.null(subset_feature)|length(subset_feature)>0,msg="`subset_feature` must be `NULL` or of length>=1.")
assertthat::assert_that(is.null(subset_sample)|length(subset_sample)>0,msg="`subset_sample` must be `NULL` or of length>=1.")
assertthat::assert_that(is.list(deseq_opts),msg="`deseq_opts` must be a list.")
assertthat::assert_that(is.list(lfc_shrink_opts),msg="`lfc_shrink_opts` must be a list.")
assertthat::assert_that(is.logical(return_dds),msg="`return_dds` must be TRUE or FALSE.")
assertthat::assert_that(methods::is(BPPARAM,"BiocParallelParam"),msg="Please provide a valid BiocParallelParam object.")
warning("It is strongly advised to provide the whole tximport object (the resulting object of 'import_dge_counts()') for DGE calling, as additional effective transcript length information can enhance result quality.")
}
if(is.list(counts)){
assertthat::assert_that(all(c("counts","length","countsFromAbundance")%in%names(counts)),msg="The provided counts object is a list but not a valid tximport object.")
eff_len<-counts$length
cfa<-counts$countsFromAbundance
counts<-counts$counts
}
if(is.null(id_col)){
assertthat::assert_that(all(rownames(pd)%in%colnames(counts)),msg="Provided id_col does not match with sample names in counts.")
assertthat::assert_that(length(cond_levels)==2,msg="More than two levels found in 'cond_col'. Please specify the two levels you want to compare in 'cond_levels'.")
message("\nComparing in '",cond_col,"': '",cond_levels[1],"' vs '",cond_levels[2],"'")
