宏基因组之基因组装

简介: 宏基因组组装,即把短的reads拼装成连续的序列contig,再根据PE等关系将contig拼装成scaffold。与单个基因组组装不同,宏基因组组装最终得到的是环境样品中全部微生物的混合scaffold。理想情况下一条scaffold对应一个物种的全基因组。但由于序列太短或者覆盖度不够,很难拼出一条完整的基因组。针对高通量测序数据,出现了多种拼接算法和软件。

背景


宏基因组组装,即把短的reads拼装成连续的序列contig,再根据PE等关系将contig拼装成scaffold。与单个基因组组装不同,宏基因组组装最终得到的是环境样品中全部微生物的混合scaffold。理想情况下一条scaffold对应一个物种的全基因组。但由于序列太短或者覆盖度不够,很难拼出一条完整的基因组。针对高通量测序数据,出现了多种拼接算法和软件。

6069a84eeb37c0ed47283f519a0be1f.png

目前拼装方法主要有两类,即针对有参考基因组的重测序数据和针对无参考基因组的从头测序数据。对于之前未被测定过基因组序列的物种,并没有参照基因组可供使用,只能采用从头拼接的方法。针对从头测序数据,目前主要有两种拼接算法,即基于重叠图的算法(overlap-layout-consensus,OLC)和基于德布鲁因图(de Bruijin Graph)的算法。基于OLC的软件主要应用于第一代测序技术产生的数据,主要有PCAP,Celera Assmble。基于德布鲁因图的算法主要用于第二代高通量测序技术产生的数据,使用的软件中主要有Soapdenovo,Metahit、idba_ud,这里主要介绍下当下文章中流行的MegihitSoapdenovo软件的使用,其他软件之后有时间再进行讨论~

Megahit


安装

conda安装


conda install -c bioconda megahit

下载二进制文件


wget https://github.com/voutcn/megahit/releases/download/v1.2.8/MEGAHIT-1.2.8-Linux-x86_64-static.tar.gz
tar zvxf MEGAHIT-1.2.8-Linux-x86_64-static.tar.gz
cd MEGAHIT-1.2.8-Linux-x86_64-static/bin/
export PATH =<路径>:$PATH

参数介绍


# 常用参数:
-1/ -2 :分别输入双端序列
-m/--memory:SdBG 图构建过程中占用的最大内存字节数(如果设置为0-1,表示硬件内存的百分比。默认0.9)
--min-contig-len 设置拼接的contig最小长度
-t/--num-cpu-threads:使用的CPU线程数, 如果GPU允许激活,最少设置为2。
-o/--out-dir <字符型> 输出的文件夹
--out-prefix <字符型> 输出文件名的前缀
--k-list <奇数形式> 逗号分隔的kmer大小的列表,每个数必须是15-255的奇数,每两个相邻的kmer之间的增幅是小于28;默认:[21,29,39,59,79,99,119,141]
# 高级参数:
--no-mercy  Mercy k-mer用于宏基因组组装低覆盖序列。对于一般数据集>=30x,megahit可以使用-no-mercy选项生成更好的结果。
--low-local-ratio 定义低的局部覆盖的比率阈值 [0.2]
--max-tip-len  删除小于此值的端点[2*k]
--no-local  关闭局部组装模式

组装策略

假设我们手头有俩个样品(a和b)经过质控和去宿主得到的Clean data,我们先使用MEGAHIT组装软件进行组装分析(Assembly Analysisi)。

宏基因组中常见有俩种组装策略:

方法一:混拼

直接将各个样品通过上述的多样品组装的方式直接拼接为Contig即可。优点是简单快速,基因基本不用冗余,增加低丰度菌测序深度提高拼接长度,但是这种方法对内存性能需求大,时间周期较长,且混拼会提高错误拼接的风险。


# 多样品组装
megahit -1 a1.fq,b1.fq -2 a2.fq,b2.fq -o final_megahit.out 
# 组装时候选组MEGHIT默认参数(MEGAHIT默认参数组装得到的contigs N50较⾼,质量较好)进⾏组装。也可以自己设置kemr的范围及最小Contig长度`--k-list 21,29,39,59,79,99,119,141 --min-contig-len 500`。

方法二:单拼

首先逐个对样本进行组装,再将所有样本中未参与组装的reads混合起来进行组装成最终的Contig。资源消耗少,错误组装的几率小,但是这种方法基因会有大量冗余,后续要进行去冗余操作。

单样品组装


# 逐个单独组装
megahit -t 2 -m 0.95 --min-contig-len 500 -1 a1_.fq -2  a2_.fq -o  a_Contig.out
megahit -t 2 -m 0.95 --min-contig-len 500 -1 b1_.fq -2  b2_.fq -o  b_Contig.out

得到未参与组装的reads(可选)

这步目的是若你不想浪费数据可以再将未参与的双端reads提取出来再用同样的参数进行组装,最后和之前的Contig混合使用。

