处理组装中的污染序列

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

但是这一次,我测的一个样本的基因组足足比我预期的大了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进行验证。

©著作权归作者所有,转载或内容合作请联系作者
  • 序言:七十年代末,一起剥皮案震惊了整个滨河市,随后出现的几起案子,更是在滨河造成了极大的恐慌,老刑警刘岩,带你破解...
    沈念sama阅读 155,348评论 4 355
  • 序言:滨河连续发生了三起死亡事件,死亡现场离奇诡异,居然都是意外死亡,警方通过查阅死者的电脑和手机,发现死者居然都...
    沈念sama阅读 65,906评论 1 283
  • 文/潘晓璐 我一进店门,熙熙楼的掌柜王于贵愁眉苦脸地迎上来,“玉大人,你说我怎么就摊上这事。” “怎么了?”我有些...
    开封第一讲书人阅读 105,191评论 0 235
  • 文/不坏的土叔 我叫张陵,是天一观的道长。 经常有香客问我,道长,这世上最难降的妖魔是什么? 我笑而不...
    开封第一讲书人阅读 43,295评论 0 201
  • 正文 为了忘掉前任,我火速办了婚礼,结果婚礼上,老公的妹妹穿的比我还像新娘。我一直安慰自己,他们只是感情好,可当我...
    茶点故事阅读 51,606评论 3 283
  • 文/花漫 我一把揭开白布。 她就那样静静地躺着,像睡着了一般。 火红的嫁衣衬着肌肤如雪。 梳的纹丝不乱的头发上,一...
    开封第一讲书人阅读 40,107评论 1 202
  • 那天,我揣着相机与录音,去河边找鬼。 笑死,一个胖子当着我的面吹牛,可吹牛的内容都是我干的。 我是一名探鬼主播,决...
    沈念sama阅读 31,507评论 2 305
  • 文/苍兰香墨 我猛地睁开眼,长吁一口气:“原来是场噩梦啊……” “哼!你这毒妇竟也来了?” 一声冷哼从身侧响起,我...
    开封第一讲书人阅读 30,217评论 0 193
  • 序言:老挝万荣一对情侣失踪,失踪者是张志新(化名)和其女友刘颖,没想到半个月后,有当地人在树林里发现了一具尸体,经...
    沈念sama阅读 33,848评论 1 234
  • 正文 独居荒郊野岭守林人离奇死亡,尸身上长有42处带血的脓包…… 初始之章·张勋 以下内容为张勋视角 年9月15日...
    茶点故事阅读 30,183评论 2 238
  • 正文 我和宋清朗相恋三年,在试婚纱的时候发现自己被绿了。 大学时的朋友给我发了我未婚夫和他白月光在一起吃饭的照片。...
    茶点故事阅读 31,734评论 1 255
  • 序言:一个原本活蹦乱跳的男人离奇死亡,死状恐怖,灵堂内的尸体忽然破棺而出,到底是诈尸还是另有隐情,我是刑警宁泽,带...
    沈念sama阅读 28,084评论 2 247
  • 正文 年R本政府宣布,位于F岛的核电站,受9级特大地震影响,放射性物质发生泄漏。R本人自食恶果不足惜,却给世界环境...
    茶点故事阅读 32,608评论 3 228
  • 文/蒙蒙 一、第九天 我趴在偏房一处隐蔽的房顶上张望。 院中可真热闹,春花似锦、人声如沸。这庄子的主人今日做“春日...
    开封第一讲书人阅读 25,921评论 0 8
  • 文/苍兰香墨 我抬头看了看天上的太阳。三九已至,却和暖如春,着一层夹袄步出监牢的瞬间,已是汗流浃背。 一阵脚步声响...
    开封第一讲书人阅读 26,627评论 0 190
  • 我被黑心中介骗来泰国打工, 没想到刚下飞机就差点儿被人妖公主榨干…… 1. 我叫王不留,地道东北人。 一个月前我还...
    沈念sama阅读 35,103评论 2 261
  • 正文 我出身青楼,却偏偏与公主长得像,于是被迫代替她去往敌国和亲。 传闻我的和亲对象是个残疾皇子,可洞房花烛夜当晚...
    茶点故事阅读 35,054评论 2 257