基因组组装程序linux,基因组组装—-SPAdes(一)

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

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

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

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

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 :left reads, 代表第几个文库。

–pe-2 :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)。

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 -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结果

文章知识点与官方知识档案匹配，可进一步学习相关知识CS入门技能树Linux入门初识Linux24810 人正在系统学习中 相关资源：电脑耗电量测量软件