-
-
Notifications
You must be signed in to change notification settings - Fork 3
/
Copy pathVCF_to_TSV.py
209 lines (192 loc) · 17.3 KB
/
VCF_to_TSV.py
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
print('''
Python3-скрипт, конвертирующий VCF4-файлы (Variant Call Format, version 4.x) в TSV (Tab-Separated Values).
Автор: Платон Быкадоров (platon.work@gmail.com), 2018.
Версия: V2.1.
Лицензия: GNU General Public License version 3.
Поддержать проект: https://money.yandex.ru/to/41001832285976
Простое руководство по установке среды разработки и запуску скриптов:
github.com/PlatonB/bioinformatic-python-scripts#Установка-среды-разработки
Спецификации формата VCF:
https://samtools.github.io/hts-specs/VCFv4.2.pdf
Работа скрипта сосредоточена на восьмом (INFO) столбце VCF.
Именно этот столбец нуждается в табификации, чтобы
информация оттуда стала парсибельнее и, возможно, нагляднее.
Скрипт также выполняет роль валидатора
правильности оформления ваших VCF-файлов.
''')
import os, sys, re
from collections import OrderedDict
src_dir_path = input('Путь к папке с исходными VCF-файлами: ')
trg_dir_path = input('\nПуть к папке для конечных файлов: ')
meta_lines_fate = input('''\nСохранять строки мета-информации в конечные файлы?
(игнорирование ввода ==> сохранять)
[yes(|y|<enter>)|no(|n)]: ''')
if meta_lines_fate != 'yes' and meta_lines_fate != 'y' and \
meta_lines_fate != 'no' and meta_lines_fate != 'n' and meta_lines_fate != '':
print(f'\n{meta_lines_fate} - недопустимая опция')
sys.exit()
#Обязательные согласно стандарту VCF заголовки
#полей (столбцов) и их порядковые номера.
#Это множество пригодится для проверки
#правильности шапки в пользовательских таблицах.
#Далее оно будет называться референсом.
fixed_fields = set(['1_#CHROM',
'2_POS',
'3_ID',
'4_REF',
'5_ALT',
'6_QUAL',
'7_FILTER',
'8_INFO'])
#Работа с исходными файлами.
src_file_names = os.listdir(src_dir_path)
for src_file_name in src_file_names:
if src_file_name.startswith('.~lock.') == True:
continue
print('Обрабатывается файл', src_file_name)
with open(os.path.join(src_dir_path, src_file_name)) as src_file_opened:
#Построение имени конечного файла.
#Открытие конечного файла на запись.
trg_file_name = '.'.join(src_file_name.split('.')[:-1]) + '_tab.tsv'
with open(os.path.join(trg_dir_path, trg_file_name), 'w') as trg_file_opened:
#Прочтение первой строки VCF-файла.
#Проверка по мягкому критерию: хотя
#бы одна строка в начале файла должна
#быть строкой мета-информации.
#Здесь и далее, в случае выявления
#несоответствий стандарту VCF, уже созданный
#текущий конечный файл будет удалён.
line = src_file_opened.readline().split('\n')[0]
if line.startswith('##') == False:
print(f'''\nОшибка в файле {src_file_name}:
VCF-файл должен содержать строки мета-информации, начинающиеся с ##''')
os.remove(os.path.join(trg_dir_path, trg_file_name))
sys.exit()
#В каждой итерации цикла while осуществляется
#как проверка считанной в прошлой итерации строки
#на наличие подстроки ##, так и прочтение новой строки.
#Строки, определённые как строки мета-информации,
#в зависимости от выбора пользователя могут быть
#прописаны или не прописаны в конечный файл.
#Проверка строк мета-информации уже по более
#жёсткому критерию: среди них должно быть не менее
#одной строки, описывающей элементы столбца INFO.
#Извлечение из строк мета-информации названий этих элементов
#и их сохранение в упорядоченный словарь в качестве ключей.
#Значения упорядоченного словаря - литералы
#пустой строки, что будет объяснено позже.
#После строк мета-информации должна идти шапка с одним #.
#Считывание первой же строки без ## приведёт к выходу из цикла.
info_structure = OrderedDict()
while line[:2] == '##':
if line.startswith('##INFO=<ID=') == True:
info_subfield_id = re.search(r'(?<=##INFO=<ID=)[^,]+', line).group()
info_structure[info_subfield_id] = ''
if meta_lines_fate == 'yes' or meta_lines_fate == 'y' or meta_lines_fate == '':
trg_file_opened.write(line + '\n')
line = src_file_opened.readline().split('\n')[0]
if info_structure == OrderedDict():
print(f'''\nОшибка в файле {src_file_name}:
VCF-файл должен содержать строки мета-информации, описывающие INFO-столбец''')
os.remove(os.path.join(trg_dir_path, trg_file_name))
sys.exit()
#Первый этап проверки, является ли строка,
#идущая после строк мета-информации, шапкой.
#Заключается он в выявлении наличия #
#в качестве первого символа строки.
if line.startswith('#') == False:
print(f'''\nОшибка в файле {src_file_name}:
VCF-файл должен содержать шапку, начинающуюся с #''')
os.remove(os.path.join(trg_dir_path, trg_file_name))
sys.exit()
#Второй, более сложный этап.
#Сверяем первые 8 элементов шапки с референсом.
#Для этого создаём множество, элементы которого
#построены идентичным образом с референсными:
#номер столбца_заголовок столбца.
#Сверяем количество элементов шапки
#пользовательского VCF и референса.
#Если оно не одинаково, то вычитаем из
#пользовательского множества референсное, выявляя
#тем самым неправильно озаглавленные столбцы.
#Сортируем и перебираем полученные неправильные элементы
#и выводим сообщения об ошибках, где указываем позиции
#ошибочных заголовков столбцов и сами заголовки.
header_row = line.split('\n')[0].split('\t')
header_check_set = set(map(lambda col_num, cell: col_num + '_' + cell, [str(i) for i in range(1, 9)][:8], header_row))
if header_check_set != fixed_fields:
for cell in sorted(list(header_check_set - fixed_fields)):
print(f'''\nОшибка в файле {src_file_name}:
{cell[2:]} - неправильный заголовок столбца {cell[:1]} VCF-файла''')
os.remove(os.path.join(trg_dir_path, trg_file_name))
sys.exit()
#Шапка прошла проверку, но пока не приняла
#свой окончательный вид, поскольку ещё не заменён
#на подзаголовки заголовок "INFO" 8-го столбца.
#В качестве подзаголовков будут ключи словаря,
#созданного ранее на основе строк мета-информации.
#Пересобираем и прописываем в конечный файл шапку.
header_row[7] = '\t'.join(info_structure.keys())
trg_file_opened.write('\t'.join(header_row) + '\n')
#Эффективная с точки зрения использования
#RAM работа с основной частью таблицы.
for line in src_file_opened:
data_row = line.split('\n')[0].split('\t')
#Перебираем элементы INFO-столбца.
#С точки зрения обработки данных, их можно
#подразделить на 2 типа: ключ-значение и флаг.
#Особенность флагов - у них не бывает значений.
#Ключи в парах ключ-значение и флаги представляют собой
#константы, а значения могут от строки к строке меняться.
#Общее свойство любых элементов - что они могут
#присутствовать, а могут отсутствовать в текущей строке.
#Отсутствующие элементы значительно осложняют табификацию,
#т.к. сложно предсказать, в какое место табулированной
#строки размещать непустой элемент, а в какое - пустой.
#Игнорирование отсутствующих элементов INFO-столбца
#приведёт к многочисленным смещениям его крайних ячеек
#влево и вправо, поэтому, если INFO-столбец - не последний,
#идущие после него данные потеряют табличный вид.
#Тут выручает упорядоченный словарь, хранящий
#заголовки столбцов, формируемых из элементов INFO.
#Если VCF составлен правильно, то и ключи элементов
#типа ключ-значение, и элементы-флаги будут совпадать
#с заголовками из этого упорядоченного словаря.
#Значения словаря - литералы пустой строки - меняются
#либо на значения пар ключ-значение, либо на флаги.
#Если в текущей INFO-строке отсутствуют некоторые
#элементы, то в качестве соответствующих значений
#ключей-заголовков останутся литералы пустой строки.
#Далее из значений упорядоченного словаря собирается
#табулированная строка, в которой каждый
#отсутствующий в INFO-столбце элемент будет
#как бы пустой ячейкой между двумя табуляциями.
for cell_element in data_row[7].split(';'):
info_key = cell_element.split('=')[0]
try:
info_val = cell_element.split('=')[1]
if info_key in info_structure:
info_structure[info_key] = info_val
else:
print(f'''\nОшибка в файле {src_file_name}:
Отсутствует или неправильно составлена строка мета-информации,
описывающая элемент {info_key} INFO-строк VCF-файла''')
os.remove(os.path.join(trg_dir_path, trg_file_name))
sys.exit()
except IndexError:
if info_key in info_structure:
info_structure[info_key] = cell_element
else:
print(f'''\nОшибка в файле {src_file_name}:
Отсутствует или неправильно составлена строка мета-информации,
описывающая элемент {info_key} INFO-строк VCF-файла''')
os.remove(os.path.join(trg_dir_path, trg_file_name))
sys.exit()
data_row[7] = '\t'.join(info_structure.values())
#Прописывание обновлённых строк
#основной части таблицы в конечный файл.
trg_file_opened.write('\t'.join(data_row) + '\n')
#Перед прочтением следующей строки необходимо привести словарь
#в исходный вид, т.е. когда значения - литералы пустой строки.
for key in info_structure.keys():
info_structure[key] = ''