Есть ли способ python для сопоставления фрагментов в пределах диапазона
У меня есть несколько генов, которые были охвачены различными чтениями. Одна проблема заключается в том, что в некоторых случаях, когда один ген покрывается несколькими прочтениями. Я хочу поймать все фрагменты прочтений, которые соответствуют длине гена (диапазону) [столбец 6 (генестарт) и столбец 7 генстоп], и добавить конечное поле, указывающее имя гена вместе с алфавитом «а», и для последовательных таких находок я хочу имя гена с добавлением b к нему.
Вход:
contig 230 3888 read1 1 8397 Gene1
contig 3343 3885 read2 1 8397 Gene1
contig 6030 8257 read4 1 8397 Gene1
contig 6030 8257 read10 1 8397 Gene1
contig 6260 8257 read7 1 8397 Gene1
contig 6030 8257 read8 1 8397 Gene1
contig 6190 6345 read10 1 8397 Gene1
contig 0 3124 read27 26 406 Gene12
contig 100 200 read30 26 406 Gene12
Желаемый результат:
contig 230 3888 read1 1 8397 Gene1 Gene1a
contig 3343 3885 read2 1 8397 Gene1 covered
contig 6030 8257 read4 1 8397 Gene1 Gene1b
contig 6030 8257 read10 1 8397 Gene1 covered
contig 6260 8257 read7 1 8397 Gene1 covered
contig 6030 8257 read8 1 8397 Gene1 covered
contig 6190 6345 read10 1 8397 Gene1 covered
contig 0 3124 read27 26 406 Gene12 Gene12a
contig 100 200 read30 26 406 Gene12 covered
Код:
with open(input_file, "r") as f_in, open(output_file, "w") as f_out:
gene_dict = {}
for line in f_in:
fields = line.strip().split()
gene_name = fields[6]
read_start = int(fields[1])
read_end = int(fields[2])
gene_start = int(fields[4])
gene_end = int(fields[5])
if gene_name not in gene_dict:
gene_dict[gene_name] = []
gene_found = False
for gene_frag in gene_dict[gene_name]:
frag_start = gene_frag[0]
frag_end = gene_frag[1]
if read_start <= frag_start <= read_end or read_start <= frag_end <= read_end:
gene_found = True
break
if gene_found:
gene_dict[gene_name].append((gene_start, gene_end))
gene_count = chr(ord('a') + len(gene_dict[gene_name]) - 1)
fields.append(f"{gene_name}{gene_count}")
else:
gene_dict[gene_name] = [(gene_start, gene_end)]
fields.append(f"{gene_name}")
f_out.write("\t".join(fields) + "\n")
Я не получаю желаемого результата. Результат, который я получил
contig 230 3888 read1 1 8397 Gene1 Gene1
contig 3343 3885 read2 1 8397 Gene1 Gene1
contig 6030 8257 read4 1 8397 Gene1 Gene1
contig 6030 8257 read10 1 8397 Gene1 Gene1
contig 6260 8257 read7 1 8397 Gene1 Gene1
contig 6030 8257 read8 1 8397 Gene1 Gene1
contig 6190 6345 read10 1 8397 Gene1 Gene1
contig 0 3124 read27 26 406 Gene12 Gene12
contig 100 200 read30 26 406 Gene12 Gene12
Файл разделен табуляцией, хотя здесь я не создал вкладку примера. Пожалуйста, укажите, где я ошибаюсь.
Кроме того, кстати, chr(ord('a')) — это действительно безумный способ написать просто 'a' (-:
Я хочу проверить, попадает ли диапазон значений в столбце 2 и столбце 3 в диапазон в столбцах 5 и 6. В моем случае несколько случаев диапазона в столбце 2 и столбце 3 попадают, поэтому я хочу добавить ген naem+алфавит, чтобы указать каждое совпадение с диапазоном. Для последовательных совпадений внутри гена я хочу написать имя гена и увеличить алфавит.
Для Gene1, read_start <= frag_start <= read_end or read_start <= frag_end <= read_end никогда не будет правдой. Пожалуйста, перепроверьте условия. Кроме того, Gene1a никогда не будет выводиться. Возможно вам нужен fields.append(f"{gene_name}a") за ген не найден. И, очевидно, ваш код не содержит никаких утверждений о covered.
Ответы 1
Приведенный ниже код имеет желаемый результат. Я взял ввод в виде многострочной строки и вывод через отпечатки для удобства с моей стороны, поэтому требуется некоторая адаптация.
Я пытался прокомментировать шаги, но спрашивайте, если не ясно. Я использовал логику range() и set(), чтобы найти совпадения.
data='''contig 230 3888 read1 1 8397 Gene1
contig 3343 3885 read2 1 8397 Gene1
contig 6030 8257 read4 1 8397 Gene1
contig 6030 8257 read10 1 8397 Gene1
contig 6260 8257 read7 1 8397 Gene1
contig 6030 8257 read8 1 8397 Gene1
contig 6190 6345 read10 1 8397 Gene1
contig 0 3124 read27 26 406 Gene12
contig 100 200 read30 26 406 Gene12 '''
#initialise variables
last_gene=''
gene_span=range(0)
total_read_span=range(0)
for contig in data.split('\n'): #iterate over contigs
fields=contig.split() #split lines into separate fields. Will need to be .split('\t) for tab-separated data
read_span=range(int(fields[1]),int(fields[2])+1)
if (gene := fields[6]) != last_gene: #check if new gene or previous gene
last_gene=gene
alphabet=iter('abcdefghijklmnopqrstuvqxyz') #set up an alphabet iterator to assign gene suffixes
gene_span=range(int(fields[4]),int(fields[5])+1)
total_read_span=range(0) #initialise total_read_span, i.e. the range of the sequence read
if len(set(read_span)&set(gene_span)): #check if read or gene spans overlap (the length of a set of intersecting values > 0)
if len(set(read_span)-set(total_read_span)): #check if the contig read has values outside of current total read span
print(f'{contig} {gene}{next(alphabet)}')
total_read_span=range(min(total_read_span.start,read_span.start),max(total_read_span.stop,read_span.stop)) #update total read span
else:
print(f'{contig} covered')
else:
print(f'{fields[3]} has no overlap with {gene}')
Спасибо получилось хорошо!! Итерация алфавитов была действительно хорошей.
хороший, что касается алфавита, в string библиотеке Python есть ascii_lowercase какая именно эта строка :) docs.python.org/3/library/string.html#string.ascii_lowercase
Это так, но реализация этого занимает примерно столько времени, сколько нужно напечатать, а потом я тоже должен это запомнить!
Тем из нас, кто не занимается биоинформатикой, трудно следить за вашим изложением. Проблема кажется несложной, но, возможно, вы захотите еще немного развернуть экспозицию.