Подсчет последовательностей ДНК с помощью python/biopython

Мой сценарий ниже подсчитывает количество вхождений последовательностей «CCCCAAAA» и «GGGGTTTT» из стандартного файла FASTA:

>contig00001  
CCCCAAAACCCCAAAACCCCAAAACCCCTAcGAaTCCCcTCATAATTGAAAGACTTAAACTTTAAAACCCTAGAAT

Скрипт подсчитывает здесь последовательность CCCCAAAA 3 раза.

CCCCAAAACCCCAAAACCCCAAAA(CCCC не в счет)

Может кто-нибудь посоветовать, как включить последовательность CCCC в конце как половинный счет, чтобы вернуть для этого значение 3,5.

Мои попытки пока безуспешны.

Мой сценарий выглядит следующим образом...

from Bio import SeqIO

input_file = open('telomer.test.fasta', 'r')
output_file = open('telomer.test1.out.tsv','w')
output_file.write('Contig\tCCCCAAAA\tGGGGTTTT\n')

for cur_record in SeqIO.parse(input_file, "fasta") :


    contig = cur_record.name
    CCCCAAAA_count = cur_record.seq.count('CCCCAAAA')
    CCCC_count = cur_record.seq.count('CCCC')

    GGGGTTTT_count = cur_record.seq.count('GGGGTTTT')
    GGGG_count = cur_record.seq.count('GGGG')
    #length = len(cur_record.seq)

    splittedContig1=contig.split(CCCCAAAA_count)

    splittedContig2=contig.split(GGGGTTTT_count)

    cnt1=len(splittedContig1)-1
    cnt2=len(splittedContig2)

  cnt1+sum([0.5 for e in splittedContig1 if e.startswith(CCCC_count)])) = CCCCAAAA_count
  cnt2+sum([0.5 for e in splittedContig2 if e.startswith(GGGG_count)])) = GGGGTTTT_count

    output_line = '%s\t%i\t%i\n' % \
    (CONTIG, CCCCAAAA_count, GGGGTTTT_count)


    output_file.write(output_line)

output_file.close()

input_file.close()

python bioinformatics biopython dna-sequence
person sheaph    schedule 01.04.2014    source источник

Ответы (1)


arrow_upward
2
arrow_downward

Вы можете использовать понимание списка split и startwith следующим образом:

contig="CCCCAAAACCCCAAAACCCCAAAACCCCTAcGAaTCCCcTCATAATTGAAAGACTTAAACTTTAAAACCCTAGAAT"
splitbase="CCCCAAAA"
halfBase="CCCC"
splittedContig=contig.split(splitbase)
cnt=len(splittedContig)-1
print cnt+sum([0.5 for e in splittedContig if e.startswith(halfBase)])

Выход:

3.5
  1. разделить строки на основе CCCCAAAA. Дал бы список, в списке будут удалены элементы CCCCAAAA
  2. длина разделенного - 1 дает количество вхождений CCCCAAAA
  3. в разделенном элементе ищите элементы, начинающиеся с CCCC. Если найдено, добавьте 0,5 к счету для каждого случая.
person venpa    schedule 01.04.2014
comment
Это работает очень хорошо само по себе, но когда я пытаюсь интегрировать в свой сценарий, я получаю сообщение об ошибке: SyntaxError: can't assign to function call Возможно ли это сделать? - person sheaph; 01.04.2014
comment
Как вы интегрировались со своим сценарием? - person venpa; 01.04.2014
comment
Я новичок в Био. Насколько я понимаю, вы должны сделать: cur_record.id.split('CCCCAAAA') - person venpa; 01.04.2014
comment
Кроме того, когда вы присваиваете результат, вы должны изменить LHS на RHS. Например, ваше назначение должно быть CCCCAAAA_count = cnt1+sum([0,5 для e в splittedContig1, если e.startswith(CCCC_count)])) - person venpa; 01.04.2014