Merge pull request #843 from maxplanck-ie/allelic_fixes

Salmon Allelic
maxplanck-ie · Dec 13, 2022 · 5ca6002 · 5ca6002
2 parents 39ee63e + 03b9ada
commit 5ca6002
Show file tree

Hide file tree

Showing 22 changed files with 546 additions and 110 deletions.
diff --git a/.ci_stuff/test_dag.sh b/.ci_stuff/test_dag.sh
@@ -184,47 +184,50 @@ if [ ${PIPESTATUS[0]} -ne 0 ] || [ $WC -ne 333 ]; then exit 1 ; fi
 
 # mRNA-seq
 WC=`mRNA-seq -i PE_input -o output --snakemakeOptions " --dryrun --conda-prefix /tmp" .ci_stuff/organism.yaml | tee >(cat 1>&2) | grep -v "Conda environment" | sed '/^\s*$/d' | wc -l`
-if [ ${PIPESTATUS[0]} -ne 0 ] || [ $WC -ne 882 ]; then exit 1 ; fi
+if [ ${PIPESTATUS[0]} -ne 0 ] || [ $WC -ne 879 ]; then exit 1 ; fi
 WC=`mRNA-seq -i PE_input -o output --sampleSheet .ci_stuff/test_sampleSheet.tsv --snakemakeOptions " --dryrun --conda-prefix /tmp" .ci_stuff/organism.yaml | tee >(cat 1>&2) | grep -v "Conda environment" | sed '/^\s*$/d' | wc -l`
-if [ ${PIPESTATUS[0]} -ne 0 ] || [ $WC -ne 892 ]; then exit 1 ; fi
+if [ ${PIPESTATUS[0]} -ne 0 ] || [ $WC -ne 889 ]; then exit 1 ; fi
 WC=`mRNA-seq -i PE_input -o output --rMats --sampleSheet .ci_stuff/test_sampleSheet.tsv --snakemakeOptions " --dryrun --conda-prefix /tmp" .ci_stuff/organism.yaml | tee >(cat 1>&2) | grep -v "Conda environment" | sed '/^\s*$/d' | wc -l`
-if [ ${PIPESTATUS[0]} -ne 0 ] || [ $WC -ne 908 ]; then exit 1 ; fi
+if [ ${PIPESTATUS[0]} -ne 0 ] || [ $WC -ne 905 ]; then exit 1 ; fi
 WC=`mRNA-seq -i PE_input -o output --sampleSheet .ci_stuff/test_sampleSheet.tsv --snakemakeOptions " --dryrun --conda-prefix /tmp" -m "alignment" .ci_stuff/organism.yaml | tee >(cat 1>&2) | grep -v "Conda environment" | sed '/^\s*$/d' | wc -l`
-if [ ${PIPESTATUS[0]} -ne 0 ] || [ $WC -ne 602 ]; then exit 1 ; fi
+if [ ${PIPESTATUS[0]} -ne 0 ] || [ $WC -ne 599 ]; then exit 1 ; fi
 WC=`mRNA-seq -i PE_input -o output --sampleSheet .ci_stuff/test_sampleSheet.tsv --snakemakeOptions " --dryrun --conda-prefix /tmp" -m "alignment,deepTools_qc" --trim .ci_stuff/organism.yaml | tee >(cat 1>&2) | grep -v "Conda environment" | sed '/^\s*$/d' | wc -l`
-if [ ${PIPESTATUS[0]} -ne 0 ] || [ $WC -ne 948 ]; then exit 1 ; fi
+if [ ${PIPESTATUS[0]} -ne 0 ] || [ $WC -ne 945 ]; then exit 1 ; fi
 WC=`mRNA-seq -i PE_input -o output --sampleSheet .ci_stuff/test_sampleSheet.tsv --snakemakeOptions " --dryrun --conda-prefix /tmp" -m "alignment-free,deepTools_qc" .ci_stuff/organism.yaml | tee >(cat 1>&2) | grep -v "Conda environment" | sed '/^\s*$/d' | wc -l`
