使用Python/Biopython计算DNA序列数目

Question

使用Python/Biopython计算DNA序列数目

pythonbioinformaticsbiopythondna-sequence

3

我的脚本在从标准FASTA文件中计算序列'CCCCAAAA'和'GGGGTTTT'的出现次数：

>contig00001  
CCCCAAAACCCCAAAACCCCAAAACCCCTAcGAaTCCCcTCATAATTGAAAGACTTAAACTTTAAAACCCTAGAAT

这个脚本在这里计数了3次CCCCAAAA序列

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()

- sheaph

1个回答

网页内容由stack overflow 提供, 点击上面的

可以查看英文原文，
原文链接

- venpa · Accepted Answer

你可以使用以下方式的分割、开始列表解析:

您可以使用split和startswith列表推导式，如下所示：

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

根据 CCCCAAAA 分割字符串。这将给出一个列表，在列表元素中，CCCCAAAA 将被删除。
分割后的长度 - 1 给出了 CCCCAAAA 的出现次数。
在分割后的元素中，查找以 CCCC 开头的元素。如果找到，则每次出现将计数增加 0.5。