美文网首页
处理组装中的污染序列

处理组装中的污染序列

作者: xuzhougeng | 来源:发表于2020-06-11 11:08 被阅读0次

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

但是这一次,我测的一个样本的基因组足足比我预期的大了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里输出内容格式如下:

  1. 'C/U': 序列是否被分类
  2. 序列ID
  3. taxonomy ID
  4. 序列长度
  5. 序列中每个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进行验证。

相关文章

  • 处理组装中的污染序列

    通常情况下,我们都不需要考虑组装中的污染问题。如果我们的样本经过一定清洗,那么测序过程中就很少会测到目标样本无关的...

  • 如何处理组装中的污染序列

    转自:https://mp.weixin.qq.com/s/8t7rNrp82tnoTW5LN_ck3g[http...

  • 评估组装中的污染

  • 宏基因组分析概述

    测序数据预处理——质控:Trimmomatic 测序数据预处理——质控统计:FastQC 序列组装与基因预测——拼...

  • 序列组装

    1.利用fastqc对模拟测序的序列进行质控分析 1.1 使用art-illumina模拟测序,生成高通量数据(a...

  • 序列组装

    序列组装是测序的下游工序,但在日益多样化的基因组学研究中,这不是必须的。尽管如此,进行序列组装无疑能够帮助我们更精...

  • 使用TransDecoder寻找转录本中的编码区

    TransDecoder能够从转录本序列中鉴定候选编码区。这些转录本序列可以来自于Trinity的从头组装,或者来...

  • 使用TransDecoder寻找转录本中的编码区

    TransDecoder能够从转录本序列中鉴定候选编码区。这些转录本序列可以来自于Trinity的从头组装,或者来...

  • 小麦三代转录组融合基因定量

    小麦三代转录组融合基因定量 三代转录组的序列普遍比较长,得到的是完整的全长的转录本,在后期处理过程中无需组装,大...

  • CheckM评估微生物基因组完整度

    测序后组装获得微生物基因组序列完不完整?污染率有多少? checkM是用于评估单菌基因组、单细胞和宏基因组的质量工...

网友评论

      本文标题:处理组装中的污染序列

      本文链接:https://www.haomeiwen.com/subject/jwavtktx.html