-if [ ${PIPESTATUS[0]} -ne 0 ] || [ $WC -ne 1005 ]; then exit 1 ; fi
+if [ ${PIPESTATUS[0]} -ne 0 ] || [ $WC -ne 1006 ]; then exit 1 ; fi
 WC=`mRNA-seq -i PE_input -o output --sampleSheet .ci_stuff/test_sampleSheet.tsv --snakemakeOptions " --dryrun --conda-prefix /tmp" -m "alignment,deepTools_qc" --bcExtract --trim .ci_stuff/organism.yaml | tee >(cat 1>&2) | grep -v "Conda environment" | sed '/^\s*$/d' | wc -l`
-if [ ${PIPESTATUS[0]} -ne 0 ] || [ $WC -ne 916 ]; then exit 1 ; fi
+if [ ${PIPESTATUS[0]} -ne 0 ] || [ $WC -ne 913 ]; then exit 1 ; fi
 WC=`mRNA-seq -i PE_input -o output --sampleSheet .ci_stuff/test_sampleSheet.tsv --snakemakeOptions " --dryrun --conda-prefix /tmp" -m "alignment,deepTools_qc" --bcExtract --UMIDedup --trim .ci_stuff/organism.yaml | tee >(cat 1>&2) | grep -v "Conda environment" | sed '/^\s*$/d' | wc -l`
-if [ ${PIPESTATUS[0]} -ne 0 ] || [ $WC -ne 966 ]; then exit 1 ; fi
+if [ ${PIPESTATUS[0]} -ne 0 ] || [ $WC -ne 963 ]; then exit 1 ; fi
 WC=`mRNA-seq -i SE_input -o output --sampleSheet .ci_stuff/test_sampleSheet.tsv --snakemakeOptions " --dryrun --conda-prefix /tmp" .ci_stuff/organism.yaml | tee >(cat 1>&2) | grep -v "Conda environment" | sed '/^\s*$/d' | wc -l`
-if [ ${PIPESTATUS[0]} -ne 0 ] || [ $WC -ne 807 ]; then exit 1 ; fi
+if [ ${PIPESTATUS[0]} -ne 0 ] || [ $WC -ne 804 ]; then exit 1 ; fi
 WC=`mRNA-seq -i SE_input -o output --sampleSheet .ci_stuff/test_sampleSheet.tsv --snakemakeOptions " --dryrun --conda-prefix /tmp" -m "alignment" .ci_stuff/organism.yaml | tee >(cat 1>&2) | grep -v "Conda environment" | sed '/^\s*$/d' | wc -l`
-if [ ${PIPESTATUS[0]} -ne 0 ] || [ $WC -ne 526 ]; then exit 1 ; fi
+if [ ${PIPESTATUS[0]} -ne 0 ] || [ $WC -ne 523 ]; then exit 1 ; fi
 WC=`mRNA-seq -i SE_input -o output --sampleSheet .ci_stuff/test_sampleSheet.tsv --snakemakeOptions " --dryrun --conda-prefix /tmp" -m "alignment,deepTools_qc" --trim .ci_stuff/organism.yaml | tee >(cat 1>&2) | grep -v "Conda environment" | sed '/^\s*$/d' | wc -l`
-if [ ${PIPESTATUS[0]} -ne 0 ] || [ $WC -ne 863 ]; then exit 1 ; fi
+if [ ${PIPESTATUS[0]} -ne 0 ] || [ $WC -ne 860 ]; then exit 1 ; fi
 WC=`mRNA-seq -i SE_input -o output --sampleSheet .ci_stuff/test_sampleSheet.tsv --snakemakeOptions " --dryrun --conda-prefix /tmp" -m "alignment-free,deepTools_qc" .ci_stuff/organism.yaml | tee >(cat 1>&2) | grep -v "Conda environment" | sed '/^\s*$/d' | wc -l`
-if [ ${PIPESTATUS[0]} -ne 0 ] || [ $WC -ne 920 ]; then exit 1 ; fi
+if [ ${PIPESTATUS[0]} -ne 0 ] || [ $WC -ne 921 ]; then exit 1 ; fi
 WC=`mRNA-seq -i SE_input -o output --sampleSheet .ci_stuff/test_sampleSheet.tsv --snakemakeOptions " --dryrun --conda-prefix /tmp" --trim --fastqc .ci_stuff/organism.yaml | tee >(cat 1>&2) | grep -v "Conda environment" | sed '/^\s*$/d' | wc -l`
