2021.3.18
持续更新中。。。
目的：比较下目前常用的四款拼接软件：Velvet、SPAdes、SOAPdenovo、ABySS
1. FastQC —— 质量评估
FastQC是一款基于Java的软件，它可以快速地对测序数据进行质量评估，并生成网页版的报告。
1.1 使用
fastqc <sequence_1.fastq.gz> <sequence_2.fastq.gz> -o <目录名> -t 线程数
重要参数
1. -o：输出文件目录
2. -t：选择程序运行的线程数
1.2 输出文件
每个原始数据文件会生成两个文件，一个为html网页文件，一个为.zip*文。打开.html可以查看序列的测序质量情况，分为三个等级：合格项为绿色√，警告是黄色的！，不合格为红色的×。（绿色越多越好）
2. Cutadapt —— 过滤
Cutadapt是一款用python写的去接头软件，可以去除任意指定的接头和序列，同时也可以用于质量过滤。（实际上原始数据过滤会过滤两种序列：去接头引物、低质量序列）
公司返回的
rawdata不含有引物序列，可能含有接头序列，需要进行过滤后进行分析
2.1 下载
conda install cutadapt -y
注意cutadapt(v3.3)的安装环境python版本不能高于3.9
2. 2 使用
cutadapt -a AGATCGGAAGAGCACACGTCTGAACTCCAGTCA -A AGATCGGAAGAGCACACGTCTGAACTCCAGTCA -o sequence_filtered_1.fastq.gz -p sequence_filtered_2.fastq.gz
-q 20,20 --max-n 0.1 -m 75 sequence_1.fastq.gz sequence_2.fastq.gz
重要参数
1. -a：paired-end测序文件中正向序列文件接头序列（文件1）
2. -A：paired-end测序文件中反向序列文件接头序列（文件2）
3. -o：文件1去接头后的结果文件
4. -p：文件2去接头后的结果文件
5. -q：序列两端碱基质量低于某一数值时被切除，用,隔开，例如：-q 20,20。
6. -m：reads1和reads2中切除接头后的序列长度最小值，低于这个数值则去除，例如：-m 75。
7. --max-n：N碱基占比一定比例时被去除，例如--max-n 0.1，表示N碱基占read比例到10%时会被去除。
8：最后是两个原始序列文件。
可选参数
1. -- pair-filter=(any|both)：any表示read1 和 read2任何一个检测到接头均舍弃；both表示 read1 和 read2 全部检测到接头才舍弃read1 和 read2。（默认any）
2. -O ：默认为3，即至少三个碱基配才认为是adapter序列。
3. -e：最大错配比例，默认是0.1。（解释：cutadapt在一条read中检测到20bp的接头序列，那么允许该20bp的接头序列有2个碱基的错配）
3. Kmergenie预测最佳kmer值
给定一组输入，KmerGenie首先计算多个K值的k-mer丰度直方图。然后，对每个K值，它预测数据集中不同基因组k-mer的数量，并返回这个数字最大化的k-mer长度。
3.1 安装
conda create -n kmer python=2.7 -y
conda install kmergenie
3.2 使用
kmergenie <fq.list> -o <result> -s 10 -t 10
重要参数：
1. fq.list：包含需要查询的文件，一行一个文件名（文件所在的绝对路径）
2. -o：结果输出的前缀名
3. -l：系统考虑的最小k值（默认：15）
4. -k：系统考虑的最大k值（默认：121）
5. -s：从最小k值到最大k值，每次增加的值（默认：10）
6. -t：线程数
3.3 输出文件
输出文件中主要看*report.html网页文件，会推算出最合适的kmer值与基因组大小。在这里，我的基因组得出的kmer=111。
话说，kmer值是越大越好还是越小越好呢？
4. 拼接 —— contigs
4.1 ⭐Velvet⭐
Velvet，基于kmer的序列拼接工具，用于基因组的de novo组装，支持多种格式，局部拼接效果好，不善于连接成scaffold。
4.1.1 安装
conda install velvet -y
如果是编译安装可以选择先修改以下配置文件参数后进行安装：
1. CATEGORIES：根据原始数据相应增减该值的小（一般设置为10）。
2. MAXKMERLENGTH：最大的kmer长度（一般设置为127）。
4.1.2 使用
分为两步执行：
步骤一：利用velveth对数据构建一个hash表
velveth <output> 111 -shortPaired -fastq -separate <sequence_filtered_1.fastq.gz> <sequence_filtered_2.fastq.gz>
重要参数：
1. output：输出文件目录
2. 111：即hash_lenghth，用来设置k-mer的大小。也可以是31,97,2的形式，指分别拼接kmer从31到97，依次增加2的序列（默认：31）
3. -shortPaired：reads的类型。
4. -fastq：输入文件的格式（默认是fasta）。
5. -separate：分开两个文件读入。
步骤二：velvetg进行序列拼接
velvetg <output> -exp_cov auto -cov_cutoff auto
重要参数：
1. <output>：velveth生成的结果目录
2. -exp_cov：期望的kmer覆盖度，设置成auto用于标准的基因组测序。该参数设置成auto后，-cov_cutoff也许设置成auto。
4.2 ⭐SPAdes⭐
SPAdes是常用的序列拼接软件之一，支持illumina、PacBio、Nanopore、Sanger、Ion Torrent等测序数据的拼接，同样适合用于混合组装来改善拼接效果（尤其适用于Ion Torrent数据的拼接）
4.2.1 安装
conda install spades -y
4.2.2 使用
spades.py --pe1-1 <sequence_filtered_1.fastq.gz> --pe1-2 <sequence_filtered_2.fastq.gz> -o <spades_out>
-k 111 -t 20 --careful
重要参数：
1. --pe<#>-1：指定输入文库，其中#表示第几个文库。例如，第一个pairend文库就可以写成--pe1-1，之后接reads1文件。
2. --pe<#>-2：和上面类似，如果是大片段的matepair文库，就使用--mp1-1等。
3. -t：线程数（默认：16）。
4. --careful：减少错误和插入确实，添加此项会消耗更多的时间。
5. -o：输出目录
可选参数：
1. -k：k-mer值列表，数字必须是小于128的奇数。注：若小片段文库(150bp×2)数据量足够(50×+)，则推荐使用-k 21,33,55,77（默认：auto）
2. --pacbio：指定pacbio数据输入。
3. --nanopore：指定nanopore数据输入。
4. -m：用于内存限制（默认：250）
5. --phred-offset：碱基质量体系，在数据纠错中会用到，现在illumina数据一般采用phred 33质量体系（默认：auto-detect）。
4.3 ⭐SOAPdenovo⭐
SOAPdenovo是华大基因开发的SOAP软件包的一部分，主要用于短序列reads拼接，尤其是illumina测序数据。从小的细菌基因组到大的动植物基因组，人基因组都适用。特点是使用起来比较简单，但是却拥有很好的拼接效果，尤其在基因组构建Scaffold方面，效果很好。
4.3.1 安装
conda install soapdenovo2 -y
4.3.2 使用
SOAPdenovo-127mer all -s config_file -K 111 -o soap -d 1 -p 20
soapdenovo里面有两个执行程序：SOAPdenovo-63mer和SOAPdenovo-127mer，根据kmer值的大小进行选择。
配置文件config_file的设置
#soapdenovo configure file
max_rd_len=150 #使用reads的最大长度，一般设置为reads的长度。
[LIB] #文库信息以此开头，如有多个文库，每个文库都需要有一个单独的[LIB]标识符
avg_ins=439 #文库平均插入长度，即测序文库大小，一般取插入文库大小的均值。
reverse_seq=0 #序列是否需要反转。如果是大片段文库，需要设置为1。
asm_flags=3 #组装的标志。1用于构建contig，2用于构建scaffold，3同时用于构建contig和scaffold，4只用于补洞，不用于拼接。（一般小片段文库设置为3，大片段设置为2）
rank=1 #文库使用的优先级。
pair_num_cutoff=3 #可选参数，确定了reads连接成contig或scaffold的阈值。
map_len=32 #可选参数，规定了在map过程中 reads和contig的比对长度必须达到该值（短片断默认：32）
q1=./reads.1.fq.gz #用于连接的数据，支持fasta和fastq格式文件，也支持压缩格式。（结尾不能有空格）
q2=./reads.2.fq.gz
重要参数：
1. -s：接配置文件
2. -o：输出文件的文件名前缀
3. -K：输入kmer的大小
4. -p：线程数（默认：8）
5. -d：去除频数不大于该值的kmer（默认：0）
6. -D： 去除频数不大于该值的由k-mer连接的边（默认：1）
7. -F：利用read对scaffold中的gap进行填补（默认：不执行）
4.4 ⭐ABySS⭐
ABySS用于从头拼接双端测序短序列或较大基因组序列。
4.4.1 安装
conda install abyss -y
4.4.2 使用
- 拼接一个paired-end文库
abyss-pe k=111 name=abyss in='sequence_filtered_1.fastq.gz sequence_filtered_2.fastq.gz'
重要参数：
1. k：kmer值
2. name：输出文件名
3. in：两个paired-end文库
- 拼接多个paired-end文库
abyss-pe k=111 name=ecoli lib='pea peb' pea='pea_1.fa pea_2.fa' peb='peb_1.fa peb_2.fa'
5. 四种拼接软件的比较
5.1 QUAST的安装和使用
conda install quast -y
quast -o <result> <velvet_contig.fa> <abyss_contig.fa> <soapdenovo_contig.fa> <spades_contig.fa> -t 15
重要参数：
1. -o：输出文件目录
2. -t：选择程序运行的线程数
3. -f：指定输入文件格式
5.2 比较结果（contigs）
综合contigs数量、总长度、N50 来看，velvet拼接contig效果最好
5.2 比较结果（scaffolds）
除了Velvet软件之外，其他三款软件在拼接contigs的同时也可以拼接成scaffolds。同样综合scaffolds数量、总长度、N50 来看，从paired-end reads数据一次性拼接成scaffolds效果最好的是ABySS。但是在使用abyss的过程中没有找到调整线程数的参数，用时较长。
6. 总结
Velvet拼接出的159个contigs已经很接近公司反馈的最终结果了，后面再用报告中提到的两款软件试试。
