diff --git a/README.md b/README.md index dbdb85e..edd5c8f 100644 --- a/README.md +++ b/README.md @@ -1,4 +1,4 @@ -# VIRTUS2 : VIRal Transcript Usage Sensor v2.0 +# VIRTUS2 : VIRal Transcript Usage Sensor v2.0.1 **!!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)** diff --git a/test/singularity.test.sh b/test/singularity.test.sh index c497c57..5b099ab 100644 --- a/test/singularity.test.sh +++ b/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! diff --git a/test/test.sh b/test/test.sh index f8900da..a747c69 100644 --- a/test/test.sh +++ b/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! \ No newline at end of file +echo Successfully completed test.sh! diff --git a/test/udocker.test.sh b/test/udocker.test.sh index 970341c..bdf7e59 100644 --- a/test/udocker.test.sh +++ b/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! diff --git a/tool/samtools/bam_filter_polyx.cwl b/tool/samtools/bam_filter_polyx.cwl index 974a338..4ed99a2 100644 --- a/tool/samtools/bam_filter_polyx.cwl +++ b/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 diff --git a/wrapper/VIRTUS_wrapper.py b/wrapper/VIRTUS_wrapper.py index 1ca6ec1..5b03ed8 100755 --- a/wrapper/VIRTUS_wrapper.py +++ b/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.') diff --git a/wrapper/scattermap.py b/wrapper/scattermap.py index 13e1f54..05ffbe7 100644 --- a/wrapper/scattermap.py +++ b/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