前言

??最近在學(xué)習(xí)單細(xì)胞方面的知識，遇到了一個小的需求就是截取需要的fastq序列。先來說一下為什么有這個需求，一般來說單細(xì)胞的測序文件有三個如test_R1.fastq.gz、test_R2.fastq.gz、test_I1.fastq.gz。其中R1、R2就是測序的read1、read2，I1則是樣本的index文件(基本上用不到)。這里需要提醒的是一般單細(xì)胞的read1是cell barcode和umi序列文件(10x的barcode一般為16nt，umi為10nt)，read2文件中的序列才是真正的基因表達(dá)量的reads。雖然測序得到雙端150bp的fastq，但read1中有很大一部分是無效的序列。為了只保留有效的序列就要去原始序列中截取有效的序列丟掉無效的數(shù)據(jù)，那么如何解決這個問題？一時半會沒有想到有什么現(xiàn)成的工具可以完成，于是乎我就自己寫了一個python腳本來完成這個需求。

結(jié)果

下面來看看我到底做了啥，比如我現(xiàn)在有如下的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<---------

現(xiàn)在需要截取5'端的26個堿基，生成新的截短的fastq文件：

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

??其實(shí)完成這個任務(wù)很簡單，用Linux的awk命令就可以解決，像這樣80多個G的單細(xì)胞數(shù)據(jù)，其實(shí)用awk來干這件事估計(jì)七八個小時的時間也能處理完。但是如果用python寫一個處理的小腳本，這樣以后就可以重復(fù)利用了，還是挺好的。于是就是寫了一段代碼，內(nèi)容如下：

#!/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

??代碼看起來還是聽簡單的，腳本使用起來也是相當(dāng)?shù)暮唵危恍枰齻€命令行參數(shù)分別為原始的fastq文件(gz壓縮和非壓縮的格式都可以)、需要截取的長度、新的fastq文件名的前綴(如示例的代碼會生成截取后的名為trim_SRX8108993.fastq.gz壓縮格式的新文件)。

最后

?emm，今天就分享到這里，有需要的朋友可以拷貝代碼使用，有什么不對的地方可以留言告訴我一下。