基因组组装----SOAPdenovo2

1.基因组组装的流程

基因组组装的大概流程如下：

(1) 测序得到raw reads序列。

(2) Reads质量评估。

(3) 原始reads质控。

(4) 选择合适的参数进行组装。

(5) Reads组装成contigs/scaffolds。

(6)组装结果评估。

file

2.reads质量控制

(1)认识你的数据。(read类型，read数量，其GC含量，可能的污染和其他问题)

(2)数据清理。(组装前清理原始数据可以使组装结果更好,因为移除了低质量和污染的reads)

(3)为组装软件提供一些组装的参数

2.1检查原始read数据的质量

对于fastq文件，建议的工具为fastqc。

我们感兴趣的:
(1)read长度:在设置组装的最大k-mer值时将很重要。

(2)质量编码类型:对于quality trimming软件很重要。(Phred33、64)

(3)GC含量:高GC的生物体往往组装得不好，并且read覆盖率分布可能不均匀。

(4)总的reads数:让你了解coverage范围。

(4)在read的开始、中间或结尾附近质量下降:确定可能的修整/清除方法和参数。

(5)出现高度重复的k-mers:说明序列可能存在污染。

(6)reads中存在大量的N:可能表明测序质量较差。你需要修剪这些read，以删除N.
下面是fastqc的命令：

