Mercurial > repos > greg > vsnp_statistics
annotate vsnp_statistics.py @ 28:e219f5c06dc0 draft default tip
"planemo upload for repository https://github.com/gregvonkuster/galaxy_tools/tree/master/tools/sequence_analysis/vsnp/vsnp_statistics commit fd6150d1cb5a9e50aa10694b9ac85cb575ed7525"
| author | greg |
|---|---|
| date | Mon, 15 Nov 2021 23:24:03 +0000 |
| parents | b908bb18008a |
| children |
| rev | line source |
|---|---|
| 0 | 1 #!/usr/bin/env python |
| 2 | |
| 3 import argparse | |
| 4 import gzip | |
| 5 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
|
6 from functools import partial |
| 0 | 7 |
| 4 | 8 import numpy |
| 9 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
|
10 from Bio import SeqIO |
| 0 | 11 |
| 12 | |
| 8 | 13 class Statistics: |
| 14 | |
| 15 def __init__(self, reference, fastq_file, file_size, total_reads, mean_read_length, mean_read_quality, reads_passing_q30): | |
| 16 self.reference = reference | |
| 17 self.fastq_file = fastq_file | |
| 18 self.file_size = file_size | |
| 19 self.total_reads = total_reads | |
| 20 self.mean_read_length = mean_read_length | |
| 21 self.mean_read_quality = mean_read_quality | |
| 22 self.reads_passing_q30 = reads_passing_q30 | |
| 23 | |
| 24 | |
| 0 | 25 def nice_size(size): |
| 26 # Returns a readably formatted string with the size | |
| 27 words = ['bytes', 'KB', 'MB', 'GB', 'TB', 'PB', 'EB'] | |
| 28 prefix = '' | |
| 29 try: | |
| 30 size = float(size) | |
| 31 if size < 0: | |
| 32 size = abs(size) | |
| 33 prefix = '-' | |
| 34 except Exception: | |
| 35 return '??? bytes' | |
| 36 for ind, word in enumerate(words): | |
| 37 step = 1024 ** (ind + 1) | |
| 38 if step > size: | |
| 39 size = size / float(1024 ** ind) | |
| 40 if word == 'bytes': # No decimals for bytes | |
| 41 return "%s%d bytes" % (prefix, size) | |
| 42 return "%s%.1f %s" % (prefix, size, word) | |
| 43 return '??? bytes' | |
| 44 | |
| 45 | |
| 8 | 46 def get_statistics(dbkey, fastq_file, gzipped): |
| 47 sampling_size = 10000 | |
| 48 # Read fastq_file into a data fram to | |
| 49 # get the phred quality scores. | |
| 50 _open = partial(gzip.open, mode='rt') if gzipped else open | |
| 51 with _open(fastq_file) as fh: | |
| 52 identifiers = [] | |
| 53 seqs = [] | |
| 54 letter_annotations = [] | |
| 55 for seq_record in SeqIO.parse(fh, "fastq"): | |
| 56 identifiers.append(seq_record.id) | |
| 57 seqs.append(seq_record.seq) | |
| 58 letter_annotations.append(seq_record.letter_annotations["phred_quality"]) | |
| 59 # Convert lists to Pandas series. | |
| 60 s1 = pandas.Series(identifiers, name='id') | |
| 61 s2 = pandas.Series(seqs, name='seq') | |
| 62 # Gather Series into a data frame. | |
| 63 fastq_df = pandas.DataFrame(dict(id=s1, seq=s2)).set_index(['id']) | |
| 64 # Starting at row 3, keep every 4 row | |
| 65 # random sample specified number of rows. | |
| 66 file_size = nice_size(os.path.getsize(fastq_file)) | |
| 25 | 67 total_reads = len(seqs) |
| 8 | 68 # Mean Read Length |
| 69 if sampling_size > total_reads: | |
| 70 sampling_size = total_reads | |
| 25 | 71 try: |
| 72 fastq_df = fastq_df.iloc[3::4].sample(sampling_size) | |
| 73 except ValueError: | |
| 74 fastq_df = fastq_df.iloc[3::4].sample(sampling_size, replace=True) | |
| 8 | 75 dict_mean = {} |
| 76 list_length = [] | |
| 77 i = 0 | |
| 78 for id, seq, in fastq_df.iterrows(): | |
| 79 dict_mean[id] = numpy.mean(letter_annotations[i]) | |
| 80 list_length.append(len(seq.array[0])) | |
| 81 i += 1 | |
| 82 mean_read_length = '%.1f' % numpy.mean(list_length) | |
| 83 # Mean Read Quality | |
| 84 df_mean = pandas.DataFrame.from_dict(dict_mean, orient='index', columns=['ave']) | |
| 85 mean_read_quality = '%.1f' % df_mean['ave'].mean() | |
| 86 # Reads Passing Q30 | |
| 87 reads_gt_q30 = len(df_mean[df_mean['ave'] >= 30]) | |
| 88 reads_passing_q30 = '{:10.2f}'.format(reads_gt_q30 / sampling_size) | |
| 89 stats = Statistics(dbkey, os.path.basename(fastq_file), file_size, total_reads, mean_read_length, | |
| 90 mean_read_quality, reads_passing_q30) | |
| 91 return stats | |
| 92 | |
| 93 | |
| 94 def accrue_statistics(dbkey, read1, read2, gzipped): | |
| 95 read1_stats = get_statistics(dbkey, read1, gzipped) | |
| 96 if read2 is None: | |
| 97 read2_stats = None | |
| 98 else: | |
| 99 read2_stats = get_statistics(dbkey, read2, gzipped) | |
| 100 return read1_stats, read2_stats | |
| 101 | |
| 102 | |
| 103 def output_statistics(read1_stats, read2_stats, idxstats_file, metrics_file, output_file): | |
| 104 paired_reads = read2_stats is not None | |
| 105 if paired_reads: | |
| 23 | 106 columns = ['Read1 FASTQ', 'File Size', 'Reads', 'Mean Read Length', 'Mean Read Quality', |
| 8 | 107 'Reads Passing Q30', 'Read2 FASTQ', 'File Size', 'Reads', 'Mean Read Length', 'Mean Read Quality', |
| 108 'Reads Passing Q30', 'Total Reads', 'All Mapped Reads', 'Unmapped Reads', | |
| 109 'Unmapped Reads Percentage of Total', 'Reference with Coverage', 'Average Depth of Coverage', | |
| 23 | 110 'Good SNP Count', 'Reference'] |
| 8 | 111 else: |
| 23 | 112 columns = ['FASTQ', 'File Size', 'Mean Read Length', 'Mean Read Quality', 'Reads Passing Q30', |
| 8 | 113 'Total Reads', 'All Mapped Reads', 'Unmapped Reads', 'Unmapped Reads Percentage of Total', |
| 23 | 114 'Reference with Coverage', 'Average Depth of Coverage', 'Good SNP Count', 'Reference'] |
| 8 | 115 with open(output_file, "w") as outfh: |
| 21 | 116 # Make sure the header starts with a # so |
| 117 # MultiQC can properly handle the output. | |
| 23 | 118 outfh.write("%s\n" % "\t".join(columns)) |
| 8 | 119 line_items = [] |
| 120 # Get the current stats and associated files. | |
| 121 # Get and output the statistics. | |
| 122 line_items.append(read1_stats.fastq_file) | |
| 123 line_items.append(read1_stats.file_size) | |
| 124 if paired_reads: | |
| 125 line_items.append(read1_stats.total_reads) | |
| 126 line_items.append(read1_stats.mean_read_length) | |
| 127 line_items.append(read1_stats.mean_read_quality) | |
| 128 line_items.append(read1_stats.reads_passing_q30) | |
| 129 if paired_reads: | |
| 130 line_items.append(read2_stats.fastq_file) | |
| 131 line_items.append(read2_stats.file_size) | |
| 132 line_items.append(read2_stats.total_reads) | |
| 133 line_items.append(read2_stats.mean_read_length) | |
| 134 line_items.append(read2_stats.mean_read_quality) | |
| 135 line_items.append(read2_stats.reads_passing_q30) | |
| 1 | 136 # Total Reads |
| 8 | 137 if paired_reads: |
| 138 total_reads = read1_stats.total_reads + read2_stats.total_reads | |
| 139 else: | |
| 140 total_reads = read1_stats.total_reads | |
| 141 line_items.append(total_reads) | |
| 1 | 142 # All Mapped Reads |
| 143 all_mapped_reads, unmapped_reads = process_idxstats_file(idxstats_file) | |
| 8 | 144 line_items.append(all_mapped_reads) |
| 145 line_items.append(unmapped_reads) | |
| 1 | 146 # Unmapped Reads Percentage of Total |
| 147 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
|
148 unmapped_reads_percentage = '{:10.2f}'.format(unmapped_reads / total_reads) |
| 0 | 149 else: |
| 1 | 150 unmapped_reads_percentage = 0 |
| 8 | 151 line_items.append(unmapped_reads_percentage) |
| 1 | 152 # Reference with Coverage |
| 153 ref_with_coverage, avg_depth_of_coverage, good_snp_count = process_metrics_file(metrics_file) | |
| 8 | 154 line_items.append(ref_with_coverage) |
| 155 line_items.append(avg_depth_of_coverage) | |
| 156 line_items.append(good_snp_count) | |
| 23 | 157 line_items.append(read1_stats.reference) |
| 8 | 158 outfh.write('%s\n' % '\t'.join(str(x) for x in line_items)) |
| 0 | 159 |
| 160 | |
| 1 | 161 def process_idxstats_file(idxstats_file): |
| 162 all_mapped_reads = 0 | |
| 163 unmapped_reads = 0 | |
| 164 with open(idxstats_file, "r") as fh: | |
| 165 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
|
166 line = line.rstrip('\r\n') |
| 1 | 167 items = line.split("\t") |
| 168 if i == 0: | |
| 169 # NC_002945.4 4349904 213570 4047 | |
| 170 all_mapped_reads = int(items[2]) | |
| 171 elif i == 1: | |
| 172 # * 0 0 82774 | |
| 173 unmapped_reads = int(items[3]) | |
| 174 return all_mapped_reads, unmapped_reads | |
| 0 | 175 |
| 176 | |
| 1 | 177 def process_metrics_file(metrics_file): |
| 178 ref_with_coverage = '0%' | |
| 179 avg_depth_of_coverage = 0 | |
| 180 good_snp_count = 0 | |
| 181 with open(metrics_file, "r") as ifh: | |
| 182 for i, line in enumerate(ifh): | |
| 183 if i == 0: | |
| 184 # Skip comments. | |
| 185 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
|
186 line = line.rstrip('\r\n') |
| 1 | 187 items = line.split("\t") |
| 188 if i == 1: | |
| 189 # MarkDuplicates 10.338671 98.74% | |
| 190 ref_with_coverage = items[3] | |
| 191 avg_depth_of_coverage = items[2] | |
| 192 elif i == 2: | |
| 193 # VCFfilter 611 | |
| 194 good_snp_count = items[1] | |
| 195 return ref_with_coverage, avg_depth_of_coverage, good_snp_count | |
| 0 | 196 |
| 197 | |
| 4 | 198 parser = argparse.ArgumentParser() |
| 0 | 199 |
| 4 | 200 parser.add_argument('--dbkey', action='store', dest='dbkey', help='Reference dbkey') |
| 201 parser.add_argument('--gzipped', action='store_true', dest='gzipped', required=False, default=False, help='Input files are gzipped') | |
| 202 parser.add_argument('--output', action='store', dest='output', help='Output Excel statistics file') | |
| 203 parser.add_argument('--read1', action='store', dest='read1', help='Required: single read') | |
| 204 parser.add_argument('--read2', action='store', dest='read2', required=False, default=None, help='Optional: paired read') | |
| 205 parser.add_argument('--samtools_idxstats', action='store', dest='samtools_idxstats', help='Output of samtools_idxstats') | |
| 8 | 206 parser.add_argument('--vsnp_azc_metrics', action='store', dest='vsnp_azc_metrics', help='Output of vsnp_add_zero_coverage') |
| 0 | 207 |
| 4 | 208 args = parser.parse_args() |
| 0 | 209 |
| 8 | 210 stats_list = [] |
| 4 | 211 idxstats_files = [] |
| 212 metrics_files = [] | |
| 213 # Accumulate inputs. | |
| 8 | 214 read1_stats, read2_stats = accrue_statistics(args.dbkey, args.read1, args.read2, args.gzipped) |
| 215 output_statistics(read1_stats, read2_stats, args.samtools_idxstats, args.vsnp_azc_metrics, args.output) |