python脚本截取fastq序列

前言

最近在学习单细胞方面的知识，遇到了一个小的需求就是截取需要的fastq序列。先来说一下为什么有这个需求，一般来说单细胞的测序文件有三个如test_R1.fastq.gz、test_R2.fastq.gz、test_I1.fastq.gz。其中R1、R2就是测序的read1、read2，I1则是样本的index文件(基本上用不到)。这里需要提醒的是一般单细胞的read1是cell barcode和umi序列文件(10x的barcode一般为16nt，umi为10nt)，read2文件中的序列才是真正的基因表达量的reads。虽然测序得到双端150bp的fastq，但read1中有很大一部分是无效的序列。为了只保留有效的序列就要去原始序列中截取有效的序列丢掉无效的数据，那么如何解决这个问题？一时半会没有想到有什么现成的工具可以完成，于是乎我就自己写了一个python脚本来完成这个需求。

结果

下面来看看我到底做了啥，比如我现在有如下的fastq文件：

@SRX8108993.1 1 length=151
NTCTCCTCACCGCTAGGGTCACGCCTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTCCCATTTTAATCTTTTTTTTCCTTCTCTAAGTCCCTTTTCTTCTGCTTTCCGATTATCACATTTTTCTTTCCCATGTTTATGTTCTCTTACCCTTC
+
#AA-FFFJ7JJFJ<7<-AFJFJAFJ-JJJJJJJJJJJJJJJJJJJJJJJJJJAFA-<-------77AA-<77------7<--777--7-7)<-<---7--7-7------7-77--7---<A--7-7----)-----7---7<------7-)
@SRX8108993.2 2 length=151
NACACAACATCTCCCATCCTTGTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTGTTTTTTTTTTTTTATTTTTTATTTTTTTTTTATAATCATATCAGATTTTTATTAGAAAGTGACCTCACACAATTTTCCAATCCTCTTTACCTAC
+
#AAAFA-F-7A-JJJFFJFFJFFJFJJJJJJJJJJJJJJJJJJJJJJJJJJJJJF<---7-A<--AF7<A-7A-<7--77<7<-<<FJ7-----77--7-7-77------7---------)-)-)-7------77-----7<---------

现在需要截取5'端的26个碱基，生成新的截短的fastq文件：

@SRX8108993.1 1 length=151
NTCTCCTCACCGCTAGGGTCACGCCT
+
#AA-FFFJ7JJFJ<7<-AFJFJAFJ-
@SRX8108993.2 2 length=151
NACACAACATCTCCCATCCTTGTTTT
+
#AAAFA-F-7A-JJJFFJFFJFFJFJ

其实完成这个任务很简单，用Linux的awk命令就可以解决，像这样80多个G的单细胞数据，其实用awk来干这件事估计七八个小时的时间也能处理完。但是如果用python写一个处理的小脚本，这样以后就可以重复利用了，还是挺好的。于是就是写了一段代码，内容如下：

#!/usr/bin/env python3
import gzip
import itertools
import sys,os
from io import BufferedReader, TextIOWrapper

#read fastq to records
def read_fastq(fastq,length):
    if fastq.endswith('gz'):
        fastq_file = TextIOWrapper(BufferedReader(gzip.open(fastq, mode='rb')))
    else:
        fastq_file = open(fastq, mode='rt')
    while True:
        read=itertools.islice(fastq_file, 4)
        head=next(read)
        length=length if length else head[head.rfind('='):]
        element = ''.join([head,next(read)[0:int(length)]+os.linesep,next(read),next(read)[0:int(length)]+os.linesep])
        if element is not '':
            yield element
        else:
            break
    fastq_file.close()
    return element

#write reads to gz file
def write_fastq(reads,out):
    out = gzip.open(out+'.fastq.gz','wb')
    for read in reads:
        out.write(read.encode())

if __name__ == "__main__":
    if len(sys.argv)!=4:
        print('Error:\n\tUsage: python3 %s fastq length outprefix' % sys.argv[0])
        sys.exit(1)
    reads=read_fastq(sys.argv[1],sys.argv[2])
    write_fastq(reads,sys.argv[3])

脚本使用方法如下：

python3 fastq_trim.py SRX8108993_1.fastq.gz 26 trim_SRX8108993

代码看起来还是听简单的，脚本使用起来也是相当的简单，只需要三个命令行参数分别为原始的fastq文件(gz压缩和非压缩的格式都可以)、需要截取的长度、新的fastq文件名的前缀(如示例的代码会生成截取后的名为trim_SRX8108993.fastq.gz压缩格式的新文件)。

最后

emm，今天就分享到这里，有需要的朋友可以拷贝代码使用，有什么不对的地方可以留言告诉我一下。