|
rliterman@0
|
1 #!/usr/bin/env python3
|
|
rliterman@0
|
2
|
|
rliterman@0
|
3 import sys
|
|
rliterman@0
|
4 import os
|
|
rliterman@0
|
5 import glob
|
|
rliterman@0
|
6 import pandas as pd
|
|
rliterman@0
|
7 from itertools import chain
|
|
rliterman@0
|
8 import scipy.stats
|
|
rliterman@0
|
9 import numpy as np
|
|
rliterman@0
|
10 import datetime
|
|
rliterman@0
|
11 import time
|
|
rliterman@0
|
12 import argparse
|
|
rliterman@0
|
13
|
|
rliterman@0
|
14 def getWarnings(df):
|
|
rliterman@0
|
15
|
|
rliterman@0
|
16 df_measures = list(set(df['Measure']))
|
|
rliterman@0
|
17 warn_list = []
|
|
rliterman@0
|
18 if 'Preserved_Diff' in df_measures:
|
|
rliterman@0
|
19 for index,row in df.iterrows():
|
|
rliterman@0
|
20 if pd.isna(row['Zscore']):
|
|
rliterman@0
|
21 warn_list.append(np.nan)
|
|
rliterman@0
|
22 elif 2.5 <= row['Zscore'] < 3:
|
|
rliterman@0
|
23 warn_list.append("Warning")
|
|
rliterman@0
|
24 elif row['Zscore'] >=3:
|
|
rliterman@0
|
25 warn_list.append("Failure")
|
|
rliterman@0
|
26 else:
|
|
rliterman@0
|
27 warn_list.append("Pass")
|
|
rliterman@0
|
28 elif 'Contig_Count' in df_measures:
|
|
rliterman@0
|
29 for index,row in df.iterrows():
|
|
rliterman@0
|
30 if pd.isna(row['Zscore']):
|
|
rliterman@0
|
31 warn_list.append(np.nan)
|
|
rliterman@0
|
32 elif row['Measure'] in ["Contig_Count","L50","L90"]:
|
|
rliterman@0
|
33 if 2.5 <= row['Zscore'] < 3:
|
|
rliterman@0
|
34 warn_list.append("Warning")
|
|
rliterman@0
|
35 elif row['Zscore'] >=3:
|
|
rliterman@0
|
36 warn_list.append("Failure")
|
|
rliterman@0
|
37 else:
|
|
rliterman@0
|
38 warn_list.append("Pass")
|
|
rliterman@0
|
39 elif row['Measure'] == "Assembly_Bases":
|
|
rliterman@0
|
40 if 2.5 <= abs(row['Zscore']) < 3:
|
|
rliterman@0
|
41 warn_list.append("Warning")
|
|
rliterman@0
|
42 elif abs(row['Zscore']) >=3:
|
|
rliterman@0
|
43 warn_list.append("Failure")
|
|
rliterman@0
|
44 else:
|
|
rliterman@0
|
45 warn_list.append("Pass")
|
|
rliterman@0
|
46 elif row['Measure'] in ["N50","N90"]:
|
|
rliterman@0
|
47 if -3 < row['Zscore'] <= -2.5:
|
|
rliterman@0
|
48 warn_list.append("Warning")
|
|
rliterman@0
|
49 elif row['Zscore'] <= -3:
|
|
rliterman@0
|
50 warn_list.append("Failure")
|
|
rliterman@0
|
51 else:
|
|
rliterman@0
|
52 warn_list.append("Pass")
|
|
rliterman@0
|
53 else:
|
|
rliterman@0
|
54 sys.exit(f"{row['Measure']}")
|
|
rliterman@0
|
55 elif ('Raw_Distance_StdDev' in df_measures) | ('Preserved_Distance_StdDev' in df_measures):
|
|
rliterman@0
|
56 for index,row in df.iterrows():
|
|
rliterman@0
|
57 if pd.isna(row['Zscore']):
|
|
rliterman@0
|
58 warn_list.append(np.nan)
|
|
rliterman@0
|
59 elif 2.5 <= row['Zscore'] < 3:
|
|
rliterman@0
|
60 warn_list.append("Warning")
|
|
rliterman@0
|
61 elif row['Zscore'] >=3:
|
|
rliterman@0
|
62 warn_list.append("Failure")
|
|
rliterman@0
|
63 else:
|
|
rliterman@0
|
64 warn_list.append("Pass")
|
|
rliterman@0
|
65 elif 'Unique_Kmers' in df_measures:
|
|
rliterman@0
|
66 for index,row in df.iterrows():
|
|
rliterman@0
|
67 if pd.isna(row['Zscore']):
|
|
rliterman@0
|
68 warn_list.append(np.nan)
|
|
rliterman@0
|
69 elif row['Measure'] in ["Align_Percent_Diff","Unique_Kmers","gIndels","Missing_Kmers"]:
|
|
rliterman@0
|
70 if 2.5 <= row['Zscore'] < 3:
|
|
rliterman@0
|
71 warn_list.append("Warning")
|
|
rliterman@0
|
72 elif row['Zscore'] >=3:
|
|
rliterman@0
|
73 warn_list.append("Failure")
|
|
rliterman@0
|
74 else:
|
|
rliterman@0
|
75 warn_list.append("Pass")
|
|
rliterman@0
|
76 elif row['Measure'] in ["Compare_Aligned","Kmer_Similarity","Self_Aligned","Median_Alignment_Length"]:
|
|
rliterman@0
|
77 if -3 < row['Zscore'] <= -2.5:
|
|
rliterman@0
|
78 warn_list.append("Warning")
|
|
rliterman@0
|
79 elif row['Zscore'] <= -3:
|
|
rliterman@0
|
80 warn_list.append("Failure")
|
|
rliterman@0
|
81 else:
|
|
rliterman@0
|
82 warn_list.append("Pass")
|
|
rliterman@0
|
83 else:
|
|
rliterman@0
|
84 sys.exit(f"{row['Measure']}")
|
|
rliterman@0
|
85
|
|
rliterman@0
|
86 elif ('SNPs_Cocalled' in df_measures) | ('Raw_SNPs_Cocalled' in df_measures) | ('Preserved_SNPs_Cocalled' in df_measures):
|
|
rliterman@0
|
87 for index,row in df.iterrows():
|
|
rliterman@0
|
88 if pd.isna(row['Zscore']):
|
|
rliterman@0
|
89 warn_list.append(np.nan)
|
|
rliterman@0
|
90 elif -3 < row['Zscore'] <= -2.5:
|
|
rliterman@0
|
91 warn_list.append("Warning")
|
|
rliterman@0
|
92 elif row['Zscore'] <= -3:
|
|
rliterman@0
|
93 warn_list.append("Failure")
|
|
rliterman@0
|
94 else:
|
|
rliterman@0
|
95 warn_list.append("Pass")
|
|
rliterman@0
|
96 else:
|
|
rliterman@0
|
97 sys.exit(f"{df_measures}")
|
|
rliterman@0
|
98
|
|
rliterman@0
|
99 return warn_list
|
|
rliterman@0
|
100
|
|
rliterman@0
|
101 start_time = time.time()
|
|
rliterman@0
|
102
|
|
rliterman@0
|
103 # Get args
|
|
rliterman@0
|
104 parser = argparse.ArgumentParser(description='CSP2 SNP Pipeline Compiler')
|
|
rliterman@0
|
105 parser.add_argument('--snp_dirs_file', type=str, help='Path to the file containing SNP directories')
|
|
rliterman@0
|
106 parser.add_argument('--output_directory', type=str, help='Path to the output directory')
|
|
rliterman@0
|
107 parser.add_argument('--isolate_data_file', type=str, help='Path to the isolate data file')
|
|
rliterman@0
|
108 parser.add_argument('--mummer_data_file', type=str, help='Path to the MUMmer data file')
|
|
rliterman@0
|
109 args = parser.parse_args()
|
|
rliterman@0
|
110
|
|
rliterman@0
|
111 snp_dirs = [line.strip() for line in open(args.snp_dirs_file, 'r')]
|
|
rliterman@0
|
112 raw_snp_distance_files = list(chain.from_iterable([glob.glob(snp_dir + '/snp_distance_pairwise.tsv') for snp_dir in snp_dirs]))
|
|
rliterman@0
|
113 screening_files = list(chain.from_iterable([glob.glob(snp_dir + '/Reference_Screening.tsv') for snp_dir in snp_dirs]))
|
|
rliterman@0
|
114
|
|
rliterman@0
|
115 # Set paths
|
|
rliterman@0
|
116 output_directory = args.output_directory
|
|
rliterman@0
|
117 log_file = f"{output_directory}/Compilation.log"
|
|
rliterman@0
|
118 mean_isolate_file = f"{output_directory}/Mean_Assembly_Stats.tsv"
|
|
rliterman@0
|
119 isolate_assembly_stats_file = f"{output_directory}/Isolate_Assembly_Stats.tsv"
|
|
rliterman@0
|
120 align_stats_file = f"{output_directory}/Isolate_Alignment_Stats.tsv"
|
|
rliterman@0
|
121 ref_mean_summary_file = f"{output_directory}/Align_Summary_by_Reference.tsv"
|
|
rliterman@0
|
122 snp_comparison_file = f"{output_directory}/SNP_Distance_Summary.tsv"
|
|
rliterman@0
|
123 qc_file = f"{output_directory}/QC_Warnings_Failures.tsv"
|
|
rliterman@0
|
124
|
|
rliterman@0
|
125 with open(log_file,"w+") as log:
|
|
rliterman@0
|
126 log.write("CSP2 SNP Pipeline Compiler\n")
|
|
rliterman@0
|
127 log.write(str(datetime.datetime.now().strftime("%Y-%m-%d %H:%M:%S"))+"\n")
|
|
rliterman@0
|
128 log.write("-------------------------------------------------------\n\n")
|
|
rliterman@0
|
129 if len(raw_snp_distance_files) == 0:
|
|
rliterman@0
|
130 log.write("\t- CSP2 SNP Pipeline Compiler cannot detected any SNP pipeline output files\n")
|
|
rliterman@0
|
131 log.write("\t- Compiler stopping...\n")
|
|
rliterman@0
|
132 sys.exit(0)
|
|
rliterman@0
|
133
|
|
rliterman@0
|
134 isolate_data = pd.read_csv(args.isolate_data_file, sep="\t")
|
|
rliterman@0
|
135 raw_isolate_count = isolate_data.shape[0]
|
|
rliterman@0
|
136
|
|
rliterman@0
|
137 mummer_data = pd.read_csv(args.mummer_data_file, sep="\t")
|
|
rliterman@0
|
138
|
|
rliterman@0
|
139 # Get reference IDs
|
|
rliterman@0
|
140 reference_ids = list(set(isolate_data[isolate_data['Isolate_Type'] == "Reference"]['Isolate_ID']))
|
|
rliterman@0
|
141 raw_ref_count = len(reference_ids)
|
|
rliterman@0
|
142
|
|
rliterman@0
|
143 raw_snp_distance_df = pd.concat([pd.read_csv(file, sep='\t').assign(Reference_ID=os.path.basename(os.path.dirname(file))) for file in raw_snp_distance_files])
|
|
rliterman@0
|
144 raw_snp_distance_df['Comparison'] = raw_snp_distance_df.apply(lambda row: ';'.join(sorted([str(row['Query_1']), str(row['Query_2'])])), axis=1)
|
|
rliterman@0
|
145
|
|
rliterman@0
|
146 # Check for preserved data
|
|
rliterman@0
|
147 preserved_snp_distance_files = list(chain.from_iterable([glob.glob(snp_dir + '/snp_distance_pairwise_preserved.tsv') for snp_dir in snp_dirs]))
|
|
rliterman@0
|
148 if len(preserved_snp_distance_files) == 0:
|
|
rliterman@0
|
149 has_preserved = False
|
|
rliterman@0
|
150 else:
|
|
rliterman@0
|
151 has_preserved = True
|
|
rliterman@0
|
152 preserved_snp_distance_df = pd.concat([pd.read_csv(file, sep='\t').assign(Reference_ID=os.path.basename(os.path.dirname(file))) for file in preserved_snp_distance_files])
|
|
rliterman@0
|
153 preserved_snp_distance_df['Comparison'] = preserved_snp_distance_df.apply(lambda row: ';'.join(sorted([str(row['Query_1']), str(row['Query_2'])])), axis=1)
|
|
rliterman@0
|
154
|
|
rliterman@0
|
155 screening_df = pd.concat([pd.read_csv(file, sep='\t').assign(Reference_ID=os.path.basename(os.path.dirname(file))) for file in screening_files])
|
|
rliterman@0
|
156
|
|
rliterman@0
|
157 snp_isolates = list(set(raw_snp_distance_df['Query_1'].tolist() + raw_snp_distance_df['Query_2'].tolist()))
|
|
rliterman@0
|
158
|
|
rliterman@0
|
159 with open(log_file,"a+") as log:
|
|
rliterman@0
|
160 log.write(f"- Detected SNP distance data for {len(snp_isolates)} isolates out of {raw_isolate_count} total isolates analyzed\n")
|
|
rliterman@0
|
161 if len(snp_isolates) <= 2:
|
|
rliterman@0
|
162 log.write("\t- CSP2 SNP Pipeline Compiler cannot do much with 2 or fewer isolates\n")
|
|
rliterman@0
|
163 log.write("\t- Compiler stopping...\n")
|
|
rliterman@0
|
164 sys.exit(0)
|
|
rliterman@0
|
165 else:
|
|
rliterman@0
|
166 failed_combinations = screening_df.loc[screening_df['Screen_Category'] != "Pass"]
|
|
rliterman@0
|
167 failed_comparisons = []
|
|
rliterman@0
|
168
|
|
rliterman@0
|
169 if failed_combinations.shape[0] > 0:
|
|
rliterman@0
|
170 reference_query_dict = failed_combinations.groupby('Reference_ID')['Query_ID'].apply(list).to_dict()
|
|
rliterman@0
|
171 for index,row in failed_combinations.iterrows():
|
|
rliterman@0
|
172 failed_comparisons.append(";".join(sorted([row['Query_ID'], row['Reference_ID']])))
|
|
rliterman@0
|
173 log.write("\n- The following query-reference alignments did not pass QC\n")
|
|
rliterman@0
|
174 for key, value in reference_query_dict.items():
|
|
rliterman@0
|
175 log.write(f"\nReference {key}\n{', '.join(map(str, value))}\n")
|
|
rliterman@0
|
176
|
|
rliterman@0
|
177 # Prune isolate data
|
|
rliterman@0
|
178 isolate_data = isolate_data.loc[isolate_data['Isolate_ID'].isin(snp_isolates)]
|
|
rliterman@0
|
179 reference_ids = [x for x in reference_ids if x in snp_isolates]
|
|
rliterman@0
|
180 ref_count = len(reference_ids)
|
|
rliterman@0
|
181 with open(log_file,"a+") as log:
|
|
rliterman@0
|
182 log.write(f"- Detected SNP distance data for {ref_count} reference isolates out of {raw_ref_count} total reference isolates analyzed\n")
|
|
rliterman@0
|
183 for ref in reference_ids:
|
|
rliterman@0
|
184 log.write(f"\t- {ref}\n")
|
|
rliterman@0
|
185
|
|
rliterman@0
|
186 # Prune MUMmer data
|
|
rliterman@0
|
187 mummer_data['Comparison'] = mummer_data.apply(lambda row: ';'.join(sorted([str(row['Query_ID']), str(row['Reference_ID'])])), axis=1)
|
|
rliterman@0
|
188 mummer_data = mummer_data.loc[~mummer_data['Comparison'].isin(failed_comparisons), ['SNPDiffs_File','Query_ID','Reference_ID','Comparison','Reference_Percent_Aligned','Query_Percent_Aligned','Median_Alignment_Length','Kmer_Similarity','Reference_Unique_Kmers','Query_Unique_Kmers','gIndels']]
|
|
rliterman@0
|
189 max_align_values = np.maximum(mummer_data['Reference_Percent_Aligned'], mummer_data['Query_Percent_Aligned'])
|
|
rliterman@0
|
190 min_align_values = np.minimum(mummer_data['Reference_Percent_Aligned'], mummer_data['Query_Percent_Aligned'])
|
|
rliterman@0
|
191 mummer_data['Align_Percent_Diff'] = 100*((max_align_values - min_align_values)/min_align_values)
|
|
rliterman@0
|
192
|
|
rliterman@0
|
193 # Run basic assembly stats
|
|
rliterman@0
|
194 isolate_stats = isolate_data.melt(id_vars=['Isolate_ID', 'Isolate_Type'], value_vars = ['Contig_Count','Assembly_Bases','N50','N90','L50','L90'],value_name='Value',var_name = "Measure")
|
|
rliterman@0
|
195 isolate_stats['Zscore'] = isolate_stats.groupby('Measure')['Value'].transform(scipy.stats.zscore).astype('float').round(3)
|
|
rliterman@0
|
196 isolate_stats['QC'] = getWarnings(isolate_stats)
|
|
rliterman@0
|
197 isolate_stats['Value'] = isolate_stats['Value'].astype('int')
|
|
rliterman@0
|
198
|
|
rliterman@0
|
199 # Reformat for final TSV
|
|
rliterman@0
|
200 isolate_stats['Min'] = np.nan
|
|
rliterman@0
|
201 isolate_stats['Max'] = np.nan
|
|
rliterman@0
|
202 isolate_stats['StdDev'] = np.nan
|
|
rliterman@0
|
203 isolate_stats = isolate_stats[['Isolate_ID','Isolate_Type','Measure','Min','Value','Max','StdDev','Zscore','QC']].rename(columns = {'Value':'Mean'})
|
|
rliterman@0
|
204 isolate_stats['Count'] = 1
|
|
rliterman@0
|
205
|
|
rliterman@0
|
206 # Get mean values
|
|
rliterman@0
|
207 isolate_mean_df = isolate_stats.groupby(by=['Measure'])['Mean'].agg(Min = 'min',Mean = "mean",Max = 'max',StdDev = 'std',Count = 'count')
|
|
rliterman@0
|
208 isolate_mean_df['Mean'] = isolate_mean_df['Mean'].astype("int")
|
|
rliterman@0
|
209 isolate_mean_df['StdDev'] = isolate_mean_df['StdDev'].astype("float").round(3)
|
|
rliterman@0
|
210 with open(log_file,"a+") as log:
|
|
rliterman@0
|
211 log.write("- Read in and processed isolate data\n\n")
|
|
rliterman@0
|
212 for index, row in isolate_mean_df.iterrows():
|
|
rliterman@0
|
213 log.write(f"{index}:\tMin: {row['Min']}\tMean: {row['Mean']}\tMax: {row['Max']}\tStdDev: {row['StdDev']}\n")
|
|
rliterman@0
|
214 log.write("\n")
|
|
rliterman@0
|
215
|
|
rliterman@0
|
216 # Run basic alignment stats
|
|
rliterman@0
|
217 isolate_mummer = pd.DataFrame(columns=['Isolate_ID', 'Compare_ID', 'Self_Aligned', 'Compare_Aligned','Align_Percent_Diff','Median_Alignment_Length', 'Kmer_Similarity', 'Unique_Kmers', 'Missing_Kmers', 'gIndels', 'SNPDiffs_File'])
|
|
rliterman@0
|
218
|
|
rliterman@0
|
219 for isolate in snp_isolates:
|
|
rliterman@0
|
220
|
|
rliterman@0
|
221 temp_mummer = mummer_data[(mummer_data['Query_ID'] == isolate) | (mummer_data['Reference_ID'] == isolate)].drop_duplicates(subset=['Comparison']).assign(Isolate_ID = isolate)
|
|
rliterman@0
|
222
|
|
rliterman@0
|
223 for index, row in temp_mummer.iterrows():
|
|
rliterman@0
|
224 if row['Query_ID'] == isolate:
|
|
rliterman@0
|
225 temp_isolate_mummer = row[['Isolate_ID', 'Reference_ID', 'Query_Percent_Aligned', 'Reference_Percent_Aligned','Align_Percent_Diff', 'Median_Alignment_Length', 'Kmer_Similarity', 'Query_Unique_Kmers', 'Reference_Unique_Kmers', 'gIndels', 'SNPDiffs_File']].to_frame().T
|
|
rliterman@0
|
226 temp_isolate_mummer.columns = ['Isolate_ID', 'Compare_ID', 'Self_Aligned', 'Compare_Aligned', 'Align_Percent_Diff','Median_Alignment_Length', 'Kmer_Similarity', 'Unique_Kmers', 'Missing_Kmers', 'gIndels', 'SNPDiffs_File']
|
|
rliterman@0
|
227 isolate_mummer = pd.concat([isolate_mummer,temp_isolate_mummer])
|
|
rliterman@0
|
228
|
|
rliterman@0
|
229 elif row['Reference_ID'] == isolate:
|
|
rliterman@0
|
230 temp_isolate_mummer = row[['Isolate_ID', 'Query_ID', 'Reference_Percent_Aligned', 'Query_Percent_Aligned', 'Align_Percent_Diff', 'Median_Alignment_Length', 'Kmer_Similarity', 'Reference_Unique_Kmers', 'Query_Unique_Kmers', 'gIndels', 'SNPDiffs_File']].to_frame().T
|
|
rliterman@0
|
231 temp_isolate_mummer.columns = ['Isolate_ID', 'Compare_ID', 'Self_Aligned', 'Compare_Aligned', 'Align_Percent_Diff', 'Median_Alignment_Length', 'Kmer_Similarity', 'Unique_Kmers', 'Missing_Kmers', 'gIndels', 'SNPDiffs_File']
|
|
rliterman@0
|
232 isolate_mummer = pd.concat([isolate_mummer,temp_isolate_mummer])
|
|
rliterman@0
|
233
|
|
rliterman@0
|
234 isolate_mummer['Isolate_Type'] = isolate_mummer['Isolate_ID'].apply(lambda x: 'Reference' if x in reference_ids else 'Query')
|
|
rliterman@0
|
235 isolate_mummer_df = isolate_mummer.melt(id_vars=['Isolate_ID','Isolate_Type'], value_vars = ['Self_Aligned','Compare_Aligned', 'Align_Percent_Diff','Median_Alignment_Length','Kmer_Similarity','Unique_Kmers','Missing_Kmers','gIndels'],value_name='Value',var_name = "Measure")
|
|
rliterman@0
|
236 isolate_mummer_df['Value'] = isolate_mummer_df['Value'].astype("float")
|
|
rliterman@0
|
237 isolate_mummer_df = isolate_mummer_df.groupby(by=['Isolate_ID','Isolate_Type','Measure'])['Value'].agg(Count = "count",Min = "min",Value = "mean",Max = "max",StdDev = 'std').reset_index()
|
|
rliterman@0
|
238 isolate_mummer_df['Value'] = isolate_mummer_df['Value'].astype("float").round(2)
|
|
rliterman@0
|
239
|
|
rliterman@0
|
240 # Get Zscores
|
|
rliterman@0
|
241 isolate_mummer_df['Zscore'] = isolate_mummer_df.groupby('Measure')['Value'].transform(scipy.stats.zscore).astype('float').round(3)
|
|
rliterman@0
|
242 isolate_mummer_df['QC'] = getWarnings(isolate_mummer_df)
|
|
rliterman@0
|
243
|
|
rliterman@0
|
244 # Reformat for final TSV
|
|
rliterman@0
|
245 align_stats = isolate_mummer_df[['Isolate_ID','Isolate_Type','Measure','Min','Value','Max','StdDev','Zscore','QC','Count']].copy().rename(columns = {"Value":"Mean"})
|
|
rliterman@0
|
246 align_stats['StdDev'] = align_stats['StdDev'].astype('float').round(3)
|
|
rliterman@0
|
247
|
|
rliterman@0
|
248 with open(log_file,"a+") as log:
|
|
rliterman@0
|
249 log.write("- Read in and processed alignment data\n")
|
|
rliterman@0
|
250
|
|
rliterman@0
|
251 # Process cocalled data
|
|
rliterman@0
|
252 raw_cocalled_df = raw_snp_distance_df[['Comparison','Query_1','Query_2','Reference_ID','SNPs_Cocalled']]
|
|
rliterman@0
|
253 isolate_cocalled_df = pd.DataFrame(columns = ['Isolate_ID','Count','Min','Mean','Max','StdDev'])
|
|
rliterman@0
|
254
|
|
rliterman@0
|
255 for isolate in snp_isolates:
|
|
rliterman@0
|
256 temp_cocalled = raw_cocalled_df[(raw_cocalled_df['Query_1'] == isolate) | (raw_cocalled_df['Query_2'] == isolate)].drop_duplicates(subset=['Comparison','Reference_ID']).assign(Isolate_ID = isolate)
|
|
rliterman@0
|
257 temp_cocalled = temp_cocalled.groupby(['Isolate_ID'])['SNPs_Cocalled'].agg(Count = "count", Min = "min", Value = "mean", Max = "max",StdDev = 'std').reset_index()
|
|
rliterman@0
|
258 isolate_cocalled_df = pd.concat([isolate_cocalled_df,temp_cocalled])
|
|
rliterman@0
|
259
|
|
rliterman@0
|
260 isolate_cocalled_df['Measure'] = 'Raw_SNPs_Cocalled'
|
|
rliterman@0
|
261 isolate_cocalled_df['Value'] = isolate_cocalled_df['Value'].astype('int')
|
|
rliterman@0
|
262 isolate_cocalled_df['Zscore'] = isolate_cocalled_df['Value'].transform(scipy.stats.zscore).astype('float').round(3)
|
|
rliterman@0
|
263 isolate_cocalled_df['QC'] = getWarnings(isolate_cocalled_df)
|
|
rliterman@0
|
264
|
|
rliterman@0
|
265 # Format for final TSV
|
|
rliterman@0
|
266 isolate_cocalled_df['Isolate_Type'] = isolate_cocalled_df['Isolate_ID'].apply(lambda x: 'Reference' if x in reference_ids else 'Query')
|
|
rliterman@0
|
267 isolate_cocalled_stats = isolate_cocalled_df[['Isolate_ID','Isolate_Type','Measure','Min','Value','Max','StdDev','Zscore','QC','Count']].copy().rename(columns={'Value':'Mean'})
|
|
rliterman@0
|
268
|
|
rliterman@0
|
269 if has_preserved:
|
|
rliterman@0
|
270 preserved_cocalled_df = preserved_snp_distance_df[['Comparison','Query_1','Query_2','Reference_ID','SNPs_Cocalled']]
|
|
rliterman@0
|
271 isolate_preserved_cocalled_df = pd.DataFrame(columns = ['Isolate_ID','Count','Min','Mean','Max','StdDev'])
|
|
rliterman@0
|
272
|
|
rliterman@0
|
273 for isolate in snp_isolates:
|
|
rliterman@0
|
274 temp_cocalled = preserved_cocalled_df[(preserved_cocalled_df['Query_1'] == isolate) | (preserved_cocalled_df['Query_2'] == isolate)].drop_duplicates(subset=['Comparison','Reference_ID']).assign(Isolate_ID = isolate)
|
|
rliterman@0
|
275 temp_cocalled = temp_cocalled.groupby(['Isolate_ID'])['SNPs_Cocalled'].agg(Count = "count", Min = "min", Value = "mean", Max = "max",StdDev = 'std').reset_index()
|
|
rliterman@0
|
276 isolate_preserved_cocalled_df = pd.concat([isolate_preserved_cocalled_df,temp_cocalled])
|
|
rliterman@0
|
277
|
|
rliterman@0
|
278 isolate_preserved_cocalled_df['Measure'] = 'Preserved_SNPs_Cocalled'
|
|
rliterman@0
|
279 isolate_preserved_cocalled_df['Value'] = isolate_preserved_cocalled_df['Value'].astype('int')
|
|
rliterman@0
|
280 isolate_preserved_cocalled_df['Zscore'] = isolate_preserved_cocalled_df['Value'].transform(scipy.stats.zscore).astype('float').round(3)
|
|
rliterman@0
|
281 isolate_preserved_cocalled_df['QC'] = getWarnings(isolate_preserved_cocalled_df)
|
|
rliterman@0
|
282
|
|
rliterman@0
|
283 # Format for final TSV
|
|
rliterman@0
|
284 isolate_preserved_cocalled_df['Isolate_Type'] = isolate_preserved_cocalled_df['Isolate_ID'].apply(lambda x: 'Reference' if x in reference_ids else 'Query')
|
|
rliterman@0
|
285 isolate_preserved_cocalled_df['StdDev'] = isolate_preserved_cocalled_df['StdDev'].astype('float').round(3)
|
|
rliterman@0
|
286 isolate_cocalled_stats = pd.concat([isolate_cocalled_stats,isolate_preserved_cocalled_df[['Isolate_ID','Isolate_Type','Measure','Min','Value','Max','StdDev','Zscore','QC','Count']].copy().rename(columns={'Value':'Mean'})])
|
|
rliterman@0
|
287
|
|
rliterman@0
|
288 with open(log_file,"a+") as log:
|
|
rliterman@0
|
289 log.write("- Processed cocalled SNP data\n")
|
|
rliterman@0
|
290
|
|
rliterman@0
|
291 if has_preserved:
|
|
rliterman@0
|
292 raw_snp_df = raw_snp_distance_df[['Comparison','Query_1','Query_2','Reference_ID','SNP_Distance']].rename(columns = {'SNP_Distance':'Raw_SNP_Distance'})
|
|
rliterman@0
|
293 preserved_snp_df = preserved_snp_distance_df[['Comparison','Query_1','Query_2','Reference_ID','SNP_Distance']].rename(columns = {'SNP_Distance':'Preserved_SNP_Distance'})
|
|
rliterman@0
|
294 snp_df = pd.merge(raw_snp_df,preserved_snp_df,how="left",on=['Comparison','Query_1','Query_2','Reference_ID'])
|
|
rliterman@0
|
295 snp_df['Preserved_Diff'] = abs(snp_df['Preserved_SNP_Distance'] - snp_df['Raw_SNP_Distance'])
|
|
rliterman@0
|
296
|
|
rliterman@0
|
297 isolate_snp_df = pd.DataFrame(columns = ['Isolate_ID','Count','Min','Value','Max','StdDev'])
|
|
rliterman@0
|
298
|
|
rliterman@0
|
299 for isolate in snp_isolates:
|
|
rliterman@0
|
300 temp_snp = snp_df[(snp_df['Query_1'] == isolate) | (snp_df['Query_2'] == isolate)].drop_duplicates(subset=['Comparison','Reference_ID']).assign(Isolate_ID = isolate)
|
|
rliterman@0
|
301 temp_snp = temp_snp.groupby(['Isolate_ID'])['Preserved_Diff'].agg(Count = "count", Min = "min", Value = "mean", Max = "max",StdDev = 'std').reset_index()
|
|
rliterman@0
|
302 isolate_snp_df = pd.concat([isolate_snp_df,temp_snp])
|
|
rliterman@0
|
303
|
|
rliterman@0
|
304 isolate_snp_df['Measure'] = 'Preserved_Diff'
|
|
rliterman@0
|
305 isolate_snp_df['Value'] = isolate_snp_df['Value'].astype("float")
|
|
rliterman@0
|
306 isolate_snp_df['Zscore'] = isolate_snp_df['Value'].transform(scipy.stats.zscore).astype('float').round(3)
|
|
rliterman@0
|
307 isolate_snp_df['Value'] = isolate_snp_df['Value'].astype('float').round(3)
|
|
rliterman@0
|
308 isolate_snp_df['QC'] = getWarnings(isolate_snp_df)
|
|
rliterman@0
|
309
|
|
rliterman@0
|
310 # Format for final TSV
|
|
rliterman@0
|
311 isolate_snp_df['Isolate_Type'] = isolate_snp_df['Isolate_ID'].apply(lambda x: 'Reference' if x in reference_ids else 'Query')
|
|
rliterman@0
|
312 isolate_snp_df['StdDev'] = isolate_snp_df['StdDev'].astype('float').round(3)
|
|
rliterman@0
|
313 isolate_snp_stats = isolate_snp_df[['Isolate_ID','Isolate_Type','Measure','Min','Value','Max','StdDev','Zscore','QC','Count']].copy().rename(columns={'Value':'Mean'})
|
|
rliterman@0
|
314 with open(log_file,"a+") as log:
|
|
rliterman@0
|
315 log.write("- Processed preserved SNP data\n")
|
|
rliterman@0
|
316 else:
|
|
rliterman@0
|
317 with open(log_file,"a+") as log:
|
|
rliterman@0
|
318 log.write("- No preserved SNP data to process\n")
|
|
rliterman@0
|
319 isolate_snp_stats = pd.DataFrame(columns = ['Isolate_ID','Isolate_Type','Measure','Min','Value','Max','StdDev','Zscore','QC','Count'])
|
|
rliterman@0
|
320
|
|
rliterman@0
|
321 # Compare SNPs across refs
|
|
rliterman@0
|
322 if len(reference_ids) == 1:
|
|
rliterman@0
|
323 isolate_stdev_stats = pd.DataFrame(columns =['Isolate_ID','Isolate_Type','Measure','Min','Mean','Max','StdDev','Zscore','QC','Count'])
|
|
rliterman@0
|
324 with open(log_file,"a+") as log:
|
|
rliterman@0
|
325 log.write("- 1 reference provided, SNP distances have no comparisons\n")
|
|
rliterman@0
|
326 else:
|
|
rliterman@0
|
327 # Get comparison stats
|
|
rliterman@0
|
328 raw_comparison_df = raw_snp_distance_df.groupby(by=['Comparison'])['SNP_Distance'].agg(Count = 'count', Min = 'min', Mean = 'mean', Max = 'max', StdDev = 'std').reset_index()
|
|
rliterman@0
|
329 raw_comparison_df['StdDev'] = raw_comparison_df['StdDev'].astype('float').round(3)
|
|
rliterman@0
|
330 raw_comparison_df['Mean'] = raw_comparison_df['Mean'].astype('int')
|
|
rliterman@0
|
331 raw_comparison_df[['Query_1', 'Query_2']] = raw_comparison_df['Comparison'].str.split(';', expand=True)
|
|
rliterman@0
|
332 raw_comparison_df['SNP_Spread'] = abs(raw_comparison_df['Max'] - raw_comparison_df['Min'])
|
|
rliterman@0
|
333
|
|
rliterman@0
|
334 comparison_df = raw_comparison_df[['Comparison','Query_1','Query_2','Mean','StdDev','Min','Max','SNP_Spread','Count']].copy()
|
|
rliterman@0
|
335
|
|
rliterman@0
|
336 # Get isolate stats
|
|
rliterman@0
|
337 isolate_stdev_df = pd.DataFrame(columns = ['Isolate_ID','Count','Min','Value','Max','StdDev'])
|
|
rliterman@0
|
338
|
|
rliterman@0
|
339 for isolate in snp_isolates:
|
|
rliterman@0
|
340 temp_compare = raw_comparison_df[(raw_comparison_df['Query_1'] == isolate) | (raw_comparison_df['Query_2'] == isolate)].drop_duplicates(subset=['Comparison']).assign(Isolate_ID = isolate)
|
|
rliterman@0
|
341 temp_compare = temp_compare.groupby(by=['Isolate_ID'])['StdDev'].agg(Count = 'count', Min = 'min', Value = 'mean', Max = 'max', StdDev = 'std').reset_index()
|
|
rliterman@0
|
342 isolate_stdev_df = pd.concat([isolate_stdev_df,temp_compare])
|
|
rliterman@0
|
343
|
|
rliterman@0
|
344 isolate_stdev_df['Measure'] = "Raw_Distance_StdDev"
|
|
rliterman@0
|
345 isolate_stdev_df['Value'] = isolate_stdev_df['Value'].astype("float")
|
|
rliterman@0
|
346 isolate_stdev_df['Zscore'] = isolate_stdev_df['Value'].transform(scipy.stats.zscore).astype('float').round(3)
|
|
rliterman@0
|
347 isolate_stdev_df['Value'] = isolate_stdev_df['Value'].astype('float').round(3)
|
|
rliterman@0
|
348 isolate_stdev_df['Isolate_Type'] = isolate_stdev_df['Isolate_ID'].apply(lambda x: 'Reference' if x in reference_ids else 'Query')
|
|
rliterman@0
|
349 isolate_stdev_df['QC'] = getWarnings(isolate_stdev_df)
|
|
rliterman@0
|
350
|
|
rliterman@0
|
351 isolate_stdev_stats = isolate_stdev_df[['Isolate_ID','Isolate_Type','Measure','Min','Value','Max','StdDev','Zscore','QC','Count']].copy().rename(columns={'Value':'Mean'})
|
|
rliterman@0
|
352
|
|
rliterman@0
|
353 if has_preserved:
|
|
rliterman@0
|
354 comparison_df.columns = ['Comparison','Query_1','Query_2','Raw_Mean','Raw_StdDev','Raw_Min','Raw_Max','Raw_SNP_Spread','Raw_Count']
|
|
rliterman@0
|
355
|
|
rliterman@0
|
356 preserved_comparison_df = preserved_snp_distance_df.groupby(by=['Comparison'])['SNP_Distance'].agg(Preserved_Count = 'count', Preserved_Min = 'min', Preserved_Mean = 'mean', Preserved_Max = 'max', Preserved_StdDev = 'std').reset_index()
|
|
rliterman@0
|
357 preserved_comparison_df['Preserved_StdDev'] = preserved_comparison_df['Preserved_StdDev'].astype('float').round(3)
|
|
rliterman@0
|
358 preserved_comparison_df['Preserved_Mean'] = preserved_comparison_df['Preserved_Mean'].astype('int')
|
|
rliterman@0
|
359 preserved_comparison_df[['Query_1', 'Query_2']] = preserved_comparison_df['Comparison'].str.split(';', expand=True)
|
|
rliterman@0
|
360 preserved_comparison_df['Preserved_SNP_Spread'] = abs(preserved_comparison_df['Preserved_Max'] - preserved_comparison_df['Preserved_Min'])
|
|
rliterman@0
|
361
|
|
rliterman@0
|
362 comparison_df = comparison_df.merge(preserved_comparison_df,how = "left", on=['Comparison','Query_1','Query_2'])
|
|
rliterman@0
|
363 comparison_df['Mean_Preserved_Diff'] = abs(comparison_df['Preserved_Mean'] - comparison_df['Raw_Mean'])
|
|
rliterman@0
|
364 comparison_df = comparison_df[['Query_1','Query_2','Raw_Mean','Preserved_Mean','Mean_Preserved_Diff','Raw_StdDev','Preserved_StdDev','Raw_SNP_Spread','Preserved_SNP_Spread','Raw_Min','Raw_Max','Preserved_Min','Preserved_Max','Raw_Count','Preserved_Count']]
|
|
rliterman@0
|
365
|
|
rliterman@0
|
366 isolate_stdev_df = pd.DataFrame(columns = ['Isolate_ID','Count','Min','Value','Max','StdDev'])
|
|
rliterman@0
|
367
|
|
rliterman@0
|
368 for isolate in snp_isolates:
|
|
rliterman@0
|
369 temp_compare = preserved_comparison_df[(preserved_comparison_df['Query_1'] == isolate) | (preserved_comparison_df['Query_2'] == isolate)].drop_duplicates(subset=['Comparison']).assign(Isolate_ID = isolate)
|
|
rliterman@0
|
370 temp_compare = temp_compare.groupby(by=['Isolate_ID'])['Preserved_StdDev'].agg(Count = 'count', Min = 'min', Value = 'mean', Max = 'max', StdDev = 'std').reset_index()
|
|
rliterman@0
|
371 isolate_stdev_df = pd.concat([isolate_stdev_df,temp_compare])
|
|
rliterman@0
|
372
|
|
rliterman@0
|
373 isolate_stdev_df['Measure'] = "Preserved_Distance_StdDev"
|
|
rliterman@0
|
374 isolate_stdev_df['Value'] = isolate_stdev_df['Value'].astype('float')
|
|
rliterman@0
|
375 isolate_stdev_df['Zscore'] = isolate_stdev_df['Value'].transform(scipy.stats.zscore).astype('float').round(3)
|
|
rliterman@0
|
376 isolate_stdev_df['Value'] = isolate_stdev_df['Value'].astype('float').round(3)
|
|
rliterman@0
|
377 isolate_stdev_df['Isolate_Type'] = isolate_stdev_df['Isolate_ID'].apply(lambda x: 'Reference' if x in reference_ids else 'Query')
|
|
rliterman@0
|
378 isolate_stdev_df['QC'] = getWarnings(isolate_stdev_df)
|
|
rliterman@0
|
379 isolate_stdev_stats = pd.concat([isolate_stdev_stats,isolate_stdev_df[['Isolate_ID','Isolate_Type','Measure','Min','Value','Max','StdDev','Zscore','QC','Count']].copy().rename(columns={'Value':'Mean'})])
|
|
rliterman@0
|
380 with open(log_file,"a+") as log:
|
|
rliterman@0
|
381 log.write("- Compared results across references\n")
|
|
rliterman@0
|
382 else:
|
|
rliterman@0
|
383 with open(log_file,"a+") as log:
|
|
rliterman@0
|
384 log.write("- Compared results across references\n")
|
|
rliterman@0
|
385 # Group by ref
|
|
rliterman@0
|
386
|
|
rliterman@0
|
387 #### Isolate ####
|
|
rliterman@0
|
388 ref_isolate_df = isolate_stats.loc[isolate_stats['Isolate_Type'] == "Reference"][['Isolate_ID','Measure','Mean','StdDev','Min','Max','Zscore','QC','Count']].rename(columns = {'Isolate_ID':'Reference_ID'})
|
|
rliterman@0
|
389
|
|
rliterman@0
|
390 #### StdDev ####
|
|
rliterman@0
|
391 ref_stdev_df = isolate_stdev_stats.loc[isolate_stdev_stats['Isolate_Type'] == "Reference"][['Isolate_ID','Measure','Mean','StdDev','Min','Max','Zscore','QC','Count']].rename(columns = {'Isolate_ID':'Reference_ID'})
|
|
rliterman@0
|
392
|
|
rliterman@0
|
393 #### MUMmer ####
|
|
rliterman@0
|
394 ref_mummer_df = pd.DataFrame(columns = ['Reference_ID','Measure','Mean','StdDev','Min','Max','Count'])
|
|
rliterman@0
|
395 for ref in reference_ids:
|
|
rliterman@0
|
396 ref_mummer = isolate_mummer[(isolate_mummer['Isolate_ID'] == ref) | (isolate_mummer['Compare_ID'] == ref)].assign(Focal_Reference = ref)
|
|
rliterman@0
|
397 ref_mummer['Comparison'] = ref_mummer.apply(lambda row: ';'.join(sorted([str(row['Isolate_ID']), str(row['Compare_ID'])])), axis=1)
|
|
rliterman@0
|
398 ref_mummer = ref_mummer.drop_duplicates(subset=['Comparison'])
|
|
rliterman@0
|
399
|
|
rliterman@0
|
400 ref_mummer = ref_mummer.melt(id_vars=['Focal_Reference','Isolate_ID','Compare_ID'], value_vars = ['Align_Percent_Diff','Median_Alignment_Length','Kmer_Similarity','gIndels'],value_name='Value',var_name = "Measure")
|
|
rliterman@0
|
401 ref_mummer['Value'] = ref_mummer['Value'].astype("float")
|
|
rliterman@0
|
402 ref_mummer = ref_mummer.groupby(by=['Measure'])['Value'].agg(Count = "count",Min = "min",Mean = "mean",Max = "max",StdDev = 'std').reset_index().assign(Reference_ID = ref)
|
|
rliterman@0
|
403 ref_mummer = ref_mummer[['Reference_ID','Measure','Mean','StdDev','Min','Max','Count']]
|
|
rliterman@0
|
404 ref_mummer_df = pd.concat([ref_mummer_df,ref_mummer])
|
|
rliterman@0
|
405 ref_mummer_df['QC'] = np.nan
|
|
rliterman@0
|
406 ref_mummer_df['Zscore'] = np.nan
|
|
rliterman@0
|
407
|
|
rliterman@0
|
408 ref_mummer_summary_df = pd.concat([ref_mummer_df[['Reference_ID','Measure','Mean','StdDev','Min','Max','Zscore','QC','Count']],align_stats.loc[(align_stats['Isolate_Type'] == "Reference") & (align_stats['Measure'].isin(['Self_Aligned','Compare_Aligned','Unique_Kmers','Missing_Kmers']))][['Isolate_ID','Measure','Mean','StdDev','Min','Max','Zscore','QC','Count']].rename(columns = {'Isolate_ID':'Reference_ID'})])
|
|
rliterman@0
|
409
|
|
rliterman@0
|
410 #### Cocalled ####
|
|
rliterman@0
|
411 ref_cocalled_summary_df = raw_cocalled_df.groupby(by=['Reference_ID'])['SNPs_Cocalled'].agg(Mean = "mean",StdDev = 'std',Min = "min",Max = "max",Count = 'count').reset_index()
|
|
rliterman@0
|
412 ref_cocalled_summary_df['Measure'] = "Raw_SNPs_Cocalled"
|
|
rliterman@0
|
413 ref_cocalled_summary_df['QC'] = np.nan
|
|
rliterman@0
|
414 ref_cocalled_summary_df['Zscore'] = np.nan
|
|
rliterman@0
|
415 ref_cocalled_summary_df = ref_cocalled_summary_df[['Reference_ID','Measure','Mean','StdDev','Min','Max','Zscore','QC','Count']]
|
|
rliterman@0
|
416
|
|
rliterman@0
|
417 if has_preserved:
|
|
rliterman@0
|
418 preserved_cocalled_summary = preserved_cocalled_df.groupby(by=['Reference_ID'])['SNPs_Cocalled'].agg(Mean = "mean",StdDev = 'std',Min = "min",Max = "max",Count = 'count').reset_index()
|
|
rliterman@0
|
419 preserved_cocalled_summary['Measure'] = "Preserved_SNPs_Cocalled"
|
|
rliterman@0
|
420 preserved_cocalled_summary['QC'] = np.nan
|
|
rliterman@0
|
421 preserved_cocalled_summary['Zscore'] = np.nan
|
|
rliterman@0
|
422 ref_cocalled_summary_df = pd.concat([ref_cocalled_summary_df,preserved_cocalled_summary[['Reference_ID','Measure','Mean','StdDev','Min','Max','Zscore','QC','Count']]])
|
|
rliterman@0
|
423
|
|
rliterman@0
|
424 #### Preserved Diff ####
|
|
rliterman@0
|
425 if has_preserved:
|
|
rliterman@0
|
426 ref_summary_preserved_df = snp_df.groupby(by=['Reference_ID'])['Preserved_Diff'].agg(Mean = "mean",StdDev = 'std',Min = "min",Max = "max",Count = 'count').reset_index()
|
|
rliterman@0
|
427 ref_summary_preserved_df['Measure'] = "Preserved_Diff"
|
|
rliterman@0
|
428 ref_summary_preserved_df['QC'] = np.nan
|
|
rliterman@0
|
429 ref_summary_preserved_df['Zscore'] = np.nan
|
|
rliterman@0
|
430 ref_summary_preserved_df = ref_summary_preserved_df[['Reference_ID','Measure','Mean','StdDev','Min','Max','Count','Zscore','QC']].copy()
|
|
rliterman@0
|
431
|
|
rliterman@0
|
432 #### Compile ####
|
|
rliterman@0
|
433 ref_summary_df = pd.concat([ref_mummer_summary_df,
|
|
rliterman@0
|
434 ref_cocalled_summary_df,
|
|
rliterman@0
|
435 ref_isolate_df,ref_stdev_df]).sort_values(by=['Measure'])
|
|
rliterman@0
|
436
|
|
rliterman@0
|
437 ref_summary_df['Mean'] = ref_summary_df['Mean'].astype("float").round(3)
|
|
rliterman@0
|
438 ref_summary_df['Min'] = ref_summary_df['Min'].astype("float").round(3)
|
|
rliterman@0
|
439 ref_summary_df['Max'] = ref_summary_df['Max'].astype("float").round(3)
|
|
rliterman@0
|
440 ref_summary_df['StdDev'] = ref_summary_df['StdDev'].astype("float").round(3)
|
|
rliterman@0
|
441
|
|
rliterman@0
|
442 # Catch warnings and failures
|
|
rliterman@0
|
443 all_isolate_stats = pd.concat([isolate_stats,align_stats,isolate_stdev_stats,isolate_cocalled_stats,isolate_snp_stats]).sort_values(by=['Zscore'])
|
|
rliterman@0
|
444
|
|
rliterman@0
|
445 warn_fail_df = all_isolate_stats.loc[all_isolate_stats['QC'].isin(['Failure','Warning'])].copy()
|
|
rliterman@0
|
446 warn_fail_df['abs_Zscore'] = warn_fail_df['Zscore'].abs()
|
|
rliterman@0
|
447 warn_fail_df = warn_fail_df.sort_values(by='abs_Zscore',ascending=False).drop('abs_Zscore',axis=1)
|
|
rliterman@0
|
448
|
|
rliterman@0
|
449 warn_fail_isolates = list(set(warn_fail_df['Isolate_ID']))
|
|
rliterman@0
|
450 if len(warn_fail_isolates) > 0:
|
|
rliterman@0
|
451 with open(log_file,"a+") as log:
|
|
rliterman@0
|
452 log.write("\n- The following samples had QC warnings or failures:\n")
|
|
rliterman@0
|
453 for isolate in warn_fail_isolates:
|
|
rliterman@0
|
454 isolate_warn_fail = warn_fail_df.loc[warn_fail_df['Isolate_ID'] == isolate]
|
|
rliterman@0
|
455
|
|
rliterman@0
|
456 if isolate in reference_ids:
|
|
rliterman@0
|
457 log.write(f"\n{isolate} (Reference):\n")
|
|
rliterman@0
|
458 else:
|
|
rliterman@0
|
459 log.write(f"\n{isolate} (Query):\n")
|
|
rliterman@0
|
460
|
|
rliterman@0
|
461 for index,row in isolate_warn_fail.iterrows():
|
|
rliterman@0
|
462 log.write(f"\t- {row['Measure']} - Mean: {row['Mean']}; Zscore: {row['Zscore']}; QC: {row['QC']}\n")
|
|
rliterman@0
|
463 else:
|
|
rliterman@0
|
464 with open(log_file,"a+") as log:
|
|
rliterman@0
|
465 log.write("-There were no QC warnings or failures\n")
|
|
rliterman@0
|
466
|
|
rliterman@0
|
467 # Output data
|
|
rliterman@0
|
468
|
|
rliterman@0
|
469 # Mean assembly stats
|
|
rliterman@0
|
470 isolate_mean_df.reset_index().to_csv(mean_isolate_file,sep='\t',index=False)
|
|
rliterman@0
|
471
|
|
rliterman@0
|
472 # Isolate assembly stats
|
|
rliterman@0
|
473 isolate_assembly_stats = isolate_stats.loc[isolate_stats['Measure'].isin(['Contig_Count','Assembly_Bases','L50','L90','N50','N90'])].drop(['Min','Max','StdDev','Count'],axis=1).rename(columns = {'Mean':'Value'})
|
|
rliterman@0
|
474 isolate_assembly_stats.to_csv(isolate_assembly_stats_file,sep='\t',index=False)
|
|
rliterman@0
|
475
|
|
rliterman@0
|
476 # Isolate alignment stats
|
|
rliterman@0
|
477 isolate_align_stats = pd.concat([align_stats,isolate_cocalled_stats,isolate_snp_stats,isolate_stdev_stats]).reset_index(drop=True)
|
|
rliterman@0
|
478 for col in ['Min', 'Mean', 'Max', 'StdDev', 'Zscore']:
|
|
rliterman@0
|
479 isolate_align_stats[col] = isolate_align_stats[col].astype("float").round(3)
|
|
rliterman@0
|
480 isolate_align_stats.to_csv(align_stats_file,sep='\t',index=False)
|
|
rliterman@0
|
481
|
|
rliterman@0
|
482 # Reference Assembly Stats
|
|
rliterman@0
|
483 ref_align_summary_df = ref_summary_df.loc[(~ref_summary_df['Measure'].isin(['Contig_Count','Assembly_Bases','L50','L90','N50','N90'])) & (~pd.isna(ref_summary_df['Zscore']))]
|
|
rliterman@0
|
484 ref_mean_summary_df = ref_summary_df.loc[(~ref_summary_df['Measure'].isin(['Contig_Count','Assembly_Bases','L50','L90','N50','N90'])) & (pd.isna(ref_summary_df['Zscore']))].drop(['Zscore','QC'],axis =1)
|
|
rliterman@0
|
485 ref_mean_summary_df['Zscore'] = np.nan
|
|
rliterman@0
|
486 ref_mean_summary_df['QC'] = np.nan
|
|
rliterman@0
|
487
|
|
rliterman@0
|
488 # Add alignment stats
|
|
rliterman@0
|
489 if has_preserved:
|
|
rliterman@0
|
490 ref_mean_summary_df = pd.concat([ref_mean_summary_df,ref_summary_preserved_df])
|
|
rliterman@0
|
491
|
|
rliterman@0
|
492 ref_isolate_align_stats = align_stats.loc[(align_stats['Isolate_Type'] == "Reference") & (align_stats['Measure'].isin(['Self_Aligned','Compare_Aligned']))].drop(['Isolate_Type'],axis=1).rename(columns = {'Isolate_ID':'Reference_ID'})[['Reference_ID','Measure','Mean','StdDev','Min','Max','Count','Zscore','QC']]
|
|
rliterman@0
|
493
|
|
rliterman@0
|
494 ref_mean_summary_stats = pd.concat([ref_mean_summary_df,ref_isolate_align_stats])
|
|
rliterman@0
|
495 ref_mean_summary_stats.to_csv(ref_mean_summary_file,sep='\t',index=False)
|
|
rliterman@0
|
496
|
|
rliterman@0
|
497 end_time = time.time()
|
|
rliterman@0
|
498
|
|
rliterman@0
|
499 with open(log_file,"a+") as log:
|
|
rliterman@0
|
500 log.write(f"\n- Completed compilation in {end_time - start_time:.2f} seconds\n")
|
|
rliterman@0
|
501 log.write(f"\t- Saved mean isolate assembly data to {mean_isolate_file}\n")
|
|
rliterman@0
|
502 log.write(f"\t- Saved raw isolate assembly data to {isolate_assembly_stats_file}\n")
|
|
rliterman@0
|
503 log.write(f"\t- Saved isolate alignment data to {align_stats_file}\n")
|
|
rliterman@0
|
504 log.write(f"\t- Saved reference summary data to {ref_mean_summary_file}\n")
|
|
rliterman@0
|
505
|
|
rliterman@0
|
506 # Comparisons if multiple refs
|
|
rliterman@0
|
507 if len(reference_ids) > 1:
|
|
rliterman@0
|
508 comparison_df.to_csv(snp_comparison_file,sep="\t",index = False)
|
|
rliterman@0
|
509 log.write(f"\t- Saved SNP distance comparisons across references to {snp_comparison_file}\n")
|
|
rliterman@0
|
510
|
|
rliterman@0
|
511 # Failures/warnings
|
|
rliterman@0
|
512 if warn_fail_df.shape[0] > 0:
|
|
rliterman@0
|
513 warn_fail_df.to_csv(qc_file,sep="\t",index=False)
|
|
rliterman@0
|
514 log.write(f"\t- Saved QC warnings/failures to {qc_file}\n") |