-if [ ${PIPESTATUS[0]} -ne 0 ] || [ $WC -ne 975 ]; then exit 1 ; fi
+if [ ${PIPESTATUS[0]} -ne 0 ] || [ $WC -ne 972 ]; then exit 1 ; fi
 WC=`mRNA-seq -i BAM_input/filtered_bam -o output --sampleSheet .ci_stuff/test_sampleSheet.tsv --snakemakeOptions " --dryrun --conda-prefix /tmp" --fromBAM .ci_stuff/organism.yaml | tee >(cat 1>&2) | grep -v "Conda environment" | sed '/^\s*$/d' | wc -l`
-if [ ${PIPESTATUS[0]} -ne 0 ] || [ $WC -ne 659 ]; then exit 1 ; fi
+if [ ${PIPESTATUS[0]} -ne 0 ] || [ $WC -ne 657 ]; then exit 1 ; fi
 #multiple comparison groups
 WC=`mRNA-seq --mode alignment,alignment-free -i PE_input -o output --rMats --sampleSheet .ci_stuff/test_sampleSheet_multiComp.tsv --snakemakeOptions " --dryrun --conda-prefix /tmp" .ci_stuff/organism.yaml | tee >(cat 1>&2) | grep -v "Conda environment" | sed '/^\s*$/d' | wc -l`
-if [ ${PIPESTATUS[0]} -ne 0 ] || [ $WC -ne 867 ]; then exit 1 ; fi
+if [ ${PIPESTATUS[0]} -ne 0 ] || [ $WC -ne 869 ]; then exit 1 ; fi
 #allelic
 WC=`mRNA-seq -m allelic-mapping,deepTools_qc -i PE_input -o output --snakemakeOptions " --dryrun --conda-prefix /tmp" --VCFfile allelic_input/file.vcf.gz --strains strain1,strain2 .ci_stuff/organism.yaml | tee >(cat 1>&2) | grep -v "Conda environment" | sed '/^\s*$/d' | wc -l`
-if [ ${PIPESTATUS[0]} -ne 0 ] || [ $WC -ne 1360 ]; then exit 1 ; fi
+if [ ${PIPESTATUS[0]} -ne 0 ] || [ $WC -ne 1357 ]; then exit 1 ; fi
+WC=`mRNA-seq -m allelic-mapping,deepTools_qc -i allelic_BAM_input/filtered_bam --fromBAM --bamExt '.filtered.bam' -o output --sampleSheet .ci_stuff/test_sampleSheet.tsv --snakemakeOptions " --dryrun --conda-prefix /tmp" --SNPfile allelic_input/snpfile.txt --NMaskedIndex allelic_input/Ngenome .ci_stuff/organism.yaml | tee >(cat 1>&2) | grep -v "Conda environment" | sed '/^\s*$/d' | wc -l`
+if [ ${PIPESTATUS[0]} -ne 0 ] || [ $WC -ne 1118 ]; then exit 1 ; fi
 WC=`mRNA-seq -m allelic-mapping,deepTools_qc -i PE_input -o output --sampleSheet .ci_stuff/test_sampleSheet.tsv --snakemakeOptions " --dryrun --conda-prefix /tmp" --VCFfile allelic_input/file.vcf.gz --strains strain1,strain2 .ci_stuff/organism.yaml | tee >(cat 1>&2) | grep -v "Conda environment" | sed '/^\s*$/d' | wc -l`
-if [ ${PIPESTATUS[0]} -ne 0 ] || [ $WC -ne 1370 ]; then exit 1 ; fi
+if [ ${PIPESTATUS[0]} -ne 0 ] || [ $WC -ne 1367 ]; then exit 1 ; fi
 WC=`mRNA-seq -m allelic-mapping,deepTools_qc -i PE_input -o output --sampleSheet .ci_stuff/test_sampleSheet.tsv --snakemakeOptions " --dryrun --conda-prefix /tmp" --SNPfile allelic_input/snpfile.txt --NMaskedIndex allelic_input/Ngenome .ci_stuff/organism.yaml | tee >(cat 1>&2) | grep -v "Conda environment" | sed '/^\s*$/d' | wc -l`
