Merge pull request #6 from yyoshiaki/v2.0.1

V2.0.1

README.md

@@ -1,4 +1,4 @@
-# VIRTUS2 : VIRal Transcript Usage Sensor v2.0 <img src="https://github.com/yyoshiaki/VIRTUS/raw/master/img/VIRTUS.jpg" width="20%" align="right" />
+# VIRTUS2 : VIRal Transcript Usage Sensor v2.0.1 <img src="https://github.com/yyoshiaki/VIRTUS/raw/master/img/VIRTUS.jpg" width="20%" align="right" />
 
 **!!Note : We updated VIRTUS to version2. In this version, we removed the gene quantification step by Salmon and single virus mode, and added coverage on viral genomes to the result to focus on virus-wide exploration. If you want to use the single virus mode, visit [https://github.com/yyoshiaki/VIRTUS](https://github.com/yyoshiaki/VIRTUS)**
 

test/singularity.test.sh

@@ -20,7 +20,7 @@ fi
 echo "null" > null.txt
 ls | grep -v -E 'test.sh' | grep -v -E 'ERR3240275' | grep -v -E 'SRR8315715' | xargs rm -r
 
-cwltool --rm-tmpdir --singularity ../workflow/createindex.cwl ../workflow/createindex.job.yaml
+cwltool --tmp-outdir-prefix=${PWD}/tmp_cwl/ --tmpdir-prefix=${PWD}/tmp_cwl/ --rm-tmpdir --singularity ../workflow/createindex.cwl ../workflow/createindex.job.yaml
 
 cd ERR3240275
 if [[ ! -e ./ERR3240275_1.fastq.gz ]]; then
@@ -29,7 +29,7 @@ if [[ ! -e ./ERR3240275_1.fastq.gz ]]; then
   wget ftp://ftp.sra.ebi.ac.uk/vol1/fastq/ERR324/005/ERR3240275/ERR3240275_1.fastq.gz
   wget ftp://ftp.sra.ebi.ac.uk/vol1/fastq/ERR324/005/ERR3240275/ERR3240275_2.fastq.gz
 fi
-cwltool --rm-tmpdir --singularity ../../workflow/VIRTUS.PE.cwl ../../workflow/VIRTUS.PE.job.yaml
+cwltool --tmp-outdir-prefix=${PWD}/tmp_cwl/ --tmpdir-prefix=${PWD}/tmp_cwl/ --rm-tmpdir --singularity ../../workflow/VIRTUS.PE.cwl ../../workflow/VIRTUS.PE.job.yaml
 cd ..
 
 
@@ -40,7 +40,7 @@ if [[ ! -e ./SRR8315715_1.fastq.gz ]]; then
   wget ftp://ftp.sra.ebi.ac.uk/vol1/fastq/SRR831/005/SRR8315715/SRR8315715_1.fastq.gz
   wget ftp://ftp.sra.ebi.ac.uk/vol1/fastq/SRR831/005/SRR8315715/SRR8315715_2.fastq.gz
 fi
-cwltool --rm-tmpdir --singularity ../../workflow/VIRTUS.SE.cwl ../../workflow/VIRTUS.SE.job.yaml
+cwltool --tmp-outdir-prefix=${PWD}/tmp_cwl/ --tmpdir-prefix=${PWD}/tmp_cwl/ --rm-tmpdir --singularity ../../workflow/VIRTUS.SE.cwl ../../workflow/VIRTUS.SE.job.yaml
 cd ..
 
 echo Successfully completed singularity.test.sh!

test/test.sh

@@ -20,7 +20,7 @@ fi
 echo "null" > null.txt
 ls | grep -v -E 'test.sh' | grep -v -E 'ERR3240275' | grep -v -E 'SRR8315715' | xargs rm -r
 
-cwltool --rm-tmpdir ../workflow/createindex.cwl ../workflow/createindex.job.yaml
+cwltool  --tmp-outdir-prefix=${PWD}/tmp_cwl/ --tmpdir-prefix=${PWD}/tmp_cwl/ --rm-tmpdir ../workflow/createindex.cwl ../workflow/createindex.job.yaml
 
 cd ERR3240275
 if [[ ! -e ./ERR3240275_1.fastq.gz ]]; then
@@ -29,7 +29,7 @@ if [[ ! -e ./ERR3240275_1.fastq.gz ]]; then
   wget ftp://ftp.sra.ebi.ac.uk/vol1/fastq/ERR324/005/ERR3240275/ERR3240275_1.fastq.gz
   wget ftp://ftp.sra.ebi.ac.uk/vol1/fastq/ERR324/005/ERR3240275/ERR3240275_2.fastq.gz
 fi
-cwltool --rm-tmpdir ../../workflow/VIRTUS.PE.cwl ../../workflow/VIRTUS.PE.job.yaml
+cwltool --tmp-outdir-prefix=${PWD}/tmp_cwl/ --tmpdir-prefix=${PWD}/tmp_cwl/ --rm-tmpdir ../../workflow/VIRTUS.PE.cwl ../../workflow/VIRTUS.PE.job.yaml
 cd ..
 
 
@@ -40,7 +40,7 @@ if [[ ! -e ./SRR8315715_1.fastq.gz ]]; then
 # wget ftp://ftp.sra.ebi.ac.uk/vol1/fastq/SRR831/005/SRR8315715/SRR8315715_1.fastq.gz
   wget ftp://ftp.sra.ebi.ac.uk/vol1/fastq/SRR831/005/SRR8315715/SRR8315715_2.fastq.gz
 fi
-cwltool --rm-tmpdir ../../workflow/VIRTUS.SE.cwl ../../workflow/VIRTUS.SE.job.yaml
+cwltool --tmp-outdir-prefix=${PWD}/tmp_cwl/ --tmpdir-prefix=${PWD}/tmp_cwl/ --rm-tmpdir ../../workflow/VIRTUS.SE.cwl ../../workflow/VIRTUS.SE.job.yaml
 cd ..
 
-echo Successfully completed test.sh!
+echo Successfully completed test.sh!

test/udocker.test.sh

@@ -20,7 +20,7 @@ fi
 echo "null" > null.txt
 ls | grep -v -E 'test.sh' | grep -v -E 'ERR3240275' | grep -v -E 'SRR8315715' | xargs rm -r
 
-cwltool --rm-tmpdir --user-space-docker-cmd=udocker ../workflow/createindex.cwl ../workflow/createindex.job.yaml
+cwltool --tmp-outdir-prefix=${PWD}/tmp_cwl/ --tmpdir-prefix=${PWD}/tmp_cwl/ --rm-tmpdir --user-space-docker-cmd=udocker ../workflow/createindex.cwl ../workflow/createindex.job.yaml
 
 cd ERR3240275
 if [[ ! -e ./ERR3240275_1.fastq.gz ]]; then
@@ -29,7 +29,7 @@ if [[ ! -e ./ERR3240275_1.fastq.gz ]]; then
   wget ftp://ftp.sra.ebi.ac.uk/vol1/fastq/ERR324/005/ERR3240275/ERR3240275_1.fastq.gz
   wget ftp://ftp.sra.ebi.ac.uk/vol1/fastq/ERR324/005/ERR3240275/ERR3240275_2.fastq.gz
 fi
-cwltool --rm-tmpdir --user-space-docker-cmd=udocker ../../workflow/VIRTUS.PE.cwl ../../workflow/VIRTUS.PE.job.yaml
+cwltool --tmp-outdir-prefix=${PWD}/tmp_cwl/ --tmpdir-prefix=${PWD}/tmp_cwl/ --rm-tmpdir --user-space-docker-cmd=udocker ../../workflow/VIRTUS.PE.cwl ../../workflow/VIRTUS.PE.job.yaml
 cd ..
 
 
@@ -40,7 +40,7 @@ if [[ ! -e ./SRR8315715_1.fastq.gz ]]; then
   wget ftp://ftp.sra.ebi.ac.uk/vol1/fastq/SRR831/005/SRR8315715/SRR8315715_1.fastq.gz
   wget ftp://ftp.sra.ebi.ac.uk/vol1/fastq/SRR831/005/SRR8315715/SRR8315715_2.fastq.gz
 fi
-cwltool --rm-tmpdir --user-space-docker-cmd=udocker ../../workflow/VIRTUS.SE.cwl ../../workflow/VIRTUS.SE.job.yaml
+cwltool --tmp-outdir-prefix=${PWD}/tmp_cwl/ --tmpdir-prefix=${PWD}/tmp_cwl/ --rm-tmpdir --user-space-docker-cmd=udocker ../../workflow/VIRTUS.SE.cwl ../../workflow/VIRTUS.SE.job.yaml
 cd ..
 
 echo Successfully completed udocker.test.sh!

tool/samtools/bam_filter_polyx.cwl

@@ -20,5 +20,5 @@ label: bam_filter_polyX
 requirements:
   - class: ShellCommandRequirement
   - class: DockerRequirement
-    dockerPull: 'yyasumizu/bam_filter_polyx:1.2'
+    dockerPull: 'yyasumizu/bam_filter_polyx:1.3'
 stdout: virusAligned.filtered.sortedByCoord.out.bam

wrapper/VIRTUS_wrapper.py

@@ -206,15 +206,31 @@
 uval = pd.Series(dtype = "float64")
 pval = pd.Series(dtype = "float64")
 
-if df["Group"].nunique() == 2:
-    print("Conducting Mann-Whitney U-test")
-    for v in summary.index:
-        u, p = stats.mannwhitneyu(summary.loc[v, ['{}_rate'.format(x) for x in df.loc[df['Group']==df['Group'].unique()[0], 'Name']]],
-        summary.loc[v, ['{}_rate'.format(x) for x in df.loc[df['Group']==df['Group'].unique()[1], 'Name']]], alternative = "two-sided")
-        uval[v] = u
-        pval[v] = p
+# print(summary)
+
+print("### Stats ###")
+print('Threshold coverage : ', th_cov)
+print('Threshold rate : ', th_rate)
+
+print('# detected viruses : ', df_res.shape[0])
+print('Max coverage : ', df_res.coverage.max())
+print('Max rate : ', df_res.rate_hit.max())
 
-    fdr = pd.Series(multipletests(pval,method = "fdr_bh")[1], index = pval.index)
+print('# retained viruses :', summary.shape[0])
+
+if df["Group"].nunique() == 2:
+    if summary.shape[0] > 0:
+        print("Conducting Mann-Whitney U-test")
+        for v in summary.index:
+            u, p = stats.mannwhitneyu(summary.loc[v, ['{}_rate'.format(x) for x in df.loc[df['Group']==df['Group'].unique()[0], 'Name']]],
+            summary.loc[v, ['{}_rate'.format(x) for x in df.loc[df['Group']==df['Group'].unique()[1], 'Name']]], alternative = "two-sided")
+            uval[v] = u
+            pval[v] = p
+
+        fdr = pd.Series(multipletests(pval,method = "fdr_bh")[1], index = pval.index)
+    else:
+        print("Skipped Mann-Whitney U-test")
+        fdr = ""
 
     summary["u-value"] = uval
     summary["p-value"] = pval
@@ -224,24 +240,25 @@
 
 
 # %% Graph drawing
-figsize = (int(args.figsize.split(',')[0]), int(args.figsize.split(',')[1]))
-# g = sns.clustermap(summary, method = "ward", metric="euclidean", figsize=figsize)
-# g.savefig("clustermap.pdf", bbox_inches='tight')
-
-
-with sns.axes_style("white"):
-    plt.figure(figsize=figsize)
-    ax = scattermap(df_rate, square=True, marker_size=df_cov, cmap='viridis_r',
-    cbar_kws={'label': 'v/h rate'})
-
-    #make a legend:
-    pws = [20, 40, 60, 80, 100]
-    for pw in pws:
-        plt.scatter([], [], s=(pw), c="k",label=str(pw))
-
-    h, l = plt.gca().get_legend_handles_labels()
-    plt.legend(h[1:], l[1:], labelspacing=.3, title="coverage(%)", borderpad=0, framealpha=0, edgecolor="w",
-                bbox_to_anchor=(1.1, -.1), ncol=1, loc='upper left', borderaxespad=0)
-    plt.savefig('scattermap.pdf' , bbox_inches='tight')
+if summary.shape[0] > 0:
+    figsize = (int(args.figsize.split(',')[0]), int(args.figsize.split(',')[1]))
+    # g = sns.clustermap(summary, method = "ward", metric="euclidean", figsize=figsize)
+    # g.savefig("clustermap.pdf", bbox_inches='tight')
+
+
+    with sns.axes_style("white"):
+        plt.figure(figsize=figsize)
+        ax = scattermap(df_rate, square=True, marker_size=df_cov, cmap='viridis_r',
+        cbar_kws={'label': 'v/h rate'})
+
+        #make a legend:
+        pws = [20, 40, 60, 80, 100]
+        for pw in pws:
+            plt.scatter([], [], s=(pw), c="k",label=str(pw))
+
+        h, l = plt.gca().get_legend_handles_labels()
+        plt.legend(h[1:], l[1:], labelspacing=.3, title="coverage(%)", borderpad=0, framealpha=0, edgecolor="w",
+                    bbox_to_anchor=(1.1, -.1), ncol=1, loc='upper left', borderaxespad=0)
+        plt.savefig('scattermap.pdf' , bbox_inches='tight')
 
 print('All processes succeeded.')

wrapper/scattermap.py

@@ -5,7 +5,7 @@
 import numpy as np
 import pandas as pd
 
-from seaborn.external.six import string_types
+# from seaborn.external.six import string_types
 
 from seaborn.utils import despine, axis_ticklabels_overlap, relative_luminance, to_utf8
 
@@ -67,12 +67,12 @@ def plot(self, ax, cax, kws):
                 cb.solids.set_rasterized(True)
 
         # Add row and column labels
-        if isinstance(self.xticks, string_types) and self.xticks == "auto":
+        if isinstance(self.xticks, str) and self.xticks == "auto":
             xticks, xticklabels = self._auto_ticks(ax, self.xticklabels, 0)
         else:
             xticks, xticklabels = self.xticks, self.xticklabels
 
-        if isinstance(self.yticks, string_types) and self.yticks == "auto":
+        if isinstance(self.yticks, str) and self.yticks == "auto":
             yticks, yticklabels = self._auto_ticks(ax, self.yticklabels, 1)
         else:
             yticks, yticklabels = self.yticks, self.yticklabels