0基础学习基因组三代测序组装（2）——数据污染评估

做完基因组三代测序数据质控之后，我们把所有reads的Q值控制在7以上，每个read的长度在1000bp以上。我们不能明确自己的测序数据是否被其他物种污染，这个时候就要用balst比对的方法确定测序数据是否被污染，以及污染的来源。

1 下载balst+工具和数据库

在之前的一篇博客中，我详细介绍了如何本地安装NCBI的blast+工具，以及下载nr/nt库，建立本地的数据库。详情点击这里。

在做数据污染评估的时候，我们还需要知道blast最佳结果对应的物种名，因此还需要下载分类数据库的以下两个子库：

# ascp工具下载大数据，wget命令下载小文件（md5校验文件）
ascp -QT -i /public/home/wlxie/miniconda3/envs/biosoft/etc/asperaweb_id_dsa.openssh -k1 -l 500m anonftp@ftp.ncbi.nlm.nih.gov:/pub/taxonomy/accession2taxid/nucl_gb.accession2taxid.gz ./

wget https://ftp.ncbi.nlm.nih.gov/pub/taxonomy/accession2taxid/nucl_gb.accession2taxid.gz.md5

ascp -QT -i /public/home/wlxie/miniconda3/envs/biosoft/etc/asperaweb_id_dsa.openssh -k1 -l 500m anonftp@ftp.ncbi.nlm.nih.gov:/pub/taxonomy/taxdump.tar.gz ./

wget https://ftp.ncbi.nlm.nih.gov/pub/taxonomy/taxdump.tar.gz.md5

md5文件校验完成之后，两个数据库分别解压。注意.gz文件用gunzip，.tar.gz文件用tar -zxvf

看看这两个数据库长什么样：

第一张图片是nucl_gb.accession2taxid，我们需要用到第二列版本信息和第三列的taxid。

第二张图片是names.dmp，我们需要用到有taxid，学名和scientific name字符串的行。

这两个数据库怎么使用后面会详细说明。分析思路来自于CSDN的博主风风是超人，遗憾的是从17年开始，NCBI不再提供gi号与blastn结果的关联，博主的本地数据库可能版本比较早，采用的是gi号分析。

我将后续的代码做了修改，下载的也都是最新的数据库。总的逻辑是利用blast结果的version号，得到nucl_gb.accession2taxid数据库中的taxid号，最后通过names.dmp中的taxid号得到学名。代码方面做了少许优化，对集群服务器可能更友好一点？

2 fq文件处理和blast

质控后的数据fq文件是“@”开头的，我们要改成fa格式也就是“>”开头。取前10000条序列，每个序列有4行，只取第一行标题和第二行序列。

# NR表示当前行，判断除以4的余数，余数1为标题行，只输出第一个元素即reads id；余数2则为序列行，输出所有元素也就是整条序列。最后替换@符号，文件名为test.fa
zcat clean_filter.fq.gz | head -n 40000 | awk '{if(NR%4==1){print $1}else if(NR%4==2){print $0}}' | sed 's/@/>/g' >test.fa

# 批量blast程序
blastn -query test.fa -out test_blastn_nt.xml -db nt -outfmt 5 -evalue 1e-5 -num_threads 20 -max_target_seqs 1

批量blast程序注意下我们输出的格式为xml格式，也就是- outfmt 5。为什么要用xml格式，因为xml格式能给出的信息最全，我们需要知道输出的版本号

evalue值根据需要设定，这里我设置1e-5

最大匹配数量注意下设置1，我们只需要知道和哪个物种相似度最高，一个输出结果就足够了（虽然设置1会有警告）。

看下blast生成的test_blastn_nt.xml这个结果文件：

虽然是第一次接触xml格式，但是感觉非常熟悉！之前做的一个微博爬虫小程序就是扒了一个类似的html格式的文件……xml格式也挺容易解析的，可以看到比对信息以标签<Iteration>开始，以</Iteration>标签结束，<Hit>标签开始表示的是比对上的结果（因为我设置了最大比对序列数量是1，所以<Hit_num>只有1）；<Hsp>标签表示某一块的比对结果（同一条序列，若干片段比对上），因此<Hsp_num>标签的数量可能不止一个。

当然，这些都可以不用关心，分析需要的信息我用红框标了出来。比较重要的是<Hit_def>标签，里面的字符串是空格隔开的，第一个元素是我们需要的物种版本号。

3 XML文件解析

前面说了解析的思路，以下是代码的实现。因为用的python语言写的程序，我的建议是在vscode一类的编程软件中写这些代码，如果有错误可以及时调试。

import sys
import re
from collections import defaultdict

# 可以不写，我是为了确保导入父目录的模块不出错
sys.path.append("./")

