annotate scripts/ReMatCh/modules/rematch_module.py @ 6:20ff3dca457f draft default tip

planemo upload commit 6857c749c21f580c828aba3543e294b69d32b662
author iss
date Mon, 23 Oct 2023 11:45:36 +0000
parents c6bab5103a14
children
Ignore whitespace changes - Everywhere: Within whitespace: At end of lines:
rev   line source
0
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
1 import os.path
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
2 import multiprocessing
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
3 import functools
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
4 import sys
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
5 import pickle
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
6
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
7 # https://chrisyeh96.github.io/2017/08/08/definitive-guide-python-imports.html#case-2-syspath-could-change
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
8 sys.path.insert(0, os.path.dirname(os.path.realpath(__file__)))
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
9 import utils
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
10
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
11
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
12 def index_fasta_samtools(fasta, region_none, region_outfile_none, print_comand_true):
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
13 command = ['samtools', 'faidx', fasta, '', '', '']
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
14 shell_true = False
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
15 if region_none is not None:
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
16 command[3] = region_none
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
17 if region_outfile_none is not None:
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
18 command[4] = '>'
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
19 command[5] = region_outfile_none
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
20 shell_true = True
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
21 run_successfully, stdout, stderr = utils.run_command_popen_communicate(command, shell_true, None, print_comand_true)
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
22 return run_successfully, stdout
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
23
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
24
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
25 # Indexing reference file using Bowtie2
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
26 def index_sequence_bowtie2(reference_file, threads):
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
27 if os.path.isfile(str(reference_file + '.1.bt2')):
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
28 run_successfully = True
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
29 else:
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
30 command = ['bowtie2-build', '--threads', str(threads), reference_file, reference_file]
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
31 run_successfully, stdout, stderr = utils.run_command_popen_communicate(command, False, None, True)
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
32 return run_successfully
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
33
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
34
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
35 # Mapping with Bowtie2
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
36 def mapping_bowtie2(fastq_files, reference_file, threads, outdir, num_map_loc,
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
37 bowtie_algorithm='--very-sensitive-local', bowtie_opt=None):
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
38 sam_file = os.path.join(outdir, str('alignment.sam'))
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
39
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
40 # Index reference file
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
41 run_successfully = index_sequence_bowtie2(reference_file, threads)
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
42
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
43 if run_successfully:
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
44 command = ['bowtie2', '', '', '-q', bowtie_algorithm, '--threads', str(threads), '-x',
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
45 reference_file, '', '--no-unal', '', '-S', sam_file]
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
46
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
47 if num_map_loc is not None and num_map_loc > 1:
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
48 command[1] = '-k'
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
49 command[2] = str(num_map_loc)
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
50
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
51 if len(fastq_files) == 1:
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
52 command[9] = '-U ' + fastq_files[0]
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
53 elif len(fastq_files) == 2:
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
54 command[9] = '-1 ' + fastq_files[0] + ' -2 ' + fastq_files[1]
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
55 else:
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
56 return False, None
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
57
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
58 if bowtie_opt is not None:
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
59 command[11] = bowtie_opt
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
60
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
61 run_successfully, stdout, stderr = utils.run_command_popen_communicate(command, False, None, True)
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
62
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
63 if not run_successfully:
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
64 sam_file = None
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
65
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
66 return run_successfully, sam_file
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
67
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
68
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
69 def split_cigar(cigar):
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
70 cigars = ['M', 'I', 'D', 'N', 'S', 'H', 'P', '=', 'X']
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
71
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
72 splited_cigars = []
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
73 numbers = ''
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
74 for char in cigar:
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
75 if char not in cigars:
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
76 numbers += char
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
77 else:
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
78 splited_cigars.append([int(numbers), char])
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
79 numbers = ''
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
80
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
81 return splited_cigars
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
82
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
83
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
84 def recode_cigar_based_on_base_quality(cigar, bases_quality, soft_clip_base_quality, mapping_position,
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
85 direct_strand_true, soft_clip_cigar_flag_recode):
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
86 cigar = split_cigar(cigar)
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
87 soft_left = []
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
88 soft_right = []
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
89 cigar_flags_for_reads_length = ('M', 'I', 'S', '=', 'X')
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
90 read_length_without_right_s = sum([cigar_part[0] for cigar_part in cigar if
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
91 cigar_part[1] in cigar_flags_for_reads_length]) - \
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
92 (cigar[len(cigar) - 1][0] if cigar[len(cigar) - 1][1] == 'S' else 0)
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
93 for x, base in enumerate(bases_quality):
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
94 if ord(base) - 33 >= soft_clip_base_quality:
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
95 if x <= cigar[0][0] - 1:
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
96 if cigar[0][1] == 'S':
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
97 soft_left.append(x)
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
98 elif x > read_length_without_right_s - 1:
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
99 if cigar[len(cigar) - 1][1] == 'S':
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
100 soft_right.append(x)
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
101
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
102 left_changed = (False, 0)
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
103 if len(soft_left) > 0:
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
104 soft_left = min(soft_left) + 1
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
105 if soft_left == 1:
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
106 cigar = [[cigar[0][0], soft_clip_cigar_flag_recode]] + cigar[1:]
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
107 left_changed = (True, cigar[0][0])
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
108 elif cigar[0][0] - soft_left > 0:
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
109 cigar = [[soft_left, 'S']] + [[cigar[0][0] - soft_left, soft_clip_cigar_flag_recode]] + cigar[1:]
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
110 left_changed = (True, cigar[0][0] - soft_left)
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
111
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
112 right_changed = (False, 0)
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
113 if len(soft_right) > 0:
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
114 soft_right = max(soft_right) + 1
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
115 cigar = cigar[:-1]
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
116 if soft_right - read_length_without_right_s > 0:
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
117 cigar.append([soft_right - read_length_without_right_s, soft_clip_cigar_flag_recode])
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
118 right_changed = (True, soft_right - read_length_without_right_s)
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
119 if len(bases_quality) - soft_right > 0:
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
120 cigar.append([len(bases_quality) - soft_right, 'S'])
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
121
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
122 if left_changed[0]:
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
123 # if direct_strand_true:
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
124 mapping_position = mapping_position - left_changed[1]
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
125 # if right_changed[0]:
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
126 # if not direct_strand_true:
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
127 # mapping_position = mapping_position + right_changed[1]
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
128
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
129 return ''.join([str(cigar_part[0]) + cigar_part[1] for cigar_part in cigar]), str(mapping_position)
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
130
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
131
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
132 def verify_is_forward_read(sam_flag_bit):
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
133 # 64 = 1000000
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
134 forward_read = False
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
135 bit = format(sam_flag_bit, 'b').zfill(7)
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
136 if bit[-7] == '1':
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
137 forward_read = True
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
138 return forward_read
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
139
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
140
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
141 def verify_mapped_direct_strand(sam_flag_bit):
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
142 # 16 = 10000 -> mapped in reverse strand
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
143 direct_strand = False
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
144 bit = format(sam_flag_bit, 'b').zfill(5)
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
145 if bit[-5] == '0':
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
146 direct_strand = True
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
147 return direct_strand
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
148
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
149
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
150 def verify_mapped_tip(reference_length, mapping_position, cigar):
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
151 tip = False
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
152 if 'S' in cigar:
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
153 cigar = split_cigar(cigar)
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
154 if cigar[0][1] == 'S':
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
155 if mapping_position - cigar[0][0] < 0:
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
156 tip = True
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
157 if cigar[len(cigar) - 1][1] == 'S':
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
158 if mapping_position + cigar[len(cigar) - 1][0] > reference_length:
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
159 tip = True
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
160 return tip
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
161
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
162
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
163 def change_sam_flag_bit_mapped_reverse_strand_2_direct_strand(sam_flag_bit):
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
164 bit = list(format(sam_flag_bit, 'b').zfill(5))
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
165 bit[-5] = '0'
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
166 return int(''.join(bit), 2)
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
167
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
168
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
169 def change_sam_flag_bit_mate_reverse_strand_2_direct_strand(sam_flag_bit):
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
170 bit = list(format(sam_flag_bit, 'b').zfill(6))
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
171 bit[-6] = '0'
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
172 return int(''.join(bit), 2)
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
173
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
174
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
175 def move_read_mapped_reverse_strand_2_direct_strand(seq, bases_quality, sam_flag_bit, cigar):
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
176 seq = utils.reverse_complement(seq)
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
177 bases_quality = ''.join(reversed(list(bases_quality)))
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
178 sam_flag_bit = change_sam_flag_bit_mapped_reverse_strand_2_direct_strand(sam_flag_bit)
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
179 cigar = ''.join([str(cigar_part[0]) + cigar_part[1] for cigar_part in reversed(split_cigar(cigar))])
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
180 return seq, bases_quality, str(sam_flag_bit), cigar
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
181
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
182
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
183 @utils.trace_unhandled_exceptions
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
184 def parallelized_recode_soft_clipping(line_collection, pickle_file, soft_clip_base_quality, sequences_length,
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
185 soft_clip_cigar_flag_recode):
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
186 lines_sam = []
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
187 for line in line_collection:
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
188 line = line.rstrip('\r\n')
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
189 if len(line) > 0:
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
190 if line.startswith('@'):
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
191 lines_sam.append(line)
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
192 else:
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
193 line = line.split('\t')
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
194 if not verify_mapped_tip(sequences_length[line[2]], int(line[3]), line[5]):
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
195 line[5], line[3] = recode_cigar_based_on_base_quality(line[5], line[10], soft_clip_base_quality,
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
196 int(line[3]),
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
197 verify_mapped_direct_strand(int(line[1])),
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
198 soft_clip_cigar_flag_recode)
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
199 lines_sam.append('\t'.join(line))
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
200 with open(pickle_file, 'wb') as writer:
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
201 pickle.dump(lines_sam, writer)
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
202
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
203
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
204 def recode_soft_clipping_from_sam(sam_file, outdir, threads, soft_clip_base_quality, reference_dict,
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
205 soft_clip_cigar_flag_recode):
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
206 pickle_files = []
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
207 sequences_length = {}
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
208 for x, seq_info in list(reference_dict.items()):
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
209 sequences_length[seq_info['header']] = seq_info['length']
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
210
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
211 with open(sam_file, 'rtU') as reader:
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
212 pool = multiprocessing.Pool(processes=threads)
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
213 line_collection = []
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
214 x = 0
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
215 for x, line in enumerate(reader):
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
216 line_collection.append(line)
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
217 if x % 10000 == 0:
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
218 pickle_file = os.path.join(outdir, 'remove_soft_clipping.' + str(x) + '.pkl')
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
219 pickle_files.append(pickle_file)
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
220 pool.apply_async(parallelized_recode_soft_clipping, args=(line_collection, pickle_file,
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
221 soft_clip_base_quality, sequences_length,
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
222 soft_clip_cigar_flag_recode,))
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
223 line_collection = []
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
224 if len(line_collection) > 0:
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
225 pickle_file = os.path.join(outdir, 'remove_soft_clipping.' + str(x) + '.pkl')
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
226 pickle_files.append(pickle_file)
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
227 pool.apply_async(parallelized_recode_soft_clipping, args=(line_collection, pickle_file,
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
228 soft_clip_base_quality, sequences_length,
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
229 soft_clip_cigar_flag_recode,))
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
230 pool.close()
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
231 pool.join()
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
232
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
233 os.remove(sam_file)
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
234
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
235 new_sam_file = os.path.join(outdir, 'alignment_with_soft_clipping_recoded.sam')
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
236 with open(new_sam_file, 'wt') as writer:
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
237 for pickle_file in pickle_files:
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
238 if os.path.isfile(pickle_file):
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
239 lines_sam = None
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
240 with open(pickle_file, 'rb') as reader:
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
241 lines_sam = pickle.load(reader)
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
242 if lines_sam is not None:
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
243 for line in lines_sam:
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
244 writer.write(line + '\n')
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
245 os.remove(pickle_file)
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
246
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
247 return new_sam_file
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
248
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
249
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
250 # Sort alignment file
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
251 def sort_alignment(alignment_file, output_file, sort_by_name_true, threads):
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
252 out_format_string = os.path.splitext(output_file)[1][1:].lower()
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
253 command = ['samtools', 'sort', '-o', output_file, '-O', out_format_string, '', '-@', str(threads), alignment_file]
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
254 if sort_by_name_true:
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
255 command[6] = '-n'
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
256 run_successfully, stdout, stderr = utils.run_command_popen_communicate(command, False, None, True)
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
257 if not run_successfully:
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
258 output_file = None
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
259 return run_successfully, output_file
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
260
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
261
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
262 # Index alignment file
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
263 def index_alignment(alignment_file):
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
264 command = ['samtools', 'index', alignment_file]
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
265 run_successfully, stdout, stderr = utils.run_command_popen_communicate(command, False, None, True)
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
266 return run_successfully
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
267
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
268
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
269 def mapping_reads(fastq_files, reference_file, threads, outdir, num_map_loc, rematch_run,
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
270 soft_clip_base_quality, soft_clip_recode_run, reference_dict, soft_clip_cigar_flag_recode,
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
271 bowtie_algorithm, bowtie_opt, clean_run=True):
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
272 # Create a symbolic link to the reference_file
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
273 if clean_run:
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
274 reference_link = os.path.join(outdir, os.path.basename(reference_file))
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
275 if os.path.islink(reference_link):
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
276 os.unlink(reference_link)
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
277 os.symlink(reference_file, reference_link)
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
278 reference_file = reference_link
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
279
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
280 bam_file = None
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
281 # Mapping reads using Bowtie2
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
282 run_successfully, sam_file = mapping_bowtie2(fastq_files=fastq_files, reference_file=reference_file,
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
283 threads=threads, outdir=outdir, num_map_loc=num_map_loc,
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
284 bowtie_algorithm=bowtie_algorithm, bowtie_opt=bowtie_opt)
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
285
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
286 if run_successfully:
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
287 # Remove soft clipping
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
288 if rematch_run == soft_clip_recode_run or soft_clip_recode_run == 'both':
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
289 print('Recoding soft clipped regions')
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
290 sam_file = recode_soft_clipping_from_sam(sam_file, outdir, threads, soft_clip_base_quality, reference_dict,
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
291 soft_clip_cigar_flag_recode)
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
292
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
293 # Convert sam to bam and sort bam
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
294 run_successfully, bam_file = sort_alignment(sam_file, str(os.path.splitext(sam_file)[0] + '.bam'), False,
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
295 threads)
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
296
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
297 if run_successfully:
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
298 os.remove(sam_file)
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
299 # Index bam
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
300 run_successfully = index_alignment(bam_file)
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
301
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
302 return run_successfully, bam_file, reference_file
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
303
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
304
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
305 def create_vcf(bam_file, sequence_to_analyse, outdir, counter, reference_file):
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
306 gene_vcf = os.path.join(outdir, 'samtools_mpileup.sequence_' + str(counter) + '.vcf')
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
307
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
308 command = ['samtools', 'mpileup', '--count-orphans', '--no-BAQ', '--min-BQ', '0', '--min-MQ', str(7), '--fasta-ref',
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
309 reference_file, '--region', sequence_to_analyse, '--output', gene_vcf, '--VCF', '--uncompressed',
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
310 '--output-tags', 'INFO/AD,AD,DP', bam_file]
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
311
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
312 run_successfully, stdout, stderr = utils.run_command_popen_communicate(command, False, None, False)
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
313 if not run_successfully:
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
314 gene_vcf = None
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
315 return run_successfully, gene_vcf
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
316
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
317
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
318 # Read vcf file
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
319 class Vcf:
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
320 def __init__(self, vcf_file, encoding=None, newline=None):
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
321 try:
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
322 self.vcf = open(vcf_file, 'rt', encoding=encoding, newline=newline)
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
323 except TypeError:
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
324 self.vcf = open(vcf_file, 'rt')
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
325 self.line_read = self.vcf.readline()
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
326 self.contigs_info_dict = {}
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
327 while self.line_read.startswith('#'):
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
328 if self.line_read.startswith('##contig=<ID='):
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
329 seq = self.line_read.split('=')[2].split(',')[0]
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
330 seq_len = self.line_read.split('=')[3].split('>')[0]
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
331 self.contigs_info_dict[seq] = int(seq_len)
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
332 self.line_read = self.vcf.readline()
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
333 self.line = self.line_read
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
334
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
335 def readline(self):
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
336 line_stored = self.line
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
337 self.line = self.vcf.readline()
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
338 return line_stored
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
339
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
340 def close(self):
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
341 self.vcf.close()
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
342
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
343 def get_contig_legth(self, contig):
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
344 return self.contigs_info_dict[contig]
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
345
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
346
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
347 def get_variants(gene_vcf, seq_name, encoding=None, newline=None):
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
348 variants = {}
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
349
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
350 vfc_file = Vcf(vcf_file=gene_vcf, encoding=encoding, newline=newline)
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
351 line = vfc_file.readline()
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
352 counter = 1
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
353 while len(line) > 0:
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
354 fields = line.rstrip('\r\n').split('\t')
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
355 if len(fields) > 0:
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
356 fields[1] = int(fields[1])
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
357
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
358 info_field = {}
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
359 try:
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
360 for i in fields[7].split(';'):
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
361 i = i.split('=')
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
362 if len(i) > 1:
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
363 info_field[i[0]] = i[1]
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
364 else:
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
365 info_field[i[0]] = None
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
366 except IndexError:
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
367 if counter > vfc_file.get_contig_legth(contig=seq_name):
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
368 break
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
369 else:
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
370 raise IndexError
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
371
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
372 format_field = {}
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
373 format_field_name = fields[8].split(':')
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
374 format_data = fields[9].split(':')
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
375
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
376 for i in range(0, len(format_data)):
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
377 format_field[format_field_name[i]] = format_data[i].split(',')
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
378
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
379 fields_to_store = {'REF': fields[3], 'ALT': fields[4].split(','), 'info': info_field,
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
380 'format': format_field}
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
381 if fields[1] in variants:
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
382 variants[fields[1]][len(variants[fields[1]])] = fields_to_store
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
383 else:
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
384 variants[fields[1]] = {0: fields_to_store}
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
385
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
386 try:
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
387 line = vfc_file.readline()
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
388 except UnicodeDecodeError:
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
389 if counter + 1 > vfc_file.get_contig_legth(contig=seq_name):
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
390 break
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
391 else:
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
392 raise UnicodeDecodeError
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
393
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
394 counter += 1
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
395 vfc_file.close()
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
396
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
397 return variants
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
398
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
399
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
400 def indel_entry(variant_position):
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
401 entry_with_indel = []
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
402 entry_with_snp = None
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
403 for i in variant_position:
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
404 keys = list(variant_position[i]['info'].keys())
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
405 if 'INDEL' in keys:
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
406 entry_with_indel.append(i)
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
407 else:
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
408 entry_with_snp = i
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
409
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
410 return entry_with_indel, entry_with_snp
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
411
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
412
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
413 def get_alt_no_matter(variant_position, indel_true):
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
414 dp = sum(map(int, variant_position['format']['AD']))
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
415 index_alleles_sorted_position = sorted(zip(list(map(int, variant_position['format']['AD'])),
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
416 list(range(0, len(variant_position['format']['AD'])))),
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
417 reverse=True)
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
418 index_dominant_allele = None
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
419 if not indel_true:
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
420 ad_idv = index_alleles_sorted_position[0][0]
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
421
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
422 if len([x for x in index_alleles_sorted_position if x[0] == ad_idv]) > 1:
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
423 index_alleles_sorted_position = sorted([x for x in index_alleles_sorted_position if x[0] == ad_idv])
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
424
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
425 index_dominant_allele = index_alleles_sorted_position[0][1]
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
426 if index_dominant_allele == 0:
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
427 alt = '.'
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
428 else:
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
429 alt = variant_position['ALT'][index_dominant_allele - 1]
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
430
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
431 else:
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
432 ad_idv = int(variant_position['info']['IDV'])
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
433
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
434 if float(ad_idv) / float(dp) >= 0.5:
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
435 if len([x for x in index_alleles_sorted_position if x[0] == index_alleles_sorted_position[0][0]]) > 1:
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
436 index_alleles_sorted_position = sorted([x for x in index_alleles_sorted_position if
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
437 x[0] == index_alleles_sorted_position[0][0]])
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
438
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
439 index_dominant_allele = index_alleles_sorted_position[0][1]
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
440 if index_dominant_allele == 0:
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
441 alt = '.'
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
442 else:
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
443 alt = variant_position['ALT'][index_dominant_allele - 1]
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
444 else:
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
445 ad_idv = int(variant_position['format']['AD'][0])
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
446 alt = '.'
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
447
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
448 return alt, dp, ad_idv, index_dominant_allele
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
449
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
450
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
451 def count_number_diferences(ref, alt):
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
452 number_diferences = 0
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
453
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
454 if len(ref) != len(alt):
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
455 number_diferences += 1
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
456
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
457 for i in range(0, min(len(ref), len(alt))):
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
458 if alt[i] != 'N' and ref[i] != alt[i]:
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
459 number_diferences += 1
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
460
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
461 return number_diferences
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
462
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
463
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
464 def get_alt_correct(variant_position, alt_no_matter, dp, ad_idv, index_dominant_allele, minimum_depth_presence,
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
465 minimum_depth_call, minimum_depth_frequency_dominant_allele):
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
466 alt = None
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
467 low_coverage = False
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
468 multiple_alleles = False
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
469
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
470 if dp >= minimum_depth_presence:
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
471 if dp < minimum_depth_call:
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
472 alt = 'N' * len(variant_position['REF'])
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
473 low_coverage = True
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
474 else:
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
475 if ad_idv < minimum_depth_call:
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
476 alt = 'N' * len(variant_position['REF'])
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
477 low_coverage = True
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
478 if float(ad_idv) / float(dp) < minimum_depth_frequency_dominant_allele:
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
479 multiple_alleles = True
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
480 else:
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
481 if float(ad_idv) / float(dp) < minimum_depth_frequency_dominant_allele:
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
482 alt = 'N' * len(variant_position['REF'])
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
483 if index_dominant_allele is not None:
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
484 variants_coverage = [int(variant_position['format']['AD'][i]) for i in
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
485 range(0, len(variant_position['ALT']) + 1) if i != index_dominant_allele]
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
486 if sum(variants_coverage) > 0:
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
487 if float(max(variants_coverage)) / float(sum(variants_coverage)) > 0.5:
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
488 multiple_alleles = True
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
489 elif float(max(variants_coverage)) / float(sum(variants_coverage)) == 0.5 and \
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
490 len(variants_coverage) > 2:
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
491 multiple_alleles = True
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
492 else:
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
493 multiple_alleles = True
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
494 else:
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
495 alt = alt_no_matter
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
496 else:
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
497 low_coverage = True
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
498
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
499 return alt, low_coverage, multiple_alleles
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
500
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
501
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
502 def get_alt_alignment(ref, alt):
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
503 if alt is None:
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
504 alt = 'N' * len(ref)
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
505 else:
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
506 if len(ref) != len(alt):
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
507 if len(alt) < len(ref):
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
508 if alt == '.':
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
509 alt = ref
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
510 alt += 'N' * (len(ref) - len(alt))
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
511 else:
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
512 if alt[:len(ref)] == ref:
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
513 alt = '.'
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
514 else:
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
515 alt = alt[:len(ref)]
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
516
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
517 return alt
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
518
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
519
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
520 def get_indel_more_likely(variant_position, indels_entry):
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
521 indel_coverage = {}
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
522 for i in indels_entry:
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
523 indel_coverage[i] = int(variant_position['info']['IDV'])
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
524 return indel_coverage.index(str(max(indel_coverage.values())))
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
525
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
526
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
527 def determine_variant(variant_position, minimum_depth_presence, minimum_depth_call,
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
528 minimum_depth_frequency_dominant_allele, indel_true):
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
529 alt_no_matter, dp, ad_idv, index_dominant_allele = get_alt_no_matter(variant_position, indel_true)
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
530
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
531 alt_correct, low_coverage, multiple_alleles = get_alt_correct(variant_position, alt_no_matter, dp, ad_idv,
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
532 index_dominant_allele, minimum_depth_presence,
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
533 minimum_depth_call,
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
534 minimum_depth_frequency_dominant_allele)
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
535
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
536 alt_alignment = get_alt_alignment(variant_position['REF'], alt_correct)
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
537
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
538 return variant_position['REF'], alt_correct, low_coverage, multiple_alleles, alt_no_matter, alt_alignment
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
539
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
540
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
541 def confirm_nucleotides_indel(ref, alt, variants, position_start_indel, minimum_depth_presence, minimum_depth_call,
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
542 minimum_depth_frequency_dominant_allele, alignment_true):
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
543 alt = list(alt)
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
544
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
545 for i in range(0, len(alt) - 1):
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
546 if len(alt) < len(ref):
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
547 new_position = position_start_indel + len(alt) - i - 1
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
548 alt_position = len(alt) - i - 1
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
549 else:
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
550 if i + 1 > len(ref):
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
551 break
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
552 new_position = position_start_indel + 1 + i
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
553 alt_position = 1 + i
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
554
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
555 if alt[alt_position] != 'N':
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
556 if new_position not in variants:
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
557 if alignment_true:
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
558 alt[alt_position] = 'N'
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
559 else:
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
560 alt = alt[: alt_position]
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
561 break
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
562
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
563 entry_with_indel, entry_with_snp = indel_entry(variants[new_position])
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
564 new_ref, alt_correct, low_coverage, multiple_alleles, alt_no_matter, alt_alignment = \
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
565 determine_variant(variants[new_position][entry_with_snp], minimum_depth_presence, minimum_depth_call,
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
566 minimum_depth_frequency_dominant_allele, False)
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
567 if alt_no_matter != '.' and alt[alt_position] != alt_no_matter:
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
568 alt[alt_position] = alt_no_matter
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
569
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
570 return ''.join(alt)
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
571
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
572
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
573 def snp_indel(variants, position, minimum_depth_presence, minimum_depth_call, minimum_depth_frequency_dominant_allele):
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
574 entry_with_indel, entry_with_snp = indel_entry(variants[position])
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
575
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
576 if len(entry_with_indel) == 0:
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
577 ref, alt_correct, low_coverage, multiple_alleles, alt_no_matter, alt_alignment = \
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
578 determine_variant(variants[position][entry_with_snp], minimum_depth_presence, minimum_depth_call,
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
579 minimum_depth_frequency_dominant_allele, False)
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
580 else:
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
581 ref_snp, alt_correct_snp, low_coverage_snp, multiple_alleles_snp, alt_no_matter_snp, alt_alignment_snp = \
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
582 determine_variant(variants[position][entry_with_snp], minimum_depth_presence, minimum_depth_call,
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
583 minimum_depth_frequency_dominant_allele, False)
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
584
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
585 indel_more_likely = entry_with_indel[0]
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
586 if len(entry_with_indel) > 1:
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
587 indel_more_likely = get_indel_more_likely(variants[position], entry_with_indel)
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
588
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
589 ref, alt_correct, low_coverage, multiple_alleles, alt_no_matter, alt_alignment = \
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
590 determine_variant(variants[position][indel_more_likely], minimum_depth_presence, minimum_depth_call,
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
591 minimum_depth_frequency_dominant_allele, True)
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
592
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
593 if alt_no_matter == '.':
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
594 ref, alt_correct, low_coverage, multiple_alleles, alt_no_matter, alt_alignment = \
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
595 ref_snp, alt_correct_snp, low_coverage_snp, multiple_alleles_snp, alt_no_matter_snp, alt_alignment_snp
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
596 else:
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
597 if alt_correct is None and alt_correct_snp is not None:
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
598 alt_correct = alt_correct_snp
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
599 elif alt_correct is not None and alt_correct_snp is not None:
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
600 if alt_correct_snp != '.' and alt_correct[0] != alt_correct_snp:
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
601 alt_correct = alt_correct_snp + alt_correct[1:] if len(alt_correct) > 1 else alt_correct_snp
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
602 if alt_no_matter_snp != '.' and alt_no_matter[0] != alt_no_matter_snp:
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
603 alt_no_matter = alt_no_matter_snp + alt_no_matter[1:] if len(alt_no_matter) > 1 else alt_no_matter_snp
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
604 if alt_alignment_snp != '.' and alt_alignment[0] != alt_alignment_snp:
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
605 alt_alignment = alt_alignment_snp + alt_alignment[1:] if len(alt_alignment) > 1 else alt_alignment_snp
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
606
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
607 # if alt_no_matter != '.':
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
608 # alt_no_matter = confirm_nucleotides_indel(ref, alt_no_matter, variants, position, minimum_depth_presence, minimum_depth_call, minimum_depth_frequency_dominant_allele, False)
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
609 # if alt_correct is not None and alt_correct != '.':
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
610 # alt_correct = confirm_nucleotides_indel(ref, alt_correct, variants, position, minimum_depth_presence, minimum_depth_call, minimum_depth_frequency_dominant_allele, False)
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
611 # if alt_alignment != '.':
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
612 # alt_alignment = confirm_nucleotides_indel(ref, alt_alignment, variants, position, minimum_depth_presence, minimum_depth_call, minimum_depth_frequency_dominant_allele, True)
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
613
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
614 return ref, alt_correct, low_coverage, multiple_alleles, alt_no_matter, alt_alignment
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
615
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
616
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
617 def get_true_variants(variants, minimum_depth_presence, minimum_depth_call, minimum_depth_frequency_dominant_allele,
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
618 sequence):
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
619 variants_correct = {}
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
620 variants_no_matter = {}
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
621 variants_alignment = {}
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
622
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
623 correct_absent_positions = {}
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
624 correct_last_absent_position = ''
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
625
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
626 no_matter_absent_positions = {}
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
627 no_matter_last_absent_position = ''
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
628
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
629 multiple_alleles_found = []
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
630
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
631 counter = 1
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
632 while counter <= len(sequence):
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
633 if counter in variants:
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
634 no_matter_last_absent_position = ''
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
635
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
636 ref, alt_correct, low_coverage, multiple_alleles, alt_no_matter, alt_alignment = \
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
637 snp_indel(variants, counter, minimum_depth_presence, minimum_depth_call,
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
638 minimum_depth_frequency_dominant_allele)
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
639
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
640 if alt_alignment != '.':
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
641 variants_alignment[counter] = {'REF': ref, 'ALT': alt_alignment}
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
642
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
643 if alt_no_matter != '.':
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
644 variants_no_matter[counter] = {'REF': ref, 'ALT': alt_no_matter}
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
645
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
646 if alt_correct is None:
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
647 if counter - len(correct_last_absent_position) in correct_absent_positions:
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
648 correct_absent_positions[counter - len(correct_last_absent_position)]['REF'] += ref
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
649 else:
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
650 correct_absent_positions[counter] = {'REF': ref, 'ALT': ''}
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
651 correct_last_absent_position += ref
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
652 else:
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
653 if alt_correct != '.':
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
654 if len(alt_correct) < len(ref):
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
655 if len(alt_correct) == 1:
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
656 correct_absent_positions[counter + 1] = {'REF': ref[1:], 'ALT': ''}
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
657 else:
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
658 correct_absent_positions[counter + 1] = {'REF': ref[1:], 'ALT': alt_correct[1:]}
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
659
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
660 correct_last_absent_position = ref[1:]
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
661 else:
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
662 variants_correct[counter] = {'REF': ref, 'ALT': alt_correct}
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
663 correct_last_absent_position = ''
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
664 else:
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
665 correct_last_absent_position = ''
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
666
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
667 if multiple_alleles:
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
668 multiple_alleles_found.append(counter)
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
669
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
670 counter += len(ref)
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
671 else:
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
672 variants_alignment[counter] = {'REF': sequence[counter - 1], 'ALT': 'N'}
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
673
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
674 if counter - len(correct_last_absent_position) in correct_absent_positions:
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
675 correct_absent_positions[counter - len(correct_last_absent_position)]['REF'] += sequence[counter - 1]
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
676 else:
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
677 correct_absent_positions[counter] = {'REF': sequence[counter - 1], 'ALT': ''}
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
678 correct_last_absent_position += sequence[counter - 1]
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
679
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
680 if counter - len(no_matter_last_absent_position) in no_matter_absent_positions:
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
681 no_matter_absent_positions[counter - len(no_matter_last_absent_position)]['REF'] += \
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
682 sequence[counter - 1]
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
683 else:
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
684 no_matter_absent_positions[counter] = {'REF': sequence[counter - 1], 'ALT': ''}
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
685 no_matter_last_absent_position += sequence[counter - 1]
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
686
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
687 counter += 1
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
688
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
689 for position in correct_absent_positions:
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
690 if position == 1:
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
691 variants_correct[position] = {'REF': correct_absent_positions[position]['REF'], 'ALT': 'N'}
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
692 else:
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
693 if position - 1 not in variants_correct:
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
694 variants_correct[position - 1] = \
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
695 {'REF': sequence[position - 2] + correct_absent_positions[position]['REF'],
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
696 'ALT': sequence[position - 2] + correct_absent_positions[position]['ALT']}
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
697 else:
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
698 variants_correct[position - 1] = \
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
699 {'REF': variants_correct[position - 1]['REF'] +
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
700 correct_absent_positions[position]['REF'][len(variants_correct[position - 1]['REF']) - 1:],
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
701 'ALT': variants_correct[position - 1]['ALT'] +
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
702 correct_absent_positions[position]['ALT'][len(variants_correct[position - 1]['ALT']) - 1 if
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
703 len(variants_correct[position - 1]['ALT']) > 0
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
704 else 0:]}
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
705
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
706 for position in no_matter_absent_positions:
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
707 if position == 1:
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
708 variants_no_matter[position] = {'REF': no_matter_absent_positions[position]['REF'], 'ALT': 'N'}
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
709 else:
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
710 if position - 1 not in variants_no_matter:
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
711 variants_no_matter[position - 1] = \
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
712 {'REF': sequence[position - 2] + no_matter_absent_positions[position]['REF'],
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
713 'ALT': sequence[position - 2] + no_matter_absent_positions[position]['ALT']}
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
714 else:
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
715 variants_no_matter[position - 1] = \
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
716 {'REF': variants_no_matter[position - 1]['REF'] +
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
717 no_matter_absent_positions[position]['REF'][len(variants_no_matter[position - 1]['REF']) -
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
718 1:],
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
719 'ALT': variants_no_matter[position - 1]['ALT'] +
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
720 no_matter_absent_positions[position]['ALT'][len(variants_no_matter[position - 1]['ALT']) -
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
721 1 if
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
722 len(variants_no_matter[position - 1]['ALT']) > 0
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
723 else 0:]}
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
724
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
725 return variants_correct, variants_no_matter, variants_alignment, multiple_alleles_found
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
726
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
727
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
728 def clean_variant_in_extra_seq_left(variant_dict, position, length_extra_seq, multiple_alleles_found,
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
729 number_multi_alleles):
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
730 number_diferences = 0
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
731
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
732 if position + len(variant_dict[position]['REF']) - 1 > length_extra_seq:
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
733 if multiple_alleles_found is not None and position in multiple_alleles_found:
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
734 number_multi_alleles += 1
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
735
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
736 temp_variant = variant_dict[position]
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
737 del variant_dict[position]
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
738 variant_dict[length_extra_seq] = {}
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
739 variant_dict[length_extra_seq]['REF'] = temp_variant['REF'][length_extra_seq - position:]
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
740 variant_dict[length_extra_seq]['ALT'] = temp_variant['ALT'][length_extra_seq - position:] if \
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
741 len(temp_variant['ALT']) > length_extra_seq - position else \
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
742 temp_variant['REF'][length_extra_seq - position]
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
743 number_diferences = count_number_diferences(variant_dict[length_extra_seq]['REF'],
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
744 variant_dict[length_extra_seq]['ALT'])
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
745 else:
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
746 del variant_dict[position]
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
747
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
748 return variant_dict, number_multi_alleles, number_diferences
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
749
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
750
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
751 def clean_variant_in_extra_seq_rigth(variant_dict, position, sequence_length, length_extra_seq):
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
752 if position + len(variant_dict[position]['REF']) - 1 > sequence_length - length_extra_seq:
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
753 variant_dict[position]['REF'] = \
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
754 variant_dict[position]['REF'][: - (position - (sequence_length - length_extra_seq)) + 1]
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
755 variant_dict[position]['ALT'] = \
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
756 variant_dict[position]['ALT'][: - (position - (sequence_length - length_extra_seq)) + 1] if \
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
757 len(variant_dict[position]['ALT']) >= - (position - (sequence_length - length_extra_seq)) + 1 else \
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
758 variant_dict[position]['ALT']
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
759
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
760 number_diferences = count_number_diferences(variant_dict[position]['REF'], variant_dict[position]['ALT'])
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
761
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
762 return variant_dict, number_diferences
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
763
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
764
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
765 def cleanning_variants_extra_seq(variants_correct, variants_no_matter, variants_alignment, multiple_alleles_found,
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
766 length_extra_seq, sequence_length):
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
767 number_multi_alleles = 0
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
768 number_diferences = 0
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
769
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
770 counter = 1
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
771 while counter <= sequence_length:
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
772 if counter <= length_extra_seq:
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
773 if counter in variants_correct:
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
774 variants_correct, number_multi_alleles, number_diferences = \
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
775 clean_variant_in_extra_seq_left(variants_correct, counter, length_extra_seq, multiple_alleles_found,
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
776 number_multi_alleles)
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
777 if counter in variants_no_matter:
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
778 variants_no_matter, ignore, ignore = \
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
779 clean_variant_in_extra_seq_left(variants_no_matter, counter, length_extra_seq, None, None)
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
780 if counter in variants_alignment:
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
781 variants_alignment, ignore, ignore = \
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
782 clean_variant_in_extra_seq_left(variants_alignment, counter, length_extra_seq, None, None)
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
783 elif sequence_length - length_extra_seq >= counter > length_extra_seq:
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
784 if counter in variants_correct:
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
785 if counter in multiple_alleles_found:
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
786 number_multi_alleles += 1
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
787 variants_correct, number_diferences_found = \
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
788 clean_variant_in_extra_seq_rigth(variants_correct, counter, sequence_length, length_extra_seq)
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
789 number_diferences += number_diferences_found
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
790 if counter in variants_no_matter:
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
791 variants_no_matter, ignore = \
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
792 clean_variant_in_extra_seq_rigth(variants_no_matter, counter, sequence_length, length_extra_seq)
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
793 if counter in variants_alignment:
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
794 variants_alignment, ignore = \
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
795 clean_variant_in_extra_seq_rigth(variants_alignment, counter, sequence_length, length_extra_seq)
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
796 else:
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
797 if counter in variants_correct:
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
798 del variants_correct[counter]
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
799 if counter in variants_no_matter:
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
800 del variants_no_matter[counter]
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
801 if counter in variants_alignment:
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
802 del variants_alignment[counter]
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
803
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
804 counter += 1
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
805
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
806 return variants_correct, variants_no_matter, variants_alignment, number_multi_alleles, number_diferences
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
807
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
808
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
809 def get_coverage(gene_coverage):
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
810 coverage = {}
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
811
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
812 with open(gene_coverage, 'rtU') as reader:
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
813 for line in reader:
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
814 line = line.rstrip('\r\n')
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
815 if len(line) > 0:
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
816 line = line.split('\t')
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
817 coverage[int(line[1])] = int(line[2])
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
818
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
819 return coverage
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
820
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
821
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
822 def get_coverage_report(coverage, sequence_length, minimum_depth_presence, minimum_depth_call, length_extra_seq):
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
823 if len(coverage) == 0:
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
824 return sequence_length - 2 * length_extra_seq, 100.0, 0.0
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
825
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
826 count_absent = 0
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
827 count_low_coverage = 0
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
828 sum_coverage = 0
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
829
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
830 counter = 1
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
831 while counter <= sequence_length:
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
832 if sequence_length - length_extra_seq >= counter > length_extra_seq:
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
833 if coverage[counter] < minimum_depth_presence:
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
834 count_absent += 1
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
835 else:
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
836 if coverage[counter] < minimum_depth_call:
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
837 count_low_coverage += 1
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
838 sum_coverage += coverage[counter]
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
839 counter += 1
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
840
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
841 mean_coverage = 0
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
842 percentage_low_coverage = 0
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
843 if sequence_length - 2 * length_extra_seq - count_absent > 0:
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
844 mean_coverage = float(sum_coverage) / float(sequence_length - 2 * length_extra_seq - count_absent)
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
845 percentage_low_coverage = \
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
846 float(count_low_coverage) / float(sequence_length - 2 * length_extra_seq - count_absent) * 100
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
847
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
848 return count_absent, percentage_low_coverage, mean_coverage
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
849
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
850
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
851 # Get genome coverage data
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
852 def compute_genome_coverage_data(alignment_file, sequence_to_analyse, outdir, counter):
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
853 genome_coverage_data_file = os.path.join(outdir, 'samtools_depth.sequence_' + str(counter) + '.tab')
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
854 command = ['samtools', 'depth', '-a', '-q', '0', '-r', sequence_to_analyse, alignment_file, '>',
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
855 genome_coverage_data_file]
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
856 run_successfully, stdout, stderr = utils.run_command_popen_communicate(command, True, None, False)
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
857 return run_successfully, genome_coverage_data_file
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
858
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
859
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
860 def write_variants_vcf(variants, outdir, sequence_to_analyse, sufix):
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
861 vcf_file = os.path.join(outdir, str(sequence_to_analyse + '.' + sufix + '.vcf'))
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
862 with open(vcf_file, 'wt') as writer:
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
863 writer.write('##fileformat=VCFv4.2' + '\n')
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
864 writer.write('#' + '\t'.join(['SEQUENCE', 'POSITION', 'ID_unused', 'REFERENCE_sequence', 'ALTERNATIVE_sequence',
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
865 'QUALITY_unused', 'FILTER_unused', 'INFO_unused', 'FORMAT_unused']) + '\n')
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
866 for i in sorted(variants.keys()):
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
867 writer.write('\t'.join([sequence_to_analyse, str(i), '.', variants[i]['REF'], variants[i]['ALT'], '.', '.',
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
868 '.', '.']) + '\n')
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
869
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
870 compressed_vcf_file = vcf_file + '.gz'
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
871 command = ['bcftools', 'convert', '-o', compressed_vcf_file, '-O', 'z', vcf_file]
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
872 run_successfully, stdout, stderr = utils.run_command_popen_communicate(command, False, None, False)
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
873 if run_successfully:
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
874 command = ['bcftools', 'index', compressed_vcf_file]
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
875 run_successfully, stdout, stderr = utils.run_command_popen_communicate(command, False, None, False)
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
876
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
877 if not run_successfully:
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
878 compressed_vcf_file = None
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
879
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
880 return run_successfully, compressed_vcf_file
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
881
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
882
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
883 def parse_fasta_in_memory(fasta_memory):
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
884 fasta_memory = fasta_memory.splitlines()
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
885 sequence_dict = {}
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
886 for line in fasta_memory:
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
887 if len(line) > 0:
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
888 if line.startswith('>'):
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
889 sequence_dict = {'header': line[1:], 'sequence': ''}
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
890 else:
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
891 sequence_dict['sequence'] += line
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
892
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
893 return sequence_dict
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
894
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
895
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
896 def compute_consensus_sequence(reference_file, sequence_to_analyse, compressed_vcf_file, outdir):
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
897 sequence_dict = None
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
898
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
899 gene_fasta = os.path.join(outdir, str(sequence_to_analyse + '.fasta'))
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
900
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
901 run_successfully, stdout = index_fasta_samtools(reference_file, sequence_to_analyse, gene_fasta, False)
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
902 if run_successfully:
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
903 command = ['bcftools', 'norm', '-c', 's', '-f', gene_fasta, '-Ov', compressed_vcf_file]
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
904 run_successfully, stdout, stderr = utils.run_command_popen_communicate(command, False, None, False)
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
905 if run_successfully:
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
906 command = ['bcftools', 'consensus', '-f', gene_fasta, compressed_vcf_file, '-H', '1']
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
907 run_successfully, stdout, stderr = utils.run_command_popen_communicate(command, False, None, False)
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
908 if run_successfully:
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
909 sequence_dict = parse_fasta_in_memory(stdout)
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
910
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
911 return run_successfully, sequence_dict
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
912
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
913
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
914 def create_sample_consensus_sequence(outdir, sequence_to_analyse, reference_file, variants, minimum_depth_presence,
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
915 minimum_depth_call, minimum_depth_frequency_dominant_allele, sequence,
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
916 length_extra_seq):
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
917 variants_correct, variants_noMatter, variants_alignment, multiple_alleles_found = \
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
918 get_true_variants(variants, minimum_depth_presence, minimum_depth_call, minimum_depth_frequency_dominant_allele,
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
919 sequence)
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
920
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
921 variants_correct, variants_noMatter, variants_alignment, number_multi_alleles, number_diferences = \
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
922 cleanning_variants_extra_seq(variants_correct, variants_noMatter, variants_alignment, multiple_alleles_found,
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
923 length_extra_seq, len(sequence))
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
924
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
925 run_successfully = False
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
926 consensus = {'correct': {}, 'noMatter': {}, 'alignment': {}}
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
927 for variant_type in ['variants_correct', 'variants_noMatter', 'variants_alignment']:
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
928 run_successfully, compressed_vcf_file = \
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
929 write_variants_vcf(eval(variant_type), outdir, sequence_to_analyse, variant_type.split('_', 1)[1])
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
930 if run_successfully:
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
931 run_successfully, sequence_dict = \
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
932 compute_consensus_sequence(reference_file, sequence_to_analyse, compressed_vcf_file, outdir)
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
933 if run_successfully:
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
934 consensus[variant_type.split('_', 1)[1]] = \
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
935 {'header': sequence_dict['header'],
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
936 'sequence': sequence_dict['sequence'][length_extra_seq:len(sequence_dict['sequence']) -
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
937 length_extra_seq]}
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
938
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
939 return run_successfully, number_multi_alleles, consensus, number_diferences
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
940
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
941
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
942 @utils.trace_unhandled_exceptions
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
943 def analyse_sequence_data(bam_file, sequence_information, outdir, counter, reference_file, length_extra_seq,
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
944 minimum_depth_presence, minimum_depth_call, minimum_depth_frequency_dominant_allele):
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
945 count_absent = None
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
946 percentage_low_coverage = None
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
947 mean_coverage = None
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
948 number_diferences = 0
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
949 number_multi_alleles = 0
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
950 consensus_sequence = {'correct': {}, 'noMatter': {}, 'alignment': {}}
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
951
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
952 # Create vcf file (for multiple alleles check)
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
953 run_successfully, gene_vcf = create_vcf(bam_file, sequence_information['header'], outdir, counter, reference_file)
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
954 if run_successfully:
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
955 # Create coverage tab file
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
956 run_successfully, gene_coverage = \
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
957 compute_genome_coverage_data(bam_file, sequence_information['header'], outdir, counter)
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
958
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
959 if run_successfully:
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
960 try:
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
961 variants = get_variants(gene_vcf=gene_vcf, seq_name=sequence_information['header'],
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
962 encoding=sys.getdefaultencoding())
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
963 except UnicodeDecodeError:
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
964 try:
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
965 print('It was found an enconding error while parsing the following VCF, but lets try forcing it to'
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
966 ' "utf_8" encoding: {}'.format(gene_vcf))
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
967 variants = get_variants(gene_vcf=gene_vcf, seq_name=sequence_information['header'],
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
968 encoding='utf_8')
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
969 except UnicodeDecodeError:
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
970 print('It was found an enconding error while parsing the following VCF, but lets try forcing it to'
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
971 ' "latin_1" encoding: {}'.format(gene_vcf))
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
972 variants = get_variants(gene_vcf=gene_vcf, seq_name=sequence_information['header'],
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
973 encoding='latin_1')
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
974
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
975 coverage = get_coverage(gene_coverage)
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
976
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
977 run_successfully, number_multi_alleles, consensus_sequence, number_diferences = \
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
978 create_sample_consensus_sequence(outdir, sequence_information['header'], reference_file, variants,
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
979 minimum_depth_presence, minimum_depth_call,
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
980 minimum_depth_frequency_dominant_allele,
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
981 sequence_information['sequence'], length_extra_seq)
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
982
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
983 try:
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
984 count_absent, percentage_low_coverage, mean_coverage = \
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
985 get_coverage_report(coverage, sequence_information['length'], minimum_depth_presence,
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
986 minimum_depth_call, length_extra_seq)
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
987 except KeyError:
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
988 print('ERROR: KeyError')
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
989 print(sequence_information)
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
990 raise KeyError
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
991
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
992 utils.save_variable_to_pickle([run_successfully, counter, number_multi_alleles, count_absent,
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
993 percentage_low_coverage, mean_coverage, consensus_sequence, number_diferences],
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
994 outdir, str('coverage_info.' + str(counter)))
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
995
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
996
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
997 def clean_header(header):
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
998 problematic_characters = ["|", " ", ",", ".", "(", ")", "'", "/", ":"]
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
999 new_header = str(header)
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
1000 if any(x in header for x in problematic_characters):
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
1001 for x in problematic_characters:
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
1002 new_header = new_header.replace(x, '_')
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
1003 return header, new_header
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
1004
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
1005
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
1006 def get_sequence_information(fasta_file, length_extra_seq):
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
1007 sequence_dict = {}
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
1008 headers = {}
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
1009 headers_changed = False
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
1010
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
1011 with open(fasta_file, 'rtU') as reader:
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
1012 blank_line_found = False
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
1013 sequence_counter = 0
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
1014 temp_sequence_dict = {}
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
1015 for line in reader:
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
1016 line = line.rstrip('\r\n')
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
1017 if len(line) > 0:
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
1018 if not blank_line_found:
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
1019 if line.startswith('>'):
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
1020 if len(temp_sequence_dict) > 0:
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
1021 if list(temp_sequence_dict.values())[0]['length'] - 2 * length_extra_seq > 0:
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
1022 sequence_dict[list(temp_sequence_dict.keys())[0]] = list(temp_sequence_dict.values())[0]
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
1023 else:
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
1024 print('{header} sequence ignored due to'
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
1025 ' length = 0'.format(header=list(temp_sequence_dict.values())[0]['header']))
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
1026 del headers[list(temp_sequence_dict.values())[0]['header']]
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
1027 temp_sequence_dict = {}
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
1028
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
1029 original_header, new_header = clean_header(line[1:])
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
1030 if new_header in headers:
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
1031 sys.exit('Found duplicated sequence'
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
1032 ' headers: {original_header}'.format(original_header=original_header))
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
1033
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
1034 sequence_counter += 1
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
1035 temp_sequence_dict[sequence_counter] = {'header': new_header, 'sequence': '', 'length': 0}
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
1036 headers[new_header] = str(original_header)
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
1037 if new_header != original_header:
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
1038 headers_changed = True
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
1039 else:
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
1040 temp_sequence_dict[sequence_counter]['sequence'] += line.replace(' ', '')
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
1041 temp_sequence_dict[sequence_counter]['length'] += len(line.replace(' ', ''))
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
1042 else:
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
1043 sys.exit('It was found a blank line between the fasta file above line ' + line)
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
1044 else:
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
1045 blank_line_found = True
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
1046
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
1047 if len(temp_sequence_dict) > 0:
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
1048 if list(temp_sequence_dict.values())[0]['length'] - 2 * length_extra_seq > 0:
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
1049 sequence_dict[list(temp_sequence_dict.keys())[0]] = list(temp_sequence_dict.values())[0]
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
1050 else:
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
1051 print('{header} sequence ignored due to'
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
1052 ' length <= 0'.format(header=list(temp_sequence_dict.values())[0]['header']))
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
1053 del headers[list(temp_sequence_dict.values())[0]['header']]
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
1054
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
1055 return sequence_dict, headers, headers_changed
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
1056
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
1057
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
1058 def sequence_data(sample, reference_file, bam_file, outdir, threads, length_extra_seq, minimum_depth_presence,
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
1059 minimum_depth_call, minimum_depth_frequency_dominant_allele, debug_mode_true, not_write_consensus):
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
1060 sequence_data_outdir = os.path.join(outdir, 'sequence_data', '')
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
1061 utils.remove_directory(sequence_data_outdir)
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
1062 os.mkdir(sequence_data_outdir)
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
1063
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
1064 sequences, headers, headers_changed = get_sequence_information(reference_file, length_extra_seq)
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
1065
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
1066 pool = multiprocessing.Pool(processes=threads)
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
1067 for sequence_counter in sequences:
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
1068 sequence_dir = os.path.join(sequence_data_outdir, str(sequence_counter), '')
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
1069 utils.remove_directory(sequence_dir)
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
1070 os.makedirs(sequence_dir)
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
1071 pool.apply_async(analyse_sequence_data, args=(bam_file, sequences[sequence_counter], sequence_dir,
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
1072 sequence_counter, reference_file, length_extra_seq,
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
1073 minimum_depth_presence, minimum_depth_call,
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
1074 minimum_depth_frequency_dominant_allele,))
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
1075 pool.close()
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
1076 pool.join()
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
1077
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
1078 run_successfully, sample_data, consensus_files, consensus_sequences = \
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
1079 gather_data_together(sample, sequence_data_outdir, sequences, outdir.rsplit('/', 2)[0], debug_mode_true,
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
1080 length_extra_seq, not_write_consensus)
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
1081
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
1082 return run_successfully, sample_data, consensus_files, consensus_sequences
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
1083
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
1084
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
1085 def chunkstring(string, length):
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
1086 return (string[0 + i:length + i] for i in range(0, len(string), length))
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
1087
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
1088
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
1089 def write_consensus(outdir, sample, consensus_sequence):
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
1090 consensus_files = {}
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
1091 for consensus_type in ['correct', 'noMatter', 'alignment']:
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
1092 consensus_files[consensus_type] = os.path.join(outdir, str(sample + '.' + consensus_type + '.fasta'))
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
1093 with open(consensus_files[consensus_type], 'at') as writer:
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
1094 writer.write('>' + consensus_sequence[consensus_type]['header'] + '\n')
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
1095 fasta_sequence_lines = chunkstring(consensus_sequence[consensus_type]['sequence'], 80)
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
1096 for line in fasta_sequence_lines:
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
1097 writer.write(line + '\n')
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
1098 return consensus_files
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
1099
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
1100
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
1101 def gather_data_together(sample, data_directory, sequences_information, outdir, debug_mode_true, length_extra_seq,
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
1102 not_write_consensus):
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
1103 run_successfully = True
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
1104 counter = 0
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
1105 sample_data = {}
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
1106
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
1107 consensus_files = None
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
1108 consensus_sequences_together = {'correct': {}, 'noMatter': {}, 'alignment': {}}
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
1109
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
1110 write_consensus_first_time = True
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
1111
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
1112 genes_directories = [d for d in os.listdir(data_directory) if
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
1113 not d.startswith('.') and
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
1114 os.path.isdir(os.path.join(data_directory, d, ''))]
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
1115 for gene_dir in genes_directories:
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
1116 gene_dir_path = os.path.join(data_directory, gene_dir, '')
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
1117
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
1118 files = [f for f in os.listdir(gene_dir_path) if
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
1119 not f.startswith('.') and
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
1120 os.path.isfile(os.path.join(gene_dir_path, f))]
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
1121 for file_found in files:
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
1122 if file_found.startswith('coverage_info.') and file_found.endswith('.pkl'):
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
1123 file_path = os.path.join(gene_dir_path, file_found)
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
1124
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
1125 if run_successfully:
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
1126 run_successfully, sequence_counter, multiple_alleles_found, count_absent, percentage_low_coverage, \
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
1127 mean_coverage, consensus_sequence, \
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
1128 number_diferences = utils.extract_variable_from_pickle(file_path)
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
1129
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
1130 if not not_write_consensus:
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
1131 for consensus_type in consensus_sequence:
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
1132 consensus_sequences_together[consensus_type][sequence_counter] = \
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
1133 {'header': consensus_sequence[consensus_type]['header'],
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
1134 'sequence': consensus_sequence[consensus_type]['sequence']}
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
1135
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
1136 if write_consensus_first_time:
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
1137 for consensus_type in ['correct', 'noMatter', 'alignment']:
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
1138 file_to_remove = os.path.join(outdir, str(sample + '.' + consensus_type + '.fasta'))
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
1139 if os.path.isfile(file_to_remove):
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
1140 os.remove(file_to_remove)
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
1141 write_consensus_first_time = False
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
1142 consensus_files = write_consensus(outdir, sample, consensus_sequence)
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
1143
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
1144 gene_identity = 0
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
1145 if sequences_information[sequence_counter]['length'] - 2 * length_extra_seq - count_absent > 0:
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
1146 gene_identity = 100 - \
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
1147 (float(number_diferences) /
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
1148 (sequences_information[sequence_counter]['length'] - 2 * length_extra_seq -
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
1149 count_absent)) * 100
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
1150
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
1151 sample_data[sequence_counter] = \
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
1152 {'header': sequences_information[sequence_counter]['header'],
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
1153 'gene_coverage': 100 - (float(count_absent) /
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
1154 (sequences_information[sequence_counter]['length'] - 2 *
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
1155 length_extra_seq)) * 100,
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
1156 'gene_low_coverage': percentage_low_coverage,
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
1157 'gene_number_positions_multiple_alleles': multiple_alleles_found,
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
1158 'gene_mean_read_coverage': mean_coverage,
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
1159 'gene_identity': gene_identity}
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
1160 counter += 1
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
1161
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
1162 if not debug_mode_true:
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
1163 utils.remove_directory(gene_dir_path)
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
1164
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
1165 if counter != len(sequences_information):
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
1166 run_successfully = False
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
1167
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
1168 return run_successfully, sample_data, consensus_files, consensus_sequences_together
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
1169
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
1170
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
1171 rematch_timer = functools.partial(utils.timer, name='ReMatCh module')
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
1172
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
1173
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
1174 @rematch_timer
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
1175 def run_rematch_module(sample, fastq_files, reference_file, threads, outdir, length_extra_seq, minimum_depth_presence,
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
1176 minimum_depth_call, minimum_depth_frequency_dominant_allele, minimum_gene_coverage,
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
1177 debug_mode_true, num_map_loc, minimum_gene_identity, rematch_run,
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
1178 soft_clip_base_quality, soft_clip_recode_run, reference_dict, soft_clip_cigar_flag_recode,
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
1179 bowtie_algorithm, bowtie_opt, gene_list_reference, not_write_consensus, clean_run=True):
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
1180 rematch_folder = os.path.join(outdir, 'rematch_module', '')
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
1181
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
1182 utils.remove_directory(rematch_folder)
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
1183 os.mkdir(rematch_folder)
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
1184
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
1185 # Map reads
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
1186 run_successfully, bam_file, reference_file = mapping_reads(fastq_files=fastq_files, reference_file=reference_file,
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
1187 threads=threads, outdir=rematch_folder,
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
1188 num_map_loc=num_map_loc, rematch_run=rematch_run,
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
1189 soft_clip_base_quality=soft_clip_base_quality,
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
1190 soft_clip_recode_run=soft_clip_recode_run,
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
1191 reference_dict=reference_dict,
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
1192 soft_clip_cigar_flag_recode=soft_clip_cigar_flag_recode,
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
1193 bowtie_algorithm=bowtie_algorithm, bowtie_opt=bowtie_opt,
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
1194 clean_run=clean_run)
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
1195 if run_successfully:
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
1196 # Index reference file
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
1197 run_successfully, stdout = index_fasta_samtools(reference_file, None, None, True)
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
1198 if run_successfully:
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
1199 print('Analysing alignment data')
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
1200 run_successfully, sample_data, consensus_files, consensus_sequences = \
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
1201 sequence_data(sample, reference_file, bam_file, rematch_folder, threads, length_extra_seq,
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
1202 minimum_depth_presence, minimum_depth_call, minimum_depth_frequency_dominant_allele,
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
1203 debug_mode_true, not_write_consensus)
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
1204
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
1205 if run_successfully:
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
1206 print('Writing report file')
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
1207 number_absent_genes = 0
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
1208 number_genes_multiple_alleles = 0
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
1209 mean_sample_coverage = 0
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
1210 with open(os.path.join(outdir, 'rematchModule_report.txt'), 'wt') as writer:
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
1211 print('\n'.join(outdir, 'rematchModule_report.txt'))
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
1212 writer.write('\t'.join(['#gene', 'percentage_gene_coverage', 'gene_mean_read_coverage',
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
1213 'percentage_gene_low_coverage', 'number_positions_multiple_alleles',
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
1214 'percentage_gene_identity']) + '\n')
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
1215 for i in range(1, len(sample_data) + 1):
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
1216 writer.write('\t'.join([gene_list_reference[sample_data[i]['header']],
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
1217 str(round(sample_data[i]['gene_coverage'], 2)),
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
1218 str(round(sample_data[i]['gene_mean_read_coverage'], 2)),
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
1219 str(round(sample_data[i]['gene_low_coverage'], 2)),
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
1220 str(sample_data[i]['gene_number_positions_multiple_alleles']),
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
1221 str(round(sample_data[i]['gene_identity'], 2))]) + '\n')
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
1222
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
1223 if sample_data[i]['gene_coverage'] < minimum_gene_coverage or \
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
1224 sample_data[i]['gene_identity'] < minimum_gene_identity:
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
1225 number_absent_genes += 1
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
1226 else:
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
1227 mean_sample_coverage += sample_data[i]['gene_mean_read_coverage']
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
1228 if sample_data[i]['gene_number_positions_multiple_alleles'] > 0:
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
1229 number_genes_multiple_alleles += 1
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
1230
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
1231 if len(sample_data) - number_absent_genes > 0:
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
1232 mean_sample_coverage = \
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
1233 float(mean_sample_coverage) / float(len(sample_data) - number_absent_genes)
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
1234 else:
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
1235 mean_sample_coverage = 0
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
1236
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
1237 writer.write('\n'.join(['#general',
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
1238 '>number_absent_genes', str(number_absent_genes),
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
1239 '>number_genes_multiple_alleles', str(number_genes_multiple_alleles),
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
1240 '>mean_sample_coverage', str(round(mean_sample_coverage, 2))]) + '\n')
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
1241
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
1242 print('\n'.join([str('number_absent_genes: ' + str(number_absent_genes)),
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
1243 str('number_genes_multiple_alleles: ' + str(number_genes_multiple_alleles)),
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
1244 str('mean_sample_coverage: ' + str(round(mean_sample_coverage, 2)))]))
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
1245
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
1246 if not debug_mode_true:
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
1247 utils.remove_directory(rematch_folder)
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
1248
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
1249 return run_successfully, sample_data if 'sample_data' in locals() else None, \
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
1250 {'number_absent_genes': number_absent_genes if 'number_absent_genes' in locals() else None,
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
1251 'number_genes_multiple_alleles': number_genes_multiple_alleles if
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
1252 'number_genes_multiple_alleles' in locals() else None,
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
1253 'mean_sample_coverage': round(mean_sample_coverage, 2) if 'mean_sample_coverage' in locals() else None}, \
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
1254 consensus_files if 'consensus_files' in locals() else None,\
c6bab5103a14 "planemo upload commit 6abf3e299d82d07e6c3cf8642bdea80e96df64c3-dirty"
iss
parents:
diff changeset
1255 consensus_sequences if 'consensus_sequences' in locals() else None