Mercurial > repos > cstrittmatter > test_eurl_vtec_wgs_pt

comparison scripts/ReMatCh/utils/strip_alignment.py @ 0:965517909457 draft

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
planemo upload commit 15239f1674081ab51ab8dd75a9a40cf1bfaa93e8
author cstrittmatter
date Wed, 22 Jan 2020 08:41:44 -0500
parents
children 0cbed1c0a762
comparison
equal deleted inserted replaced
-1:000000000000 0:965517909457
1 #!/usr/bin/env python
2
3 # -*- coding: utf-8 -*-
4
5 """
6 strip_alignment.py - Strip alignment positions containing gaps,
7 missing data and invariable positions
8 <https://github.com/B-UMMI/ReMatCh/>
9
10 Copyright (C) 2017 Miguel Machado <mpmachado@medicina.ulisboa.pt>
11
12 Last modified: March 20, 2017
13
14 This program is free software: you can redistribute it and/or modify
15 it under the terms of the GNU General Public License as published by
16 the Free Software Foundation, either version 3 of the License, or
17 (at your option) any later version.
18
19 This program is distributed in the hope that it will be useful,
20 but WITHOUT ANY WARRANTY; without even the implied warranty of
21 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
22 GNU General Public License for more details.
23
24 You should have received a copy of the GNU General Public License
25 along with this program. If not, see <http://www.gnu.org/licenses/>.
26 """
27
28 from Bio import SeqIO
29 import os
30 import argparse
31 import sys
32
33
34 version = '0.1'
35
36
37 def get_sequences(infile):
38 print 'Getting sequences'
39 sequences_SeqIO = list(SeqIO.parse(infile, 'fasta'))
40
41 sequence_length = None
42 sequences_dict = {}
43 all_executed_printed = False
44 for x, sequence in enumerate(sequences_SeqIO):
45 if sequence_length is None:
46 sequence_length = len(sequence.seq)
47 if sequence_length != len(sequence.seq):
48 sys.exit('Sequences with different length!')
49 sequences_dict[sequence.id] = list(sequence.seq)
50
51 if (x + 1) % 10 == 0:
52 print '\n' + str(round((float(x + 1) / len(sequences_SeqIO)) * 100, 2)) + '% of sequences already processed (getting sequences)'
53 if x + 1 == len(sequences_SeqIO):
54 all_executed_printed = True
55 if not all_executed_printed:
56 print '\n' + str(round((float(x + 1) / len(sequences_SeqIO)) * 100, 2)) + '% of sequences already processed (getting sequences)'
57
58 return sequences_dict, sequence_length
59
60
61 def positions_type(sequences_dict, sequence_length, notGAPs, notMissing, notInvariable):
62 print 'Determining positions type'
63 positions_2_keep = []
64 invariable = []
65 missing = []
66 gaps = []
67 gaps_missing = 0
68 all_executed_printed = False
69 for i in range(0, sequence_length):
70 data = []
71 for sample in sequences_dict:
72 data.append(sequences_dict[sample][i])
73 possibilities = set(data)
74 if len(possibilities) == 1:
75 invariable.append(i)
76 if len(possibilities.intersection(set(['N']))) > 0:
77 missing.append(i)
78 if len(possibilities.intersection(set(['-']))) > 0:
79 gaps.append(i)
80 if len(possibilities.intersection(set(['N', '-']))) > 0:
81 gaps_missing += 1
82 if len(possibilities) > 1 and len(possibilities.intersection(set(['N', '-']))) == 0:
83 positions_2_keep.append(i)
84
85 if (i + 1) % 10000 == 0:
86 print '\n' + str(round((float(i + 1) / sequence_length) * 100, 2)) + '% of positions already processed (determining positions type)'
87 if i + 1 == len(sequences_dict):
88 all_executed_printed = True
89 if not all_executed_printed:
90 print '\n' + str(round((float(i + 1) / sequence_length) * 100, 2)) + '% of positions already processed (determining positions type)'
91
92 print 'Positions to keep (no matter): ' + str(len(positions_2_keep))
93 print 'Invariable sites: ' + str(len(invariable))
94 print 'Positions with missing data ("N"): ' + str(len(missing))
95 print 'Positions with GAPs ("-"): ' + str(len(gaps))
96 print 'Positions with GAPs or missing data: ' + str(gaps_missing)
97
98 if notGAPs:
99 positions_2_keep.extend(gaps)
100 if notMissing:
101 positions_2_keep.extend(missing)
102 if notInvariable:
103 positions_2_keep.extend(invariable)
104
105 positions_2_keep = sorted(set(positions_2_keep))
106
107 print 'Positions to keep (final): ' + str(len(positions_2_keep))
108
109 return positions_2_keep
110
111
112 def chunkstring(string, length):
113 return (string[0 + i:length + i] for i in range(0, len(string), length))
114
115
116 def write_fasta(sequences_dict, positions_2_keep, outfile):
117 print 'Writing stripped sequences'
118 all_executed_printed = False
119 with open(outfile, 'wt') as writer:
120 for x, sample in enumerate(sequences_dict):
121 writer.write('>' + sample + '\n')
122 fasta_sequence_lines = chunkstring(''.join([sequences_dict[sample][i] for i in positions_2_keep]), 80)
123 for line in fasta_sequence_lines:
124 writer.write(line + '\n')
125
126 if (x + 1) % 100 == 0:
127 print '\n' + str(round((float(x + 1) / len(sequences_dict)) * 100, 2)) + '% of sequences already processed (writing stripped sequences)'
128 if x + 1 == len(sequences_dict):
129 all_executed_printed = True
130 if not all_executed_printed:
131 print '\n' + str(round((float(x + 1) / len(sequences_dict)) * 100, 2)) + '% of sequences already processed (writing stripped sequences)'
132
133
134 def strip_alignment(args):
135 outdir = os.path.dirname(os.path.abspath(args.outfile))
136 if not os.path.isdir(outdir):
137 os.makedirs(outdir)
138
139 outfile = os.path.abspath(args.outfile)
140
141 infile = os.path.abspath(args.infile.name)
142
143 sequences_dict, sequence_length = get_sequences(infile)
144 positions_2_keep = positions_type(sequences_dict, sequence_length, args.notGAPs, args.notMissing, args.notInvariable)
145 write_fasta(sequences_dict, positions_2_keep, outfile)
146
147
148 def main():
149 parser = argparse.ArgumentParser(prog='strip_alignment.py', description='Strip alignment positions containing gaps, missing data and invariable positions', formatter_class=argparse.ArgumentDefaultsHelpFormatter)
150 parser.add_argument('--version', help='Version information', action='version', version=str('%(prog)s v' + version))
151
152 parser_required = parser.add_argument_group('Required options')
153 parser_required.add_argument('-i', '--infile', type=argparse.FileType('r'), metavar='/path/to/aligned/input/file.fasta', help='Path to the aligned fasta file', required=True)
154 parser_required.add_argument('-o', '--outfile', type=str, metavar='/path/to/stripped/output/file.fasta', help='Stripped output fasta file', required=True, default='alignment_stripped.fasta')
155
156 parser_optional_general = parser.add_argument_group('General facultative options')
157 parser_optional_general.add_argument('--notGAPs', action='store_true', help='Not strip positions with GAPs')
158 parser_optional_general.add_argument('--notMissing', action='store_true', help='Not strip positions with missing data')
159 parser_optional_general.add_argument('--notInvariable', action='store_true', help='Not strip invariable sites')
160
161 args = parser.parse_args()
162
163 strip_alignment(args)
164
165
166 if __name__ == "__main__":
167 main()