如何处理组装中的污染序列
转自：https://mp.weixin.qq.com/s/8t7rNrp82tnoTW5LN_ck3g
通常情况下，我们都不需要考虑组装中的污染问题。如果我们的样本经过一定清洗，那么测序过程中就很少会测到目标样本无关的序列，因此在组装过程中，这些序列由于信号过低将不会被组装出来，也就不会出现在最终结果中。
但是这一次，我测的一个样本的基因组足足比我预期的大了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
然后将/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
这里使用--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
最终的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版本为1.10）
后续可以使用awk进行过滤，然后和kraken2的分析结果进行比较，或者使用BLASTN进行验证。
这是我目前想到的两种方法，我也自己查了查，发现似乎还有软件专门做这种事情，但是我感觉好像都不符合我的需求。
精选留言
- 试过一次从NCBI下载某植物的genome survey数据进行组装注释，在浏览注释出来的基因时发现了昆虫的基因，多看了几个类似的情况之后差不多可以确定是来自烟粉虱。猜想是叶片上有烟粉虱若虫而取样者未发现。因为数据量不大，我下载了烟粉虱基因组，将reads比对之后过滤，过滤掉大约5% reads，然后用干净数据重新进行了组装。有趣的是同一个组上传的RNAseq数据并没有类似污染。
回复：可能是因为昆虫的RNA抽提比较困难，跟植物的相比
2.前段时间组装一个五十多m的真菌，最后组装出来七十m，不过好在当时测了Hic，用hic数据绘制了热图，发现大概十几m的序列在热图上是一片空白，怀疑是污染序列就去掉了
作者回复： 估计是HiC和三代是独立的两次实验，所以污染源不同，也相当于过滤了污染。
3.我发现可以利用下ANI，比如FastANI过滤？先利用Diamond 比对到Nr库，找出所有比对到奇怪的物种的序列，用FastANI计算是哪种物种序列。
4.同一个物种基因组，没有组装到染色体上的序列，GC含量差别太大，有的50%，染色体34%，但是Blastn就是一个物种的序列。我可以很明确，GC不行
