使用SHOREmap做mapping-by-sequencing

简介

SHOREmap可以用来分析传统作图群体（自然系natural strains和分化系,diverged accession杂交，或outcrossing)或近等作图群体（isogenic mapping population, 诱变后代与未诱变亲本进行杂交，即会交,backcrossing)所产生的重测序数据。根据作图群体构建方式的不同，SHOREmap的outcross或backcross采用不同基于滑窗(sliding)方式对等位基因频率进行分析。

SHOREmap的backcross和outcross都需要从突变重组库中获得的一致的碱基识别信息

安装

前置安装

SHOREmap需要DISLIN科学库进行数据可视化

但是在安装DISLIN之前还需要保证存在/usr/lib/libXm.so*和/usr/lib/libXm.so*，这两者的安全需要root权限，所以要么联系管理员，要么想办法绕开(这个办法，我还没有想到).

sudo apt-get update
sudo apt-get install libmotif4
sudo apt-get install libxt-dev

开始安装doslin库

cd /path/to/src
# 下载
wget ftp://ftp.gwdg.de/pub/grafik/dislin/linux/i586_64/dislin-11.0.linux.i586_64.tar.gz
# 解压缩
tar -zxvf dislin-11.0.linux.i586_64.tar.gz
cd dislin-11.0
# 加入系统路径
mkdir -p $HOME/biosoft/dislin
DISLIN=$HOME/biosoft/dislin
export DISLIN
# 安装
./INSTALL
# 复制dislin_d.h 到dislin的文件下
cp ./example/dislin_d.h $DISLIN
# 删除安装文件(可选)
rm -rf dislin-11.0

安装SHOREmap v3.x

我这次安装的是当前最新的3.4版本，其他版本估计换汤不换药。

