Подсчитать определенный шаблон в seqID python

Question

У меня на самом деле есть огромный файл multitifasta seq, такой:

>Seq_1_0035_0035
ATTGGAT
>Seq_2_0042_0035
ATTGAGGA
>EOGWX56TR_0035_0042 (busco)
ATGGAGAT
>EOGWX56TR_0042_0042 (busco)
ATGGATGG
>Seq6_035_0042
ATGGGAATAG
>EOG55FTG_0035_0042 (busco)
AATGGATA
>EOG5GFFTA_0042_0042 (busco)
ATGGAGATA
>Seq56_0035_0042
ATGGAGATAT
>EOGATTT_0035_0042  (busco)
AAATGAGATA
>EOGATTT_0042_0042  (busco)
ATGGAAT
>EOGATTA_0042_0042  (busco)
ATAGGAGAT

, и я действительно хочу посчитать, сколько генов Busco у меня есть в моем файле (все они начинаются с имени >EOG).поэтому у меня есть такой скрипт:

count=1
for record in SeqIO.parse("concatenate_with_busco_names_0035_0042_aa.fa", "fasta"):
    count+=1
print(count)

set_of_labels = set()

with open("concatenate_with_busco_names_0035_0042_aa.fa") as f:
  for line in f:
    if line.startswith('>EOG'):
      label = line[4:].split('_')[0]
      set_of_labels.add(label)

print("Total number of Busco genes: " + str(len(set_of_labels)))

Но мне также хотелось бы узнать, сколько у меня генов между каждым соответствующим компортиментом.Я объясняю это лучше:

Как вы можете видеть, в каждом seqID есть два числа such _number_number Эти числа являются конкретными, и первое _number соответствует виду, к которому относится последовательность, а второе _numberэто конкретный номер.В любом случае мне бы хотелось, если бы можно было посчитать, сколько я получу, сколько разных генов Busco я получу для seq с первым номером _0035 и _0042 И также сколько для seq ID:

_0035_0042
_0035_0042
_0042_0042
_0042_0035

В приведенном выше примере это будет:

Total busco: 5 (I count only once if the >busco is present even if _number are different)
Total busco for the specie _0035 (_0035_0042 and _0035_0035) : 3
Total busco for the specie _0042 (_0042_0042 and _0042_0035) : 4
Total busco for the specific specie  _0035_0042 : 3
Total busco for the specific specie  _0042_0035 : 0
Total busco for the specific specie  _0042_0042 : 4
Total busco for the specific specie  _0035_0035 : 0

Привет, надеюсь, это понятно, ведь первая часть (total busco:) уже выполнена моим сценарием, мне нужно только посчитать 7 других способов.

вот реальные данные данные

Answer 1

Это тривиально для класса Counter из стандартной библиотеки Python:

from collections import Counter
from io import StringIO

label_counter = Counter()
specy_counter = Counter()
specific_specy_counter = Counter()

# replace this with an open() on your real file 
finput = StringIO(""">Seq_1_0035_0035
ATTGGAT
>Seq_2_0042_0035
ATTGAGGA
>EOGWX56TR_0035_0042 (busco)
ATGGAGAT
>EOGWX56TR_0042_0042 (busco)
ATGGATGG
>Seq6_035_0042
ATGGGAATAG
>EOG55FTG_0035_0042 (busco)
AATGGATA
>EOG5GFFTA_0042_0042 (busco)
ATGGAGATA
>Seq56_0035_0042
ATGGAGATAT
>EOGATTT_0035_0042  (busco)
AAATGAGATA
>EOGATTT_0042_0042  (busco)
ATGGAAT
>EOGATTA_0042_0042  (busco)
ATAGGAGAT""")



for line in finput:
    try:
        if line.startswith('>EOG'):
            label, specy, specific = line[4:].replace(" (busco)", "").strip().split('_')
            label_counter[label] += 1
            specy_counter[specy] += 1
            specific_specy_counter[(specy, specific)] += 1
    except ValueError:
        print("Invalid line:", line)


print("Total busco:", len(label_counter))
for specy, count in specy_counter.items():
    print("Total busco for the specie {} : {}".format(specy, count))
for (specy, specific), count in specific_specy_counter.items():
    print("Total busco for the specific specy {}_{} : {}".format(specy, specific, count))

Обратите внимание, что виды или особенности с 0 значениями не будут отображаться.

Answer 2

Вы можете использовать несколько счетчиков, чтобы получить отдельные счетчики для видов и определенных видов в дополнение к счетчику busco , например:

import collections

busco = collections.defaultdict(int)  # busco counter
species = collections.defaultdict(int)  # species counter
specific_species = collections.defaultdict(int)  # specific species counter

with open("concatenate_with_busco_names_0035_0042_aa.fa", "r") as f:
    for line in f:
        if line[:4] == ">EOG":
            entry = line.split()[0][4:].split('_')
            busco[entry[0]] += 1
            species[entry[1]] += 1
            specific_species[entry[1] + "_" + entry[2]] += 1

print("Total busco: {}".format(len(busco)))
for specie, total in species.items():
    print("Total busco for the specie {}: {}".format(specie, total))
for specie, total in specific_species.items():
    print("Total busco for the specific specie {}: {}".format(specie, total))

Что должно дать:

Total busco: 5
Total busco for the specie 0035: 3
Total busco for the specie 0042: 4
Total busco for the specific specie 0035_0042: 3
Total busco for the specific specie 0042_0042: 4

Внесенные в список (конкретные) виды не будут отображаться, но если вы действительно хотите распечатать их, вы можете объединить их со счетчика species и распечатать ихзначения (по умолчанию 0):

import itertools

print("Total busco: {}".format(len(busco)))
for specie, total in species.items():
    print("Total busco for the specie {}: {}".format(specie, total))
for specie in itertools.product(species, species):
    s = "_".join(specie)
    print("Total busco for the specific specie {}: {}".format(s, specific_species[s]))

Что дает:

Total busco: 5
Total busco for the specie 0035: 3
Total busco for the specie 0042: 4
Total busco for the specific specie 0035_0035: 0
Total busco for the specific specie 0035_0042: 3
Total busco for the specific specie 0042_0035: 0
Total busco for the specific specie 0042_0042: 4

ОБНОВЛЕНИЕ : если вы хотите получить уникальные значения для busco, затем вам нужно инвертировать счет в индекс для specie / конкретный образец и собрать busco значения в set в качестве их значения.Тогда все, что вам нужно, это получить длину каждого набора, что-то вроде:

import collections
import itertools

busco = set()
species = collections.defaultdict(set)
specific_species = collections.defaultdict(set)

with open("concatenate_with_busco_names_0035_0042_aa.fa", "r") as f:
    for line in f:
        if line[:4] == ">EOG":
            entry = line.split()[0][4:].split('_')
            busco.add(entry[0])
            species[entry[1]].add(entry[0])
            specific_species[entry[1] + "_" + entry[2]].add(entry[0])

print("Total busco: {}".format(len(busco)))
for specie, buscos in species.items():
    print("Total busco for the specie {}: {}".format(specie, len(buscos)))
for specie in itertools.product(species, species):
    s = "_".join(specie)
    print("Total busco for the specific specie {}: {}".format(s, len(specific_species[s])))

Что для печати ваших полных данных:

Total busco: 421
Total busco for the specie 0035: 402
Total busco for the specie 0042: 397
Total busco for the specific specie 0035_0035: 392
Total busco for the specific specie 0035_0042: 262
Total busco for the specific specie 0042_0035: 305
Total busco for the specific specie 0042_0042: 383