通过前面两篇文章对MutMap原理和分析流程的介绍，我们对MutMap的分析已经有了比较深入的了解，是时候来点儿数据实战一下了！

MutMap的发明者已经搭建好了一套Pepline，并写好了MutMap pipeline framework供大家免费使用。我们只需要将数据按照要求进行放置和命名，并通过非常简单的几个命令，将任务提交给服务器，经过三五天的分析（以水稻为例），就可以拿到结果了。

---

1.分析环境的搭建

首先，我们需要有一台Linux服务器。由于分析过程中需要进行大量的运算，个人PC的配置根本吃不消。所以我们需要一台服务器，当然配置越高越好。在服务器上安装Linux系统，我们使用的是Redhat，当然也可以选择免费的Ubuntu或CentOS。

然后，在服务器安装如下软件：

Perl (v5.8.8) 5.22.2

R (version 2.15.0)3.0.1

BWA (version 0.5.9-r16)本版

SAMtools (0.1.8 or before) 本版

FASTX-Toolkit 使用conda安装0.0.14版

(http://hannonlab.cshl.edu/fastx_toolkit/download.html)

（括号中的为推荐版本，后面的黑体为我自己使用的版本）

关于FASTX-Toolkit版本的选择：

一定要安装0.0.13及以上的版本，因为从0.0.13开始FASTX-Toolkit才开始支持Illumina 1.8+ Phred+33的质量打分规则，否则会报错。就这一个报错，花费了我一天的时间，直到使用conda安装了0.0.14版本才解决。（2018.01.24记）。

2.MutMap pipeline framework的下载和安装

先下载MutMap pepline程序,解压到Linux系统的工作目录下（此处以myhome为例）。

#切换到安装目录

$ cd /myhome

#将framework压缩文件拷贝到相应文件夹

$ cp MutMap_framework1.4.4.tar.gz

#解压

$ tar zxvf MutMap_framework1.4.4.tar.gz

#改名

$ mv MutMap_framework1.4.4 MutMap_test

3.将测序数据链接到相应的文件夹中

将数据上传到服务器中专门的目录下，通过链接的方式，将原始文件链接到相应的操作目录下。

3.1 通过MD5值验证文件的完整性

（1）Windows中生成MD5值

参考链接：https://jingyan.baidu.com/article/f7ff0bfc18ff302e26bb1382.html

（2）在Linux中通过验证MD5值检验文件是否完整

参考链接：https://www.cnblogs.com/zhuxiaohou110908/p/5786893.html

3.2 将数据链接到野生型、分离群体混池文件夹中

$ cd /myhome/MutMap_test/1.qualify_read/anyname(野生型亲本)

$ ln –s /fastq/ZZZ/ZZZ_L003_R1_001.fastq.gz p_3_1_sequence.txt.gz

$ ln –s /fastq/ZZZ/ZZZ_L003_R2_001.fastq.gz p_3_2_sequence.txt.gz

$ ln –s /fastq/ZZZ/ZZZ_L004_R1_001.fastq.gz p_4_1_sequence.txt.gz

$ ln –s /fastq/ZZZ/ZZZ_L004_R2_001.fastq.gz p_4_2_sequence.txt.gz

$ cd /myhome/MutMap_test/1.qualify_read/mybulk(突变体混池)

$ ln –s /fastq/XXX/XXX_L005_R1_001.fastq.gz p_5_1_sequence.txt.gz

$ ln –s /fastq/XXX/XXX_L005_R2_001.fastq.gz p_5_2_sequence.txt.gz

$ ln –s /fastq/XXX/XXX_L006_R1_001.fastq.gz p_6_1_sequence.txt.gz

$ ln –s /fastq/XXX/XXX_L006_R2_001.fastq.gz p_6_2_sequence.txt.gz

4.放置参考基因组fasta文件

#将参考基因组fasta格式文件放在下面的文件夹中

/myhome/MutMap_test/downloaded_fasta/IRGSP-1.0_genome.fasta/

#为fasta格式参考基因组创建.fai文件

samtools faidx IRGSP-1.0_genome.fasta

5. 配置common.fnc文件

5.1编辑config.tx文件

在第12、15、31、52、115行做出如下修改：

#在第12行修改突变体混池名称

Key1_Bulked_sample_name="mybulk"

#在第15行修野生型亲本名称

Key1_My_cultivar_sample="anyname"

#第31行指定参考基因组的路径

Key2_Path_public_reference_FASTA="./downloaded_fasta/IRGSP1.0_genome.fasta/IRGSP-1.0_genome.fasta"

#在第52行配置BWA的线程数（实验室的服务器最大32，设置为20）

Key2_BWA_CPU=20

#第115行设置突变体混池的个体数

Key3_Individuals=20

其他设置一般不需要修改，具体设置参考MutMap Protocol

5.2 加载配置

#切换进入MutMap top目录下

$ cd /myhome/MutMap_test

#运行Bat_make_common.fnc.sh

$ ./Bat_make_common.fnc.sh

# 上面的命令会做出如下的操作

mv 1.qualify_read/mybulk 1.qualify_read/XXX

mv 1.qualify_read/anyname 1.qualify_read/ZZZ

mkdir -p 1.qualify_read/XXX/q30p90/sep_pair

mkdir -p 1.qualify_read/ZZZ/q30p90/sep_pair

6.分析数据

6.1 对测序数据进行质控和过滤

（1）对野生型亲本测序数据进行质控和过滤

# cd进入1.qualify_read文件夹下

$ cd /myhome/MutMap_test/1.qualify_read

#对野生型亲本测序数据进行质控和过滤

$ nohup ./Run_all_Bats.sh 9 &

#查看日志文件nohup.log，自动追加最新工作结果并打印到屏幕

tail -f nohup.log

（2）对突变体混池测序数据进行质控和过滤

#对突变体混池测序数据进行质控和过滤

$ nohup ./Run_all_Bats.sh 0 &

#查看日志文件nohup.log，自动追加最新工作结果并打印到屏幕

tail -f nohup.log

6.2 构建参考基因组

（1）Run “Run_all_Bats.sh”

将野生型亲本测序reads比对到参考基因组上，并替换SNP，构建自己的参考基因组：

#cd进入2.make_consensus文件夹下

$ cd /myhome/MutMap_test/2.make_consensus/

#运行程序

$nohup ./Run_all_Bats.sh &

（2）Run “Run_all_Bats.sh”

将野生型亲本的reads重新比对到新构建的参考基因组上，发现错配造成的SNP，用于下一步的排除和过滤：

#cd进入下面的文件夹

$ cd /myhome/MutMap_test/2.make_consensus/90.align_to_this_fasta/

#运行如下程序

$ nohup ./Run_all_Bats.sh &

6.3突变体混池数据分析及作图

#切换到3.analysis文件夹下

$ cd /myhome/MutMap_test/3.analysis/

#运行如下命令，进行数据分析和作图

$ nohup ./Run_all_Bats.sh &

7.数据存储位置

（1）构建的野生型参考基因组

数据存放在2.make_consensus/50.make_consensus中，格式为fasta

（2）突变体混池和野生型亲本比对鉴定出的snp信息

存储在3.analysis/70.awk_custom中

（3）SNP-index图

3.analysis/90.slidingwindow/pngs

（4）SNP低密度区域

每个窗口中的SNP数目小于10:

3.analysis/90.slidingwindow/pngs/Z1MZ2K/mut_index_X/mask10

8.使用samtools查看突变位点

使用方法：

samtools tview align.bam ref.fasta

samtools tview命令说明：

按？查看帮助，q或Esc退出

空格键向右移动一屏，Backspace键向左移动一屏

.意为和正链匹配，，意为和反向链匹配

大写字母意为，正链替换为该碱基；小写字母意为反向炼替换为该碱基

按g可跳转至指定的碱基，如：8:1256321

也可以使用IGV查看突变位点。

9. 突变位点的注释

Pepline的输出结果中并没有直接给出突变位点的注释信息，不利于突变位点的筛选，得到所有的SNP后，我们后续可以通过snpEFF、ANNOVAR、AnnTools等软件，对突变进行注释。

10. 结果展示

此处展示一张我用MutMap pipeline分析出的结果，找到了突变基因的连锁区间和候选突变位点。