cd $HOME/biosoft
wget http://bioinfo.mpipz.mpg.de/shoremap/SHOREmap_v3.4.tar.gz
# 替换SHOREmap下的dislin的一些文件
tar -zxvf SHOREmap_v3,4
rm dislin/*dislin_d.*
cp $DISLIN/*dislin_d.* dislin
# 编辑/etc/profile或.bashrc
vi .bashrc
export LD_LIBRARY_PATH=$HOME/src/SHOREmap_v3.4/dislin
# 退出保存.bashrc: Esc+:wq
source .bashrc
# 到之前安装的文件夹下
cd & cd src/SHOREmap_v3.4
(可选，如果没有g++)sudo apt-get install build-essential
make
# 将编译文件拷贝到习惯的文件夹中，然后添加执行路径
cp SHOREmap ../../biosoft/SHOREmap_v3.4
echo "export $HOME/bisoft/SHOREmap_v3.4" >> ~/.bashrc

最后，可以重新启动一下bash验证

官方网站提供的两个常见问题的解答

Note 1: if the compilation complains like "/usr/bin/ld: cannot find -lXt" (or "/usr/bin/ld: cannot find -ldislin_d"), please open the makefile with the command

vi makefile

Press keys 'esc' and 'i' on the keyboard to edit makefile; move the cursor with arrow keys to the position before -lXt, and edit -L/path/to/libXt.so/; if '-ldislin_d' is not found, edit -L/path/to/dislin_d/ before -ldislin_d). After that, press keys 'esc', type in :wq, and press enter to save editing and quit vi. ('-L' tells the linker where to find the library given by -l)

Note 2: if '/usr/lib/ld: warning: libXm.so.3, needed by ./dislin/libdislin_d.so, not found (try using -rpath or -rpath-link)' occurs, and you have installed libmotif4, do the following:

cp /usr/lib/libXm.so.4 /usr/lib/libXm.so.3

We can make SHOREmap avaiable for general use by inserting the following command into /etc/profile

export PATH=$PATH:/path/to/SHOREmap_v3.x

and

source /etc/profile

总体流程

OUTCROSS

outcross的基本步骤 描述
SHOREmap extract 提取与SNP突变相关的重测序的一致的识别
SHOREmap create 根据背景/亲本系的重测序质量创建SNP标记列表
SHOREmap outcross 进行等位基因频率分析并定义mapping interval（也就是找到突变所在的大致区域)
SHOREmap annotate 对mapping interval中的突变基因效应进行注释

BACKCROSS

backcross的基本步骤 描述
SHOREmap extract 提取与SNP突变相关的重测序的一致的识别
SHOREmap backcross 进行等位基因分析
SHOREmap annotate 对mapping interval中的突变基因效应进行注

具体步骤

下载数据

只安装软件，却没有数据，我们也只能干瞪眼。

oucross分析所需数据

# OCF2 
wget -4 -q http://bioinfo.mpipz.mpg.de/shoremap/data/software/OC.fg.reads1.fq.gz &
wget -4 -q http://bioinfo.mpipz.mpg.de/shoremap/data/software/OC.fg.reads2.fq.gz &
# Ler 
wget -4 -q http://bioinfo.mpipz.mpg.de/shoremap/data/software/OC.bg.reads1.fq.gz &
wget -4 -qh ttp://bioinfo.mpipz.mpg.de/shoremap/data/software/OC.bg.reads2.fq.gz &

backcross分析所需数据

# BCF2
wget -4 -q http://bioinfo.mpipz.mpg.de/shoremap/data/software/BC.fg.reads1.fq.gz &
wget -4 -q http://bioinfo.mpipz.mpg.de/shoremap/data/software/BC.fg.reads2.fq.gz &
# mir159a (Col)
wget -4 -q http://bioinfo.mpipz.mpg.de/shoremap/data/software/BC.bg.reads1.fq.gz
wget -4 -q http://bioinfo.mpipz.mpg.de/shoremap/data/software/BC.bg.reads2.fq.gz

其他数据

除了最基本的测序数据外，我们可能还需要参考基因组，已有的注释数据等

# 参考基因组
wget -4 -q http://bioinfo.mpipz.mpg.de/shoremap/data/software/TAIR10_chr_all.fas &
# 基因注释
wget -4 -q http://bioinfo.mpipz.mpg.de/shoremap/data/software/TAIR10_GFF3_genes.gff &
# SHORE操作的结果数据
wget -4 -q http://bioinfo.mpipz.mpg.de/shoremap/data/software/scoring_matrix_het.txt &

重测序

首先使用bwa,bowtie2等read比对工具将得到的数据比对到参考基因组上。
假设你当前处在MBS文件夹下，该文件下有如下文件

# 混池测序结果，双端
BC.fg.reads1.fq.gz
BC.fg.reads2.fq.gz
# 背景信息测序结果
BC.bg.reads1.fq.gz
BC.bg.reads2.fa.gz
# 拟南芥参考基因组
TAIR10_chr_all.fas
# 拟南芥注释信息
TAIR10_GFF3_genes.gff

以下操作都是基于上述文件进行。

第一步: 序列比对，产生SAM文件

mkdir index
# 创建比对所需索引
bowtie2-build TAIR10_chr_all.fas index/TAIR10
# 序列比对
bowtie2 -x index/TAIR10 -1 BC.fg.reads1.fq.gz -2 BC.fg.reads2.fq.gz -S FG.sam
bowtie2 -x index/TAIR10 -1 BC.bg.reads1.fq.gz -2 BC.bg.reads2.fq.gz -S BG.sam

第二步： SAMtools预测突变位点

为了加快运算速度，可以先转换格式，并排序

samtools view -b -o BG.bam BG.sam
samtools view -b -o FG.bam FG.sam
samtools sort -o BG.sorted.bam BG.bam 
samtools sort -o FG.sorted.bam FG.bam 
samtools index BG.sorted.bam
samtools index FG.sorted.bam

consensus-calling program 寻找可能的变异位点
由于官方说明中使用samtools版本为0.1.19，先要解释一下参数 mpileup -uD -f，

# -f:faidx indexed reference sequence file 前后版本一致
# -u:generate uncompress BCF output 前后版本一致
# -D:output per-sample DP in BCF (require -g/-u).与输出格式有关，目前改为-t

因此，对于1.4版本的samtools，相应的参数为

mpileup -u -t DP -f

官方说明的bcftools也是0.1.19，参数为bcftools view -vcg 旧版本的view在当前的版本用于过滤，功能被call替代

# -v Output variant sites only (force -c)
# -c属于Call variants using Bayesian inference. This option automatically invokes option -e.When -v is in use, skip loci where the fraction of samples covered by reads is below FLOAT. [0],目前被m取代
# -g Call per-sample genotypes at variant sites (force -c)，这个没有找到合适的替代参数

综上，推荐使用如下命令行

samtools mpileup -u -t DP -f ../../../index/TAIR10_chr_all.fa ../../align/bwa/default/BG.sorted.bam | bcftools call -vm -Ov > BG.vcf &
samtools mpileup -u -t DP -f ../../../index/TAIR10_chr_all.fa ../../align/bwa/default/FG.sorted.bam | bcftools call -vm -Ov > FG.vcf &

额外步骤：VCF格式转换

由于vcftools工具版本，所以最后的文件版本是4.2，而SHOREmap要求4.1。通过biostar找到高人写的降级工具（其实就是把一些字符替换一下，但是不了解vcf不同版本的差异话，是不知道怎么写）
把下面的代码存为vcf_dowgrade.sh

# If you are trying to view VCF 4.2 files in IGV - you may run into issues. This function might help you.
# This script will:
# 1. Rename the file as version 4.1
# 2. Replace parentheses in the INFO lines (IGV doesn't like these!)

function vcf_downgrade() {
  outfile=${1/.bcf/}
  outfile=${outfile/.gz/}
  outfile=${outfile/.vcf/}
  bcftools view --max-alleles 2 -O v $1 | \
  sed "s/##fileformat=VCFv4.2/##fileformat=VCFv4.1/" | \
  sed "s/(//" | \
  sed "s/)//" | \
  sed "s/,Version=\"3\">/>/" | \
  bcftools view -O z > ${outfile}.dg.vcf.gz
  tabix ${outfile}.dg.vcf.gz
}

其实对于单个文件而言，可以直接用以下命令

infile=BG.vcf
outfile=BG.vcf
bcftools  view --max-alleles 2 -O v ${infile} | \
sed "s/##fileformat=VCFv4.2/##fileformat=VCFv4.1/" | \
  sed "s/(//" | \
  sed "s/)//" | \
  sed "s/,Version=\"3\">/>/" | \
  bcftools view -O z > ${outfile}.dg.vcf.gz

使用SHOREmap寻找突变所在区

第一步：需要把bcf文件通过SHOREmap convert转换成SHOREmap能认识的格式

SHOREmap convert --marker samtools.vcf --folder path/to/folder -runid int

会生成三个文件3_converted_consen.txt, 3_converted_variant.txt and 3_converted_reference.txt.

第二步：提取候选分子标记的consensus information(mapping pool)

SHOREmap extract --chrsizes chromsize.txt --folder ../SHOREmap_analysis --marker 11_converted_variant.txt --consen 11_converted_consen.txt -verbose

第三步：
然后使用SHOREmap backcross分析。SHOREmap backcross可用来分析回交作图群体所得到重组后代混池数据。相对于传统作图群体，只有诱变剂产生的突变会分离，也只有这些才会用于突变定位。

SHOREmap backcross会尝试过滤出所有参考基因组和测序池之间不同部分用于找到突变点特异部分。为了保证不是自然变异或者是测序错误，测序池选择的部分要多次出现在亲本或背景中。然后根据前景和/或背景的(识别碱基,base calls,质量/覆盖率/等位基因)信息，确定是否把保留的SNP位点作为分子标记。在正确的筛选后（拟南芥大概有上百个标记），SHOREmap backcross就能在分析marker的AF后识别大致的峰。进一步对变异注释后，就能找到目标性状的候选基因了。

SHOREmap backcross所需的输入文件如下：

  • 染色体大小文件，--chrsizes。分为两行，一行是染色体位置，一行是染色体大小。scaffold同理
  • 候选marker文件。也就是使用SHOREmap convert通过vcf生成的converted_variant.txt，每一列的含义如下。
    Column Description
    1 Project name
    2 Identity of chromosome
    3 Position of the SNP-marker
    4 Reference base
    5 Alternative base (or mutant base)
    6 Quality of the alternative base (ranging from 0 to 40)
    7 Number of reads supporting the predicted base
    8 Ratio of reads supporting the predicted base to total coverage
SHOREmap backcross --chrsizes chromsize.txt --marker ../convert/11_converted_variant.txt --consen extracted_consensus_0.txt --folder ../BC_analysis -plot-bc --marker-score 40 --marker-freq 0.0 --min-converage 10 --max-coverage 80 -bg ../convert/12_converted_variant.txt  --bg-cov 1 --bg-freq 0.0 --bg-score 1 -non-EMS --cluster 1 --marker-hit 1 -verbose

第四步：对结果进行注释

SHOREmap annotate --chrsizes chromsize.txt --folder ../BC_analysis/ann --snp ../convert/11_converted_variant.txt --chrom 2 --start 1 --end 4000000 --genome ../../TAIR10_chr_all.fas --gff ../../TAIR10_GFF3_genes.gff

注意点：

Alignment：不同的比对软件和参数设计会对结果有多大影响
SNP calling: samtools+vcftools 或者是GATK产生的vcf文件所包含的结果对结果有什么影响
SHOREmap： 如何调整阈值，如何辨认结果。

目标：

比对基本上是高通量的数据分析的必经之路，必须要对比对有足够的理解。
vcf格式描述了

最后编辑于 ：
©著作权归作者所有,转载或内容合作请联系作者
  • 人面猴
    序言：七十年代末，一起剥皮案震惊了整个滨河市，随后出现的几起案子，更是在滨河造成了极大的恐慌，老刑警刘岩，带你破解...
    沈念sama阅读 201,681评论 5赞 474
  • 死咒
    序言：滨河连续发生了三起死亡事件，死亡现场离奇诡异，居然都是意外死亡，警方通过查阅死者的电脑和手机，发现死者居然都...
    沈念sama阅读 84,710评论 2赞 377
  • 救了他两次的神仙让他今天三更去死
    文/潘晓璐 我一进店门，熙熙楼的掌柜王于贵愁眉苦脸地迎上来，“玉大人，你说我怎么就摊上这事。” “怎么了？”我有些...
    开封第一讲书人阅读 148,623评论 0赞 334
  • 道士缉凶录：失踪的卖姜人
    文/不坏的土叔 我叫张陵，是天一观的道长。 经常有香客问我，道长，这世上最难降的妖魔是什么？ 我笑而不...
    开封第一讲书人阅读 54,202评论 1赞 272
  • 港岛之恋（遗憾婚礼）
    正文 为了忘掉前任，我火速办了婚礼，结果婚礼上，老公的妹妹穿的比我还像新娘。我一直安慰自己，他们只是感情好，可当我...
    茶点故事阅读 63,232评论 5赞 363
  • 恶毒庶女顶嫁案：这布局不是一般人想出来的
    文/花漫 我一把揭开白布。 她就那样静静地躺着，像睡着了一般。 火红的嫁衣衬着肌肤如雪。 梳的纹丝不乱的头发上，一...
    开封第一讲书人阅读 48,368评论 1赞 281
  • 城市分裂传说
    那天，我揣着相机与录音，去河边找鬼。 笑死，一个胖子当着我的面吹牛，可吹牛的内容都是我干的。 我是一名探鬼主播，决...
    沈念sama阅读 37,795评论 3赞 393
  • 双鸳鸯连环套：你想象不到人心有多黑
    文/苍兰香墨 我猛地睁开眼，长吁一口气：“原来是场噩梦啊……” “哼！你这毒妇竟也来了？” 一声冷哼从身侧响起，我...
    开封第一讲书人阅读 36,461评论 0赞 256
  • 万荣杀人案实录
    序言：老挝万荣一对情侣失踪，失踪者是张志新（化名）和其女友刘颖，没想到半个月后，有当地人在树林里发现了一具尸体，经...
    沈念sama阅读 40,647评论 1赞 295
  • 护林员之死
    正文 独居荒郊野岭守林人离奇死亡，尸身上长有42处带血的脓包…… 初始之章·张勋 以下内容为张勋视角 年9月15日...
    茶点故事阅读 35,476评论 2赞 317
  • 白月光启示录
    正文 我和宋清朗相恋三年，在试婚纱的时候发现自己被绿了。 大学时的朋友给我发了我未婚夫和他白月光在一起吃饭的照片。...
    茶点故事阅读 37,525评论 1赞 329
  • 活死人
    序言：一个原本活蹦乱跳的男人离奇死亡，死状恐怖，灵堂内的尸体忽然破棺而出，到底是诈尸还是另有隐情，我是刑警宁泽，带...
    沈念sama阅读 33,226评论 3赞 318
  • 日本核电站爆炸内幕
    正文 年R本政府宣布，位于F岛的核电站，受9级特大地震影响，放射性物质发生泄漏。R本人自食恶果不足惜，却给世界环境...
    茶点故事阅读 38,785评论 3赞 303
  • 男人毒药：我在死后第九天来索命
    文/蒙蒙 一、第九天 我趴在偏房一处隐蔽的房顶上张望。 院中可真热闹，春花似锦、人声如沸。这庄子的主人今日做“春日...
    开封第一讲书人阅读 29,857评论 0赞 19
  • 一桩弑父案，背后竟有这般阴谋
    文/苍兰香墨 我抬头看了看天上的太阳。三九已至，却和暖如春，着一层夹袄步出监牢的瞬间，已是汗流浃背。 一阵脚步声响...
    开封第一讲书人阅读 31,090评论 1赞 258
  • 情欲美人皮
    我被黑心中介骗来泰国打工， 没想到刚下飞机就差点儿被人妖公主榨干…… 1. 我叫王不留，地道东北人。 一个月前我还...
    沈念sama阅读 42,647评论 2赞 348
  • 代替公主和亲
    正文 我出身青楼，却偏偏与公主长得像，于是被迫代替她去往敌国和亲。 传闻我的和亲对象是个残疾皇子，可洞房花烛夜当晚...
    茶点故事阅读 42,215评论 2赞 341

推荐阅读更多精彩内容

  • Spring Cloud
    Spring Cloud为开发人员提供了快速构建分布式系统中一些常见模式的工具（例如配置管理，服务发现，断路器，智...
    卡卡罗2017阅读 134,570评论 18赞 139
  • [儿童几何艺术搭建]用吸管玩游戏造故事
    吸管只能用来喝水?孩子们的脑洞那么大，怎么可能放过任何好玩的东西，在游戏中发挥创意，在故事里实现梦想，在生活中感悟...
    设计狮妈咪阅读 1,451评论 0赞 0
  • 刚毕业的我们，何去何从？
    如果你在一个企业学不到自己想要的知识，工作经验半年，你会跳槽吗？ 被问到这个问题，感觉回到了当初的自己。那时候，我...
    左手清香阅读 165评论 0赞 0
  • 1107
    梦到经历了很多很多的事情，最后父亲用电瓶车载我回家，醒来后鼻子一酸，真是对不起自己的家人，对不起自己逝去的岁月
    268fb13a2916阅读 259评论 0赞 0
  • 记得你，爱别人。
    “嘭”的一声，橙子先生的车门被关上，接着就看到一个清瘦的身影渐渐走远了，没有回头。 谁也不知道安妮心中此时在想些...
    等风来cll阅读 316评论 0赞 2