changeset 3:6853676d2bae draft

"planemo upload for repository https://github.com/galaxyproject/tools-iuc/tree/master/tools/vsnp commit 92f46d4bb55b582f05ac3c4b094307f114cbf98f"
author iuc
date Fri, 27 Aug 2021 11:45:39 +0000
parents 9b4011d5e8cb
children e3016c6c5994
files test-data/13-1941-6_S4_L001_R1_600000.fastq.gz test-data/13-1941-6_S4_L001_R2_600000.fastq.gz test-data/add_zc_metrics.tabular test-data/add_zc_metrics1.tabular test-data/add_zc_metrics2.tabular test-data/add_zc_metrics3.tabular test-data/add_zc_metrics4.tabular test-data/add_zc_metrics5.tabular test-data/samtools_idxstats.tabular test-data/samtools_idxstats1.tabular test-data/samtools_idxstats2.tabular test-data/samtools_idxstats3.tabular test-data/samtools_idxstats4.tabular test-data/samtools_idxstats5.tabular test-data/vsnp_statistics1.tabular test-data/vsnp_statistics2.tabular test-data/vsnp_statistics3.tabular test-data/vsnp_statistics4.tabular vsnp_determine_ref_from_data.xml vsnp_statistics.py
diffstat 20 files changed, 237 insertions(+), 3 deletions(-) [+]
line wrap: on
line diff
Binary file test-data/13-1941-6_S4_L001_R1_600000.fastq.gz has changed
Binary file test-data/13-1941-6_S4_L001_R2_600000.fastq.gz has changed
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/test-data/add_zc_metrics.tabular	Fri Aug 27 11:45:39 2021 +0000
@@ -0,0 +1,3 @@
+# File	Number of Good SNPs	Average Coverage	Genome Coverage
+MarkDuplicates on data 4_ MarkDuplicates BAM output		10.338671	98.74%
+VCFfilter_ on data 7	611		
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/test-data/add_zc_metrics1.tabular	Fri Aug 27 11:45:39 2021 +0000
@@ -0,0 +1,3 @@
+# File	Number of Good SNPs	Average Coverage	Genome Coverage
+Mcap_Deer_DE_SRR650221_fastq_gz		0.439436	8.27%
+Mcap_Deer_DE_SRR650221_fastq_gz	36		
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/test-data/add_zc_metrics2.tabular	Fri Aug 27 11:45:39 2021 +0000
@@ -0,0 +1,3 @@
+# File	Number of Good SNPs	Average Coverage	Genome Coverage
+MarkDuplicates on data 4_ MarkDuplicates BAM output		10.338671	98.74%
+VCFfilter_ on data 7	611		
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/test-data/add_zc_metrics3.tabular	Fri Aug 27 11:45:39 2021 +0000
@@ -0,0 +1,3 @@
+# File	Number of Good SNPs	Average Coverage	Genome Coverage
+13-1941-6_S4_L001_R1_600000_fastq_gz		0.001252	0.13%
+13-1941-6_S4_L001_R1_600000_fastq_gz	0		
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/test-data/add_zc_metrics4.tabular	Fri Aug 27 11:45:39 2021 +0000
@@ -0,0 +1,3 @@
+# File	Number of Good SNPs	Average Coverage	Genome Coverage
+Mcap_Deer_DE_SRR650221_fastq_gz		0.439436	8.27%
+Mcap_Deer_DE_SRR650221_fastq_gz	36		
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/test-data/add_zc_metrics5.tabular	Fri Aug 27 11:45:39 2021 +0000
@@ -0,0 +1,3 @@
+# File	Number of Good SNPs	Average Coverage	Genome Coverage
+13-1941-6_S4_L001_600000_fastq		0.002146	0.16%
+13-1941-6_S4_L001_600000_fastq	0		
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/test-data/samtools_idxstats.tabular	Fri Aug 27 11:45:39 2021 +0000
@@ -0,0 +1,2 @@
+NC_002945.4	4349904	45	4047
+*	0	0	5
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/test-data/samtools_idxstats1.tabular	Fri Aug 27 11:45:39 2021 +0000
@@ -0,0 +1,2 @@
+NC_002945.4	4349904	17063	0
+*	0	0	223
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/test-data/samtools_idxstats2.tabular	Fri Aug 27 11:45:39 2021 +0000
@@ -0,0 +1,2 @@
+NC_002945.4	4349904	45	4047
+*	0	0	5
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/test-data/samtools_idxstats3.tabular	Fri Aug 27 11:45:39 2021 +0000
@@ -0,0 +1,2 @@
+NC_002945.4	4349904	24	0
+*	0	0	2
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/test-data/samtools_idxstats4.tabular	Fri Aug 27 11:45:39 2021 +0000
@@ -0,0 +1,2 @@
+NC_002945.4	4349904	17063	0
+*	0	0	223
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/test-data/samtools_idxstats5.tabular	Fri Aug 27 11:45:39 2021 +0000
@@ -0,0 +1,2 @@
+NC_002945.4	4349904	46	2
+*	0	0	4
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/test-data/vsnp_statistics1.tabular	Fri Aug 27 11:45:39 2021 +0000
@@ -0,0 +1,2 @@
+	Reference	File Size	Mean Read Length	Mean Read Quality	Reads Passing Q30	Total Reads	All Mapped Reads	Unmapped Reads	Unmapped Reads Percentage of Total	Reference with Coverage	Average Depth of Coverage	Good SNP Count
+Mcap_Deer_DE_SRR650221_fastq_gz	89	1.6 MB	121.0	29.7	      0.53	4317	17063	223	      0.05	8.27%	0.439436	36
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/test-data/vsnp_statistics2.tabular	Fri Aug 27 11:45:39 2021 +0000
@@ -0,0 +1,3 @@
+	Reference	File Size	Mean Read Length	Mean Read Quality	Reads Passing Q30	Total Reads	All Mapped Reads	Unmapped Reads	Unmapped Reads Percentage of Total	Reference with Coverage	Average Depth of Coverage	Good SNP Count
+13-1941-6_S4_L001_R1_600000_fastq_gz	89	8.7 KB	100.0	65.7	      1.00	25	45	5	      0.20	98.74%	10.338671	611
+13-1941-6_S4_L001_R2_600000_fastq_gz	89	8.5 KB	100.0	66.3	      1.00	25	45	5	      0.20	98.74%	10.338671	611
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/test-data/vsnp_statistics3.tabular	Fri Aug 27 11:45:39 2021 +0000
@@ -0,0 +1,3 @@
+	Reference	File Size	Mean Read Length	Mean Read Quality	Reads Passing Q30	Total Reads	All Mapped Reads	Unmapped Reads	Unmapped Reads Percentage of Total	Reference with Coverage	Average Depth of Coverage	Good SNP Count
+13-1941-6_S4_L001_R1_600000_fastq_gz	89	8.7 KB	100.0	65.7	      1.00	25	24	2	      0.08	0.13%	0.001252	0
+Mcap_Deer_DE_SRR650221_fastq_gz	89	1.6 MB	121.0	29.7	      0.53	4317	17063	223	      0.05	8.27%	0.439436	36
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/test-data/vsnp_statistics4.tabular	Fri Aug 27 11:45:39 2021 +0000
@@ -0,0 +1,3 @@
+	Reference	File Size	Mean Read Length	Mean Read Quality	Reads Passing Q30	Total Reads	All Mapped Reads	Unmapped Reads	Unmapped Reads Percentage of Total	Reference with Coverage	Average Depth of Coverage	Good SNP Count
+13-1941-6_S4_L001_R1_600000_fastq_gz	89	8.7 KB	100.0	65.7	      1.00	25	46	4	      0.16	0.16%	0.002146	0
+13-1941-6_S4_L001_R2_600000_fastq_gz	89	8.5 KB	100.0	66.3	      1.00	25	46	4	      0.16	0.16%	0.002146	0
--- a/vsnp_determine_ref_from_data.xml	Thu Jan 14 09:11:43 2021 +0000
+++ b/vsnp_determine_ref_from_data.xml	Fri Aug 27 11:45:39 2021 +0000
@@ -1,11 +1,11 @@
-<tool id="vsnp_determine_ref_from_data" name="vSNP: determine reference" version="@WRAPPER_VERSION@.1" profile="@PROFILE@">
+<tool id="vsnp_determine_ref_from_data" name="vSNP: determine reference" version="@WRAPPER_VERSION@.1+galaxy0" profile="@PROFILE@">
     <description>from input data</description>
     <macros>
         <import>macros.xml</import>
     </macros>
     <requirements>
