在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()
首先是处理速度非常慢，然后需要非常大的内存空间，代码功底还是差很多