这里有两种方式:

  1. 一种是将各样品质控和去宿主后的 Clean Data 采⽤ Bowtie2 软件⽐对⾄各样品组装后的contigs上,后续根据samtools软件获取未被利⽤上的 PE reads(时间较长)。
  2. 使用Soapaligner短序列比对软件,这种方法速度非常快(推荐使用),而且还可直接生成未比对的序列unmapping.fa(未比对上的reads)。流程如下:

构建contig索引


# 这里以样品a为例,首先根据它的Congtig序列构建索引
2bwt-builder final.contigs.fa  # 在a_Contig.out目录下

比对


# 将质控后的序列比对至contig上
soap -p 6 -r 2 -m 200 -x 400 -M 4 -a a1_.fq -b a2_.fq -D final.contigs.fa.index -o a_PE.soap -2 a_SE.soap -u a_UN.fasta
# 参数介绍
-p  处理器数目,默认为1;
-r  比对到多个位置的序列如何输出,1表示不输出,2表示任意输出一次,3表示全部输出,默认为1;
-m  最小插入片段长度,默认400,single-end不需要设置该参数;
-x  最大插入片段长度,默认600,single-end不需要设置该参数;
-a  输入文件,pair-end时,输入其中一端文件;
-b  single-end不需要设置,pair-end输入另一端文件;
-u  没有比对上的reads;
-M  <int>   vi为每个read或read的seed部分匹配模式,不应该超过俩个误配。

ok,我们看一下生成的文件a_UN.fasta内容,可以看到序列的ID后面分别以1和2标识双端的reads,所以我们可以通过一些软件或者自编脚本分别将双端序列提取出来。

28dbb30fd01837e7fbc63107f10c313.png

这里提供两种简单解决思路:

方式一: 先提取出序列id,再根据id取序列,推荐使用软件seqtk工具(seqkit也可).


1. 可以分别先将标识1 和 2的序列ID取出来,然后再根据id取序列
less a_UN.fasta |grep "1$" |sed 's/>//'> a1_unmapping.id
less a_UN.fasta |grep "2$" | sed 's/>//'> a2_unmapping.id
2. 根据id取对应的序列
seqtk subseq a_UN.fasta a1_unmapping.id  > a1_unmapping.fa
seqtk subseq a_UN.fasta a2_unmapping.id   > a2_unmapping.fa

方式二:

提供一个简单python脚本,使用:python split_unmaping.py a1_unmapping.fa a2_unmapping.fa


####split_unmaping.py 
import sys
seq_file = {}
with open(sys.argv[1],"r") as old_fas: 
    for line in old_fas:
        line = line.strip()
        if line[0] == '>':
            seq_id = line
            seq_file[seq_id] = ''
        else:
            seq_file[seq_id] += line
new_fas1 = open(sys.argv[2],'w')
new_fas2 = open(sys.argv[3],'w')
for key,value in seq_file.items():
    if  '/1' in key:
        print(key + '\n' +value,file  = new_fas1)
    else:
        print(key + '\n'+ value, file = new_fas2)

混合组装所有样品未参与组装的reads

之后可以将两个样品未被利⽤上的 reads 放在⼀起,使⽤MEGAHIT进⾏混合组装,得到所有样品混合组装的contigs;组装参数与单样本组装相同,得到的Contigs在于之前各样品的Contigs进行组合得到最终的Contig。


# 同样的参数进行组装
megahit --min-contig-len 500 -1 a1_unmapping.fa,b1_unmapping.fa -2 a1_unmapping.fa,b2_unmapping.fa -o unmapping_megahit.out

Soapdenovo


安装


# 下载文件
wget https://github.com/aquaskyline/SOAPdenovo2/archive/r241.tar.gz
tar xzvf r241.tar.gz
cd SOAPdenovo2-r241/
make
# conda安装也可以

编译成功后路径下会生成3个文件:

  • SOAPdenovo-63mer
  • SOAPdenovo-127mer
  • SOAPdenovo-fusion。

前2个可执行文件用于组装, 63mer代表支持的kmer最大长度为63,127mer代表支持的kmer最大长度为127,除了支持的kmer长度不同外,其他用法完全相同。注意:该软件不像Megahit可以指定kmer list一个范围,每次组装只能设置一个kmer,我们使用时候一般要由低到高设置Kmer值比对组装效果(N50),不同类型样品所选择的kmer范围不同,大家请根据样品或者参考文献自行选择~~

准备config文件

Soapdenovo组装软件的使用需要提前准备一个配置文件config,将reads长度,插入片段,文库等信息写入其中。


##### 配置文件config
max_rd_len=151
[LIB]
avg_ins=379 #输出插入长度
reverse_seq=0
asm_flags=3
rd_len_cutoff=151
rank=1
q1=./a1.fq
q2=./a2.fq
########**配置文件说明**  
max_rd_len    表示read的最大长度;
[LIB]    表示文库信息标签;
avg_ins    文库中插入片段平均长度
reverse_seq    序列是否需要被反转,0(不反转),1(反转),一般插入片段大于等于2k文库,在建库是会将插入片段进行环化,此时须设置该参数为1:
asm_flags    表示reads用于组装哪个部分,可设为1,2,3, 1表示reads仅用于contig组装,2表示reads仅用于scaffold组装,3表示reads同时用于contig和scaffold组装;
rank    构建scaffold时,不同文库中reads的使用顺序,文库中reads序列越短,级别越高;
q1,q2    用于组装的双端fq文件。

