Python初体验
2024-08-25 09:37:46
今天开始所有的工作脚本全都从perl转变到python,开发速度明显降低了不少,相信以后随着熟练度提升会好起来。贴一下今天一个工作代码,由于之前去一家小公司测序时,序列长度竟然都没有达到要求,为了之后的索赔事宜,写了个脚本统计所有序列的结果,主要包括总的reads数,bases数,和达到测序策略要求长度的reads数(双端),bases数,高质量(Q30)bases数,高质量reads数(双端)等等......多个文件的统计工作一般会写一个单独处理一个文件的脚本,然后再写一个脚本用来生成多个文件的处理的shell脚本,然后想办法并行处理这个shell就可以,效率会快很多。由于测序数据往往比较大,IO操作时,逐行读取是上策。
统计单个文件测序数据情况脚本:
from __future__ import division
from Bio import SeqIO as fq
import os
import sys
import re
read1_gzfile = sys.argv[1]
read2_gzfile = sys.argv[2]
gz_handle1 = os.popen( 'gunzip -cd %s' % read1_gzfile )
gz_handle2 = os.popen( 'gunzip -cd %s' % read2_gzfile )
basename1 = os.path.basename(read1_gzfile)
basename1 = re.match('(\S+)_R1_001\.fastq\.gz',basename1).group(1)
basename2 = os.path.basename(read2_gzfile)
basename2 = re.match('(\S+)_R2_001\.fastq\.gz',basename2).group(1)
if basename1 != basename2:
raise 'Two Read are not mapped!'
cwd = os.getcwd()
out_handle = open('%s/%s.stat'%(cwd,basename1),'w')
out_handle.write('AllReadsNum\tRead1_PE300_ReadsNum\tRead2_PE300_ReadsNum\tUseful_ReadsNum(Read1>=300 and Read2>=300)\tAll_Bases\tRead1_Q30_Bases(PE300)\tRead2_Q30_Bases(PE300)\tQ30_PE_Reads(Q30>50%)\tUseful_Bases(All)\tUseful_Ratio\n') AllReadsNum = 0
AllBases = 0
Read1_PE300_ReadsNum = 0
Read2_PE300_ReadsNum = 0
Useful_ReadsNum = 0
Read1_Q30_Bases = 0
Read2_Q30_Bases = 0
Q30_PE_Reads = 0
Useful_Bases = 0 def PE300(seq):
if len(seq) >= 300:
return True
else:
return False def Q30(qual_list):
num = 0
for qual in qual_list:
if qual >= 30:
num += 1
return num reads2 = fq.parse(gz_handle2,'fastq')
for read1 in fq.parse(gz_handle1,'fastq'):
read2 = reads2.next()
seq1 = read1.seq
qual1 = read1.letter_annotations['phred_quality']
seq2 = read2.seq
qual2 = read2.letter_annotations['phred_quality']
AllReadsNum += 1
AllBases += len(seq1)
AllBases += len(seq2)
R1_300 = PE300(seq1)
R2_300 = PE300(seq2)
if R1_300 and R2_300:
Useful_ReadsNum +=1
R1_Q30 = Q30(qual1)
R2_Q30 = Q30(qual2)
Read1_Q30_Bases += R1_Q30
Read2_Q30_Bases += R2_Q30
if ( R1_Q30 / len(seq1) >= 0.5 ) and ( R2_Q30 / len(seq2) >= 0.5 ):
Q30_PE_Reads += 1
Useful_Bases += R1_Q30
Useful_Bases += R2_Q30
elif R1_300:
Read1_PE300_ReadsNum += 1
elif R2_300:
Read2_PE300_ReadsNum += 1 Useful_Ratio = Useful_Bases / AllBases
out_handle.write('%i\t%i\t%i\t%i\t%i\t%i\t%i\t%i\t%i\t%f\n'%(AllReadsNum,Read1_PE300_ReadsNum,Read2_PE300_ReadsNum,Useful_ReadsNum,AllBases,Read1_Q30_Bases,Read2_Q30_Bases,Q30_PE_Reads,Useful_Bases,Useful_Ratio))
summary.py
生成脚本:
import os
out = open('summary.sh','w')
cwd = os.getcwd()
with open('templist') as gzfiles:
for gzfile1 in gzfiles:
gzfile2 = gzfiles.next()
out.write('python %s/summary.py %s %s\n'%(cwd,gzfile1.strip(),gzfile2.strip()))
run_summary.py
使用qsub_sge方法,并行投递生成的summary.sh就完成了
最新文章
- Word 宏命令大全
- autoLayout 纯代码
- 测试JdbcTemplate执行SQL语句和存储过程
- BZOJ 2054 疯狂的馒头
- hadoop job执行完的统计信息
- Response.End(); 用HttpContext.Current.ApplicationInstance.CompleteRequest 代替
- html标签及用法小结
- Android进阶(二十一)创建Android虚拟机
- python3 生成器初识 NLP第五条
- SimMechanics/Second Generation倒立摆模型建立及初步仿真学习
- layui---form表单模块
- WebApi上传文件
- 转:C# 对委托的BeginInvoke,EndInvoke 及Control 的BeginInvoke,EndInvoke 的理解
- ef-codefirst方式配置实体类,生成数据库
- LisView控件
- javascript 高级程序设计 十一
- Javascript中怎样获取统一管理的Java提示语
- /etc/grub.conf
- Navicat MySql乱码解决
- 11.15 Daily Scrum
热门文章
- CDH集群搭建部署
- Protocol Buffer序列化对比Java序列化.
- Android Toolbar 标题居中及字体样式自定义
- Python2/3中的urllib库
- struts2摘抄
- Javascript实现简单的下拉二级菜单
- HDU 1043 Eight (BFS&;#183;八数码&;#183;康托展开)
- HDU 4960 Another OCD Patient(记忆化搜索)
- 时光轴三之 ExpandableListView版时光轴效果
- JS中数组的迭代方法和归并方法