|||
通过前面两篇文章对MutMap原理和分析流程的介绍,我们对MutMap的分析已经有了比较深入的了解,是时候来点儿数据实战一下了!
MutMap的发明者已经搭建好了一套Pepline,并写好了MutMap pipeline framework供大家免费使用。我们只需要将数据按照要求进行放置和命名,并通过非常简单的几个命令,将任务提交给服务器,经过三五天的分析(以水稻为例),就可以拿到结果了。
首先,我们需要有一台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记)。
先下载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
将数据上传到服务器中专门的目录下,通过链接的方式,将原始文件链接到相应的操作目录下。
参考链接:https://jingyan.baidu.com/article/f7ff0bfc18ff302e26bb1382.html
参考链接:https://www.cnblogs.com/zhuxiaohou110908/p/5786893.html
$ 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
#将参考基因组fasta格式文件放在下面的文件夹中 /myhome/MutMap_test/downloaded_fasta/IRGSP-1.0_genome.fasta/ #为fasta格式参考基因组创建.fai文件 samtools faidx IRGSP-1.0_genome.fasta
在第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
#切换进入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
# cd进入1.qualify_read文件夹下 $ cd /myhome/MutMap_test/1.qualify_read #对野生型亲本测序数据进行质控和过滤 $ nohup ./Run_all_Bats.sh 9 & #查看日志文件nohup.log,自动追加最新工作结果并打印到屏幕 tail -f nohup.log
#对突变体混池测序数据进行质控和过滤 $ nohup ./Run_all_Bats.sh 0 & #查看日志文件nohup.log,自动追加最新工作结果并打印到屏幕 tail -f nohup.log
将野生型亲本测序reads比对到参考基因组上,并替换SNP,构建自己的参考基因组:
#cd进入2.make_consensus文件夹下 $ cd /myhome/MutMap_test/2.make_consensus/ #运行程序 $nohup ./Run_all_Bats.sh &
将野生型亲本的reads重新比对到新构建的参考基因组上,发现错配造成的SNP,用于下一步的排除和过滤:
#cd进入下面的文件夹 $ cd /myhome/MutMap_test/2.make_consensus/90.align_to_this_fasta/ #运行如下程序 $ nohup ./Run_all_Bats.sh &
#切换到3.analysis文件夹下 $ cd /myhome/MutMap_test/3.analysis/ #运行如下命令,进行数据分析和作图 $ nohup ./Run_all_Bats.sh &
数据存放在2.make_consensus/50.make_consensus
中,格式为fasta
存储在3.analysis/70.awk_custom
中
3.analysis/90.slidingwindow/pngs
每个窗口中的SNP数目小于10:3.analysis/90.slidingwindow/pngs/Z1MZ2K/mut_index_X/mask10
使用方法:
samtools tview align.bam ref.fasta
samtools tview命令说明:
按?
查看帮助,q
或Esc
退出
空格键
向右移动一屏,Backspace
键向左移动一屏
.
意为和正链匹配,,
意为和反向链匹配
大写字母意为,正链替换为该碱基;小写字母意为反向炼替换为该碱基
按g可跳转至指定的碱基,如:8:1256321
也可以使用IGV查看突变位点。
Pepline的输出结果中并没有直接给出突变位点的注释信息,不利于突变位点的筛选,得到所有的SNP后,我们后续可以通过snpEFF、ANNOVAR、AnnTools等软件,对突变进行注释。
此处展示一张我用MutMap pipeline分析出的结果,找到了突变基因的连锁区间和候选突变位点。
Archiver|手机版|科学网 ( 京ICP备07017567号-12 )
GMT+8, 2024-9-1 09:00
Powered by ScienceNet.cn
Copyright © 2007- 中国科学报社