单样品组装


SOAPdenovo-127mer all -s config -K 69 -M 3 -R -F -d 1 -u -o a.contig
## 参数
-s     config配置文件;
-K    k-mer的长度;
-o    输出文件前缀;
-d    [INT], kmerFreqCutoff, 去除频数小于等于该值的kmers,默认为0;
-M   [INT], mergeLevel连接contigs时, 合并相似序列的等级,默认为1,最小值为0,最大值为3,
-u    构建scaffold前屏蔽过高或过低覆盖度contigs,默认屏蔽;
-F     利用reads对scaffolds的gap进行填补,默认不执行;
-p    需要使用的cpu数目,默认8;

注意:这里我们演示的是单样品组装,若混拼只需在config文件q1,q2下方接着增添其他样品(如q1=./b1.fq  q2=./b2.fq)路径信息即可。

未参与组装的reads进行混装 (可选)

这里我们同样可以将未参与组装的reads利用起来,还是根据最开始提到Soapaligner比对的方法先得到个样品未mapping上的的双端reads,然后在配置文件输入进去:


max_rd_len=151
[LIB]
avg_ins=379 #输出插入长度
reverse_seq=0
asm_flags=3
rd_len_cutoff=151
rank=1
q1=./a1_unmapping.fa
q2=./a2_unmapping.fa
q1=./b1_unmapping.fa
q2=./b2_unmapping.fa

最后


若采用单拼之后再混合的策略,以上两种软件进行组装的最后步骤就是将合并单样品和混合组装⽣成最终的contigs,并进⾏统计分析和后续基因预测去冗余,建议大家可以尝试不同组装软件进行测试评估,比对N50,N90,Contig大小等指标差异,一般选择最大N50的软件即可。



作者:凯凯何_Boy

链接:https://www.jianshu.com/p/a53a43a0c9f5

来源:简书

著作权归作者所有。商业转载请联系作者获得授权,非商业转载请注明出处。

相关文章
|
算法 芯片
DNA测序原理:illumina和Pacbio对比介绍
DNA测序原理:illumina和Pacbio对比介绍
|
机器学习/深度学习 安全 数据挖掘
用于 DNA 测序的机器学习模型,理论上可以解码任何测序读数中所编码的数据值
用于 DNA 测序的机器学习模型,理论上可以解码任何测序读数中所编码的数据值
171 0
用于 DNA 测序的机器学习模型,理论上可以解码任何测序读数中所编码的数据值
|
数据采集
宏基因组数据 数据预处理
宏基因组下载数据过程
285 2
|
数据库 网络架构 索引
宏基因组之基因丰度计算
目前有两种方式可计算宏基因组基因的丰度,一种是基于比对,比如bwa,bowtie,soapaligner等主流的比对软件,另一种是不比对快速估计基因丰度,可以用近俩年来流行的salmon软件,这个软件在转录组的数据比对中也经常用到,可以直接计算出原始的Counts值和标准化的TPM值,此外由于是基于非比对,计算的速度得到很大的提升,同时也节省了部分的内存(少了庞大的sam/bam文件),可以说是多快好省,但是目前的高分文章中也还是不少用基于比对的方法去计算宏基因组的基因丰度的,下面我就分别简单介绍一下基于比对的soapaligner和不比对快速估计的samlon俩个软件的操作流程!!
786 0
|
Shell Python Perl
宏基因组之基因预测
12年有篇BMC的文献对几款预测的软件做了评估,其实参考大多数的文献中最常见的俩个软件也就是Prodigal和Metagenemark这俩个软件,分析过程中我这俩个软件都感受一下,现在将过程记录一下~~有兴趣的话可以看看这篇文献哦。
255 0
|
数据可视化 数据库 Python
scRNA分析 | 解决可能的报错,从0开始教你完成细胞通讯分析-cellphoneDB
scRNA分析 | 解决可能的报错,从0开始教你完成细胞通讯分析-cellphoneDB
642 0
|
存储 算法 数据可视化
利用TCseq包进行基因表达趋势分析
TCseq包提供了一个统一的套件去处理不同时序类型的数据分析,可以应用于转录组或者像ATAC-seq,Chip-seq的表观基因组时序型数据分析。该包主要的集中于不同时间点的差异分析,时间趋势分析及可视化作图。
436 0
|
数据采集 算法
测序质控和基因组组装原理
测序质控和基因组组装原理
|
数据采集 设计模式 存储
全基因组重测序流程【超细致!!】
全基因组重测序流程【超细致!!】
|
机器学习/深度学习 算法 数据可视化