根据sam文件计算reads的GC含量

最新推荐文章于 2023-11-20 10:53:19 发布

原创

最新推荐文章于 2023-11-20 10:53:19 发布 · 1.1k 阅读

1 ·

CC 4.0 BY-SA版权

本文介绍如何根据sam文件计算DNA序列reads的GC含量。内容包括输入sam文件的格式，计算目标，输出结果文件的结构，以及在代码实现过程中如何处理可能出现的错误，如使用try语句和上下文管理器确保程序的稳定性。

根据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 = row[9]
		gc = cal_gc(seq)
		#print([read_id, seq, gc]) #检查
		out_handle.write('%s\t%s\t%s\n' % (read_id,seq,gc))   # %s占位符
	cnt += 1

in_handle.close()
out_handle.close()