唠唠叨叨：
因为处理ATAC-seq的原始fastq数据太大（已知GEO号在ncbi上下载），所以被老师要求用python编程分割fastq脚本，其实用软件seqkit就可以进行分割，命令简单还快速，想分多少份就分多少份，但老师说我不练怎么会编程呢，所以没办法，我这刚入门的小白，就在网上各种查资料，发现都是处理fasta格式的，因为fastq格式和fasta格式并不一样，不知道可以借鉴不，日夜研究别人的代码，以及翻找网络，读到一篇将fastq文件写入和读出的代码，理解并做出了尝试，然后又绞尽脑汁如何进行分割，后来想到按行分配，只是需要每次都查看一下文件有多少行，然后又搜了一下代码，想把两者结合起来，中间报了各种error,就反复调试，7.1的那一天上午成功了，但是我取的是10000行，不知道生成的文件为什么是7000多行，没找到原因，没有打开生成的文件看看，下午老师过来看我发现是没有加换行符的问题，实际上输出是对的
分割fastq脚本：
需要引入时间函数（此处参考网上），可以看脚本分割用了多久
gz文件打开需要引入gzip函数
from datetime import datetime
import re
import sys
def Main():#定义函数，读入fastq文件
source_dir = sys.argv[1] #传入的参数1：目标文件夹所在目录
read_number=int(sys.argv[2]) #传入的参数2：需要分割的行数，需要是4的倍数，脚本里需要验证一下
target_dir = sys.argv[3]#传入的参数3：分割后文件所在的目录
#print(target_dir)
fas={}
id=None
name = 1
N=0 #读取数据的行数数值，初始值为0
if read_number%4 != 0:
print("参数错误，程序退出")
sys.exit()
else:
print("开始。。。。。")
print(datetime.now().strftime('%Y-%m-%d %H:%M:%S'))
with gzip.open(source_dir,'rb') as f:
for line in f:
N+=1
if line.decode()[0]=="@":
id=line.decode().strip().split()[0][1:]#对二进制文件进行解码
fas[id]=[] #将该id对应的初始值为空
else:
seq=fas[id].append(line.decode())
#print(fas)
if N==read_number:
# with gzip.open(target_dir+'/57.%04d.fq.gz'%(name),'wb') as out:
# print(target_dir+'/57.%04d.fq.gz'%(name))
if re.findall('r1', source_dir, flags=re.IGNORECASE): #正则匹配可以不区分大小写 或直接判断R1是否在source_dir里，需大写：if R in source_dir:
with gzip.open(target_dir+'.R1.%04d.fq.gz'%(name),'wb') as out:
print(target_dir+'.R1.%04d.fq.gz'%(name))
for id,seq in fas.items():
res = '@'+id+"/1"+"\n"+seq[0]+seq[1]+seq[2]
out.write(res.encode())
name+=1
N=0
fas={}
else:
with gzip.open(target_dir+'.R2.%04d.fq.gz'%(name),'wb') as out:
print(target_dir+'.R2.%04d.fq.gz'%(name))
for id,seq in fas.items():
res = '@'+id+"/2"+"\n"+seq[0]+seq[1]+seq[2]
out.write(res.encode())
name+=1
N=0
fas={}
else:
continue
print("完成。。。。。")
print(datetime.now().strftime('%Y-%m-%d %H:%M:%S'))
if __name__ == "__main__":
Main()
优化：此脚本由我的导师提出几个优化步骤又进行改进，以上是完整版，改进为：一传入参数；二对传入的文件名进行判断是含R1还是R2，如果含R1在文件的read_name后加上/1，反之R2加上/2；三将输出的文件名的数固定住
对于我的脚本有一些缺陷：如果文件生于的行数怎么办，虽然此脚本有考虑到，但我没运行成功，二就是分割的时间太慢
7月的第一天，我一个小白能写出脚本感觉异常欣慰，二是班长通知六级报名，而我去年六级竟然裸考过了，终于不用再去奔赴六级的考场了，哈哈，开心，感谢我读过的英文文献
老师版：
第二天我们老师在我旁边又教我写出另一个版本，即按read数进行分割fastq文件，不得不说我们老师真厉害，代码比我简洁多了,而且教我写出只用了一个小时，时间上也会慢些，应该是文件较大，但是不用考虑剩余行的问题，只需要提前计算好分割的reads数
import gzip
from datetime import datetime
import re
import sys
def Main():
test_R1=False
test_R2=False
source_dir = sys.argv[1]
read_number=int(sys.argv[2])
target_dir = sys.argv[3]
#print(target_dir)
fas={}
id=None
name = 1
N=0 #读取数据的行数数值，初始值为0
print("开始。。。。。")
print(datetime.now().strftime('%Y-%m-%d %H:%M:%S'))
with gzip.open(source_dir,'rb') as f:
for line in f:
if N%(4*read_number)==0:
if re.findall('r1', source_dir, flags=re.IGNORECASE):
f_out = gzip.open(target_dir+'.R1.%04d.fq.gz'%(name),'wb')
print(target_dir+'.R1.%04d.fq.gz'%(name))
test_R1 = True
suffix="/1"
elif re.findall('r2', source_dir, flags=re.IGNORECASE):
f_out = gzip.open(target_dir+'.R2.%04d.fq.gz'%(name),'wb')
print(target_dir+'.R2.%04d.fq.gz'%(name))
test_R2= True
suffix="/2"
else:
print("Error: R1 or R2 in filename not detected，exit")
sys.exit()
name+=1
N+=1
if N%4==1:
if test_R1:
if re.findall('/1', line.decode(), flags=re.IGNORECASE):
f_out.write(line)
else:
f_out.write((line.decode().strip().split()[0]+suffix+"\n").encode())
if test_R2:
if re.findall('/2', line.decode(), flags=re.IGNORECASE):
f_out.write(line)
else:
f_out.write((line.decode().strip().split()[0]+suffix+"\n").encode())
elif N%4==3:
f_out.write(("+\n").encode())
else:
f_out.write(line)
print("完成。。。。。")
print(datetime.now().strftime('%Y-%m-%d %H:%M:%S'))
if __name__ == "__main__":
Main()
这两段代码都已运行成功
网友评论