统计fastq 中GC含量、Q20、Q30

话不多说，直接上代码

import sys
from Bio import SeqIO
import gzip

fq1,fq2 = sys.argv[1:]  #输入文件

### 判断输入文件为fastq还是fastq.gz
def trans2fastq(fin) :
      fout = fin.replace('.gz','')
      g = gzip.GzipFile(mode='rb', fileobj=open(fin,'rb'))
      open(fout,'wb').write(g.read())
      return(fout)

### 将fastq.gz 转为fastq
if '.gz' in fq1 : fq1 = trans2fastq(fq1)
if '.gz' in fq2 : fq2 = trans2fastq(fq2)

#### 计算
for f in [fq1, fq2] :
        reads_count = 0
        gc_base = 0
        all_base = 0
        all_qual = 0
        q20_qual = 0
        q30_qual = 0
        record = SeqIO.parse(f,'fastq')
        for lane in record :
                ## 序列数
                reads_count+=1
                ### 碱基统计
                gc_base+=len([i for i in lane.seq if i=='G' or i=='C'])
                all_base+=(len(lane.seq))
                ### 质量值统计
                qual = lane.letter_annotations['phred_quality']
                all_qual+=(len(qual))
                q20_qual+=len([q for q in qual if q >=20])
                q30_qual+=len([q for q in qual if q >=30])
        ### 计算比例
        q20 = round(q20_qual/all_qual,4) * 100
        q30 = round(q30_qual/all_qual,4) * 100
        gc = round(gc_base/all_base,4) * 100
        print(reads_count,all_base,gc,q20,q30)