python中计算dna序列的GC含量
001、对G、C计数进行统计
[root@pc1 test01]# ls a.fa test.py [root@pc1 test01]# cat a.fa ## 测试DNA序列 >Rosalind_6404 CCTGCGGAAGATCGGCACTAGAATAGCCAGAACCGTTTCTCTGAGGCTTCCGGCCTTCCC TCCCACTAATAATTCTGAGG >Rosalind_5959 CCATCGGTAGCGCATCCTTAGTCCAATTAAGTCCCTATCCAGGCGCTCCGCCGAAGGTCT ATATCCATTTGTCAGCAGACACGC >Rosalind_0808 CCACCCTCGTGGTATGGCTAGGCATTCAGGAACCGGAGAACGCTTCAGACCAGCCCGGAC TGGGAACCTGCGGGCAGTAGGTGGAAT [root@pc1 test01]# cat test.py ## 统计程序 #!/usr/bin/env python # -*- coding: utf-8 -*- in_file = open("a.fa", "r") dict1 = dict() for i in in_file: i = i.strip() if i.startswith(">"): temp = i dict1[temp] = str() else: dict1[temp] += i in_file.close() dict2 = dict() for i,j in dict1.items(): count = 0 for k in j: if k == "C" or k == "G": count += 1 dict2[i] = count/len(j) for i,j in dict2.items(): print("%s %.2f" % (i, j * 100)) max_value = max(dict2.values()) for i,j in dict2.items(): if j == max_value: print(i.replace(">", "")) print("%.6f" % (j * 100)) [root@pc1 test01]# python3 test.py ## 输出结果 >Rosalind_6404 53.75 >Rosalind_5959 53.57 >Rosalind_0808 60.92 Rosalind_0808 60.919540
002、利用函数结构实现
[root@pc1 test01]# ls a.fa test.py [root@pc1 test01]# cat a.fa ## 测试fasta文件 >Rosalind_6404 CCTGCGGAAGATCGGCACTAGAATAGCCAGAACCGTTTCTCTGAGGCTTCCGGCCTTCCC TCCCACTAATAATTCTGAGG >Rosalind_5959 CCATCGGTAGCGCATCCTTAGTCCAATTAAGTCCCTATCCAGGCGCTCCGCCGAAGGTCT ATATCCATTTGTCAGCAGACACGC >Rosalind_0808 CCACCCTCGTGGTATGGCTAGGCATTCAGGAACCGGAGAACGCTTCAGACCAGCCCGGAC TGGGAACCTGCGGGCAGTAGGTGGAAT [root@pc1 test01]# cat test.py ## 统计程序 #!/usr/bin/env python # -*- coding: utf-8 -*- in_file = open("a.fa", "r") dict1 = dict() for i in in_file: i = i.strip() if i[0] == ">": temp = i dict1[temp] = "" else: dict1[temp] += i def gc(sequence): return (sequence.count("C") + sequence.count("G")) * 100/len(sequence) def max_gc(fa): dict2 = {} for i,j in fa.items(): dict2[i] = gc(j) return max(dict2.items(), key = lambda x:x[1]) print(max_gc(dict1)) [root@pc1 test01]# python3 test.py ## 统计结果 ('>Rosalind_0808', 60.91954022988506)
003、改进
[root@pc1 test01]# ls a.fa test.py [root@pc1 test01]# cat a.fa ## 测试DNA序列 >Rosalind_6404 CCTGCGGAAGATCGGCACTAGAATAGCCAGAACCGTTTCTCTGAGGCTTCCGGCCTTCCC TCCCACTAATAATTCTGAGG >Rosalind_5959 CCATCGGTAGCGCATCCTTAGTCCAATTAAGTCCCTATCCAGGCGCTCCGCCGAAGGTCT ATATCCATTTGTCAGCAGACACGC >Rosalind_0808 CCACCCTCGTGGTATGGCTAGGCATTCAGGAACCGGAGAACGCTTCAGACCAGCCCGGAC TGGGAACCTGCGGGCAGTAGGTGGAAT [root@pc1 test01]# cat test.py ## 计算程序 #!/usr/bin/env python # -*- coding: utf-8 -*- in_file = open("a.fa", "r") dict1 = dict() for i in in_file: i = i.strip() if i[0] == ">": temp = i dict1[temp] = "" else: dict1[temp] += i def gc(sequence): return (sequence.count("C") + sequence.count("G")) * 100/len(sequence) def max_gc(fa): dict2 = {} for i,j in fa.items(): dict2[i] = gc(j) tmp = max(dict2.items(), key = lambda x:x[1]) return (tmp[0].lstrip(">") + "\n" + str(tmp[1])) print(max_gc(dict1)) [root@pc1 test01]# python3 test.py ## 计算结果 Rosalind_0808 60.91954022988506
。