Как получить количество последовательностей (в фасте) с условиями использования Python? - PullRequest
Новая
       21

Как получить количество последовательностей (в фасте) с условиями использования Python?

0 голосов
/ 18 апреля 2019

У меня есть файл fasta (fasta - это файл, в котором строка заголовка начинается с >, за которым следует строка последовательности, соответствующая этому заголовку).Я хочу получить количество последовательностей, соответствующих TRINITY, и общее количество последовательностей, которое начинается с >K после каждой >TRINITY последовательностей.Я смог получить счет для >TRINITY последовательностей, но не знал, как получить счет для >K для соответствующей группы >TRINITY последовательностей.Как я могу сделать это в Python?

myfasta.fasta: 

>TRINITY_DN12824_c0_g1_i1
TGGTGACCTGAATGGTCACCACGTCCATACAGA
>K00363:119:HTJ23BBXX:1:1212:18730:9403 1:N:0:CGATGTAT
CACTATTACAATTCTGATGTTTTAATTACTGAGACAT
>K00363:119:HTJ23BBXX:1:2228:9678:46223_(reversed) 1:N:0:CGATGTAT
TAGATTTAAAATAGACGCTTCCATAGA
>TRINITY_DN12824_c0_g1_i1
TGGTGACCTGAATGGTCACCACGTCCATACAGA
>K00363:119:HTJ23BBXX:1:1212:18730:9403 1:N:0:CGATGTAT
CACTATTACAATTCTGATGTTTTAATTACTGAGACAT
>TRINITY_DN555_c0_g1_i1
>K00363:119:HTJ23BBXX:1:2228:9658:46188_(reversed) 1:N:0:CGATGTAT
CGATGCTAGATTTAAAATAGACG
>K00363:119:HTJ23BBXX:1:2106:15260:10387_(reversed) 1:N:0:CGATGTAT
TTAAAATAGACGCTTCCATAGAGA

Результат, который я хочу: 

reference   reference_counts    Corresponding_K_sequences
>TRINITY_DN12824_c0_g1_i1   2   3
>TRINITY_DN555_c0_g1_i1 1   2

Вот код, который я написал, который толькоучитывает >TRINITY счетчиков последовательностей, но не может распространить его на бит, где он также будет считать соответствующие >K последовательности, поэтому любая помощь будет принята.Для запуска: python code.py myfasta.fasta output.txt

import sys
import os
from Bio import SeqIO
from collections import defaultdict
filename = sys.argv[1]
outfile = sys.argv[2]
dedup_records = defaultdict(list)

for record in SeqIO.parse(filename, "fasta"):
    #print(record)
    #print(record.id)
    if record.id.startswith('TRINITY'):
        #print(record.id)
    # Use the sequence as the key and then have a list of id's as the value
        dedup_records[str(record.seq)].append(record.id)
        #print(dedup_records)
with open(outfile, 'w') as output:
#   # to get the counts of duplicated TRINITY ids (sorted order)
    for seq, ids in sorted(dedup_records.items(), key = lambda t: len(t[1]), reverse=True):
        #output.write("{}   {}\n".format(ids,len(ids)))
        print(ids, len(ids))

1 Ответ

2 голосов
/ 18 апреля 2019

У вас правильное мышление, но вам нужно отслеживать последний заголовок, который начинается с «ТРИНИТИ», и немного изменить свою структуру: 

from Bio import SeqIO
from collections import defaultdict

TRIN, d = None, defaultdict(lambda: [0,0])

for r in SeqIO.parse('myfasta.fasta', 'fasta'):
    if r.id.startswith('TRINITY'):
        TRIN = r.id
        d[TRIN][0] += 1
    elif r.id.startswith('K'):
        if TRIN:
            d[TRIN][1] += 1

print('reference\treference_counts\tCorresponding_K_sequences')
for k,v in d.items():
    print('{}\t{}\t{}'.format(k,v[0],v[1]))

Выходы: 

reference   reference_counts    Corresponding_K_sequences
TRINITY_DN12824_c0_g1_i1    2   3
TRINITY_DN555_c0_g1_i1  1   2
...