使用MutMap快速定位突变基因:实战篇

通过前面两篇文章对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分析出的结果,找到了突变基因的连锁区间和候选突变位点。

©著作权归作者所有,转载或内容合作请联系作者
  • 序言:七十年代末,一起剥皮案震惊了整个滨河市,随后出现的几起案子,更是在滨河造成了极大的恐慌,老刑警刘岩,带你破解...
    沈念sama阅读 214,922评论 6 497
  • 序言:滨河连续发生了三起死亡事件,死亡现场离奇诡异,居然都是意外死亡,警方通过查阅死者的电脑和手机,发现死者居然都...
    沈念sama阅读 91,591评论 3 389
  • 文/潘晓璐 我一进店门,熙熙楼的掌柜王于贵愁眉苦脸地迎上来,“玉大人,你说我怎么就摊上这事。” “怎么了?”我有些...
    开封第一讲书人阅读 160,546评论 0 350
  • 文/不坏的土叔 我叫张陵,是天一观的道长。 经常有香客问我,道长,这世上最难降的妖魔是什么? 我笑而不...
    开封第一讲书人阅读 57,467评论 1 288
  • 正文 为了忘掉前任,我火速办了婚礼,结果婚礼上,老公的妹妹穿的比我还像新娘。我一直安慰自己,他们只是感情好,可当我...
    茶点故事阅读 66,553评论 6 386
  • 文/花漫 我一把揭开白布。 她就那样静静地躺着,像睡着了一般。 火红的嫁衣衬着肌肤如雪。 梳的纹丝不乱的头发上,一...
    开封第一讲书人阅读 50,580评论 1 293
  • 那天,我揣着相机与录音,去河边找鬼。 笑死,一个胖子当着我的面吹牛,可吹牛的内容都是我干的。 我是一名探鬼主播,决...
    沈念sama阅读 39,588评论 3 414
  • 文/苍兰香墨 我猛地睁开眼,长吁一口气:“原来是场噩梦啊……” “哼!你这毒妇竟也来了?” 一声冷哼从身侧响起,我...
    开封第一讲书人阅读 38,334评论 0 270
  • 序言:老挝万荣一对情侣失踪,失踪者是张志新(化名)和其女友刘颖,没想到半个月后,有当地人在树林里发现了一具尸体,经...
    沈念sama阅读 44,780评论 1 307
  • 正文 独居荒郊野岭守林人离奇死亡,尸身上长有42处带血的脓包…… 初始之章·张勋 以下内容为张勋视角 年9月15日...
    茶点故事阅读 37,092评论 2 330
  • 正文 我和宋清朗相恋三年,在试婚纱的时候发现自己被绿了。 大学时的朋友给我发了我未婚夫和他白月光在一起吃饭的照片。...
    茶点故事阅读 39,270评论 1 344
  • 序言:一个原本活蹦乱跳的男人离奇死亡,死状恐怖,灵堂内的尸体忽然破棺而出,到底是诈尸还是另有隐情,我是刑警宁泽,带...
    沈念sama阅读 34,925评论 5 338
  • 正文 年R本政府宣布,位于F岛的核电站,受9级特大地震影响,放射性物质发生泄漏。R本人自食恶果不足惜,却给世界环境...
    茶点故事阅读 40,573评论 3 322
  • 文/蒙蒙 一、第九天 我趴在偏房一处隐蔽的房顶上张望。 院中可真热闹,春花似锦、人声如沸。这庄子的主人今日做“春日...
    开封第一讲书人阅读 31,194评论 0 21
  • 文/苍兰香墨 我抬头看了看天上的太阳。三九已至,却和暖如春,着一层夹袄步出监牢的瞬间,已是汗流浃背。 一阵脚步声响...
    开封第一讲书人阅读 32,437评论 1 268
  • 我被黑心中介骗来泰国打工, 没想到刚下飞机就差点儿被人妖公主榨干…… 1. 我叫王不留,地道东北人。 一个月前我还...
    沈念sama阅读 47,154评论 2 366
  • 正文 我出身青楼,却偏偏与公主长得像,于是被迫代替她去往敌国和亲。 传闻我的和亲对象是个残疾皇子,可洞房花烛夜当晚...
    茶点故事阅读 44,127评论 2 352

推荐阅读更多精彩内容

  • 郭芳艳 焦点网络五期 焦作 坚持原创分享第25天 今天一大早一到单位,校长和同事都跟我说了一个坏消息——...
    冰山蓝鹰阅读 166评论 0 0
  • 最近看到两句话觉得很有感悟,第一句是人们宁可去关心一个蹩脚的演员的吃喝拉撒,也不愿去关注一个平凡人波澜壮阔的内...
    粉红嬷嬷阅读 175评论 0 0
  • 组别 301期 利他一组 【日精进打卡第165天 【知~学习】 诵读《六项精进大纲》3遍,累计261遍;诵读《大...
    J0hn先生阅读 178评论 0 0
  • 依然是昨天那个时候醒来,不知为什么近日常常失眠了。撩开窗帘,刷的透过一股清新的空气,吸到肚里,格外舒服。深冬大雪节...
    宏波_阅读 551评论 0 2
  • 冰冰冰_阅读 137评论 0 0