通过 blast 结果查看 测序数据fastq是否被污染，以及污染reads所属物种、所占比例

由于NGS测序时混库的原因，我们得到的数据有被其他物种污染的可能，所以当我们拿到测序数据时，如何确定我们的数据是否被污染呢？
下面我们通过blast的结果，来检查测序数据fastq 是否被污染，以及污染来源、所占比例。

1.下载blast

版本地址：ftp://ftp.ncbi.nlm.nih.gov/blast/executables/LATEST/
blast使用手册《BLAST Command Line Applications User Manual》

wget ftp://ftp.ncbi.nlm.nih.gov/blast/executables/blast+/LATEST/ncbi-blast-2.10.0+-x64-linux.tar.gz

下载完成后解压

2.下载blast 所需数据库

下载blast的数据库，nt为核酸数据库，nr为蛋白质数据库

2.1 在NCBI官网下载

进入NCBI官网：https://www.ncbi.nlm.nih.gov/



这里这些 nt.xx.tar.gz 就是需要下载的 nt 库啦

2.2 直接使用blast的脚本下载

https://www.ncbi.nlm.nih.gov/books/NBK537770/

update_blastdb.pl --decompress nr [*]
update_blastdb.pl --decompress nt [*]

不过第二种方式下载的比较慢，因为他是串联下载的，还容易断掉，所以个人感觉还是用第一种方法直接wget的好，只需要把 nt.xx.tar.gz 中间的xx编号更换，就可以批量下载了。

wget ftp://ftp.ncbi.nlm.nih.gov/blast/db/nt.00.tar.gz

下载完成后，将 nt.xx.tar.gz 解压即可

tar zxvf nt.xx.tar.gz

3. NCBI Taxonomy 数据库下载

从Taxonomy :数据库中，我们可以知道所有生物的分类和命名。

wget -c https://ftp.ncbi.nlm.nih.gov/pub/taxonomy/accession2taxid/nucl_gb.accession2taxid.gz
wget -c https://ftp.ncbi.nlm.nih.gov/pub/taxonomy/taxdump.tar.gz

下载完成后解压。

nucl_gb.accession2taxid 数据库格式如下：

第一列Accession ： 序列标识码
第一列Accession.version : 带版本号的序列标识码
第三列： 序列的taxid 号，即物种分类号。如 Homo sapiens 的是9606.
第四列：序列的gi号

这个文件比较大，我们后边只用三四列信息，所以可以把这个库处理一下，只留下第3,4列。

cut -f 3,4 nucl_gb.accession2taxid > cutted_nucl_gb.accession2taxid

taxdump.tar.gz 解压后会有7个库，我们用 names.dmp ，其中包含物种的taxid号与物种学名。对应的格式如下：

第1列为 物种的taxid号。
第2列为物种名称。
我们后面会选择scientific name 对应的物种学名。

ok,到此，我们所需的所有数据库（nt ，nucl_gb.accession2taxid ，names.dmp ）都准备好了。下面需要得到blast结果~

4 fastq文件前处理

从fastq中抽取序列,保存为fasta格式。我们这里从fq文件中提取20000行，也就是5000条reads。

zcat myfile.fastq.gz | head -n 20000 | awk '{if(NR%4==1){print ">"$1}else if(NR%4==2){print $0}}'|sed 's/@//g' > myfile.fa

5 对抽取序列进行blast

./blastn -query myfile.fa -out out.xml -max_target_seqs 1 -outfmt 5 -db database_dir/nt -num_threads 2 -evalue 1e-5

需要注意的是-db 不能只写到存放nt库的目录这一级，后面需要加上库的类型，这里使用的nt库，所以加上nt。 输出结果文件类型选择的 -outfmt 5 ， 也就是XML格式的结果文件，因为这种结果信息比较全。XML文件的格式内容有时间再补充。

6 从XML结果文件提取gi号

提取的信息包括：

Iteration_query-def：reads id
Hit_id : 匹配序列的 id ,信息中包括gi号
Hit_def ： 匹配序列物种信息

