生信编程2.hg38基因组序列的一些探究
有一些五六年前的学生们都成长为了各个生物信息学相关公司的小领导，而且他们都有了自己的公众号，知乎号，也算是一番人物。最近他们跟我反馈面试找不到或者说很难直接考核筛选到认真干活的生信工程师，挺有意思的。让我想起来了早在生信技能树论坛创立之初我为了引流，而规划的200个生信工程师面试题。值得继续分享：
统计hg38基因组中的GC含量
统计N和GC含量==> 分别统计ATGCN的数量，然后使用（N数量/总数，G数量+C数量/总数）
首先要构建数据结构，也就是说用这种数据结构来统计碱基数量 python中的字典
{
'chr1':{
'A':n,
'T':n,
'C':n,
'G':n,
'N':n
},
'chr2':{
'A':n,
'T':n,
'C':n,
'G':n,
'N':n
},
}
初级脚本
hg38 = open('test.fa','r')
sum_atgc = {}
for line in hg38:
if line.startswith('>'):
chr_id = line[1:]
sum_atgc[chr_id] = {}
sum_atgc[chr_id]['A'] = 0
sum_atgc[chr_id]['T'] = 0
sum_atgc[chr_id]['G'] = 0
sum_atgc[chr_id]['C'] = 0
sum_atgc[chr_id]['N'] = 0
else:
line = line.upper()
sum_atgc[chr_id]['A'] += line.count('A')
sum_atgc[chr_id]['T'] += line.count('T')
sum_atgc[chr_id]['G'] += line.count('G')
sum_atgc[chr_id]['C'] += line.count('C')
sum_atgc[chr_id]['N'] += line.count('N')
for chr_id,atgc_count in sum_atgc.items():
print(chr_id)
GC = int(atgc_count['G']) + int(atgc_count['C'])
total = sum(atgc_count.values())
print("A : %s" % (atgc_count['A']))
print( "T : %s" % (atgc_count['T']))
print( "G : %s" % (atgc_count['G']))
print( "C : %s" % (atgc_count['C']))
print( "N : %s" % (atgc_count['N']))
print( "total : %s" % (sum(atgc_count.values())))
print("CG_per : %s" % (GC/total))
print("N_per : %s" % (int(atgc_count['N'])/total))
虽然逻辑上符合，而且也能够计算出目的的结果，但是重复代码太多，需要进行优化 ==注意== 由于基因组中可能同时存在大小写不同的碱基，所以在进行统计的时候需要统一大小写line.upper()
减少代码重复问题
#!/usr/bin/python
import sys
args = sys.argv
hg38 = open(args[1],'r')
bases = ['A','T','G','C','N']
sum_atgc = {}
for line in hg38:
if line.startswith('>'):
chr_id = line[1:]
sum_atgc[chr_id] = {}
for base in bases:
sum_atgc[chr_id][base] = 0
else:
line = line.upper()
for base in bases:
sum_atgc[chr_id][base] += line.count(base)
for chr_id,atgc_count in sum_atgc.items():
print(chr_id)
GC = int(atgc_count['G']) + int(atgc_count['C'])
total = sum(atgc_count.values())
for base in bases:
print("%s : %s" % (base, atgc_count[base]))
print( "total : %s" % (total))
print("CG_per : %s" % (GC/total))
print("N_per : %s" % (int(atgc_count['N'])/total))
==注意== 将需要进行处理的碱基放到一个列表中，然后对列表进行循环，这样就能避免重复写相同的代码
从1开始写python
0.5 开始课程 1 基础语法 2 编程技巧（打开文件，参数处理等） 3 效率（程序的运行时间）
提高效率
#!/usr/bin/python
import sys
import time
import re
start = time.process_time()
args = sys.argv
sum_atgc = {}
bases = ['A','T','G','C','N']
def get_chr(buf):
(buf, tmp) = buf.split('>',1)
if '>' in tmp:
get_chr(tmp)
return buf,tmp
with open(args[1],'r') as f:
tmp = f.readline()
chr_id = re.split(r'\s',tmp)[0][1:]
sum_atgc[chr_id] = {}
for base in bases:
sum_atgc[chr_id][base] = 0
while 1:
buffer = f.read(1024 * 1024)
if not buffer:
break
if '>' in buffer:
#print(len(buffer.split('>')))
(buffer,tmp) = get_chr(buffer)
buffer = buffer.upper()
for base in bases:
sum_atgc[chr_id][base] += buffer.count(base)
(tmp, buffer) = tmp.split('\n', 1)
chr_id = re.split(r'\s', tmp)[0][1:]
sum_atgc[chr_id] = {}
for base in bases:
sum_atgc[chr_id][base] = 0
buffer = buffer.upper()
for base in bases:
sum_atgc[chr_id][base] += buffer.count(base)
else:
buffer = buffer.upper()
for base in bases:
sum_atgc[chr_id][base] += buffer.count(base)
for chr_id,atgc_count in sum_atgc.items():
print(chr_id)
GC = int(atgc_count['G']) + int(atgc_count['C'])
total = sum(atgc_count.values())
for base in bases:
print("%s : %s" % (base, atgc_count[base]))
print( "total : %s" % (total))
print("CG_per : %s" % (GC/total))
print("N_per : %s" % (int(atgc_count['N'])/total))
end = time.process_time()
print("time : %s" % (int(end-start)))
==解析1== 使用f.read()读取文件比较慢，通过使用buffer读取可以提升速度，每次读取的字节数为（1024*1024）。对于fasta格式的基因组文件来说，需要区分'>'开头的序列名和序列。首先要读第一行的序列名称tmp=f.readline()，后面的内容就可以按照buffer的要求来读取。一个buffer读入的文件有三种形式：
ATCGATGACTGA\n >chr1 \nGCTGATGAC
包含序列，中间混入一个序列名称 ATCGATGCA\n> chrx \nAGCTGATCGATCGATCGATGATCGATGCATGCTA\n >chry \nATGCTAGCTAGTC
对于某些比较短的染色体来说，可能存在序列中混入了两个序列名称的情况 ATCGATGCGGCTAGCTAGCTGATCGAGGCCGATGACG
只有序列信息，没有序列名称
==解析2== 使用while循环读入buffer时，一定要判断是否有buffer，否则文件读完之后，也不会退出循环
==解析3== 对buffer中含有两个甚至多个以'>'开头的序列名的情况进行处理：最好的解决方式是递归
def get_chr(buf):
(buf, tmp) = buf.split('>',1)
if '>' in tmp:
get_chr(tmp)
return buf,tmp