玩转基因组拼图：初次探秘SPAdes工具

关于基因组组装的流程可以参考我的简书:https://www.jianshu.com/p/e02ecd537c83

前面介绍了SOAPdenovo2的使用，但是由于结果整体不理想，尝试使用另外一个软件SPAdes。

该软件目前的版本支持paired-end reads, mate-pairs及unpaired reads。SPAdes一般用于小基因组组装，不适合大型基因组组装(例:哺乳动物)。官方文档可知第一步和第二步耗内存(cpu)，第三步耗磁盘空间(实际是缓存，属于高I/O操作)。

SPAdes包含几个单独的模块：

(1) BayesHammer:Illumina reads的read纠错。在单细胞和标准数据集上效果也很好。

(2) IonHammer:IonTorrent(半导体测序，通过半导体芯片直接把化学信号转为数字信号)数据的read纠错。

(3) SPAdes: 迭代short reads基因组组装模块。K的值是根据read长度和数据集类型自动选择的，也可以自己选择。

(4) MismatchCorrector:一种工具，可改善所产生的contig和scaffold的错配和较短的插入缺失。默认是关的，建议打开(对于小型基因组而言)。

运行SPAdes时建议运行BayesHammer或 IonHammer以获得高质量的组装(默认打开)，如果你使用别的read纠错工具，该模块可以关闭。

1.组装前准备

(1)文件夹组织形式:

clean_data:存放测序得到的reads。
fastqc_results:存放两次fastqc结果。(看reads质量)
kmergenie_results:存放kmergenie结果。(基因组调查)
QC_results:质控后的结果。
quast_results:quast结果。(组装结果评价)
reference:物种参考基因组。(用于quast，如果没有参考基因组可以不用)
spades_results:spades组装结果。
tools:组装所需工具。

(2)所有软件写进环境变量。

vim ~/.bashrc
#注:把所需工具写入到环境变量中。(fastqc,trimmomatic,kmergenie,SPAdes,quast)
source ~/.bashrc

2.read质量控制

所用软件为fastqc和trimmomatic(软件安装略)。

