Mercurial > repos > greg > vsnp_statistics
annotate vsnp_statistics.py @ 7:de2af65c4633 draft
Uploaded
| author | greg |
|---|---|
| date | Thu, 22 Jul 2021 18:05:22 +0000 |
| parents | d0fbdeaaa488 |
| children | 1becb6606626 |
| rev | line source |
|---|---|
| 0 | 1 #!/usr/bin/env python |
| 2 | |
| 3 import argparse | |
|
5
d0fbdeaaa488
"planemo upload for repository https://github.com/gregvonkuster/galaxy_tools/tree/master/tools/sequence_analysis/vsnp/vsnp_statistics commit 770e89322a15829580ed9577a853660f63233f32"
greg
parents:
4
diff
changeset
|
4 import csv |
| 0 | 5 import gzip |
| 6 import os | |
|
5
d0fbdeaaa488
"planemo upload for repository https://github.com/gregvonkuster/galaxy_tools/tree/master/tools/sequence_analysis/vsnp/vsnp_statistics commit 770e89322a15829580ed9577a853660f63233f32"
greg
parents:
4
diff
changeset
|
7 from functools import partial |
| 0 | 8 |
| 4 | 9 import numpy |
| 10 import pandas | |
|
5
d0fbdeaaa488
"planemo upload for repository https://github.com/gregvonkuster/galaxy_tools/tree/master/tools/sequence_analysis/vsnp/vsnp_statistics commit 770e89322a15829580ed9577a853660f63233f32"
greg
parents:
4
diff
changeset
|
11 from Bio import SeqIO |
| 0 | 12 |
| 13 | |
| 14 def nice_size(size): | |
| 15 # Returns a readably formatted string with the size | |
| 16 words = ['bytes', 'KB', 'MB', 'GB', 'TB', 'PB', 'EB'] | |
| 17 prefix = '' | |
| 18 try: | |
| 19 size = float(size) | |
| 20 if size < 0: | |
| 21 size = abs(size) | |
| 22 prefix = '-' | |
| 23 except Exception: | |
| 24 return '??? bytes' | |
| 25 for ind, word in enumerate(words): | |
| 26 step = 1024 ** (ind + 1) | |
| 27 if step > size: | |
| 28 size = size / float(1024 ** ind) | |
| 29 if word == 'bytes': # No decimals for bytes | |
| 30 return "%s%d bytes" % (prefix, size) | |
| 31 return "%s%.1f %s" % (prefix, size, word) | |
| 32 return '??? bytes' | |
| 33 | |
| 34 | |
| 4 | 35 def output_statistics(fastq_files, idxstats_files, metrics_files, output_file, gzipped, dbkey): |
| 1 | 36 # Produce an Excel spreadsheet that |
| 37 # contains a row for each sample. | |
| 38 columns = ['Reference', 'File Size', 'Mean Read Length', 'Mean Read Quality', 'Reads Passing Q30', | |
| 39 'Total Reads', 'All Mapped Reads', 'Unmapped Reads', 'Unmapped Reads Percentage of Total', | |
| 40 'Reference with Coverage', 'Average Depth of Coverage', 'Good SNP Count'] | |
| 41 data_frames = [] | |
| 4 | 42 for i, fastq_file in enumerate(fastq_files): |
| 1 | 43 idxstats_file = idxstats_files[i] |
| 44 metrics_file = metrics_files[i] | |
| 45 file_name_base = os.path.basename(fastq_file) | |
| 46 # Read fastq_file into a data frame. | |
|
5
d0fbdeaaa488
"planemo upload for repository https://github.com/gregvonkuster/galaxy_tools/tree/master/tools/sequence_analysis/vsnp/vsnp_statistics commit 770e89322a15829580ed9577a853660f63233f32"
greg
parents:
4
diff
changeset
|
47 _open = partial(gzip.open, mode='rt') if gzipped else open |
|
d0fbdeaaa488
"planemo upload for repository https://github.com/gregvonkuster/galaxy_tools/tree/master/tools/sequence_analysis/vsnp/vsnp_statistics commit 770e89322a15829580ed9577a853660f63233f32"
greg
parents:
4
diff
changeset
|
48 with _open(fastq_file) as fh: |
|
d0fbdeaaa488
"planemo upload for repository https://github.com/gregvonkuster/galaxy_tools/tree/master/tools/sequence_analysis/vsnp/vsnp_statistics commit 770e89322a15829580ed9577a853660f63233f32"
greg
parents:
4
diff
changeset
|
49 identifiers = [] |
|
d0fbdeaaa488
"planemo upload for repository https://github.com/gregvonkuster/galaxy_tools/tree/master/tools/sequence_analysis/vsnp/vsnp_statistics commit 770e89322a15829580ed9577a853660f63233f32"
greg
parents:
4
diff
changeset
|
50 seqs = [] |
|
d0fbdeaaa488
"planemo upload for repository https://github.com/gregvonkuster/galaxy_tools/tree/master/tools/sequence_analysis/vsnp/vsnp_statistics commit 770e89322a15829580ed9577a853660f63233f32"
greg
parents:
4
diff
changeset
|
51 letter_annotations = [] |
|
d0fbdeaaa488
"planemo upload for repository https://github.com/gregvonkuster/galaxy_tools/tree/master/tools/sequence_analysis/vsnp/vsnp_statistics commit 770e89322a15829580ed9577a853660f63233f32"
greg
parents:
4
diff
changeset
|
52 for seq_record in SeqIO.parse(fh, "fastq"): |
|
d0fbdeaaa488
"planemo upload for repository https://github.com/gregvonkuster/galaxy_tools/tree/master/tools/sequence_analysis/vsnp/vsnp_statistics commit 770e89322a15829580ed9577a853660f63233f32"
greg
parents:
4
diff
changeset
|
53 identifiers.append(seq_record.id) |
|
d0fbdeaaa488
"planemo upload for repository https://github.com/gregvonkuster/galaxy_tools/tree/master/tools/sequence_analysis/vsnp/vsnp_statistics commit 770e89322a15829580ed9577a853660f63233f32"
greg
parents:
4
diff
changeset
|
54 seqs.append(seq_record.seq) |
|
d0fbdeaaa488
"planemo upload for repository https://github.com/gregvonkuster/galaxy_tools/tree/master/tools/sequence_analysis/vsnp/vsnp_statistics commit 770e89322a15829580ed9577a853660f63233f32"
greg
parents:
4
diff
changeset
|
55 letter_annotations.append(seq_record.letter_annotations["phred_quality"]) |
|
d0fbdeaaa488
"planemo upload for repository https://github.com/gregvonkuster/galaxy_tools/tree/master/tools/sequence_analysis/vsnp/vsnp_statistics commit 770e89322a15829580ed9577a853660f63233f32"
greg
parents:
4
diff
changeset
|
56 # Convert lists to Pandas series. |
|
d0fbdeaaa488
"planemo upload for repository https://github.com/gregvonkuster/galaxy_tools/tree/master/tools/sequence_analysis/vsnp/vsnp_statistics commit 770e89322a15829580ed9577a853660f63233f32"
greg
parents:
4
diff
changeset
|
57 s1 = pandas.Series(identifiers, name='id') |
|
d0fbdeaaa488
"planemo upload for repository https://github.com/gregvonkuster/galaxy_tools/tree/master/tools/sequence_analysis/vsnp/vsnp_statistics commit 770e89322a15829580ed9577a853660f63233f32"
greg
parents:
4
diff
changeset
|
58 s2 = pandas.Series(seqs, name='seq') |
|
d0fbdeaaa488
"planemo upload for repository https://github.com/gregvonkuster/galaxy_tools/tree/master/tools/sequence_analysis/vsnp/vsnp_statistics commit 770e89322a15829580ed9577a853660f63233f32"
greg
parents:
4
diff
changeset
|
59 # Gather Series into a data frame. |
|
d0fbdeaaa488
"planemo upload for repository https://github.com/gregvonkuster/galaxy_tools/tree/master/tools/sequence_analysis/vsnp/vsnp_statistics commit 770e89322a15829580ed9577a853660f63233f32"
greg
parents:
4
diff
changeset
|
60 fastq_df = pandas.DataFrame(dict(id=s1, seq=s2)).set_index(['id']) |
| 1 | 61 total_reads = int(len(fastq_df.index) / 4) |
| 62 current_sample_df = pandas.DataFrame(index=[file_name_base], columns=columns) | |
| 63 # Reference | |
| 64 current_sample_df.at[file_name_base, 'Reference'] = dbkey | |
| 65 # File Size | |
| 66 current_sample_df.at[file_name_base, 'File Size'] = nice_size(os.path.getsize(fastq_file)) | |
| 67 # Mean Read Length | |
| 68 sampling_size = 10000 | |
| 69 if sampling_size > total_reads: | |
| 70 sampling_size = total_reads | |
| 71 fastq_df = fastq_df.iloc[3::4].sample(sampling_size) | |
| 72 dict_mean = {} | |
| 73 list_length = [] | |
|
5
d0fbdeaaa488
"planemo upload for repository https://github.com/gregvonkuster/galaxy_tools/tree/master/tools/sequence_analysis/vsnp/vsnp_statistics commit 770e89322a15829580ed9577a853660f63233f32"
greg
parents:
4
diff
changeset
|
74 i = 0 |
|
d0fbdeaaa488
"planemo upload for repository https://github.com/gregvonkuster/galaxy_tools/tree/master/tools/sequence_analysis/vsnp/vsnp_statistics commit 770e89322a15829580ed9577a853660f63233f32"
greg
parents:
4
diff
changeset
|
75 for id, seq, in fastq_df.iterrows(): |
|
d0fbdeaaa488
"planemo upload for repository https://github.com/gregvonkuster/galaxy_tools/tree/master/tools/sequence_analysis/vsnp/vsnp_statistics commit 770e89322a15829580ed9577a853660f63233f32"
greg
parents:
4
diff
changeset
|
76 dict_mean[id] = numpy.mean(letter_annotations[i]) |
|
d0fbdeaaa488
"planemo upload for repository https://github.com/gregvonkuster/galaxy_tools/tree/master/tools/sequence_analysis/vsnp/vsnp_statistics commit 770e89322a15829580ed9577a853660f63233f32"
greg
parents:
4
diff
changeset
|
77 list_length.append(len(seq.array[0])) |
|
d0fbdeaaa488
"planemo upload for repository https://github.com/gregvonkuster/galaxy_tools/tree/master/tools/sequence_analysis/vsnp/vsnp_statistics commit 770e89322a15829580ed9577a853660f63233f32"
greg
parents:
4
diff
changeset
|
78 i += 1 |
|
d0fbdeaaa488
"planemo upload for repository https://github.com/gregvonkuster/galaxy_tools/tree/master/tools/sequence_analysis/vsnp/vsnp_statistics commit 770e89322a15829580ed9577a853660f63233f32"
greg
parents:
4
diff
changeset
|
79 current_sample_df.at[file_name_base, 'Mean Read Length'] = '%.1f' % numpy.mean(list_length) |
| 1 | 80 # Mean Read Quality |
| 81 df_mean = pandas.DataFrame.from_dict(dict_mean, orient='index', columns=['ave']) | |
|
5
d0fbdeaaa488
"planemo upload for repository https://github.com/gregvonkuster/galaxy_tools/tree/master/tools/sequence_analysis/vsnp/vsnp_statistics commit 770e89322a15829580ed9577a853660f63233f32"
greg
parents:
4
diff
changeset
|
82 current_sample_df.at[file_name_base, 'Mean Read Quality'] = '%.1f' % df_mean['ave'].mean() |
| 1 | 83 # Reads Passing Q30 |
| 84 reads_gt_q30 = len(df_mean[df_mean['ave'] >= 30]) | |
|
5
d0fbdeaaa488
"planemo upload for repository https://github.com/gregvonkuster/galaxy_tools/tree/master/tools/sequence_analysis/vsnp/vsnp_statistics commit 770e89322a15829580ed9577a853660f63233f32"
greg
parents:
4
diff
changeset
|
85 reads_passing_q30 = '{:10.2f}'.format(reads_gt_q30 / sampling_size) |
| 1 | 86 current_sample_df.at[file_name_base, 'Reads Passing Q30'] = reads_passing_q30 |
| 87 # Total Reads | |
| 88 current_sample_df.at[file_name_base, 'Total Reads'] = total_reads | |
| 89 # All Mapped Reads | |
| 90 all_mapped_reads, unmapped_reads = process_idxstats_file(idxstats_file) | |
| 91 current_sample_df.at[file_name_base, 'All Mapped Reads'] = all_mapped_reads | |
| 92 # Unmapped Reads | |
| 93 current_sample_df.at[file_name_base, 'Unmapped Reads'] = unmapped_reads | |
| 94 # Unmapped Reads Percentage of Total | |
| 95 if unmapped_reads > 0: | |
|
5
d0fbdeaaa488
"planemo upload for repository https://github.com/gregvonkuster/galaxy_tools/tree/master/tools/sequence_analysis/vsnp/vsnp_statistics commit 770e89322a15829580ed9577a853660f63233f32"
greg
parents:
4
diff
changeset
|
96 unmapped_reads_percentage = '{:10.2f}'.format(unmapped_reads / total_reads) |
| 0 | 97 else: |
| 1 | 98 unmapped_reads_percentage = 0 |
| 99 current_sample_df.at[file_name_base, 'Unmapped Reads Percentage of Total'] = unmapped_reads_percentage | |
| 100 # Reference with Coverage | |
| 101 ref_with_coverage, avg_depth_of_coverage, good_snp_count = process_metrics_file(metrics_file) | |
| 102 current_sample_df.at[file_name_base, 'Reference with Coverage'] = ref_with_coverage | |
| 103 # Average Depth of Coverage | |
| 104 current_sample_df.at[file_name_base, 'Average Depth of Coverage'] = avg_depth_of_coverage | |
| 105 # Good SNP Count | |
| 106 current_sample_df.at[file_name_base, 'Good SNP Count'] = good_snp_count | |
| 107 data_frames.append(current_sample_df) | |
|
5
d0fbdeaaa488
"planemo upload for repository https://github.com/gregvonkuster/galaxy_tools/tree/master/tools/sequence_analysis/vsnp/vsnp_statistics commit 770e89322a15829580ed9577a853660f63233f32"
greg
parents:
4
diff
changeset
|
108 output_df = pandas.concat(data_frames) |
|
d0fbdeaaa488
"planemo upload for repository https://github.com/gregvonkuster/galaxy_tools/tree/master/tools/sequence_analysis/vsnp/vsnp_statistics commit 770e89322a15829580ed9577a853660f63233f32"
greg
parents:
4
diff
changeset
|
109 output_df.to_csv(output_file, sep='\t', quoting=csv.QUOTE_NONE, escapechar='\\') |
| 0 | 110 |
| 111 | |
| 1 | 112 def process_idxstats_file(idxstats_file): |
| 113 all_mapped_reads = 0 | |
| 114 unmapped_reads = 0 | |
| 115 with open(idxstats_file, "r") as fh: | |
| 116 for i, line in enumerate(fh): | |
|
5
d0fbdeaaa488
"planemo upload for repository https://github.com/gregvonkuster/galaxy_tools/tree/master/tools/sequence_analysis/vsnp/vsnp_statistics commit 770e89322a15829580ed9577a853660f63233f32"
greg
parents:
4
diff
changeset
|
117 line = line.rstrip('\r\n') |
| 1 | 118 items = line.split("\t") |
| 119 if i == 0: | |
| 120 # NC_002945.4 4349904 213570 4047 | |
| 121 all_mapped_reads = int(items[2]) | |
| 122 elif i == 1: | |
| 123 # * 0 0 82774 | |
| 124 unmapped_reads = int(items[3]) | |
| 125 return all_mapped_reads, unmapped_reads | |
| 0 | 126 |
| 127 | |
| 1 | 128 def process_metrics_file(metrics_file): |
| 129 ref_with_coverage = '0%' | |
| 130 avg_depth_of_coverage = 0 | |
| 131 good_snp_count = 0 | |
| 132 with open(metrics_file, "r") as ifh: | |
| 133 for i, line in enumerate(ifh): | |
| 134 if i == 0: | |
| 135 # Skip comments. | |
| 136 continue | |
|
5
d0fbdeaaa488
"planemo upload for repository https://github.com/gregvonkuster/galaxy_tools/tree/master/tools/sequence_analysis/vsnp/vsnp_statistics commit 770e89322a15829580ed9577a853660f63233f32"
greg
parents:
4
diff
changeset
|
137 line = line.rstrip('\r\n') |
| 1 | 138 items = line.split("\t") |
| 139 if i == 1: | |
| 140 # MarkDuplicates 10.338671 98.74% | |
| 141 ref_with_coverage = items[3] | |
| 142 avg_depth_of_coverage = items[2] | |
| 143 elif i == 2: | |
| 144 # VCFfilter 611 | |
| 145 good_snp_count = items[1] | |
| 146 return ref_with_coverage, avg_depth_of_coverage, good_snp_count | |
| 0 | 147 |
| 148 | |
| 4 | 149 parser = argparse.ArgumentParser() |
| 0 | 150 |
| 4 | 151 parser.add_argument('--dbkey', action='store', dest='dbkey', help='Reference dbkey') |
| 152 parser.add_argument('--gzipped', action='store_true', dest='gzipped', required=False, default=False, help='Input files are gzipped') | |
| 153 parser.add_argument('--input_idxstats_dir', action='store', dest='input_idxstats_dir', required=False, default=None, help='Samtools idxstats input directory') | |
| 154 parser.add_argument('--input_metrics_dir', action='store', dest='input_metrics_dir', required=False, default=None, help='vSNP add zero coverage metrics input directory') | |
| 155 parser.add_argument('--input_reads_dir', action='store', dest='input_reads_dir', required=False, default=None, help='Samples input directory') | |
| 156 parser.add_argument('--list_paired', action='store_true', dest='list_paired', required=False, default=False, help='Input samples is a list of paired reads') | |
| 157 parser.add_argument('--output', action='store', dest='output', help='Output Excel statistics file') | |
| 158 parser.add_argument('--read1', action='store', dest='read1', help='Required: single read') | |
| 159 parser.add_argument('--read2', action='store', dest='read2', required=False, default=None, help='Optional: paired read') | |
| 160 parser.add_argument('--samtools_idxstats', action='store', dest='samtools_idxstats', help='Output of samtools_idxstats') | |
| 161 parser.add_argument('--vsnp_azc', action='store', dest='vsnp_azc', help='Output of vsnp_add_zero_coverage') | |
| 0 | 162 |
| 4 | 163 args = parser.parse_args() |
| 0 | 164 |
| 4 | 165 fastq_files = [] |
| 166 idxstats_files = [] | |
| 167 metrics_files = [] | |
| 168 # Accumulate inputs. | |
| 169 if args.read1 is not None: | |
| 170 # The inputs are not dataset collections, so | |
| 171 # read1, read2 (possibly) and vsnp_azc will also | |
| 172 # not be None. | |
| 173 fastq_files.append(args.read1) | |
| 174 idxstats_files.append(args.samtools_idxstats) | |
| 175 metrics_files.append(args.vsnp_azc) | |
| 176 if args.read2 is not None: | |
| 177 fastq_files.append(args.read2) | |
| 0 | 178 idxstats_files.append(args.samtools_idxstats) |
| 179 metrics_files.append(args.vsnp_azc) | |
| 4 | 180 else: |
| 181 for file_name in sorted(os.listdir(args.input_reads_dir)): | |
| 182 fastq_files.append(os.path.join(args.input_reads_dir, file_name)) | |
| 183 for file_name in sorted(os.listdir(args.input_idxstats_dir)): | |
| 184 idxstats_files.append(os.path.join(args.input_idxstats_dir, file_name)) | |
| 185 if args.list_paired: | |
| 186 # Add the idxstats file for reverse. | |
| 187 idxstats_files.append(os.path.join(args.input_idxstats_dir, file_name)) | |
| 188 for file_name in sorted(os.listdir(args.input_metrics_dir)): | |
| 189 metrics_files.append(os.path.join(args.input_metrics_dir, file_name)) | |
| 190 if args.list_paired: | |
| 191 # Add the metrics file for reverse. | |
| 192 metrics_files.append(os.path.join(args.input_metrics_dir, file_name)) | |
| 193 output_statistics(fastq_files, idxstats_files, metrics_files, args.output, args.gzipped, args.dbkey) |