-if [ ${PIPESTATUS[0]} -ne 0 ] || [ $WC -ne 1353 ]; then exit 1 ; fi
+if [ ${PIPESTATUS[0]} -ne 0 ] || [ $WC -ne 1350 ]; then exit 1 ; fi
 WC=`mRNA-seq -m allelic-mapping,deepTools_qc -i PE_input -o output --sampleSheet .ci_stuff/test_sampleSheet.tsv --snakemakeOptions " --dryrun --conda-prefix /tmp" --VCFfile allelic_input/file.vcf.gz --strains strain1 .ci_stuff/organism.yaml | tee >(cat 1>&2) | grep -v "Conda environment" | sed '/^\s*$/d' | wc -l`
-if [ ${PIPESTATUS[0]} -ne 0 ] || [ $WC -ne 1370 ]; then exit 1 ; fi
+if [ ${PIPESTATUS[0]} -ne 0 ] || [ $WC -ne 1367 ]; then exit 1 ; fi
+WC=`mRNA-seq -m allelic-mapping,deepTools_qc,alignment-free -i PE_input -o output --sampleSheet .ci_stuff/test_sampleSheet.tsv --snakemakeOptions " --dryrun --conda-prefix /tmp" --VCFfile allelic_input/file.vcf.gz --strains strain1 .ci_stuff/organism.yaml | tee >(cat 1>&2) | grep -v "Conda environment" | sed '/^\s*$/d' | wc -l`
+if [ ${PIPESTATUS[0]} -ne 0 ] || [ $WC -ne 1803 ]; then exit 1 ; fi
 
-# noncoding-RNA-seq
 WC=`noncoding-RNA-seq -i PE_input -o output --snakemakeOptions " --dryrun --conda-prefix /tmp" .ci_stuff/organism.yaml | tee >(cat 1>&2) | grep -v "Conda environment" | sed '/^\s*$/d' | wc -l`
 if [ ${PIPESTATUS[0]} -ne 0 ] || [ $WC -ne 717 ]; then exit 1 ; fi
 WC=`noncoding-RNA-seq -i PE_input -o output --sampleSheet .ci_stuff/test_sampleSheet.tsv --snakemakeOptions " --dryrun --conda-prefix /tmp" .ci_stuff/organism.yaml | tee >(cat 1>&2) | grep -v "Conda environment" | sed '/^\s*$/d' | wc -l`

diff --git a/conda-recipe/meta.yaml b/conda-recipe/meta.yaml
@@ -1,6 +1,6 @@
 package:
   name: snakepipes
-  version: 2.6.1
+  version: 2.7.0
 
 source:
   path: ../

diff --git a/docs/content/News.rst b/docs/content/News.rst
@@ -1,11 +1,19 @@
 snakePipes News
 ===============
 
+
+snakePipes 2.7.0
+----------------
+
+* Added the allelic version of Salmon-based transcript quantitation to mRNA-seq workflow. Will be run if *both* 'allelic-mapping' and 'alignment-free' modes are specified. An allelic version of sleuth will be run, if sample sheet is provided, as well as DESeq2 on allelic Salmon counts.
+
+
 snakePipes 2.6.1
 ----------------
 
 * Capped tabulate version as 0.9.0 breaks snakemake
 
+
 snakePipes 2.6.0
 ----------------
 
@@ -19,6 +27,7 @@ snakePipes 2.6.0
 * Fixed a couple of issues in the ATAC-seq workflow after sofware versions update.
 * Fixed genome size conversion to string.
 
+
 snakePipes 2.5.4
 ----------------
 

