Mercurial > repos > cstrittmatter > test_eurl_vtec_wgs_pt

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

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
planemo upload commit 15239f1674081ab51ab8dd75a9a40cf1bfaa93e8
author cstrittmatter
date Wed, 22 Jan 2020 08:41:44 -0500
parents
children 0cbed1c0a762
comparison
equal deleted inserted replaced
-1:000000000000 0:965517909457
1 import functools
2 import os
3 import sys
4
5 import utils
6
7 def remove_alignment(alignment_file):
8 directory = os.path.dirname(alignment_file)
9 files = [f for f in os.listdir(directory) if not f.startswith('.') and os.path.isfile(os.path.join(directory, f))]
10 for file_found in files:
11 if file_found.startswith(os.path.splitext(os.path.basename(alignment_file))[0]):
12 file_found = os.path.join(directory, file_found)
13 os.remove(file_found)
14
15
16 def remove_reference_stuff(outdir, reference_file):
17 files = [f for f in os.listdir(outdir) if not f.startswith('.') and os.path.isfile(os.path.join(outdir, f))]
18 for file_found in files:
19 if file_found.startswith(os.path.splitext(os.path.basename(reference_file))[0]):
20 file_found = os.path.join(outdir, file_found)
21 os.remove(file_found)
22
23
24 def clean_rematch_folder(consensus_files, bam_file, reference_file, outdir, doNotRemoveConsensus, debug_mode_true):
25 if not debug_mode_true:
26 if not doNotRemoveConsensus:
27 for consensus_type, file_path in consensus_files.items():
28 if os.path.isfile(file_path):
29 os.remove(file_path)
30 if bam_file is not None:
31 remove_alignment(bam_file)
32 remove_reference_stuff(outdir, reference_file)
33
34
35 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):
36 sequence_data_outdir = os.path.join(outdir, 'sequence_data', '')
37 utils.removeDirectory(sequence_data_outdir)
38 os.mkdir(sequence_data_outdir)
39
40 sequences, headers = utils.get_sequence_information(reference_file, length_extra_seq)
41
42 threads_2_use = rematch.determine_threads_2_use(len(sequences), threads)
43
44 import multiprocessing
45
46 pool = multiprocessing.Pool(processes=threads)
47 for sequence_counter in sequences:
48 sequence_dir = os.path.join(sequence_data_outdir, str(sequence_counter), '')
49 utils.removeDirectory(sequence_dir)
50 os.makedirs(sequence_dir)
51 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, threads_2_use,))
52 pool.close()
53 pool.join()
54
55 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)
56
57 return run_successfully, sample_data, consensus_files, consensus_sequences
58
59
60 def write_report(outdir, sample_data, minimum_gene_coverage, minimum_gene_identity):
61 print 'Writing report file'
62 number_absent_genes = 0
63 number_genes_multiple_alleles = 0
64 mean_sample_coverage = 0
65 with open(os.path.join(outdir, 'rematchModule_report.txt'), 'wt') as writer:
66 writer.write('\t'.join(['#gene', 'percentage_gene_coverage', 'gene_mean_read_coverage', 'percentage_gene_low_coverage', 'number_positions_multiple_alleles', 'percentage_gene_identity']) + '\n')
67 for i in range(1, len(sample_data) + 1):
68 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')
69
70 if sample_data[i]['gene_coverage'] < minimum_gene_coverage or sample_data[i]['gene_identity'] < minimum_gene_identity:
71 number_absent_genes += 1
72 else:
73 mean_sample_coverage += sample_data[i]['gene_mean_read_coverage']
74 if sample_data[i]['gene_number_positions_multiple_alleles'] > 0:
75 number_genes_multiple_alleles += 1
76
77 if len(sample_data) - number_absent_genes > 0:
78 mean_sample_coverage = float(mean_sample_coverage) / float(len(sample_data) - number_absent_genes)
79 else:
80 mean_sample_coverage = 0
81
82 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')
83
84 print '\n'.join([str('number_absent_genes: ' + str(number_absent_genes)), str('number_genes_multiple_alleles: ' + str(number_genes_multiple_alleles)), str('mean_sample_coverage: ' + str(round(mean_sample_coverage, 2)))]) + '\n'
85
86 return number_absent_genes, number_genes_multiple_alleles, mean_sample_coverage
87
88
89 module_timer = functools.partial(utils.timer, name='Module ReMatCh')
90
91
92 @module_timer
93 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):
94 module_dir = os.path.join(outdir, 'rematch', '')
95 utils.removeDirectory(module_dir)
96 os.makedirs(module_dir)
97
98 sys.path.append(os.path.join(os.path.dirname(rematch), 'modules', ''))
99 import rematch_module as rematch
100
101 print 'Analysing alignment data'
102 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)
103
104 if run_successfully:
105 number_absent_genes, number_genes_multiple_alleles, mean_sample_coverage = write_report(outdir, sample_data, minimum_gene_coverage, minimum_gene_identity)
106
107 if not debug_mode_true:
108 utils.removeDirectory(module_dir)
109
110 clean_rematch_folder(consensus_files, bam_file, reference_file, outdir, doNotRemoveConsensus, debug_mode_true)
111
112 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