fq测序数据损坏怎么办

崔耳又又

前两天一个师兄和我说，他有一个转录组数据有问题，做质控的时候就卡住了。当时猜测可能是复制数据的过程中出现了问题。

当时我百度了一下如何修复测序数据，参考曾老师的方法（https://zhuanlan.zhihu.com/p/603304179）对数据进行了修复。我这里也把我的修复过程还原一下，整体思路是和曾老师的方法是一致的。

1 将损坏数据中的未损坏部分提取出来

这一步其实非常简单，直接将gz格式的fq数据解压，解压成功的部分即是未损坏的部分。注意在解压时候如果直接用 <span class="ne-text">gunzip sample_1.fq.gz</span>，会出现以下错误。

所以使用如下代码进行解压：

zcat sample_1.fq.gz > sample_1.fq

此时也会报错，但是没关系，未损坏的部分已经提取出来了。

一般来说末尾还残存一些损坏的数据，我先来简单说一下完整的fq数据格式是什么样子的，方便大家将损坏的部分分辨出来。

fq数据中4行为1个单位，第一行是@开头的序列ID，第二行是序列reads，第三行为“+”，第四行是对应的测序质量ASCII码（质控主要就是用软件对第四行的信息进行分析）。

下图中是正确的fq文件格式，为了让大家易于理解，每个单位之间用红色虚线进行了分隔。

我们输入 <span class="ne-text">tail -500 sample_1.fq|less -SN</span>，发现在后500行的第285行开始为不完整的fq格式，即整个文件的后“500-284=216（行）”为不完整的数据。

我们通过以下代码将文件的后216行删除。

sed -i  ':a;$d;N;2,216ba;P;D' sample_1.fq  #216可替换为您需要的行数

输入 <span class="ne-text">tail -12 sample_.fq</span>检查一下，发现已经将损害的部分剔除。

2 将双端数据统一

同一个样品的双端测序数据的两个fq文件里面需要保证reads的一致性，所以第二步我们需要将双端数据统一。首先将损坏数据中未损坏数据的ID提取出来，我这边是“_1”端的数据损坏，大家如果是“_2”端数据损坏自行灵活替换。输入以下代码将ID提取到“id.lis”中：

grep '@' sample_1.fq > id.lis

然后根据我自己写的python脚本将“_2”端的fq数据进行提取

python getfq2.py --fq sample_2.fq

python脚本如下，写的比较简单，各位大佬勿喷。

#coding:UTF8
import argparse

parser=argparse.ArgumentParser(description='根据ID提取fq数据')
parser.add_argument('--fq', type=str,help='fq文件', required=True)
args=parser.parse_args()
fq=args.fq

save=open(fq+'_new','w')
dic={}
for id in open('id.lis','r'):
    dic[id.split('/')[0]]=id.split('/')[0]
fq=open(fq,'r').readlines()
for num in range(0,len(fq),4):
    if fq[num].split('/')[0] in dic:
        save.write(fq[num]+fq[num+1]+fq[num+2]+fq[num+3])
print ('succeed')

现在我们就拥有了一对新的数据，但是我们可以发现拯救出来的只有原始大小的二十分之一，只能够勉勉强强用，所以大家平时还是要提前做好备份，做好md5检测文件，以备不时之需。

微信扫一扫分享文章

微生信

保存数据最关键。这种情况下，只能问公司要，或者从bam提取也好的。

崔耳又又

微生信 发表于 2024-8-13 12:00
保存数据最关键。这种情况下，只能问公司要，或者从bam提取也好的。

确实确实，主要是数据太久远了