Patch 1.0.5 (#206)

* simplify results output @vyepez88 * add MAE error messages for empty output @mumichae * more robust dependency graph generation @nickhsmith @mumichae * Bump version: 1.0.4 → 1.0.5 Co-authored-by: Vicente <[email protected]> Co-authored-by: nickhsmith <[email protected]>
gagneurlab · Jun 11, 2021 · 1527d24 · 1527d24
1 parent 9069bad
commit 1527d24
Show file tree

Hide file tree

Showing 13 changed files with 68 additions and 45 deletions.
diff --git a/docs/source/conf.py b/docs/source/conf.py
@@ -23,7 +23,7 @@
 author = 'Michaela Müller'
 
 # The full version, including alpha/beta/rc tags
-release_ = '1.0.4'
+release_ = '1.0.5'
 
 
 # -- General configuration ---------------------------------------------------

diff --git a/drop/__init__.py b/drop/__init__.py
@@ -4,4 +4,4 @@
 from . import utils
 from . import demo
 
-__version__ = "1.0.4"
+__version__ = "1.0.5"
diff --git a/drop/cli.py b/drop/cli.py
@@ -17,7 +17,7 @@
 
 @click.group()
 @click_log.simple_verbosity_option(logger)
-@click.version_option('1.0.4',prog_name='drop')
+@click.version_option('1.0.5',prog_name='drop')
 def main():
     pass
 

diff --git a/drop/modules/aberrant-expression-pipeline/OUTRIDER/Summary.R b/drop/modules/aberrant-expression-pipeline/OUTRIDER/Summary.R
@@ -15,6 +15,8 @@
 #'  output:
 #'   - wBhtml: '`sm config["htmlOutputPath"] + 
 #'              "/AberrantExpression/Outrider/{annotation}/Summary_{dataset}.html"`'
+#'   - res_html: '`sm config["htmlOutputPath"] + 
+#'              "/AberrantExpression/Outrider/{annotation}/OUTRIDER_results_{dataset}.tsv"`'
 #'  type: noindex
 #' output:
 #'  html_document:
@@ -125,19 +127,13 @@ if (nrow(res) > 0) {
 
 
 #' ## Results table
-#+echo=F
+
+## Save results table in the html folder and provide link to download
+file <- snakemake@output$res_html
+fwrite(res, file, sep = '\t', quote = F)
+#+ echo=FALSE, results='asis'
+cat(paste0("<a href='./", basename(file), "'>Download OUTRIDER results table</a>"))
+
 res[, pValue := format(pValue, scientific = T, digits = 2)]
 res[, padjust := format(padjust, scientific = T, digits = 2)]
 DT::datatable(res, caption = "OUTRIDER results", style = 'bootstrap', filter = 'top')
-
-#' ### Download results table
-web_dir <- snakemake@config$webDir
-if (!is.null(web_dir)) {
-  results_link <- paste0(web_dir, 
-                         "/aberrant_expression/results/", snakemake@wildcards$annotation,
-                         "/outrider/", snakemake@wildcards$dataset ,"/OUTRIDER_results.tsv")
-} else {
-  results_link <- snakemake@input$results
-}
-#' [Download OUTRIDER results table](`r results_link`)
-
diff --git a/drop/modules/aberrant-expression-pipeline/OUTRIDER/results.R b/drop/modules/aberrant-expression-pipeline/OUTRIDER/results.R
@@ -64,11 +64,3 @@ if(!is.null(sa$HPO_TERMS) & nrow(res) > 0){
 
 # Save results 
 fwrite(res, snakemake@output$results, sep = "\t", quote = F)
-
-web_dir <- snakemake@config$webDir
-if (!is.null(web_dir)) {
-    pub_res <- paste0(web_dir, 
-                      "/aberrant_expression/results/",{snakemake@wildcards$annotation},"/outrider/",
-                      {snakemake@wildcards$dataset},"/OUTRIDER_results.tsv")
-    fwrite(res, pub_res, sep = "\t", quote = F)
-}
diff --git a/drop/modules/aberrant-expression-pipeline/Snakefile b/drop/modules/aberrant-expression-pipeline/Snakefile
@@ -13,7 +13,12 @@ rule aberrantExpression:
 
 rule aberrantExpression_dependency:
     output: AE_graph_file
-    shell: "snakemake --nolock --rulegraph {AE_index_output} | dot -Tsvg -Grankdir=TB > {output}"
+    shell:
+        """
+        snakemake --nolock --rulegraph {AE_index_output} | \
+            sed -ne '/digraph snakemake_dag/,/}}/p' | \
+            dot -Tsvg -Grankdir=TB > {output}
+        """
 
 rule aberrantExpression_bamStats:
     input: 

diff --git a/drop/modules/aberrant-splicing-pipeline/FRASER/Summary.R b/drop/modules/aberrant-splicing-pipeline/FRASER/Summary.R
@@ -14,6 +14,8 @@
 #'  output:
 #'   - wBhtml: '`sm config["htmlOutputPath"] +
 #'               "/AberrantSplicing/{dataset}--{annotation}_summary.html"`'
+#'   - res_html: '`sm config["htmlOutputPath"] +
+#'               "/AberrantSplicing/FRASER_results_{dataset}--{annotation}.tsv"`'
 #'  type: noindex
 #'---
 
@@ -93,15 +95,12 @@ for(type in psiTypes){
   after
 }
 
-#' # Results
+#' ## Results
 res <- fread(snakemake@input$results)
-file <- gsub(".html$", ".tsv", snakemake@output$wBhtml)
-write_tsv(res, file=file)
-
-#'
-#' The results table can also be downloaded with the link below.
+file <- snakemake@output$res_html
+write_tsv(res, file = file)
 #+ echo=FALSE, results='asis'
-cat(paste0("<a href='./", basename(file), "'>Download results table</a>"))
+cat(paste0("<a href='./", basename(file), "'>Download FRASER results table</a>"))
 
 # round numbers
 if(nrow(res) > 0){
@@ -114,8 +113,5 @@ if(nrow(res) > 0){
   res[, padjustGene := signif(padjustGene, 2)]
 }
 
-#' ## Results table
 DT::datatable(res, options=list(scrollX=TRUE), escape=FALSE, filter = 'top')
 
-#' ## Samples table
-DT::datatable(as.data.table(colData(fds)), options=list(scrollX=TRUE))
diff --git a/drop/modules/aberrant-splicing-pipeline/Snakefile b/drop/modules/aberrant-splicing-pipeline/Snakefile
@@ -13,4 +13,9 @@ rule aberrantSplicing:
 
 rule aberrantSplicing_dependency:
     output: AS_graph_file
-    shell: "snakemake --nolock --rulegraph {AS_index_output} | dot -Tsvg -Grankdir=TB > {output}"
+    shell:
+        """
+        snakemake --nolock --rulegraph {AS_index_output} | \
+            sed -ne '/digraph snakemake_dag/,/}}/p' | \
+            dot -Tsvg -Grankdir=TB > {output}
+        """
diff --git a/drop/modules/mae-pipeline/MAE/ASEReadCounter.sh b/drop/modules/mae-pipeline/MAE/ASEReadCounter.sh
@@ -31,8 +31,8 @@ header+="lowBaseQDepth\trawDepth\totherBases\timproperPairs"
 echo -e $header >> $tmp
 
 # get chr format
-bam_chr=$($samtools idxstats ${bam_file} | grep -vP "\t0\t0" | cut -f1) # only chr with coverage
-vcf_chr=$($bcftools view ${vcf_file} | cut -f1 | grep -v '#' | uniq)
+bam_chr=$($samtools idxstats ${bam_file} | grep -vP "\t0\t0" | cut -f1 | sort -u) # only chr with coverage
+vcf_chr=$($bcftools view ${vcf_file} | cut -f1 | grep -v '#' | sort -u)
 if [ "$(echo ${vcf_chr} | grep -c 'chr')" -eq 0 ]; then
   echo "use NCBI format"
   canonical=$ncbi2ucsc
@@ -42,11 +42,10 @@ else
 fi
 
 # subset to standard chromosomes
-chr_subset=$(comm -12 <(cut -f1 -d" " ${canonical} | sort -u) <(echo "${vcf_chr}" | sort -u))
-chr_subset=$(comm -12 <(echo "${bam_chr}" | sort -u) <(echo "${chr_subset}") | uniq)
+chr_subset=$(comm -12 <(cut -f1 -d" " ${canonical} | sort -u) <(echo "${vcf_chr}"))
+chr_subset=$(comm -12 <(echo "${bam_chr}") <(echo "${chr_subset}") | uniq)
 
 for chr in $chr_subset; do
-  echo $chr
   $gatk ASEReadCounter \
     -R ${fasta} \
     -I ${bam_file} \
@@ -58,10 +57,23 @@ for chr in $chr_subset; do
     tail -n+2 >>$tmp
 done
 
-echo $mae_id
 cat $tmp | awk -v id="${mae_id}" \
   -F $'\t' 'BEGIN {OFS = FS} NR==1{print $0, "ID"} NR>1{print $0, id}' |
   bgzip >${output}
 rm ${tmp}
 
+num_out=$(zcat "${output}" | wc -l )
+if [ "${num_out}" -lt 2 ]
+then
+  printf  "%s\n" "" "ERROR: No allele-specific counts" \
+    "  Make sure that the chromosome styles of the FASTA reference and BAM file match." \
+    "  If that isn't the issue, check that your VCF and BAM files are correctly formatted." \
+    "  If this problem persists and if this is your only sample causing issues, consider removing it from your analysis, as a last resort." \
+    "" "  MAE ID: ${mae_id}" \
+    "  VCF file: ${vcf_file}" \
+    "  BAM file: ${bam_file}" \
+    "  FASTA file: ${fasta}"
+  exit 1
+fi
+
 zcat ${output} | head
diff --git a/drop/modules/mae-pipeline/MAE/filterSNVs.sh b/drop/modules/mae-pipeline/MAE/filterSNVs.sh
@@ -51,5 +51,17 @@ else  # VCF and BAM have same chromosome format
     rm ${tmp}.tbi
 fi
 
+num_out=$(zcat "${output}" | grep -vc '#' )
+if [ "${num_out}" -eq 0 ]
+then
+  printf  "%s\n" "" "ERROR: No entries after filtering for SNVs" \
+  "  Make sure that the VCF file is correctly formatted and contains heterozygous variants." \
+  "  This analysis is independent per sample, so consider removing the sample from your analysis as a last resort." \
+  "" "  VCF ID: ${vcf_id}" \
+  "  VCF file: ${vcf_file}" \
+  "  BAM file: ${bam_file}"
+  exit 1
+fi
+
 $bcftools index -t ${output}
 
diff --git a/drop/modules/mae-pipeline/Snakefile b/drop/modules/mae-pipeline/Snakefile
@@ -13,7 +13,12 @@ rule mae:
 
 rule mae_dependency:
     output: MAE_graph_file
-    shell: "snakemake --nolock --rulegraph {MAE_index_output} | dot -Tsvg -Grankdir=TB > {output}"
+    shell:
+        """
+        snakemake --nolock --rulegraph {MAE_index_output} | \
+            sed -ne '/digraph snakemake_dag/,/}}/p' | \
+            dot -Tsvg -Grankdir=TB > {output}
+        """
 
 rule sampleQC:
     input: cfg.getHtmlFromScript(MAE_WORKDIR / "QC" / "Datasets.R")

diff --git a/setup.cfg b/setup.cfg
@@ -1,5 +1,5 @@
 [bumpversion]
-current_version = 1.0.4
+current_version = 1.0.5
 commit = True
 
 [bumpversion:file:setup.py]

diff --git a/setup.py b/setup.py
@@ -21,7 +21,7 @@
 
 setuptools.setup(
     name="drop",
-    version="1.0.4",
+    version="1.0.5",
     author="Michaela Müller, Daniela Klaproth-Andrade, Vicente Yépez, Christian Mertes",
     author_email="[email protected]",
     description="Detection of RNA Outlier Pipeline",