xmlfile = open("/public/home/wlxie/luobuma/luobuma/sample_1_rawdata/Third_generation_sequencing/test_blastn_nt.xml","r")
outfile = open("/public/home/wlxie/luobuma/luobuma/sample_1_rawdata/Third_generation_sequencing/filted_accession_version.txt","w")

# 定义一个空的字典，提取有queryID和subjectID的行
dict1 = defaultdict(list)
for lines in xmlfile:
    line = lines.strip()
    read_id = re.match('<Iteration_query-def>.*</Iteration_query-def>',line)
    Hit_def = re.match('<Hit_def>.*</Hit_def>',line)

    # 解析queryID
    if read_id != None:     
        read_id = read_id.group()
        read_id = read_id.split("<")[1].split(">")[1]
        key=read_id
    
    # 字典key值和value值赋值
    elif Hit_def !=None:        
        Hit_def = Hit_def.group()
        Hit_def = Hit_def.split("<")[1].split(">")[1]
        dict1[key].append(Hit_def)

# 写入文件，制表符分割
for key in dict1:
    outfile.write(key + "\t" + "\t".join(dict1[key])+"\n")

看下运行结束后解析得到的文件：

一共是两列，第一列是reads的queryID，第二列subjectID就是比对上的序列信息。可以看到第二列可以以空格为分隔符，提取第一个元素也就是物种版本号，后面会说。

为什么要把物种版本号提出来而不直接用这段内容里的物种名呢？因为不同的物种名字段数量和位置不一样，无法用统一的命令直接提取，精确的版本号可以对应唯一一个taxid，从而被精准地注释上物种学名。

4 匹配物种学名

这里需要注意一个问题，blast用的nt库还有物种分类用到的两个数据库，他们的更新时间是不一致的。也就是说，物种版本号不一定能完全匹配上taxid，而taxid也不一定能匹配上学名。

而python语言写的程序，用到字典类型数据的时候，如果没有对应的key值匹配是会报错的，不会继续执行下去。一会儿解释，代码如下：

import sys
import gc

sys.path.append("./")

accession = open("/public/home/wlxie/luobuma/luobuma/sample_1_rawdata/Third_generation_sequencing/filted_accession_version.txt","r")
accession2taxid = open("/public/home/wlxie/luobuma/luobuma/sample_1_rawdata/Third_generation_sequencing/nucl_gb.accession2taxid","r")
taxid2name = open("/public/home/wlxie/luobuma/luobuma/sample_1_rawdata/Third_generation_sequencing/names.dmp","r")
final_res = open("/public/home/wlxie/luobuma/luobuma/sample_1_rawdata/Third_generation_sequencing/final_res.txt","w")

# 从names.dmp提取taxid和学名，匹配有scientific name的行
taxid_name_dict = {}        
for lines in taxid2name:
    if "scientific name" in lines:
        line = lines.strip().split("|")
        taxid = line[0].strip()
        name = line[1].strip()
        taxid_name_dict[taxid] = name

# 从nucl_gb.accession2taxid提取taxid和版本号
accession_taxid_dict = {}       
for lines in accession2taxid:
    line = lines.strip().split("\t")
    TAXID = line[2]
    VERSION = line[1]
    accession_taxid_dict[VERSION] = TAXID

# 添加两个判断条件，版本号匹配不上taxid和taxid匹配不上学名的情况。gc.collect()释放内存。
for lines in accession:
    line = lines.strip().split("\t")
    version = line[1].split()[0]
    if version in accession_taxid_dict:
        taxid = accession_taxid_dict[version]
        if taxid in taxid_name_dict:
            final_res.write("\t".join(line)+"\t"+taxid_name_dict[taxid]+"\n")
        else:
            final_res.write("\t".join(line)+"\t"+"INVALID TAXID"+"\n")
    else:
        final_res.write("\t".join(line)+"\t"+"INVALID ACCESSION VERSION"+"\n")
    gc.collect()

如果不加最后两个判断条件，程序会在报错的那行read序列终止。

通过比较两个输出结果文件行数是否一致来判断匹配是否完全。

两个文件输出结果一致说明匹配完成，为什么这里是9338而不是我们一开始blast的10000条序列呢？那是因为有662条序列balst结果的E值大于1e-5，没有在nt中比对上合适的序列

5 输出物种注释分布结果

到这一步就有很多种处理方法了，可以把结果文件直接用excel打开，统计reads在nt库的分布情况和比对上的物种分布。也可以直接写个python脚本做个数据统计。

统计前我们先检查一下是否存在上一步匹配失败的reads。

1
2
3

# 统计匹配失败的reads
cat final_res.txt | grep "INVALID TAXID"
cat final_res.txt | grep "INVALID ACCESSION VERSION"

提示有8条reads的物种版本号比对不上taxid，且都是Pyrus x bretschneideri这个物种，说明这个物种还未在nucl_gb.accession2taxid这个NCBI官方数据库中更新。在结果文件中将其替换掉。

