使用python以滑动窗口方式统计基因组fasta文件中各序列的各碱基(如GC碱基)平均含量
001、
(base) root@PC1:/home/test2# ls a.fasta test.py (base) root@PC1:/home/test2# head a.fasta ## 测试数据 >scaffold_1 CCCGGGTAAAACGGGTCTTCAAGAAAACGCTCCTCCGTTAATGCCGGCCGATTCAAATAA CCTCTGGCAACACCCGCTCCGGCAATGTATAGTTCACCGATACATCCAACAGGCAGCATC CGCTGATTCTGATTCAGGATATACAATCTGACATGATGAACAGGTTTTCCAATTGGAATC CGTTCAAGTTTTTCTTGCGGCGGACAATCAAAGAATGCAGCTTCTACGGTTGCTTCCGTT GGCCCATAGGAATTGGTTATTGAAACATTTGGAAGCAACACGTGAAATCGGGAGACAAGA TGGGTCCCCAGCTGTTCTCCCCCAGAAAACACTCGCTTGAGTCTGTTTGTCTTAATCGGT ACAGAGCGATATTTTATATGTTCTAAAAATGCATGGAGCATTGAAGGCACAAAATGCATA GCTGTGATCTTTTGTTCTTCTATGGCCTTCGCGATCACTTCAGGTTCTTTTTCGCCTCCC TGAGGAAGCAGATAAACAGAAGCTCCGGCATAAGGCCACCAAAACAGCTCCCATATTGAA (base) root@PC1:/home/test2# cat test.py ## 测试程序 #!/usr/bin/python import re in_file = open("a.fasta", "r") out_file = open("result.txt", "w") dict1 = dict() for i in in_file: i = i.strip() if i[0] == ">": key = i dict1[key] = "" else: dict1[key] += i step = 5000 print("id", "start", "end", "a", "c", "g", "t", "cg", file = out_file, sep = "\t") for i,j in dict1.items(): length = len(j) start = 0 if length < step and start == 0: block = j[::] a = len(re.findall("[aA]", block)) c = len(re.findall("[cC]", block)) g = len(re.findall("[gG]", block)) t = len(re.findall("[tT]", block)) cg = len(re.findall("[cCgG]", block)) len_block = len(block) print(i, 1, len(j), a/len_block, c/len_block, g/len_block, t/len_block, cg/len_block, file = out_file, sep = "\t") while length > step: end = start + step block = j[start:end] a = len(re.findall("[aA]", block)) c = len(re.findall("[cC]", block)) g = len(re.findall("[gG]", block)) t = len(re.findall("[tT]", block)) cg = len(re.findall("[cCgG]", block)) len_block = len(block) print(i, start + 1, end, a/len_block, c/len_block, g/len_block, t/len_block, cg/len_block, file = out_file, sep = "\t") start += step length -= step if length < step: block = j[end::] a = len(re.findall("[aA]", block)) c = len(re.findall("[cC]", block)) g = len(re.findall("[gG]", block)) t = len(re.findall("[tT]", block)) cg = len(re.findall("[cCgG]", block)) len_block = len(block) print(i, end, len(j), a/len_block, c/len_block, g/len_block, t/len_block, cg/len_block, file = out_file, sep = "\t") in_file.close() out_file.close() (base) root@PC1:/home/test2# python test.py ## 运行程序 (base) root@PC1:/home/test2# ls a.fasta result.txt test.py (base) root@PC1:/home/test2# head result.txt ## 运行结果 id start end a c g t cg >scaffold_1 1 5000 0.2562 0.2254 0.2242 0.2942 0.4496 >scaffold_1 5001 10000 0.2416 0.233 0.2286 0.2968 0.4616 >scaffold_1 10001 15000 0.249 0.2544 0.221 0.2756 0.4754 >scaffold_1 15001 20000 0.2556 0.2422 0.2158 0.2864 0.458 >scaffold_1 20001 25000 0.2604 0.2352 0.2104 0.294 0.4456 >scaffold_1 25001 30000 0.2524 0.2404 0.219 0.2882 0.4594 >scaffold_1 30001 35000 0.3002 0.219 0.1934 0.2874 0.4124 >scaffold_1 35001 40000 0.287 0.219 0.2126 0.2814 0.4316 >scaffold_1 40001 45000 0.2424 0.25 0.2096 0.298 0.4596
参考:https://mp.weixin.qq.com/s?__biz=MzIxNzc1Mzk3NQ==&mid=2247491497&idx=1&sn=5759613bc25c175b71aaf5a5091d0a51&chksm=97f5afb1a08226a7b9eac750c7df99fdd6b5ec110f2f21e79ee0ef542caf93c8105d9b5358c4&scene=178&cur_album_id=2403674812188688386#rd