diff --git a/snakePipes/__init__.py b/snakePipes/__init__.py
@@ -1 +1 @@
-__version__ = '2.6.1'
+__version__ = '2.7.0'
diff --git a/snakePipes/common_functions.py b/snakePipes/common_functions.py
@@ -22,6 +22,7 @@ def set_env_yamls():
     return {'CONDA_SHARED_ENV': 'envs/shared.yaml',
             'CONDA_CREATE_INDEX_ENV': 'envs/createIndices.yaml',
             'CONDA_RNASEQ_ENV': 'envs/rna_seq.yaml',
+            'CONDA_SLEUTH_ENV': 'envs/sleuth.yaml',
             'CONDA_RMATS_ENV': 'envs/rMats.yaml',
             'CONDA_scRNASEQ_ENV': 'envs/sc_rna_seq.yaml',
             'CONDA_seurat_ENV': 'envs/sc_rna_seq_seurat.yaml',

diff --git a/snakePipes/shared/rscripts/DESeq2.R b/snakePipes/shared/rscripts/DESeq2.R
@@ -79,7 +79,7 @@ if(isTRUE(tximport)) {
   tx2gene <- read.delim(tx2gene_file, header = FALSE)
   tx2gene <- tx2gene[c(1,2)]
   # check setup table and import
-  countdata <- checktable(sampleSheet = sampleInfo, salmon_dir = dirname(countFilePath), tx2gene_annot = tx2gene)
+  countdata <- checktable(sampleSheet = sampleInfo, salmon_dir = dirname(countFilePath), tx2gene_annot = tx2gene, alleleSpecific = allelic_info)
 } else {
   sampleInfo$name <- make.names(sampleInfo$name)
   rownames(sampleInfo)<-sampleInfo$name
@@ -93,19 +93,23 @@ if(isTRUE(tximport)) {
 ## ~~~~~~~ 3. run DESeq wrapper ~~~~~~~~
 #in case of the allelic-specific workflow, allow for 1 condition and skip deseq2 basic in this case 
 if(length(unique(sampleInfo$condition))>1){
-    seqout <- DESeq_basic(countdata, coldata = sampleInfo, fdr = fdr, alleleSpecific = allelic_info, from_salmon = tximport)
+    if(tximport & allelic_info){
+        message("Detected allelic Salmon counts. Skipping DESeq_basic.")
+    }else{
+        seqout <- DESeq_basic(countdata, coldata = sampleInfo, fdr = fdr, alleleSpecific = allelic_info, from_salmon = tximport)
 
-    DESeq_writeOutput(DEseqout = seqout,
+        DESeq_writeOutput(DEseqout = seqout,
                 fdr = fdr, outprefix = "DEseq_basic",
                 geneNamesFile = geneNamesFilePath)
+        }
 }
 #DESeq_downstream(DEseqout = seqout, countdata, sampleInfo,
 #             fdr = fdr, outprefix = "DEseq_basic", heatmap_topN = topN,
 #             geneNamesFile = geneNamesFilePath)
 
 ## Run allele-sepecific DESeq wrapper (if asked for)
 if (isTRUE(allelic_info)) {
-    seqout_allelic <- DESeq_allelic(countdata, coldata = sampleInfo, fdr = fdr)
+    seqout_allelic <- DESeq_allelic(countdata, coldata = sampleInfo, fdr = fdr, from_salmon=tximport)
 
     DESeq_writeOutput(DEseqout = seqout_allelic,
                  fdr = fdr, outprefix = "DEseq_allelic",
@@ -132,19 +136,21 @@ write.bibtex(bib, file = 'citations.bib')
 file.copy(rmdTemplate, to = 'DESeq2_report_basic.Rmd')
 
 if(length(unique(sampleInfo$condition))>1){
-  outprefix = "DEseq_basic"
-  cite_options(citation_format = "text",style = "html",cite.style = "numeric",hyperlink = TRUE)
-  render('DESeq2_report_basic.Rmd',
+  if(!tximport & !allelic_info){
+      outprefix = "DEseq_basic"
+      cite_options(citation_format = "text",style = "html",cite.style = "numeric",hyperlink = TRUE)
+      render('DESeq2_report_basic.Rmd',
               output_format = "html_document",
               clean = TRUE,
               params = list(
                   DEseqoutRdata = paste0(outprefix, "_DESeq.Rdata"),
                   ddr.df = paste0(outprefix, "_DEresults.tsv"),
                   countdata = countFilePath,
                   coldata = sampleInfo,
-                  fdr = 0.05,
+                  fdr = fdr,
                   heatmap_topN = 20,
                   geneNamesFile = geneNamesFilePath))
+    }
 }
 
 if (isTRUE(allelic_info)) {
@@ -158,7 +164,7 @@ if (isTRUE(allelic_info)) {
                       ddr.df = paste0(outprefix, "_DEresults.tsv"),
                       countdata = countFilePath,
                       coldata = sampleInfo,
-                      fdr = 0.05,
+                      fdr = fdr,
                       heatmap_topN = 20,
                       geneNamesFile = geneNamesFilePath))
 }

diff --git a/snakePipes/shared/rscripts/DESeq2Report.Rmd b/snakePipes/shared/rscripts/DESeq2Report.Rmd
@@ -81,7 +81,7 @@ load(params$DEseqoutRdata) ## dds and ddr
 ddr.df <- params$ddr.df  #_DEresults.tsv file
 countdata <- params$countdata
 coldata <- params$coldata
-fdr <- params$fdr
+fdr <- as.numeric(params$fdr)
 heatmap_topN <- params$heatmap_topN
 geneNamesFile <- params$geneNamesFile
 

diff --git a/snakePipes/shared/rscripts/DE_functions.R b/snakePipes/shared/rscripts/DE_functions.R
@@ -31,8 +31,13 @@ checktable <- function(countdata = NA, sampleSheet = NA, alleleSpecific = FALSE,
   ## check files
   if(!is.na(salmon_dir)) {
 
-    # mode = Salmon : check whether salmon output files exist in Salmon dir
-    files <- paste0(salmon_dir,"/",sampleSheet$name,".quant.sf")
+    #mode = Salmon : check whether salmon output files exist in Salmon dir
+    if(alleleSpecific){
+        files <- c(paste0(salmon_dir,"/",sampleSheet$name,".genome1.quant.sf"),paste0(salmon_dir,"/",sampleSheet$name,".genome2.quant.sf"))
+    }else{
+        files <- paste0(salmon_dir,"/",sampleSheet$name,".quant.sf")
+    }
+    #files<-dir(salmon_dir,pattern=".quant.sf",full.names=TRUE)
     filecheck <- file.exists(files)
     if(!(all(filecheck == TRUE))) {
       cat("Error! The following File(s) don't exist : ")
@@ -63,7 +68,7 @@ checktable <- function(countdata = NA, sampleSheet = NA, alleleSpecific = FALSE,
         }
     }
   }
-  if(all(is.integer(countdata))){
+  if(all(is.integer(countdata)) | !is.na(salmon_dir)){
         print("All countdata is integer.")
     }else{
         print("Non-integer counts detected. The data will be rounded, as this is well within the expected sampling variation of a technical replicate.")
@@ -88,25 +93,26 @@ DESeq_basic <- function(countdata, coldata, fdr, alleleSpecific = FALSE, from_sa
     cnames.sub<-unique(colnames(coldata)[2:which(colnames(coldata) %in% "condition")])
     d<-as.formula(noquote(paste0("~",paste(cnames.sub,collapse="+"))))
 
+
     # Normal DESeq
     print("Performing basic DESeq: test vs control")
     if(isTRUE(from_salmon)) {
       print("Using input from tximport")
         dds <- DESeq2::DESeqDataSetFromTximport(countdata,
                                   colData = coldata, design =d)
 
-    } else {
-      print("Using input from count table")
-      if(isTRUE(alleleSpecific)) {
-          rnasamp <- dplyr::select(countdata, dplyr::ends_with("_all"))
-          rownames(coldata)<-colnames(rnasamp)
-          dds <- DESeq2::DESeqDataSetFromMatrix(countData = rnasamp,
-                                    colData = coldata, design =d)
       } else {
-          dds <- DESeq2::DESeqDataSetFromMatrix(countData = countdata,
+          if(isTRUE(alleleSpecific)) {
+            rnasamp <- dplyr::select(countdata, dplyr::ends_with("_all"))
+            rownames(coldata)<-colnames(rnasamp)
+            countdata<-rnasamp
+           }
+
+           dds <- DESeq2::DESeqDataSetFromMatrix(countData = countdata,
                                     colData = coldata, design =d)
+
       }
-    }
+
     if(length(size_factors) > 1) {
         print("applying size factors")
         print(size_factors)
@@ -135,33 +141,50 @@ DESeq_basic <- function(countdata, coldata, fdr, alleleSpecific = FALSE, from_sa
 #'
 #'
 
-DESeq_allelic <- function(countdata, coldata, fdr) {
+DESeq_allelic <- function(countdata, coldata, fdr, from_salmon=FALSE) {
 
     # AlleleSpecific DEseq
     print("Performing Allele-specific DESeq using Interaction design : Genome2 vs Genome1")
+    if(isTRUE(from_salmon)) {
+        # create alleleSpecific design matrix
+        coldata_allelic <- data.frame(name = colnames(as.data.frame(countdata$counts)),
+                   allele = rep(c("genome1", "genome2"), nrow(coldata)),
+                   condition = rep(coldata$condition, each = 2) )
+        rownames(coldata_allelic)<-colnames(as.data.frame(countdata$counts))
+        coldata_allelic$allele<-factor(coldata_allelic$allele,levels=c("genome1","genome2"))
+        coldata_allelic$condition<-factor(coldata_allelic$condition,levels=unique(coldata_allelic$condition))
+        print("Using input from tximport")
+        dds <- DESeq2::DESeqDataSetFromTximport(countdata,
+                                  colData = coldata_allelic, design =~1)
+
+    } else {
+      print("Using input from count table")
     rnasamp <- dplyr::select(countdata, -dplyr::ends_with("_all"))
 
     # create alleleSpecific design matrix
-    design <- data.frame(name = colnames(rnasamp),
+    coldata_allelic <- data.frame(name = colnames(rnasamp),
                    allele = rep(c("genome1", "genome2"), nrow(coldata)),
                    condition = rep(coldata$condition, each = 2) )
-    rownames(design)<-colnames(rnasamp)
+    rownames(coldata_allelic)<-colnames(rnasamp)
+    coldata_allelic$allele<-factor(coldata_allelic$allele,levels=c("genome1","genome2"))
+    coldata_allelic$condition<-factor(coldata_allelic$condition,levels=unique(coldata_allelic$condition))
+    dds <- DESeq2::DESeqDataSetFromMatrix(rnasamp, colData = coldata_allelic,
+                              design = ~1)
+    rownames(dds) <- rownames(rnasamp)
 
+    }
+
     # Run DESeq
-    if(length(unique(design$condition))>1){
-      dds <- DESeq2::DESeqDataSetFromMatrix(rnasamp, colData = design,
-                              design = ~allele + condition + allele:condition)
-      rownames(dds) <- rownames(rnasamp)
+    if(length(unique(coldata_allelic$condition))>1){
+      DESeq2::design(dds) <- formula(~allele + condition + allele:condition)
       dds <- DESeq2::DESeq(dds,betaPrior = FALSE)
       ddr <- DESeq2::results(dds, name=paste0("allelegenome2.condition",unique(coldata$condition)[2]))
       ddr_shrunk <- DESeq2::lfcShrink(dds,coef=paste0("allelegenome2.condition",unique(coldata$condition)[2]),type="apeglm",res=ddr)
     } else {
-    dds <- DESeq2::DESeqDataSetFromMatrix(rnasamp, colData = design,
-                              design = ~allele)
-      rownames(dds) <- rownames(rnasamp)
-      dds <- DESeq2::DESeq(dds,betaPrior = FALSE)
-      ddr <- DESeq2::results(dds, name="allele_genome2_vs_genome1")
-      ddr_shrunk <- DESeq2::lfcShrink(dds,coef="allele_genome2_vs_genome1",type="apeglm",res=ddr)
+        DESeq2::design(dds) <- formula(~allele)
+        dds <- DESeq2::DESeq(dds,betaPrior = FALSE)
+        ddr <- DESeq2::results(dds, name="allele_genome2_vs_genome1")
+        ddr_shrunk <- DESeq2::lfcShrink(dds,coef="allele_genome2_vs_genome1",type="apeglm",res=ddr)
     }
     output <- list(dds = dds, ddr = ddr, ddr_shrunk=ddr_shrunk)
     return(output)