根据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 += 1in_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 = 0for 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 += 1finally: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 = 0for 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

根据sam文件计算reads的GC含量相关推荐

python学习——把计算GC含量的代码封装成函数
把代码封装成函数的好处是可以重复使用该段代码,并且会使代码结构清晰例如要计算chr1以及chr2染色体的GC含量,代码如下: 1 # 将代码封装为函数并重复使用,例如计算染色体的GC含量 2 chr ...
下载FATSQ，读取10条序列并计算每条序列的长度和GC含量
FASTQ 这是目前存储测序数据最普遍.最公认的一个数据格式,另一个是uBam格式,但这篇文章中不打算对其进行介绍.上面所讲的FASTA文件,它所存的都是已经排列好的序列(如参考序列),FASTQ存的 ...
推荐一个SAM文件中flag含义解释工具--转载
SAM是Sequence Alignment/Map 的缩写.像bwa等软件序列比对结果都会输出这样的文件.samtools网站上有专门的文档介绍SAM文件.具体地址:http://samtools. ...
31、SAM文件中flag含义解释工具--转载
转载:http://www.cnblogs.com/nkwy2012/p/6362996.html SAM是Sequence Alignment/Map 的缩写.像bwa等软件序列比对结果都会输出这样 ...
html文件下的flag,推荐一个SAM文件中flag含义解释工具
SAM是Sequence Alignment/Map 的缩写.像bwa等软件序列比对结果都会输出这样的文件.samtools网站上有专门的文档介绍SAM文件.具体地址:http://samtools. ...
【bioinfo】sam文件可选区域字段(Optional Feild)含义
写在前面 SAM文件格式说明 SAM文件可选字段说明上面两个链接是关于sam文件格式的官方说明文档,第一个包含所有字段,但对可选字段的说明比较简略,第二个是对可选字段的详细说明.这里根据官方文档,说 ...
bowtie结果sam文件解读
sam文件解读 @HD VN:1.0 SO:unsorted @SQ SN:chr1 LN:249250621 @SQ SN:chr2 LN:243199373 @ ...
攻破Administrator权限－－破解SAM文件法
攻破Administrator权限--破解SAM文件法事情是这样子的: 朋友来找我,说: 现在有一台PC,配置很一般,但是又常用,可是机子的Administrator用户被别人控制了,也就是他没有密 ...
windows 账户SAM文件损坏的解决办法
SAM文件损坏的解决办法<?xml:namespace prefix = o ns = "urn:schemas-microsoft-com:office:office" / ...
WINPE的DOS工具箱删除sam文件方式破解xp开机密码
返校后发现机子密码忘记了.各种方法无果.格机有点代价高昂.试试DOS. winpe启动. 里头有个DOS工具箱. 工具箱里头以后如下的工具: 搜寻sam文件.我的xp系统盘是C盘,第一个盘. 关于这个 ...

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

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

输入文件

目标

输出文件

代码

如果担心中间过程出错

根据sam文件计算reads的GC含量相关推荐

最新文章

热门文章