根据sam文件计算reads的GC含量
禁转
输入文件
DNA序列的sam文件
第一列,序列名;第十列,序列;分割 tab
目标
计算每个read的GC含量(只考虑DNA序列由ATGC组成的情况),并输出结果到文件
输出文件
第一列read id,第二列read的序列,第三列GC含量
代码
def cal_gc(seq):
"""
计算GC含量
"""
seq = seq.upper() #全部转换为大写
return round((seq.count('C') + seq.count('G')) / float(len(seq)),2)
#round 用于保留两位小数,float将序列长度转换为浮点型
infile = 'data/alignment.sam'
outfile = 'output/gc.txt'
in_handle = open(infile) #open(文件路径,打开方式默认为r)
out_handle = open(outfile, 'w')
out_handle.write('read_id\tseq\tgc\n') #表头
cnt = 0
for line in in_handle:
if cnt != 0: #用于有列名时去掉第一行列名
row = line.split('\t') #返回列表
read_id = row[0]
seq = ro