for fq_file in raw_data/*                          #所有的原始data都在raw_data文件夹中
do
        fastqc -o ./fastqc_results/first $fq_file  #-o后面代表的是用于存放输出结果的文件夹。
done

对于每个fq文件，会生成一个网页报告，如果质量好的话就不需要进行下面一步。

file

2.2reads质量控制

既然您已经对原始数据有了一定的了解，那么在组装之前使用这些信息来清理和修剪reads的操作对提高其整体质量是很重要的。这里推荐的软件是Trimmomatic(它可以处理paired-end数据)

time java -jar ${trimmomatic} PE -threads 80\ #time计算一个命令(程序)执行时间
                                              #${trimmomatic}软件所在的路径
                                              #threads线程数，cat /proc/sys/kernel/threads-max可查看最大
        raw_data/Z4-G5_FDMS190672649-1a_1.fq.gz raw_data/Z4-G5_FDMS190672649-1a_2.fq.gz \                     #PE测序的产生的read1、read2文件。
        #四个质控后的文件
        QC_results/Z4-G5_FDMS190672649-1a.paired.1.fq.gz QC_results/Z4-G5_FDMS190672649-1a.unpaired.1.fq.gz \ #read1的输出文件名
        QC_results/Z4-G5_FDMS190672649-1a.paired.2.fq.gz QC_results/Z4-G5_FDMS190672649-1a.unpaired.2.fq.gz \ #read2的输出文件名
        #质控参数
        ILLUMINACLIP:./tools/Trimmomatic-0.38/adapters/TruSeq3-PE-2.fa:2:30:10:8:True \ 
        SLIDINGWINDOW:5:15 LEADING:5 TRAILING:5 MINLEN:50 && echo "** fq QC done **"

3组装

3.1软件下载与安装

(1)组装软件1采用的是SOAPdenovo2,可在github上下载，下载完成后完成后传到服务器。
https://github.com/aquaskyline/SOAPdenovo2

#注:不要去下release里面的(有bug，最大只支持63mer)，把整个clone下来即可。
unzip SOAPdenovo2-master.zip
cd SOAPdenovo2-master/
make

编译完成后生成三个可执行文件。
第一个支持最大k-mer为63，第二个文件支持的最大kmer为127，第三个文件用于单细胞测序和宏基因组测序数据。

file

3.2 基因组调查

对于较为复杂的基因组，通常会在正式组装前进行kmer分析，以了解以下信息。

(1)基因组杂合度。

(2)重复序列比例。

(3)预估基因组大小。

(4)测序深度。

所用的软件为kmergenie，该软件可在http://kmergenie.bx.psu.edu/ 下载，下载完成后传到服务器。

#安装(使用时，确保服务器装有Pthon与R)
tar -xzvf kmergenie-1.7051.tar.gz && cd kmergenie-1.7051
python setup.py install
 
#添加可执行权限
chmod -R 755 *

接下来介绍怎么用，参考：http://wap.sciencenet.cn/blog-3406804-1159967.html

touch fq_501.txt                                 #存储fastq文件的路径。
vim fq_501.txt                                   #编辑fq.txt，输入fastq文件路径。
mkdir kmergenie_result     
./tools/kmergenie-1.7051/kmergenie fq_501.list -o ./kmergenie_result/PB_501 -k 105 -
l 15 -s 10 -t 12                               
#-o:输出文件存放的地址。
#-k:最大的k值
#-l:最小的k值
#-s:从最小的k----最大的k，每次增加的k。
#-t:线程数。
#注:刚开始s可以设大点，可以根据判断出的最佳k值，缩小s再进行判断((第一次:I:15,k:105,s:10 第二次I:85,k:140,s:6)。

在结果文件中会生成report。报告会以折线图的形式给出每种长度的kmer下，估计的基因组大小，同时给出了最佳的kmer(评估基因组总大小最高)

file

在理想条件下，该图应该为一条平滑的曲线，有明确的最大值，如下图：

file

然而在一些情况下，图往往与上图相似，有多个局部最大值。 如下图

file

这些情况表明，对于某些k值，Kmergenie中的统计模型并不总是正确地拟合输入数据。因此，由Kmergenie预测的最佳k值可能不是最佳的。在随后的分析中，还尝试使用比Kmergenie预测的k大的k(图中绿色)

在其他情况下，基因组k-mers的数量不会随着k值的升高而下降。这是一个高覆盖率(high coverage)数据集的迹象。
什么是覆盖率呢？
覆盖率(coverage)：指的是测序过程中读取单个碱基的平均次数。如果覆盖度为100×，说明平均每个碱基测序了100次。对碱基测序的频率越高，数据的质量越高。

file

在这种情况下，建议以最大k值（大于默认值120重新使用Kmergenie。 异常高的k值(>100)可能会产生良好的结果。

kmer频数分布图:

file

3.3SOAPdenovo2使用

根据Kmergenie得到的预测的最佳k值后就可以进行组装。SOAPdenovo2使用过程中最重要的一步是配置文件中参数的设置。下面以一个配置文件为例进行介绍。

#创建空的配置文件
touch config_file
vim config_file                  

#输入配置文件参数
max_rd_len=150                  #最大的read长度，可以根据fastqc得到。
[LIB]
avg_ins=350                     #文库平均插入长度(或取插入大小分布图中的峰值位置)，根据建库情况设置。              
reverse_seq=0                   #序列是否需要反转,PE(插入片段小于1kb)测序不需要，设为0，SE(插入片段在2kb-5kb)测序需要，设为1.              
asm_flags=3                     #1:表示只组装contig，2:表示只组装scaffold，3:表示同时组装contig和scaffold，4:表示只补gap。
rd_len_cutoff=150               #作用跟max_rd_len相同，大于该长度的序列会被切除至该长度，一般该参数不用。
rank=1                          #如果建库时候建了多个不同大小的插入片段文库，需设不同rank跑。对于短片段设为1。
pair_num_cutoff=3               #至少有几对reads支持才可连接contig构建scaffold,小片段为3、大片段为5。
map_len=32                      #对于pair-end，默认为32，对于mate-pair，默认为35。
#一对fastq文件
q1=/genome_assembly/clean_data/PB-501_BDSW192002279-1a_1.clean.fq.gz            
q2=/genome_assembly/clean_data/PB-501_BDSW192002279-1a_2.clean.fq.gz           

#正式组装(kmer选103,113,123,127分别跑)
#SOAPdenovo2既可以一步组装也可以分四步，如果基因组大且复杂建议分四步。
 nohup ./tools/SOAPdenovo2-master/SOAPdenovo-127mer all -s config_file -d 1 -R -F -K 113 -p 60 -o PB_501_113 > PB_501_113.log &  #注:加nohup和&是为了让它在后台运行(防止远程连接断开程序停止)

#分四步(在k=103、113、123、及127(支持的最大kmer)跑了后发现127在同样参数下contig N50最大)，下面是为了优化结果跑的。
#pregraph(d=1,2,3,4)
nohup ./tools/SOAPdenovo2-master/SOAPdenovo-127mer pregraph -s config_file -o PB_501_127_d_3 -R -K 127 -p 60 -d 2 >PB_501_127_d_3_pre.log &
nohup ./tools/SOAPdenovo2-master/SOAPdenovo-127mer contig -g PB_501_127_d_3 -R -p 60 >PB_501_127_d_3_con.log &
nohup ./tools/SOAPdenovo2-master/SOAPdenovo-127mer map -s config_file -g PB_501_127_d_3 -p 60 -k 127 >PB_501_127_d_3_map.log &  
nohup ./tools/SOAPdenovo2-master/SOAPdenovo-127mer scaff -g PB_501_127_d_3 -F -p 60 >PB_501_127_d_3_scaf.log &

3.4组装结果统计

file

4.组装结果的评估

QUAST软件下载:https://sourceforge.net/projects/quast/files/quast-5.0.2.tar.gz/download

下载完成后传到服务器。

#解压安装
tar -zxvf quast-5.0.2.tar.gz && cd quast-5.0.2
python  setup.py install_full

4.1参考基因组下载

由于我跑的是水稻数据，下载的是日本晴的参考基因组:https://phytozome.jgi.doe.gov/pz/portal.html

4.2QUAST使用

./tools/quast-5.0.2/quast.py -o ./quast_results -R ./reference/IRGSP-1.0_genome.fasta.gz -t 70 ./PB_501_results/PB_501_127_d_1/PB_501_127_d_1.scafSeq ./PB_501_results/PB_501_127_d_2/PB_501_127_d_2.scafSeq ./PB_501_results/PB_501_127_d_3/PB_501_127_d_3.scafSeq ./PB_501_results/PB_501_127_d_6/PB_501_127_d_6.scafSeq
#-o输出目录
#-R参考基因组序列
#-t线程数
#后面为要比较的contig/scaffold所在目录。

4.3结果

4.3.1统计结果

将contig/scaffold比对到参考基因组上，下面是比对得到的结果。

file file

累计长度:若与灰色线重合越好，组装结果越好。

file

错误组装情况：

file

GC含量:

file

4.3.2contig/scaffold大小可视化

file

4.3.3比对到每条染色体情况

file

比对到每条染色体可视化情况:

file

本文由博客一文多发平台 OpenWrite 发布！