Python初体验

今天开始所有的工作脚本全都从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	Read1_PE300_ReadsNum	Read2_PE300_ReadsNum	Useful_ReadsNum(Read1>=300 and Read2>=300)	All_Bases	Read1_Q30_Bases(PE300)	Read2_Q30_Bases(PE300)	Q30_PE_Reads(Q30>50%)	Useful_Bases(All)	Useful_Ratio
')

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	%i	%i	%i	%i	%i	%i	%i	%i	%f
'%(AllReadsNum,Read1_PE300_ReadsNum,Read2_PE300_ReadsNum,Useful_ReadsNum,AllBases,Read1_Q30_Bases,Read2_Q30_Bases,Q30_PE_Reads,Useful_Bases,Useful_Ratio))

 生成脚本：

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
'%(cwd,gzfile1.strip(),gzfile2.strip()))

使用qsub_sge方法，并行投递生成的summary.sh就完成了