-        <requirement type="package" version="1.76">biopython</requirement>
-        <requirement type="package" version="5.3">pyyaml</requirement>
+        <requirement type="package" version="1.79">biopython</requirement>
+        <requirement type="package" version="5.4.1">pyyaml</requirement>
     </requirements>
     <command detect_errors="exit_code"><![CDATA[
 #import re
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/vsnp_statistics.py	Fri Aug 27 11:45:39 2021 +0000
@@ -0,0 +1,193 @@
+#!/usr/bin/env python
+
+import argparse
+import csv
+import gzip
+import os
+from functools import partial
+
+import numpy
+import pandas
+from Bio import SeqIO
+
+
+def nice_size(size):
+    # Returns a readably formatted string with the size
+    words = ['bytes', 'KB', 'MB', 'GB', 'TB', 'PB', 'EB']
+    prefix = ''
+    try:
+        size = float(size)
+        if size < 0:
+            size = abs(size)
+            prefix = '-'
+    except Exception:
+        return '??? bytes'
+    for ind, word in enumerate(words):
+        step = 1024 ** (ind + 1)
+        if step > size:
+            size = size / float(1024 ** ind)
+            if word == 'bytes':  # No decimals for bytes
+                return "%s%d bytes" % (prefix, size)
+            return "%s%.1f %s" % (prefix, size, word)
+    return '??? bytes'
+
+
+def output_statistics(fastq_files, idxstats_files, metrics_files, output_file, gzipped, dbkey):
+    # Produce an Excel spreadsheet that
+    # contains a row for each sample.
+    columns = ['Reference', 'File Size', 'Mean Read Length', 'Mean Read Quality', 'Reads Passing Q30',
+               'Total Reads', 'All Mapped Reads', 'Unmapped Reads', 'Unmapped Reads Percentage of Total',
+               'Reference with Coverage', 'Average Depth of Coverage', 'Good SNP Count']
+    data_frames = []
+    for i, fastq_file in enumerate(fastq_files):
+        idxstats_file = idxstats_files[i]
+        metrics_file = metrics_files[i]
+        file_name_base = os.path.basename(fastq_file)
+        # Read fastq_file into a data frame.
+        _open = partial(gzip.open, mode='rt') if gzipped else open
+        with _open(fastq_file) as fh:
+            identifiers = []
+            seqs = []
+            letter_annotations = []
+            for seq_record in SeqIO.parse(fh, "fastq"):
+                identifiers.append(seq_record.id)
+                seqs.append(seq_record.seq)
+                letter_annotations.append(seq_record.letter_annotations["phred_quality"])
+        # Convert lists to Pandas series.
+        s1 = pandas.Series(identifiers, name='id')
+        s2 = pandas.Series(seqs, name='seq')
+        # Gather Series into a data frame.
+        fastq_df = pandas.DataFrame(dict(id=s1, seq=s2)).set_index(['id'])
+        total_reads = int(len(fastq_df.index) / 4)
+        current_sample_df = pandas.DataFrame(index=[file_name_base], columns=columns)
+        # Reference
+        current_sample_df.at[file_name_base, 'Reference'] = dbkey
+        # File Size
+        current_sample_df.at[file_name_base, 'File Size'] = nice_size(os.path.getsize(fastq_file))
+        # Mean Read Length
+        sampling_size = 10000
+        if sampling_size > total_reads:
+            sampling_size = total_reads
+        fastq_df = fastq_df.iloc[3::4].sample(sampling_size)
+        dict_mean = {}
+        list_length = []
+        i = 0
+        for id, seq, in fastq_df.iterrows():
+            dict_mean[id] = numpy.mean(letter_annotations[i])
+            list_length.append(len(seq.array[0]))
+            i += 1
+        current_sample_df.at[file_name_base, 'Mean Read Length'] = '%.1f' % numpy.mean(list_length)
+        # Mean Read Quality
+        df_mean = pandas.DataFrame.from_dict(dict_mean, orient='index', columns=['ave'])
+        current_sample_df.at[file_name_base, 'Mean Read Quality'] = '%.1f' % df_mean['ave'].mean()
+        # Reads Passing Q30
+        reads_gt_q30 = len(df_mean[df_mean['ave'] >= 30])
+        reads_passing_q30 = '{:10.2f}'.format(reads_gt_q30 / sampling_size)
+        current_sample_df.at[file_name_base, 'Reads Passing Q30'] = reads_passing_q30
+        # Total Reads
+        current_sample_df.at[file_name_base, 'Total Reads'] = total_reads
+        # All Mapped Reads
+        all_mapped_reads, unmapped_reads = process_idxstats_file(idxstats_file)
+        current_sample_df.at[file_name_base, 'All Mapped Reads'] = all_mapped_reads
+        # Unmapped Reads
+        current_sample_df.at[file_name_base, 'Unmapped Reads'] = unmapped_reads
+        # Unmapped Reads Percentage of Total
+        if unmapped_reads > 0:
+            unmapped_reads_percentage = '{:10.2f}'.format(unmapped_reads / total_reads)
+        else:
+            unmapped_reads_percentage = 0
+        current_sample_df.at[file_name_base, 'Unmapped Reads Percentage of Total'] = unmapped_reads_percentage
+        # Reference with Coverage
+        ref_with_coverage, avg_depth_of_coverage, good_snp_count = process_metrics_file(metrics_file)
+        current_sample_df.at[file_name_base, 'Reference with Coverage'] = ref_with_coverage
+        # Average Depth of Coverage
+        current_sample_df.at[file_name_base, 'Average Depth of Coverage'] = avg_depth_of_coverage
+        # Good SNP Count
+        current_sample_df.at[file_name_base, 'Good SNP Count'] = good_snp_count
+        data_frames.append(current_sample_df)
+    output_df = pandas.concat(data_frames)
+    output_df.to_csv(output_file, sep='\t', quoting=csv.QUOTE_NONE, escapechar='\\')
+
+
+def process_idxstats_file(idxstats_file):
+    all_mapped_reads = 0
+    unmapped_reads = 0
+    with open(idxstats_file, "r") as fh:
+        for i, line in enumerate(fh):
+            line = line.rstrip('\r\n')
+            items = line.split("\t")
+            if i == 0:
+                # NC_002945.4 4349904 213570 4047
+                all_mapped_reads = int(items[2])
+            elif i == 1:
+                # * 0 0 82774
+                unmapped_reads = int(items[3])
+    return all_mapped_reads, unmapped_reads
+
+
+def process_metrics_file(metrics_file):
+    ref_with_coverage = '0%'
+    avg_depth_of_coverage = 0
+    good_snp_count = 0
+    with open(metrics_file, "r") as ifh:
+        for i, line in enumerate(ifh):
+            if i == 0:
+                # Skip comments.
+                continue
+            line = line.rstrip('\r\n')
+            items = line.split("\t")
+            if i == 1:
+                # MarkDuplicates 10.338671 98.74%
+                ref_with_coverage = items[3]
+                avg_depth_of_coverage = items[2]
+            elif i == 2:
+                # VCFfilter 611
+                good_snp_count = items[1]
+    return ref_with_coverage, avg_depth_of_coverage, good_snp_count
+
+
+parser = argparse.ArgumentParser()
+
+parser.add_argument('--dbkey', action='store', dest='dbkey', help='Reference dbkey')
+parser.add_argument('--gzipped', action='store_true', dest='gzipped', required=False, default=False, help='Input files are gzipped')
+parser.add_argument('--input_idxstats_dir', action='store', dest='input_idxstats_dir', required=False, default=None, help='Samtools idxstats input directory')
+parser.add_argument('--input_metrics_dir', action='store', dest='input_metrics_dir', required=False, default=None, help='vSNP add zero coverage metrics input directory')
+parser.add_argument('--input_reads_dir', action='store', dest='input_reads_dir', required=False, default=None, help='Samples input directory')
+parser.add_argument('--list_paired', action='store_true', dest='list_paired', required=False, default=False, help='Input samples is a list of paired reads')
+parser.add_argument('--output', action='store', dest='output', help='Output Excel statistics file')
+parser.add_argument('--read1', action='store', dest='read1', help='Required: single read')
+parser.add_argument('--read2', action='store', dest='read2', required=False, default=None, help='Optional: paired read')
+parser.add_argument('--samtools_idxstats', action='store', dest='samtools_idxstats', help='Output of samtools_idxstats')
+parser.add_argument('--vsnp_azc', action='store', dest='vsnp_azc', help='Output of vsnp_add_zero_coverage')
+
+args = parser.parse_args()
+
+fastq_files = []
+idxstats_files = []
+metrics_files = []
+# Accumulate inputs.
+if args.read1 is not None:
+    # The inputs are not dataset collections, so
+    # read1, read2 (possibly) and vsnp_azc will also
+    # not be None.
+    fastq_files.append(args.read1)
+    idxstats_files.append(args.samtools_idxstats)
+    metrics_files.append(args.vsnp_azc)
+    if args.read2 is not None:
+        fastq_files.append(args.read2)
+        idxstats_files.append(args.samtools_idxstats)
+        metrics_files.append(args.vsnp_azc)
+else:
+    for file_name in sorted(os.listdir(args.input_reads_dir)):
+        fastq_files.append(os.path.join(args.input_reads_dir, file_name))
+    for file_name in sorted(os.listdir(args.input_idxstats_dir)):
+        idxstats_files.append(os.path.join(args.input_idxstats_dir, file_name))
+        if args.list_paired:
+            # Add the idxstats file for reverse.
+            idxstats_files.append(os.path.join(args.input_idxstats_dir, file_name))
+    for file_name in sorted(os.listdir(args.input_metrics_dir)):
+        metrics_files.append(os.path.join(args.input_metrics_dir, file_name))
+        if args.list_paired:
+            # Add the metrics file for reverse.
+            metrics_files.append(os.path.join(args.input_metrics_dir, file_name))
+output_statistics(fastq_files, idxstats_files, metrics_files, args.output, args.gzipped, args.dbkey)