cd .../genome_assembly/spades                                   #注:...是省略了前面的路径,后面也一样。
touch fq_trim.sh
vim fq_trim.sh
#以下均在fq_trim.sh这个shell脚本中。
trimmomatic=.../genome_assembly/tools/Trimmomatic-0.38/trimmomatic-0.38.jar     #所安装的trimmomatic所在路径
#1.first fastqc
for fq_file in clean_data/*                                            
do
fastqc -o ./fastqc_results/first $fq_file                            
done
echo "** first fastqc down! **"
cd ./clean_data
for id in $(seq 501 503)                                        #注:这里根据你自己data的文件名进行设置。       
do
    file1="PB-${id}_1.clean.fq.gz"                              #测序得到的两个reads文件的文件名。                   
    file2="PB-${id}_2.clean.fq.gz"
    sample_id="PB-${id}"                                   
    if [ -f "$file1" ]                                          #判断file1是否存在。
    then
    cd ..
    fq1=./clean_data/$file1
    fq2=./clean_data/$file2
    outdir=.../genome_assembly/spades/QC_results                #trim后存放结果的文件夹。
    sample=$sample_id
    #2.trimmomatic(QC)
    time java -jar ${trimmomatic} PE -threads 80 -phred33\
           $fq1 $fq2 \
           $outdir/${sample}.paired.1.fq.gz $outdir/${sample}.unpaired.1.fq.gz \
           $outdir/${sample}.paired.2.fq.gz $outdir/${sample}.unpaired.2.fq.gz \
           ILLUMINACLIP:.../genome_assembly/tools/Trimmomatic-0.38/adapters/TruSeq3-PE-2.fa:2:30:10 \
           SLIDINGWINDOW:5:20 LEADING:5 TRAILING:5 MINLEN:50 && echo "** fq QC done **"
    cd ./clean_data
    fi
done
#3.second fastqc
cd ..
for fq_pair_file in $outdir/*
do
    fastqc -o ./fastqc_results/second $fq_pair_file
done
echo "** second fastqc down! **"
end=`date +%s`
dif=$[ end - start ]
echo 'the shell script run time is below:'
echo $dif

#运行
nohup ./fq_trim.sh >fq_trim.log &                                #放服务器后台运行

关于fastqc结果解读请参考简书:https://www.jianshu.com/p/dacedb7f6e2f

3基因组调查

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

(1)基因组杂合度。

(2)重复序列比例。

(3)预估基因组大小。

(4)测序深度。

3.1kmergenie评估基因组

KmerGenie可以进行k-mer分析及基因组大小评估。KmerGenie的最大优点在于可以实现在多个预设k-mer下的自动分析，除了进行常规的k-mer频数统计之外，还能够基于不同k-mer自动计算基因组大小，并为基因组组装评估一个最佳组装k-mer数值作为备选。

touch fq_501.txt                                 #存储fastq文件的路径。
vim fq_501.txt                                   #编辑fq.txt，输入PB_501经过trim后fastq文件路径(上一篇SOAPdenovo2的使用过程中我未进行trim，这里的话我前面第一步增加了trim)。
touch fq_503.txt                                 
vim fq_503.txt   
#写kmergenie的脚本
touch kmergenie.sh
vim kmergenie.sh
cd /sevzone/home/ytzheng/genome_assembly/spades
kmergenie fq_501.list -o kmergenie_result/PB_501 -k 140 -l 71 -s 6 -t 12
kmergenie fq_503.list -o kmergenie_result/PB_503 -k 150 -l 100 -s 6 -t 12
#-o:输出文件存放的目录。
#-k:最大的k值
#-l:最小的k值
#-s:从最小的k----最大的k，每次增加的k。
#-t:线程数。
#注:刚开始s可以设大点，可以根据判断出的最佳k值，缩小s再进行判断。(如果你设置的k范围内没有最佳值，请调整 k范围)

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

kmer频数分布图:

3.2jellyfish && GenomeScope评估基因组

这里请移步到:https://mp.weixin.qq.com/s/3tC4w9Z1LCpPueS6llcK4w

4基因组组装

4.1SPAdes软件下载

github：http://cab.spbu.ru/files/release3.14.0/SPAdes-3.14.0-Linux.tar.gz

#解压
tar -xzf SPAdes-3.14.0-Linux.tar.gz
cd SPAdes-3.14.0-Linux/bin/
ls
#添加到环境变量(前面部分已经添加)
#测试
spades.py --test

(1)可执行文件

(2)test结果(说明每个文件的代表什么)

4.2组装

注:为了运行read纠错，reads必须是FASTQ或BAM格式。

4.2.1参数介绍

spades.py        #查看参数(怎么用具体可以看官网说明，我仅介绍其中一些参数。)
#Basic options
-o:指定输出文件目录(必需)。
其他的参数:根据你是什么数据进行设置(如果你的数据是宏基因组、质粒、单细胞、RNA-seq、 IonTorrent等数据，可以指定相应参数)。
#Pipeline options
--only-error-correction:仅进行read纠错。
--only-assembler:仅进行组装。
--careful:减少错配和短插入失。 让程序运行MismatchCorrector----不建议用于中型和大型基因组。官方文档中可以看到MismatchCorrector占用的disk space非常多。
注:默认情况下执行read纠错、组装。
--continue:简单理解就是你如果跑其中一个k时程序停止，SPAdes会继续运行该k。
#Input data(仅介绍paired-end)
--pe<#>-1 <file_name>:left reads, <#>代表第几个文库。
--pe<#>-2 <file_name>:right reads, <#>代表第几个文库。
note:只有一个文库，指定#等于1即可。
#Advanced options
-t:线程数，根据服务器有多少线程设置。
-m:memory限制(Gb),达到该值程序终止，默认是250Gb，当然不指定它在运行时候也会给出一个值，如果因为该值太低跑不下去，你也可以自己指定。
-k:使用的k-mer用逗号分隔(所有值必须是奇数，小于128并且按升序列出)。--sc：设置默认值为21,33,55。
--cov-cutoff:read coverge的截止值，可以是float值、auto、off。默认是off。
--phred-offset:碱基质量体系，33或64，可以自己指定，程序也会自动检测。

4.2.2正式组装

#1.root用户下设置打开文件的最大数量(对于高I/O,占用的是缓存,不占用CPU,为了提高速度,可通过线程数增加来提高，但这样会增加文件打开数，所以需要设高一点)。
vi /etc/security/limits.conf
xx soft nofile 6000         #这里xx代表的是你用的linux用户名。
xx hard nofile 6000

#2.非root用户组装脚本
touch spades.sh
vim spades.sh
spades.py --pe1-1 QC_results/PB-501.paired.1.fq.gz --pe1-2 QC_results/PB-501.paired.2.fq.gz -t 200 -k 97,107,117,127 -m 600 --careful --phred-offset 33 -o spades_results/PB_501_trim
#--pe1-1、--pe1-2:pair-end测序得到的两个文件名，这里我用的是经过trim后的数据。
#-t:线程数(对于CPU密集型的程序，应该使用少一点的线程，对于IO密集型任务，应该使用多一点的线程)。
#-k:kmer选择的值，一般不需要设置，默认的话21,33,55,77，但是对于我的data(水稻)组装结果明显不好，所以我决定以前面kmergenie预测出来的结果为参考设置k值。
#-m:内存。---跑不下去加到跑的下去，不能超。单位为Gb，注:free -m可以查看最大的内存。
#--careful:运行BayesHammer(read纠错)、SPAdes、MismatchCorrector。
#--phred-offset:碱基质量体系，33或64。
#-o:输出文件目录。
#note:spades还有一个优势，如果你有其他组装工具产生的contigs的话可以，可以使用--trusted-contigs或 --untrusted-contigs选项。前者是当获得高质量组装中使用的，这些contigs将用于graph的构建、间隙填补、处理重复序列。第二个选项用于相对不太可靠的contigs，这些contigs用于间隙填补、处理重复序列。

#3.一些问题。注:由于我用的SPAdes跑的基因组与官方文档的相比相对较大，在运行过程中出现过以下问题。

(1).前两步运行时间相对还好，占用内存峰值约200G(相对而言是CPU密集型)，第三步占用的缓存过多，导致服务器的缓存几乎达到缓存峰值，这里如果缓存占用太多的话可以用以下脚本释放缓存(root用户下)，但是这样释放缓存会导致程序跑的慢，因此如果服务器内存较大的话还是建议自己手动释放。
touch release_memory.sh
vim release_memory.sh
#以下内容在shell脚本中
sleeptime=1800                      #这里可以设置多少秒清除一次缓存。
for i in {1..1000}
do
   sync
   echo "sync done"
   echo 3 > /proc/sys/vm/drop_caches
   echo "release done"
   sleep ${sleeptime}
done
#运行
nohup ./release_memory.sh >release_memory.log &

(2).第三步为I/O密集型，我的数据跑了很久，为了提高速度，我把线程数调高，但是调高的同时导致我前两步报错超过了文件打开的数量(默认是1024，但还好没蹦，只是在最后结果给出warning)，这里为了防止程序crash掉，我是通过root用户来增加打开文件的数量(见前面1)。

(3).内存不足的问题(out of memory)。
free -m  #查看服务器的内存情况，如果还要很多空闲的内存，把组装脚本中的-m参数调高。

4.3结果文件

<output_dir>/corrected/:包含用BayesHammer纠错reads产生*.fastq.gz文件。
<output_dir>/scaffolds.fasta:包含产生的scaffolds。
<output_dir>/contigs.fasta:包含产生的contigs。
<output_dir>/assembly_graph.gfa:包含SPAdes组装图和scaffolds的路径。
<output_dir>/assembly_graph.fastg:包含SPAdes组装图。(FASTG格式)
<output_dir>/contigs.paths:组装图中包含与contigs.fasta对应的路径。
<output_dir>/scaffolds.paths:组装图中包含与scaffolds.fasta对应的路径。
#note:一般只要scaffolds和contigs就可以了。

5.评价组装结果

QUAST软件下载略。

5.1参考基因组下载

如果你组装的物种有参考基因组，可以去下载相应的参考基因组。

5.2QUAST使用

touch quast.sh
vim quast.sh
#以下均在shell脚本中。
cd .../genome_assembly/spades
quast.py -o ./quast_results/PB_501_trim_quast -R ./reference/xxxx.fasta.gz -t 70 ./spades_results/PB_501_trim/scaffolds.fasta ./spades_results/PB_503_trim/scaffolds.fasta
#-o输出目录
#-R参考基因组序列
#-t线程数
#后面为要比较的contig/scaffold所在目录。

5.3QUAST结果

具体结果解读可以参考:http://blog.sciencenet.cn/blog-3406804-1163959.html

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