Mercurial > repos > iuc > vsnp_add_zero_coverage

comparison vsnp_add_zero_coverage.py @ 9:40b97055bb99 draft default tip

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
planemo upload for repository https://github.com/galaxyproject/tools-iuc/tree/master/tools/vsnp commit c38fd63f7980c70390d104a73ba4c72b266444c3
author iuc
date Fri, 10 Jun 2022 06:08:02 +0000
parents 2e863710a2f0
children
comparison
equal deleted inserted replaced
8:18b59c38017e 9:40b97055bb99
31 31
32 32
33 def get_zero_df(reference): 33 def get_zero_df(reference):
34 # Create a zero coverage dictionary. 34 # Create a zero coverage dictionary.
35 zero_dict = {} 35 zero_dict = {}
36 reference_length = 0
36 for record in SeqIO.parse(reference, "fasta"): 37 for record in SeqIO.parse(reference, "fasta"):
37 chrom = record.id 38 chrom = record.id
38 total_len = len(record.seq) 39 total_len = len(record.seq)
40 reference_length = reference_length + len(record.seq)
39 for pos in list(range(1, total_len + 1)): 41 for pos in list(range(1, total_len + 1)):
40 zero_dict["%s-%s" % (str(chrom), str(pos))] = 0 42 zero_dict["%s-%s" % (str(chrom), str(pos))] = 0
41 # Convert it to a data frame with depth_x 43 # Convert it to a data frame with depth_x
42 # and depth_y columns - index is NaN. 44 # and depth_y columns - index is NaN.
43 zero_df = pandas.DataFrame.from_dict(zero_dict, orient='index', columns=["depth"]) 45 zero_df = pandas.DataFrame.from_dict(zero_dict, orient='index', columns=["depth"])
44 return zero_df 46 return zero_df, reference_length
45 47
46 48
47 def output_zc_vcf_file(base_file_name, vcf_file, zero_df, total_zero_coverage, output_vcf): 49 def output_zc_vcf_file(base_file_name, vcf_file, zero_df, total_zero_coverage, output_vcf):
48 column_names = ["CHROM", "POS", "ID", "REF", "ALT", "QUAL", "FILTER", "INFO", "FORMAT", "Sample"] 50 column_names = ["CHROM", "POS", "ID", "REF", "ALT", "QUAL", "FILTER", "INFO", "FORMAT", "Sample"]
49 vcf_df = pandas.read_csv(vcf_file, sep='\t', header=None, names=column_names, comment='#') 51 vcf_df = pandas.read_csv(vcf_file, sep='\t', header=None, names=column_names, comment='#')
80 else: 82 else:
81 shutil.move(vcf_file, output_vcf) 83 shutil.move(vcf_file, output_vcf)
82 return good_snp_count 84 return good_snp_count
83 85
84 86
85 def output_metrics_file(base_file_name, average_coverage, genome_coverage, good_snp_count, output_metrics): 87 def output_metrics_file(base_file_name, average_coverage, genome_coverage, good_snp_count, dbkey, reference,
86 bam_metrics = [base_file_name, "", "%4f" % average_coverage, genome_coverage] 88 reference_length, total_zero_coverage, percent_ref_with_zero_coverage, output_metrics):
87 vcf_metrics = [base_file_name, str(good_snp_count), "", ""] 89 columns = ["BAM File", "Reference", "reference Length", "Genome with Coverage", "Average Depth",
88 metrics_columns = ["File", "Number of Good SNPs", "Average Coverage", "Genome Coverage"] 90 "No Coverage Bases", "Percent Ref with Zero Coverage", "Quality SNPs"]
91 values = [base_file_name, dbkey, str(reference_length), genome_coverage, average_coverage,
92 str(total_zero_coverage), percent_ref_with_zero_coverage, str(good_snp_count)]
89 with open(output_metrics, "w") as fh: 93 with open(output_metrics, "w") as fh:
90 fh.write("# %s\n" % "\t".join(metrics_columns)) 94 fh.write("# %s\n" % "\t".join(columns))
91 fh.write("%s\n" % "\t".join(bam_metrics)) 95 fh.write("%s\n" % "\t".join(values))
92 fh.write("%s\n" % "\t".join(vcf_metrics))
93 96
94 97
95 def output_files(vcf_file, total_zero_coverage, zero_df, output_vcf, average_coverage, genome_coverage, output_metrics): 98 def output_files(vcf_file, total_zero_coverage, percent_ref_with_zero_coverage, zero_df, output_vcf,
99 average_coverage, genome_coverage, output_metrics, reference, reference_length, dbkey):
96 base_file_name = get_sample_name(vcf_file) 100 base_file_name = get_sample_name(vcf_file)
97 good_snp_count = output_zc_vcf_file(base_file_name, vcf_file, zero_df, total_zero_coverage, output_vcf) 101 good_snp_count = output_zc_vcf_file(base_file_name, vcf_file, zero_df, total_zero_coverage, output_vcf)
98 output_metrics_file(base_file_name, average_coverage, genome_coverage, good_snp_count, output_metrics) 102 output_metrics_file(base_file_name, average_coverage, genome_coverage, good_snp_count, dbkey, reference,
103 reference_length, total_zero_coverage, percent_ref_with_zero_coverage, output_metrics)
99 104
100 105
101 def get_coverage_and_snp_count(bam_file, vcf_file, reference, output_metrics, output_vcf): 106 def get_coverage_and_snp_count(bam_file, vcf_file, dbkey, reference, output_metrics, output_vcf):
102 coverage_df = get_coverage_df(bam_file) 107 coverage_df = get_coverage_df(bam_file)
103 zero_df = get_zero_df(reference) 108 zero_df, reference_length = get_zero_df(reference)
104 coverage_df = zero_df.merge(coverage_df, left_index=True, right_index=True, how='outer') 109 coverage_df = zero_df.merge(coverage_df, left_index=True, right_index=True, how='outer')
105 # depth_x "0" column no longer needed. 110 # depth_x "0" column no longer needed.
106 coverage_df = coverage_df.drop(columns=['depth_x']) 111 coverage_df = coverage_df.drop(columns=['depth_x'])
107 coverage_df = coverage_df.rename(columns={'depth_y': 'depth'}) 112 coverage_df = coverage_df.rename(columns={'depth_y': 'depth'})
108 # Covert the NaN to 0 coverage and get some metrics. 113 # Covert the NaN to 0 coverage and get some metrics.
109 coverage_df = coverage_df.fillna(0) 114 coverage_df = coverage_df.fillna(0)
110 coverage_df['depth'] = coverage_df['depth'].apply(int) 115 coverage_df['depth'] = coverage_df['depth'].apply(int)
111 total_length = len(coverage_df) 116 total_length = len(coverage_df)
112 average_coverage = coverage_df['depth'].mean() 117 average_coverage = "{:.2f}".format(coverage_df['depth'].mean())
113 zero_df = coverage_df[coverage_df['depth'] == 0] 118 zero_df = coverage_df[coverage_df['depth'] == 0]
114 total_zero_coverage = len(zero_df) 119 total_zero_coverage = len(zero_df)
120 percent_ref_with_zero_coverage = "{:.6%}".format(total_zero_coverage / reference_length * 100)
115 total_coverage = total_length - total_zero_coverage 121 total_coverage = total_length - total_zero_coverage
116 genome_coverage = "{:.2%}".format(total_coverage / total_length) 122 genome_coverage = "{:.2%}".format(total_coverage / total_length)
117 # Output a zero-coverage vcf fil and the metrics file. 123 # Output a zero-coverage vcf fil and the metrics file.
118 output_files(vcf_file, total_zero_coverage, zero_df, output_vcf, average_coverage, genome_coverage, output_metrics) 124 output_files(vcf_file, total_zero_coverage, percent_ref_with_zero_coverage, zero_df, output_vcf,
125 average_coverage, genome_coverage, output_metrics, reference, reference_length, dbkey)
119 126
120 127
121 if __name__ == '__main__': 128 if __name__ == '__main__':
122 parser = argparse.ArgumentParser() 129 parser = argparse.ArgumentParser()
123 130
124 parser.add_argument('--bam_input', action='store', dest='bam_input', help='bam input file') 131 parser.add_argument('--bam_input', action='store', dest='bam_input', help='bam input file')
132 parser.add_argument('--dbkey', action='store', dest='dbkey', help='bam input dbkey')
125 parser.add_argument('--output_metrics', action='store', dest='output_metrics', required=False, default=None, help='Output metrics text file') 133 parser.add_argument('--output_metrics', action='store', dest='output_metrics', required=False, default=None, help='Output metrics text file')
126 parser.add_argument('--output_vcf', action='store', dest='output_vcf', required=False, default=None, help='Output VCF file') 134 parser.add_argument('--output_vcf', action='store', dest='output_vcf', required=False, default=None, help='Output VCF file')
127 parser.add_argument('--reference', action='store', dest='reference', help='Reference dataset') 135 parser.add_argument('--reference', action='store', dest='reference', help='Reference dataset')
128 parser.add_argument('--vcf_input', action='store', dest='vcf_input', help='vcf input file') 136 parser.add_argument('--vcf_input', action='store', dest='vcf_input', help='vcf input file')
129 137
130 args = parser.parse_args() 138 args = parser.parse_args()
131 139
132 get_coverage_and_snp_count(args.bam_input, args.vcf_input, args.reference, args.output_metrics, args.output_vcf) 140 get_coverage_and_snp_count(args.bam_input, args.vcf_input, args.dbkey, args.reference, args.output_metrics, args.output_vcf)