如何处理NGS数据中的污染？

---

本次文章和大家讨论一个大家可能胡遇到很常见的一个问题，在测序中我们很难避免引入一些微生物污染或者人类的污染，例如，我想测序拟南芥，其中由于实验员的操作不够干净，很容易引入一些人类的DNA，又或者该拟南芥的叶子上也混杂着细菌真菌等其他的DNA。当你进行组装或者做一些后续分析，这些污染会造成一些可想不到的不良影响。

这篇文章会和大家讨论一下，处理NGS数据中的微生物污染？

清理原始数据

显而易见的，在一开始解决问题是最好的，最精确的。一般来说在组装和任何其他下游分析之前识别和移除那些污染的reads，可以很好的消除组装出来的基因组含有污染物序列。

在这里给大家推荐一款比较好用的工具DeconSeq (http://deconseq./)，它通过尝试将.fasta或.fastq文件与其人类，微生物，病毒和其内部数据库进行比对，自动检测并清除基因组和宏基因组数据集中的序列污染物。

DeconSeq提供可视化交互式Web界面，容易操作上手（限制为600MB文件大小限制）运行，也可以在本地下载和安装（但这会占用了一定的空间）。

要是你比较清楚你的污染物的来源，你也可以自己建立一个.fasta 文件来含括污染物的集合，然后通过简单的比对例如bowtie2 对含有污染物的reads进行识别。

首先，您需要对包含所有污染物序列的文件进行index（contaminants.fasta）：

bowtie-build contaminants.fasta IDX/contaminants_idx


然后，您可以使用--un的flag将无法与此集合比对的样本（raw1.fq）中的所有原始reads写入新的无污染的fastq文件（nocontaminants.fq）：

#这里给出一个single end read的例子，pair end 也是可以使用类似操作的。
bowtie --un nocontaminants.fq IDX/contaminants_idx raw1.fq

接着给大家安利一个比较好用的脚本去自动化这个流程，其过程包括相同的污染物去除方法以及adapter的去除和其他质量控制措施。这个脚本叫2-ScrubReads，可在github中查阅 (https://github.com/CGRL-QB3-UCBerkeley/denovoTargetCapturePopGen/blob/master/2-ScrubReads),

下面以去除E.coli为例，简单介绍该工具的使用：

2-ScrubReads -f ~/pre-clean/ -o ~/cleaned_data/ -a Adapters.fasta -b libInfo.txt -t /path/to/trimmomatic-0.32.jar -c e_coli_K12.fasta -e 200 -m 15 -z

#-f 存储初始数据的文件夹
#-o 输出的文件夹
#-a 剪除adapters的序列
#-b 文库的信息
#-t trimmomatic的路径
#-c 污染序列文件
#-e read lib的平均长度
#-m reads开始时的碱基数（nt）以推断重复序列
#-z 是否启用Fastqc去做质控


组装好后再去除污染物

在原始数据中去除污染物，是最好的策略，但在某些情况下（例如你想把一些比对不上的序列进行组装），你希望在组装后识别并清除污这些染物的contigs。你可以将组装好的cotnigs.fasta文件用blast进行去污染，因为blastn会搜索您感兴趣的样本（sample1），以查询来自潜在污染序列的.fasta的自定义blastn数据库（污染物）（contaminants.fasta） ）：

#建数据库
makeblastdb -in contaminants.fasta -parse_seqids -dbtype nucl
#进行blastn搜索
blastn -db contaminants -query sample1.fasta -outfmt 10 -num_alignments 5 -max_hsps 1 -out blast-output.m6

更加好的替代方法是对整个NCBI GenBank核苷酸数据库（nt）的本地副本进行搜索，这样你可以确保找出所有污染物(因为你所选的或者制作的contamination.fasta会不够完整：

blastn -db nt -query sample1.fasta -outfmt 10 -num_alignments 5 -max_hsps 1 -out blast-output.m6


这会产生一个具有best hits信息的文件，然后通过该hits的名称，或者其特定的GI序列标识符，你可以对污染物进行过滤。比如你想过滤掉含有非脊椎动物的contigs，你需要写上一个脚本，对不是脊椎动物的GI进行过滤。这里有一个这样的脚本的例子(https://github.com/calacademy-research/GItaxidIsVert)。当你过滤掉污染物，筛选出你想要的contigs的ID，你需要进一步将其从原来得含有污染物的cotnigs中取出，一般可以使用工具seqtk subseq (https://github.com/lh3/seqtk)。

从污染物的cotngis中去除SNP

最后，在某些情况下，你可能会发现自己需要从SNP矩阵中去除污染物位点。如果你的SNP matrix是通过将原始数据与可靠的参考基因组比对而生成的，该基因组本身没有污染物，那么应该没有问题。但是，如果你将reads与从一个或几个样本的从头组装生成的“伪参考序列”进行比对，这在使用非模型生物体时很常见，那么reads总是有可能与此伪造的含有污染物的参考序列比对上。

一般来说通过对深度进行截止过滤和其他过滤步骤将消除这些和污染物比对上的位点。但是，如果你想彻底的减少这周含有污染物伪参考序列的影响，那么最好将这些contigs留在它们比对的期间，以减少错配的可能性。

例子，
你可以使用samtools / bioawk，删除一组.bam文件中含有的可疑序列中的所有位点（C7961234; C7963448; C8091874）。

samtools view -h sorted-piledup.bam | awk '{if($3 != "C7961234" && $3 != "C7963448" && $3 != "C8091874" && $2 != "SN:C7961234" && $2 != "SN:C7963448" && $2 != "SN:C8091874"){print $0}}' | samtools view -Sb - > output-NoContamNoMito.bam

或者在一个vcf文件中，使用vcftools 然后使用 --not-chr 这个option：

vcftools --vcf input.vcf --not-chr C7961234 --not-chr C7963448 --not-chr C8091874 --recode --recode-INFO-all --out nocontamination


小结

这是我总结了几个比较常用对NGS数据去污染的方法，当然我相信很多公司或者研究组都会有他们自己的流程去进行污染过滤。欢迎在评论中补充其他的方法。

---

还有更多文章，请移步公众号阅读