1 2	# sed命令在原文件进行全局替换 sed -i 's/INVALID ACCESSION VERSION/Pyrus x bretschneideri/g' final_res.txt

修改完成，检查无误后，用以下python脚本统计物种注释分布：

from collections import Counter     # 引入counter模块
import sys

sys.path.append("./")

name_file = open("/public/home/wlxie/luobuma/luobuma/sample_1_rawdata/Third_generation_sequencing/final_res.txt","r")
res_stastics = open("/public/home/wlxie/luobuma/luobuma/sample_1_rawdata/Third_generation_sequencing/stastics.txt","w")

name_list = []
for lines in name_file:
    line = lines.strip().split("\t")
    name = line[-1]     # 取最后一列
    name_list.append(name)

# Counter()函数统计词频
count_result = Counter(name_list)
count_list = list(count_result.items())     # 注意需要创建一个list
count_list.append(('Unmap',662))        # 注意手动添加blast失败的序列条数到list中
count_list.sort(key=lambda x:x[1],reverse=True) # 以第二维数据值，即统计的物种学名出现次数排序

res_stastics.write("Name\tHit_reads\tpercentage\n")
for i in count_list:
    name = i[0]
    number = i[1]
    reads_num = 10000
    percentage ="%.2f%%"%(100*float(number)/float(reads_num))       # 浮点两位小数的百分比
    res_stastics.write(name+"\t"+str(number)+"\t"+str(percentage)+"\n")

需要注意手动添加blast失败的序列条数，方便最后一起统计。打开生成的stastics.txt文件：

这个数据是制表符分割的，可以用excel做一个分布统计表，或者用R做一个柱状图，底下展示结果

可以看到，10000条序列比对结果占比最高的前两条序列（橘黄色的）是细菌的核酸序列，总数达到3437条。992条序列比对上罗布麻，662条序列未比对上nt库。可以认为这个序测序数据被细菌污染，可以和测序公司battle要求重新测一遍了……

6 补充说明

2022/12/4 更新

秉着科学严谨的态度，再更新一些内容查漏补缺。

质控过滤后的reads有183万条，而我只取了前1万条。考虑到测序开头的低质量reads可能会对分析结果产生干扰（比如开头的电信号不稳定），我写了个python脚本对过滤后的数据随机取10000条reads，这样就只有随机误差影响分析结果了。

在第2步fq文件处理部分，为了python调用方便，先解压clean_filter.fq.gz文件

1	gunzip clean_filter.fq.gz

读取解压后的文件需要49G内存，我只能在集群上处理，接着运行如下python脚本

import random       # 调用random模块产生随机数
import linecache    # 调用linecache模块读入指定行
import gc

output_file = open("/public/home/wlxie/luobuma/luobuma/sample_1_rawdata/Third_generation_sequencing/random_test.fa","w")

reads_list = list(range(1,1834926))     # 共有1834925条reads
line = random.sample(reads_list,10000)
for i in line:
    text1 = linecache.getline("/public/home/wlxie/luobuma/luobuma/sample_1_rawdata/Third_generation_sequencing/clean_filter.fq",4*i-3)
    text2 = linecache.getline("/public/home/wlxie/luobuma/luobuma/sample_1_rawdata/Third_generation_sequencing/clean_filter.fq",4*i-2)
    query_id = text1.split()[0]
    query_id_fa = query_id.replace("@", ">")
    output_file.write(query_id_fa + "\n" + text2)
    gc.collect()

处理方式比之前多了几步，我运行了两次脚本，发现两次产生的文件大小都在135M左右，也就是随机取10000条reads产生的文件比取前10000条reads产生的文件大了40M。证明三代测序开头测得序列质量不太行（短序列不一定质量不好，但是质量不好的序列一定是短序列），拿到随机产生的10000条reads做blast，后续步骤都一样。

2023/1/13更新

与测序公司技术人员沟通后，纠正一些误区：

三代测序数据是长读长片段，直接使用长读长序列进行blastn比对，能比对上数据库的概率会大很多（总有些区域能比对上一些同源序列，一条reads也可能比对上不同的物种）。
真核生物基因组含有大量的重复序列，对三代测序数据一般要进行随机打断成250-500bp大小，然后随机取一条进行blastn比对，如果这条序列刚好是重复序列片段，就会出现比对不上的结果（重复区域测序准确度也相对较低）。

因此，以上的方法更适合二代测序数据的污染评估，三代数据污染评估需要在代码上考虑实现随机打断成短片段（250-500bp），再进行blastn比对nt数据库。

总的来说，这种随机取测序read进行blastn来确定数据污染的方法，只是简单粗暴的看一下样本污染情况，如果污染比例低（通常小于1%），说明数据可以使用，还需要结合GC-depth做具体分析（只有一个红色区域说明数据无污染）。blastn比对这部分内容一般不会在文章中做展示。