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

根据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()

如果担心中间过程出错

使用try

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 = None
out_handle = None
try:
	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

finally:
	in_handle.close()
	out_handle.close()
	#不论try中是否有异常，finally中的内容都会执行

使用上下文管理器

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'

with open(infile) as in_handle, open(outfile,'w') as out_handle:
	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