在上一篇文章《使用MutMap快速定位突变基因：原理篇》中，我对MutMap的原理进行的简单的介绍。在本篇教程中，我会对MutMap分析的流程及分析过程中使用到的软件及其参数的设置进行逐一的分解，以求对MutMap的分析流程有更深入的理解。

---

1.对重测序数据进行质控，挑出高质量的reads

（1）使用软件：FASTX-Toolkit

① fastq_quality_filter：过滤掉低质量的reads

使用说明：

fastq_quality_filter [-h] [-v] [-q N] [-p N] [-z] [-i INFILE] [-o OUTFILE]过滤低质量序列

[-q N]       = 最小的需要留下的质量值

[-p N]       = 每个reads中最少有百分之多少的碱基需要有-q的质量值

[-z]         =压缩输出

[-v]       =详细-报告序列编号，如果使用了-o则报告会直接在STDOUT，如果没有则输入到STDERR

在MutMap中，使用q30p90对数据进行质控。

【特别说明】：对于 Illumina 1.8+ Phred+33规则进行质量打分的测序数据，需要添加-Q33选项，否则程序不能识别打分字符，会报错。

② fastx_quality_stats：序列的基本信息，如GC含量等

使用说明：

fastx_quality_stats [-h] [-N] [-i INFILE] [-o OUTFILE]

[-h] = 获得帮助信息

[-i] = FASTA/Q格式的输入文件

[-o] = 输出路径，如果不设置会直接输出到屏幕

（2）结果输出

① 过滤后的高质量reads

结果输出到每个样本下的q30p90/sep_pair/*1.RnFq30p90.1或2.gz，这些reads为质控后的高质量的reads，用于构建参考基因组和后续的比对分析。

② 过滤后reads的统计信息

统计结果输出到每个样本下的q30p90/sep_pair/*1.RnFq30p90.1或2.gz_stats.txt，这些txt格式的文件中存储着过滤后的reads的统计信息。

2.构建亲本参考基因组

（1）为要分析的数据创建链接

运行程序./Bat_make_symbolic_link_of_qualified_fastq.sh，在2.make_consensus文件夹中，为该步骤中要处理的数据创建链接，链接到1.qualify_read文件夹下的q30p90/sep_pair/*1.RnFq30p90.1或2.gz数据，这些数据为第一步质控过的高质量数据。

（2）将过滤后的亲本reads比对到基因组上

① 为参考基因组建立index

使用软件：BWA

使用命令：bwa index -p 输出文件名 -a reference.fa

② 寻找SA coordinates

确定reads在参考基因组上的坐标。

使用软件：BWA

使用命令： bwa aln

bwa aln -n 0.04 -o 1 -k 2 -t 20 reference.fa  leftRead.fastq > leftRead.sai

bwa aln -n 0.04 -o 1 -k 2 -t 20 reference.fa  rightRead.fastq > rightRead.sai

参数说明：

-n  ：NUM    Maximum edit distance if the value is INT, or the fraction of missing alignments given 2% uniform base error rate if FLOAT. In the latter case, the maximum edit distance is automatically chosen for different read lengths. [0.04]

-o：Maximum number of gap opens

-k：Maximum edit distance in the seed

-t：Number of threads (multi-threading mode)

③ 将SA coordinates输出文件转换为sam格式

将上一步的SA coordinates输出的fai格式的文件转换为sam格式的文件。

使用软件：BWA

使用命令：bwa sample

bwa sample -a 1000 -n 3 -N 10 reference.fa leftRead.sai leftRead.fastq  > leftRead.sam

bwa sample -a 1000 -n 3 -N 10 reference.fa  rightRead.sai  rightread.fastq > rightRead.sam

参数说明：

-a:Maximum insert size for a read pair to be considered being mapped properly. Since 0.4.5, this option is only used when there are not enough good alignment to infer the distribution of insert sizes.

-n:Maximum number of alignments to output in the XA tag for reads paired properly. If a read has more than INT hits, the XA tag will not be written.

-N:Maximum number of alignments to output in the XA tag for disconcordant read pairs (excluding singletons). If a read has more than INT hits, the XA tag will not be written.

④ 将sam文件转换为bam文件

使用软件：samtools

使用命令：samtools view

samtools view -Sb leftRead.sam > leftRead.bam

samtools view -Sb rightRead.sam > rightRead.bam

参数说明：

-S: input is SAM

默认下输入是 BAM 文件，若是输入是 SAM 文件，则最好加该参数，否则有时候会报错。

-b: output BAM

默认下输出是 SAM 格式文件，该参数设置输出 BAM 格式

⑤ 将leftRead.bam和rightRead.bam合并

使用软件：samtools

使用命令：samtools merge

samtools merge left_and_right.bam leftRead.bam rightRead.bam

⑥ 对merge后的bam文件排序

使用软件：samtools

使用命令：samtools sort

samtools sort left_and_right.bam left_and_right_sorted.bam

⑦ 移除潜在的PCR重复

使用软件：samtools

使用命令：samtools rmdup

samtools rmdup myBAMums_sorted.bam myBAMums_MSR.bam

⑧ 对bam文件建立索引

使用软件：samtools

使用命令：samtools index

samtools index myBAMums_MSR.bam

（3）对bam文件的比对质量进行再次优化

使用软件：coval refine

将Indel附近有错配的reads重新排列，直到错配碱基数最少

去除错配数目超过设置值的reads

去除含有3个以上Indel、1个Indel和一个soft-clipped end、或者2个soft-clipped end的reads。

文件输出：2.make_consensus\20.coval_refine文件夹下

（4） variant calling

使用软件：coval call

寻找亲本重测序reads和Nip之间的snp位点和Indel，用于下一步构建亲本参考基因组。

文件输出：2.make_consensus\30.coval_call文件夹下

（5）将找到的变异位点替换Nip，构建亲本参考基因组

使用脚本：

/Bat_RYKMSWBDHV_to_ACGT.pl.sh

/Bat_make_consensus.sh

文件输出：2.make_consensus\50.make_consensus文件夹下

(6) 将野生型亲本的reads重新比对到新构建的野生型参考基因组上

找由于错配导致的变异位点，用于后续的排除和过滤

使用软件：bwa、cover refine、cover call

文件输出： 2.make_consensus\90.align_to_this_fasta\30.coval_call文件夹下

3. 将突变体混池测序数据比对到亲本参考基因组上，寻找变异位点、计算SNP-index并作图

分析步骤：

将混池reds比对到亲本参考基因组上

使用coval refine对比对结果进行过滤

使用coval call检测变异位点

排除由于错配导致的snp位点

计算将找到的可信的snp位点的深度和snp-index

计算snp的置信区间

画SNP-index的slidingwindow图

4.结果输出：

可信的变异位点位于文件夹：3.analysis\70.awk_custom

SNP-index图位于文件夹3.analysis\90.slidingwindow