# -*- coding:utf-8 -*-
import re
from collections import defaultdict

xmlfile=open("blast_result.xml","r")
outfile=open("tiqu_gi.txt","w")

dict1=defaultdict(list)
for lines in xmlfile:
	line=lines.strip()
	read_id = re.match('<Iteration_query-def>.*</Iteration_query-def>',line)
	Hit_id = re.match('<Hit_id>.*</Hit_id>',line)
	Hit_def = re.match('<Hit_def>.*</Hit_def>',line)
	
	if read_id !=None:
		read_id=read_id.group()
		read_id = read_id.split("<")[1].split(">")[1]
		key=read_id

	elif Hit_id !=None:
		Hit_id = Hit_id.group()
		Hit_id = Hit_id.split("<")[1].split(">")[1]
		dict1[key].append(Hit_id)

	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")

这样得到的处理文件如下：
共3列，TAB分割，第一列reads id,；第二例 gi号信息；第三列 物种信息。

此时的物种信息列，字段是不整齐的，如Homo sapiens ，虽然都是Homo sapiens，但是字段很不一致，不利于统计，所以需要进行学名统一。
逻辑就是从blast结果中得到gi号，通过gi号得到taxid ,通过taxid 得到物种学名。

# -*- coding:utf-8 -*-

tiqu_gi = open("tiqu_gi.txt","r")
gi2taxid = open("cutted_nucl_gb.accession2taxid","r")
taxid2name = open("names.dmp","r")
get_name = open("scientific_name.txt","w")

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
		
tiqu_dict=defaultdict(list)
for lines in tiqu_gi:
	line = lines.strip().split("\t")
	gi = line[1].split("|")[1]
	tiqu_dict[gi].append("\t".join(line))


gi_taxid_dict={}
for lines in gi2taxid:
	line = lines.strip().split("\t")
	GI = line[1]
	taxid = line[0]
	gi_taxid_dict[GI]=taxid

jiaoji=set(tiqu_dict.keys())&set(gi_taxid_dict.keys())

tax_list=taxid_name_dict.keys()

#tiqu_gi = open("result.txt","r")
tiqu_gi = open("tiqu_gi.txt","r")
for lines in tiqu_gi:
	line = lines.strip().split("\t")
	gi = line[1].split("|")[1]
	if gi in jiaoji:
		taxid=gi_taxid_dict[gi]
		if taxid in tax_list:
			get_name.write("\t".join(line)+"\t"+taxid_name_dict[taxid]+"\n")

结果文件共4列，最后一列为匹配得到的物种学名。

下面在进行各个物种reads在所有抽取reads中所占比例即可。

# -*- coding:utf-8 -*-
from collections import Counter

scientific_name=open("scientific_name.txt","r")
final_result =open("final_result.txt","w")

name_list_all=[]
for lines in scientific_name:
	line = lines.strip().split("\t")
	name = line[-1]
	name_list_all.append(name)

count_result = Counter(name_list_all)
count_list = count_result.items()
count_list.sort(key=lambda x:x[1],reverse=True)

final_result.write("Name\tHit_reads\tpercentage1\tpercentage2\n")
for i in count_list:
	name = i[0]
	number = i[1]
	reads_num = 5000
	percentage1 = "%.2f%%"%(100*float(number)/float(reads_num))
	percentage2 ="%.2f%%"%(100*float(number)/float(len(name_list_all)))
	final_result.write(name+"\t"+str(number)+"\t"+str(percentage1)+"\t"+str(percentage2)+"\n")

结果文件如下，排在前三个的 Hit 分别是 Homo sapiens ； eukaryotic synthetic construct ； Gorilla gorilla gorilla . 以及Hit在5000条reads所占比例（percentage1）, 在所有Hit中所占比例（percentage2）。

这样我们就可以知道我们数据匹配物种情况啦~
限于个人编程水平，脚本肯定有更好的方法，如有错误还请拍砖~