Mercurial > repos > iss > eurl_vtec_wgs_pt

comparison scripts/modules/run_rematch.py @ 0:c6bab5103a14 draft

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
"planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
author iss
date Mon, 21 Mar 2022 15:23:09 +0000
parents
children
comparison
equal deleted inserted replaced
-1:000000000000 0:c6bab5103a14
1 import functools
2 import os
3 import sys
4 import multiprocessing
5
6 try:
7 import modules.utils as utils
8 except ImportError:
9 from pathotyping.modules import utils as utils
10
11
12 # {'noMatter': '/home/ubuntu/NGStools/patho_typing/mpmachado_stuff.out_test/rematch/sample.noMatter.fasta', 'correct': '/home/ubuntu/NGStools/patho_typing/mpmachado_stuff.out_test/rematch/sample.correct.fasta', 'alignment': '/home/ubuntu/NGStools/patho_typing/mpmachado_stuff.out_test/rematch/sample.alignment.fasta'}
13
14
15 def remove_alignment(alignment_file):
16 directory = os.path.dirname(alignment_file)
17 files = [f for f in os.listdir(directory) if not f.startswith('.') and os.path.isfile(os.path.join(directory, f))]
18 for file_found in files:
19 if file_found.startswith(os.path.splitext(os.path.basename(alignment_file))[0]):
20 file_found = os.path.join(directory, file_found)
21 os.remove(file_found)
22
23
24 def remove_reference_stuff(outdir, reference_file):
25 files = [f for f in os.listdir(outdir) if not f.startswith('.') and os.path.isfile(os.path.join(outdir, f))]
26 for file_found in files:
27 if file_found.startswith(os.path.splitext(os.path.basename(reference_file))[0]):
28 file_found = os.path.join(outdir, file_found)
29 os.remove(file_found)
30
31
32 def clean_rematch_folder(consensus_files, bam_file, reference_file, outdir, doNotRemoveConsensus, debug_mode_true):
33 if not debug_mode_true:
34 if not doNotRemoveConsensus:
35 for consensus_type, file_path in list(consensus_files.items()):
36 if os.path.isfile(file_path):
37 os.remove(file_path)
38 if bam_file is not None:
39 remove_alignment(bam_file)
40 remove_reference_stuff(outdir, reference_file)
41
42
43 def sequence_data(sample, reference_file, bam_file, outdir, threads, length_extra_seq, minimum_depth_presence, minimum_depth_call, minimum_depth_frequency_dominant_allele, debug_mode_true, rematch):
44 sequence_data_outdir = os.path.join(outdir, 'sequence_data', '')
45 utils.removeDirectory(sequence_data_outdir)
46 os.mkdir(sequence_data_outdir)
47
48 sequences, headers = utils.get_sequence_information(reference_file, length_extra_seq)
49
50 pool = multiprocessing.Pool(processes=threads)
51 for sequence_counter in sequences:
52 sequence_dir = os.path.join(sequence_data_outdir, str(sequence_counter), '')
53 utils.removeDirectory(sequence_dir)
54 os.makedirs(sequence_dir)
55 pool.apply_async(rematch.analyse_sequence_data, args=(bam_file, sequences[sequence_counter], sequence_dir, sequence_counter, reference_file, length_extra_seq, minimum_depth_presence, minimum_depth_call, minimum_depth_frequency_dominant_allele,))
56 pool.close()
57 pool.join()
58
59 run_successfully, sample_data, consensus_files, consensus_sequences = rematch.gather_data_together(sample, sequence_data_outdir, sequences, outdir.rsplit('/', 2)[0], debug_mode_true, length_extra_seq, False)
60
61 return run_successfully, sample_data, consensus_files, consensus_sequences
62
63
64 def determine_general_statistics(sample_data, minimum_gene_coverage, minimum_gene_identity):
65 print('Writing report file')
66 number_absent_genes = 0
67 number_genes_multiple_alleles = 0
68 mean_sample_coverage = 0
69
70 with open('output_dir/rematch/rematchModule_report.txt', 'wt') as writer:
71 writer.write('\t'.join(['#gene', 'percentage_gene_coverage', 'gene_mean_read_coverage', 'percentage_gene_low_coverage', 'number_positions_multiple_alleles', 'percentage_gene_identity']) + '\n')
72 for i in range(1, len(sample_data) + 1):
73 writer.write('\t'.join([sample_data[i]['header'], str(round(sample_data[i]['gene_coverage'], 2)), str(round(sample_data[i]['gene_mean_read_coverage'], 2)), str(round(sample_data[i]['gene_low_coverage'], 2)), str(sample_data[i]['gene_number_positions_multiple_alleles']), str(round(sample_data[i]['gene_identity'], 2))]) + '\n')
74
75 if sample_data[i]['gene_coverage'] < minimum_gene_coverage or sample_data[i]['gene_identity'] < minimum_gene_identity:
76 number_absent_genes += 1
77 else:
78 mean_sample_coverage += sample_data[i]['gene_mean_read_coverage']
79 if sample_data[i]['gene_number_positions_multiple_alleles'] > 0:
80 number_genes_multiple_alleles += 1
81
82 if len(sample_data) - number_absent_genes > 0:
83 mean_sample_coverage = float(mean_sample_coverage) / float(len(sample_data) - number_absent_genes)
84 else:
85 mean_sample_coverage = 0
86
87 writer.write('\n'.join(['#general', '>number_absent_genes', str(number_absent_genes), '>number_genes_multiple_alleles', str(number_genes_multiple_alleles), '>mean_sample_coverage', str(round(mean_sample_coverage, 2))]) + '\n')
88
89 print('\n'.join([str('number_absent_genes: ' + str(number_absent_genes)),
90 str('number_genes_multiple_alleles: ' + str(number_genes_multiple_alleles)),
91 str('mean_sample_coverage: ' + str(round(mean_sample_coverage, 2)))]) + '\n')
92
93 return number_absent_genes, number_genes_multiple_alleles, mean_sample_coverage
94
95 module_timer = functools.partial(utils.timer, name='Module ReMatCh')
96
97
98 @module_timer
99 def run_rematch(rematch, outdir, reference_file, bam_file, threads, length_extra_seq, minimum_depth_presence, minimum_depth_call, minimum_depth_frequency_dominant_allele, minimum_gene_coverage, minimum_gene_identity, debug_mode_true, doNotRemoveConsensus):
100 module_dir = os.path.join(outdir, 'rematch', '')
101 utils.removeDirectory(module_dir)
102 os.makedirs(module_dir)
103
104 sys.path.append(os.path.join(os.path.dirname(rematch), 'modules'))
105 import rematch_module as rematch
106
107 print('Analysing alignment data')
108 run_successfully, sample_data, consensus_files, consensus_sequences = sequence_data('sample', reference_file, bam_file, module_dir, threads, length_extra_seq, minimum_depth_presence, minimum_depth_call, minimum_depth_frequency_dominant_allele, debug_mode_true, rematch)
109
110 if run_successfully:
111 number_absent_genes, number_genes_multiple_alleles, mean_sample_coverage = \
112 determine_general_statistics(sample_data=sample_data, minimum_gene_coverage=minimum_gene_coverage,
113 minimum_gene_identity=minimum_gene_identity)
114
115 if not debug_mode_true:
116 utils.removeDirectory(module_dir)
117
118 clean_rematch_folder(consensus_files, bam_file, reference_file, outdir, doNotRemoveConsensus, debug_mode_true)
119
120 return run_successfully, {'number_absent_genes': number_absent_genes if 'number_absent_genes' in locals() else None, 'number_genes_multiple_alleles': number_genes_multiple_alleles if 'number_genes_multiple_alleles' in locals() else None, 'mean_sample_coverage': round(mean_sample_coverage, 2) if 'mean_sample_coverage' in locals() else None}, sample_data if 'sample_data' in locals() else None