拆分fastq文件

在NCBI下载的转录组数据 本来是双端测序数据，但是不知道为啥read1 和 read2是在一个文件里，拆分的话可以使用seqkit这个工具

参考链接

https://bioinf.shenwei.me/seqkit/usage/#grep

seqkit grep -n -r -p 1$ SRR16509471.fastq.gz -o SRR16509471_1.fastq.gz
seqkit grep -n -r -p 2$ SRR16509471.fastq.gz -o SRR16509471_2.fastq.gz

-n 是根据序列的名字来

• r 是使用正则表达式

-p 是正则表达式的内容

1$是指序列的名字最后一个字符是1 同理2$是指最后一个字符是2

处理速度还相对可以

image.png

压缩文件是8个G

但是有一个问题

正常的序列 1和 2如下

image.png image.png

/1 /2 前面的数字read1和read2是一样的

但是在同一个文件里拆分后数字不一样，不知道有没有啥影响

image.png image.png

可以做个比对试一下

如果想要修改，还是使用seqkit

https://bioinf.shenwei.me/seqkit/usage/#replace

zcat SRR16509471_1.fastq.gz | head -n 12 | seqkit replace -p '.\d+ \d+/' -r ".{nr} {nr}/"

image.png image.png

也尝试了自己写python脚本

import gzip
import time
import argparse
import pandas as pd
from multiprocessing import Pool
from Bio.SeqIO.QualityIO import FastqGeneralIterator




def split_fq(fastqIterator):
    read01 = {'seq_id':[],
            'seq':[],
            'third_line':[],
            'quality':[]}
    read02 = {'seq_id':[],
            'seq':[],
            'third_line':[],
            'quality':[]}
    
    i = 0
    
    for a in fastqIterator:
        i += 1
        
        if a[0].endswith('1'):
            read01['seq_id'].append('@' + a[0])
            read01['seq'].append(a[1])
            read01['third_line'].append('+')
            read01['quality'].append(a[2])
            
            if i%10000000 == 0:
                print('[{0}] {1} {2} reads'.format(time.ctime(),str(i/1000),'w'))
        
        
        elif a[0].endswith('2'):
            read02['seq_id'].append('@' + a[0])
            read02['seq'].append(a[1])
            read02['third_line'].append('+')
            read02['quality'].append(a[2])
            
            if i%10000000 == 0:
                print('[{0}] {1} {2} reads'.format(time.ctime(),str(i/1000),'w'))
            
    
    return [read01,read02]
    

def final_run():
    parser = argparse.ArgumentParser(
            formatter_class=argparse.RawDescriptionHelpFormatter,
            description="split fq",
            epilog='''
            @author: MingYan
            '''
        )
        
    parser.add_argument('-i','--input',required=True)
    parser.add_argument('-r1','--output-read1',required=True)
    parser.add_argument('-r2','--output-read2',required=True)
        
    args = parser.parse_args()
        
    in_file = args.input
    #number_threads = args.number_threads

    r1 = args.output_read1
    r2 = args.output_read2

    print('[{0}] Ready!'.format(time.ctime()))

    obj_fastq = FastqGeneralIterator(gzip.open(in_file,'rt'))
    
    reads = split_fq(obj_fastq)

    pd.DataFrame.from_dict(reads[0]).to_csv(r1,sep="\n",index=False,header=False)
    pd.DataFrame.from_dict(reads[1]).to_csv(r2,sep="\n",index=False,header=False)

    print('[{0}] Congratulations!'.format(time.ctime()))
    
    
if __name__ == '__main__':
    final_run()

首先是处理速度非常慢，然后需要非常大的内存空间，代码功底还是差很多