基因组组装程序linux,基因组组装----SPAdes(一)
前面系统性地介绍了SOAPdenovo2的使用方法，并因结果整体效果不尽如人意而转而使用另一个软件SPAdes进行分析
该软件目前的版本能够涵盖paired-end reads、mate-pairs以及unpaired reads等多种读取类型。SPAdes主要针对较小基因组进行组装，在处理较大基因组如哺乳动物系统时表现不佳。官方文档详细说明了第一步和第二步涉及CPU资源消耗（即计算资源），而第三步则涉及存储空间（实际上是缓存操作），属于高输入/输出类型
file
SPAdes包含几个单独的模块：
(1) BayesHammer:对Illumina reads进行读码校正，在单细胞和标准数据集上的表现良好
IonHammer:IonTorrent(半导体测序,利用半导体芯片将化学信号直接转换为数字信号)的数据读取与纠错处理
SPAdes: 基于迭代算法的短读基因组组装模块。K值的选择依据是基于read长度以及所处理的数据集类型自动进行，并且还可以自行设定。
Mismatch Corrector Tool:一种专有术语,旨在修复由错配或较短插入缺失引起的contig和scaffold问题。该工具处于关闭状态,推荐在小型基因组分析中启用。
运行SPAdes时推荐启用BayesHammer或IonHammer以显著提升组装质量（此功能已开启），若需采用其他读取错误校正工具则可关闭本模块
1.组装前准备
(1)文件夹组织形式:
file
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
file
2.read质量控制
所用软件为fastqc和trimmomatic(软件安装略)。
cd .../genome_assembly/spades #注:...是省略了前面的路径,后面也一样。
touch fq_trim.sh
vim fq_trim.sh
#以下均在fq_trim.sh这个shell脚本中。
trimmomatic等于所安装的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 & #放服务器后台运行
3基因组调查
在基因组具有较高复杂度的情况下，一般会在组装过程中执行k-mer分析以获取以下数据。
(1)基因组杂合度。
(2)重复序列比例。
(3)预估基因组大小。
(4)测序深度。
3.1kmergenie评估基因组
该程序支持执行k-mers分析以及对染色体长度的估算。其显著优势在于能够实现根据不同预设参数下的自动化处理过程，在此基础之上除了提供传统的频率统计结果之外还能根据不同的选择参数自动生成适合基因组拼接的最佳窗口设置建议
touch fq_501.txt #存储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的取值范围)
该报告将被创建于result文件中，并将以折线图的形式展示不同k值对应的估计基因组大小曲线图。此外该报告还将提供最佳kmer值（基于评估基因组总大小的最大化）
file
file
kmer频数分布图:
file
3.2jellyfish && GenomeScope评估基因组
4基因组组装
4.1SPAdes软件下载
#解压
tar -xzf SPAdes-3.14.0-Linux.tar.gz
cd SPAdes-3.14.0-Linux/bin/
ls
#添加到环境变量(前面部分已经添加)
#测试
spades.py --test
(1)可执行文件
file
(2)test结果(说明每个文件的代表什么)
file
4.2组装
注:为了运行read纠错，reads必须是FASTQ或BAM格式。
4.2.1参数介绍
spades.py #查看参数(怎么用具体可以看官网说明，我仅介绍其中一些参数。)
#Basic options
-o:指定输出文件目录(必需)。
其他参数设置项；基于你所处的数据类型来配置；例如：如果是宏基因组学数据、质粒载体、单细胞测序数据、RNA测序（如RNA-seq）或IonTorrent等类型的生化数据；则应对应设置相应的参数选项。
#Pipeline options
--only-error-correction:仅进行read纠错。
--only-assembler:仅进行组装。
--careful:避免发生因错配而导致的小插入错误。
确保程序正常运行MismatchCorrector----不建议在处理规模较大的基因组时使用。
官方文档明确指出该工具对磁盘空间的需求非常大。
注:默认情况下执行read纠错、组装。
--continue标识意味着：当程序运行到某个特定的k值时会终止；SPAdes会继续执行该k
#Input data(仅介绍paired-end)
--pe-1 :left reads, 代表第几个文库。
--pe-2 :right reads, 代表第几个文库。
note:只有一个文库，指定#等于1即可。
#Advanced options
-t:线程数，根据服务器有多少线程设置。
-m: memory限制 (Gb)，当程序运行至此内存限制时会终止。默认设置为 250 GB，在程序默认情况下应用此内存限制。无需指定此内存限制，在程序运行过程中会自动调整。若因设置过低导致无法正常运行，则可以选择自行设定。
-m: memory限制 (Gb)，当程序运行至此内存限制时会终止。默认设置为 250 GB，在程序默认情况下应用此内存限制。无需指定此内存限制，在程序运行过程中会自动调整。若因设置过低导致无法正常运行，则可以选择自行设定。
-k: 设置使用的k-mer以逗号分隔（所有值必须是奇数且小于128，并且按升序排列）。--sc选项：用于指定默认值的选项（例如：21、33、55）。
-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 QualityControlResults/PB-501.paired.1.fq.gz --pe1-2 QualityControlResults/PB-501.paired.2.fq.gz -t 200 -k min_k=97,max_k=127 -m max_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等数值。然而，在我的实验中使用水稻数据时发现组装效果并不理想，在此我决定以KMERGENIE软件预测的结果作为参考来设定k值参数
使用#-m标记表示内存占用情况。当程序运行缓慢时无法增加其运行速度或资源投入（即上限）。其中单位为GB，在注释信息中提到：使用free -m命令可以查询系统中使用的最大内存值。
#--careful:运行BayesHammer(read纠错)、SPAdes、MismatchCorrector。
#--phred-offset:碱基质量体系，33或64。
#-o:输出文件目录。
#note: Spades offers another benefit. If you have contigs generated from other assembly tools available, you can utilize them by using either the --trusted-contigs or the --untrusted-contigs option. The former is applicable when constructing high-quality assembles, where these contigs will be employed for graph construction, gap filling, and handling repetitive sequences. The latter is intended for less reliable contigs, which will also be used for gap filling and managing repetitive sequences.
第3部分有一些问题。因为我在使用SPAdes进行基因组比对时所处理的基因组规模与官方文档相比相对较大，在运行过程中出现过以下问题
(1) 前两步运行的时间较为理想。所占内存峰值约为200G（相对于属于CPU密集型运算）。第三步骤消耗的缓存过大导致服务器存储容量几乎被耗尽（即接近存储容量极限）。若要清理过多缓存空间，则可参考以下命令（针对root用户）：free -h /dev/shm等操作可能会降低程序运行效率但这样做可能会降低程序运行效率因此在服务器内存较为充足的场景下仍推荐手动清理空间
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 &
第三阶段属于I/O密集型。我的数据运行时间过长，在尝试加快运行速度时将线程数量提升了一级。然而这样做却引发了错误信息：文件打开数量已超限（默认值为1024）。为了避免程序崩溃，在切换到root账户以获取权限，并重新设置文件打开限制后问题得以解决。
file
(3).内存不足的问题(out of memory)。
免费使用-free命令查看可用内存占用情况：#在命令行界面输入free -m以查看当前服务器的内存使用情况；如果发现仍有大量未使用的内存空间（即剩余swap空间较多），建议适当增加组装脚本中的-m参数以提高内存分配效率。
file
4.3结果文件
/corrected/:包含用BayesHammer纠错reads产生*.fastq.gz文件。
/scaffolds.fasta:包含产生的scaffolds。
/contigs.fasta:包含产生的contigs。
/assembly_graph.gfa:包含SPAdes组装图和scaffolds的路径。
/assembly_graph.fastg:包含SPAdes组装图。(FASTG格式)
/contigs.paths:组装图中包含与contigs.fasta对应的路径。
/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 --output ./quast_results/PB_501_trim_quast --reference ./reference/xxxx.fasta.gz --trim_level 70 --input_dir ./spades_results/PB_501_trim/ ./spades_results/PB_503_trim/
#-o输出目录
#-R参考基因组序列
#-t线程数
#后面为要比较的contig/scaffold所在目录。
5.3QUAST结果
本文由博客一文多发平台 OpenWrite 发布！