细菌基因组拼接最终流程（三）
2021.4.3
持续更新中。。。
目的：从Illumina原始下机数据拼接细菌基因组草图。
软件：FastQC、Cutadapt、Velvet、SSPACE、gapfiller
1. FastQC(v0.11.9) —— 质控
fastqc <sequence_1.fastq.gz> <sequence_2.fastq.gz> -o <目录名> -t 线程数
重要参数：
1. -o：输出文件目录（需要提前新建一个目录）
2. -t：选择程序运行的线程数
3. -f：指定输入文件格式
主要结果文件：
*.html
2. Cutadapt(v3.3) —— 过滤
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 <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. Velvet —— 拼接contigs
步骤一：利用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. SSPACE —— 拼接scaffolds
4.1 下载安装
1. 下载bowtie(或者bowtie2)
conda install bowtie2
2. 从github下载最新版本加压缩进入后进入加压缩目录。（最新的免费版本是v2.1.1）
① 最新的版本可以用bwa直接处理压缩包，可惜作者并不是免费提供的。
② 解压缩后的SSPACE_Basic.pl即是主要的执行命令
4.2 使用
步骤一：写配置文件library.txt
#中间以空格符隔开
#1 2 3 4 5 6
#文库 正向序列 反向序列 插入文库大小 偏差 reas方向(pairend测序和matepaire不一样)
lib1 FP822_filtered_R1.fastq FP822_filtered_R2.fastq 500 0.25 FR
步骤二：scaffold
perl SSPACE_Basic_v2.0.pl -l <library.txt> -s <velvet_contigs.fa> -T 20 -b standard_out
重要参数：
1. -l：后接配置文
2. -s：要连接的contig序列
3. -T：-T：线程数（默认：1）
4. -b：输出文件前缀名
可选参数：
① -m：利用reads对contig进行衍生时，最小overlap的长度（默认：32bp)
② -o：利用reads对contig进行衍生时，最小reads覆盖的数量（默认：20）
③ -k：连接两条contig连接时，最小支持reads对数（默认：5对）
④ -a：连接两条contig连接时，最大连接的比率（默认：0.7）
⑤ -n：连接两条contig连接时，最小需要的overlap长度（默认：15bp）
⑥ -z：用于连接scaffold的最小contig长度（默认：0）
⑦ -g：比对过程中，允许的gap数（默认：0）
⑧-p：是否输出dot文件用于图形展示（默认：0，不输出）
5. GapFiller —— 拼接scaffolds
由于无法获得该软件，因此补洞的时候选择了其他的替代软件：soapdenovo2-gapcloser。其配置文件的书写方法总体同soapdenovo2，需要将asm_flags的参数设置为4即可。
- 配置文件
max_rd_len=150
[LIB]
avg_ins=439
reverse_seq=0
asm_flags=4
rank=1
pair_num_cutoff=3
map_len=32
q1=FP822_R1.fastq.gz
q2=FP822_R2.fastq.gz
- 使用方法
GapCloser -a <scaffolds.fasta> -b <config_file> -o <gapcloser_scaffolds.fasta> -l 150 -t 10
重要参数：
-a：输入scaffold文件
-b：输入配置文件
-o：输出文件名
-l：read的最大程度（默认100）
-t：线程数（默认1）
6. 结果
结果比较
7. 总结
- 根据目前所获得信息（拼接知识和软件），该套流程拼接结果和公司较为接近。
- 影响最终拼接效果的因素有很多，其中kmer取值较为关键。
- 如果有更好的软件，后续会继续更新，也希望有高手指出不足之处。
