处理组装中的污染序列

通常情况下，我们都不需要考虑组装中的污染问题。如果我们的样本经过一定清洗，那么测序过程中就很少会测到目标样本无关的序列，因此在组装过程中，这些序列由于信号过低将不会被组装出来，也就不会出现在最终结果中。

但是这一次，我测的一个样本的基因组足足比我预期的大了50MB，并且这绝不可能是因为杂合导致的冗余，这就迫使我开始去寻找这背后的原因。最后，我发现这可能是因为测序数据存在大量的无关污染，导致我组装出了大量的无关序列。

这一篇文章主要介绍如何找出其中的污染序列。

方法1: 序列分类

最简单粗暴的方法，就是将我们的序列和NCBI上的NT进行BLAST，然后分析每个序列的最佳匹配对应的物种是否符合预期。使用BLAST的速度会比较慢，但是比较可靠，如果追求速度，可以考虑Kraken2。

Kraken2安装要求只要有200GB的可用空间（NT库），运行内存在64G左右，同时它依赖于sed, wget, perl, g++(支持C++11)和ncbi-blast。更为重要的一点事，网络环境比较良好，因为需要从NCBI的FTP服务器下载数据。

安装方法如下

git clone https://github.com/DerrickWood/kraken2.git
cd kraken2
./install_kraken2.sh /opt/biosoft/kraken2

然后将/opt/biosoft/kraken2/加入到环境变量PATH中，保证自己能够调用kraken2,kraken2-build,kraken2-inspect

Kraken2数据库至少包括三个文件

hash.k2d: 一个哈希字典，记录taxo和minimizer的对应关系
opts.k2d: 记录构建数据框用的参数
taxo.k2d: 构建数据库所需要的taxonomy信息

除了这三个文件以外的其他文件都是构建过程中的临时文件，例如library和taxonomy目录，删除他们对后续分析没有影响。

kraken2能够下载"archaea", "bacteria", "plasmid", "viral", "human", "fungi", "plant", "protozoa", "nr", "nt", "UniVec", "UniVec_Core"数据，然后构建数据库。我们可以使用kraken2的标准数据库，包括"archaea", "bacteria","viral", "human"和"UniVec_Core"，也可以按照自己的需求安装所需要的数据库。

由于我们只是分析污染而已，我不再介绍标准库的安装策略，只介绍nt库构建。

# 下载taxonomy
kraken2-build --use-ftp --download-taxonomy --db /data/database/kraken2
# 下载
kraken2-build --use-ftp --threads 20 --download-library nt --db /data/database/kraken2
kraken2-build --use-ftp --threads 20 --download-library UniVec_Core --db /data/database/kraken2
# 构建
kraken2-build --build --threads 20 --db /data/database/kraken2

这里使用--use-ftp，是因为我遇到了rsync: failed to connect to ftp.ncbi.nlm.nih.gov这类报错。

最终对序列的分类就变得异常简单了，只需要调用kraken2设置好数据位置即可。

kraken2 --report report.txt --use-names --db /data/database/kraken2 seqs.fa > classified.txt
# --report: 增加报告信息
# --use-names: 增加拉丁名

最终的classified.txt里输出内容格式如下:

'C/U': 序列是否被分类
序列ID
taxonomy ID
序列长度
序列中每个k-mer的LCA匹配位置，例如562:13 561:4 A:31 0:1 562:3含义如下

第一个13k-mer和#562匹配
第二个4 k-mer和#561偶偶恶奴
第三个31k-mer是模糊碱基
第四个k-mer不在互数据库
第五个3-kmers和#562匹配

方法2: 覆盖度和GC含量

这里的覆盖度(coverage)既包括三代测序回帖的覆盖度，也包括二代重测序的覆盖度。我一开始认为，由于细菌并不是测序的主体，那么它们的覆盖度应该比较低，同时细菌的GC含量应该和我们物种的GC含量存在差异，说不定可以根据覆盖度和GC含量进行作图来鉴定污染。但实际操作的时候发现，这次污染实在是有点严重，之前通过kraken2鉴定的细菌序列对应的覆盖度和主要测序样本的覆盖度并没有明显的差异。

后来我想到了可以原本用于基因组polish的二代测序进行覆盖度分析。三代测序和二代测序是两次独立的建库流程，在三代建库时可能会引入污染A，二代建库可能会引入污染B，两次的污染源应该不同，那么将二代数据回帖到三代组装结果中时，污染A对应基因组就没有回帖上的read，那么覆盖度就接近0，我们就可以将这部分序列过滤掉。

计算覆盖度的代码如下

bwa mem -t 20 ref.fasta R1.fastq.gz R2.fastq.gz | samtools sort -@ 20 > bwa_aln.bam
samtools index -@ 20 bwa_aln.bam
samtools coverage bwa_aln.bam > coverage.txt
samtools coverage -m -A bwa_aln.bam > hist.txt

后续可以使用awk进行过滤，然后和kraken2的分析结果进行比较，或者使用BLASTN进行验证。