搜索
查看: 3094|回复: 1

SHOREmap练习MBS

[复制链接]

18

主题

56

帖子

403

积分

管理员

Rank: 9Rank: 9Rank: 9

积分
403
发表于 2017-3-17 20:35:47 | 显示全部楼层 |阅读模式
本帖最后由 hoptop 于 2017-3-17 20:54 编辑

简介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权限,所以要么联系管理员,要么想办法绕开(这个办法,我还没有想到).
[Shell] 纯文本查看 复制代码
sudo apt-get update
sudo apt-get install libmotif4
sudo apt-get install libxt-dev
开始安装doslin库
[Shell] 纯文本查看 复制代码
[/size][/color]
cd /path/to/src# 下载
wget ftp:[color=#93a1a1]//ftp.gwdg.de/pub/grafik/dislin/linux/i586_64/dislin-11.0.linux.i586_64.tar.gz[/color]
# 解压缩
tar -zxvf dislin[color=#2aa198]-11.0[/color].linux.i586_64.tar.gzcd dislin[color=#2aa198]-11.0[/color]# 加入系统路径mkdir -p $HOME/biosoft/dislinDISLIN=$HOME/biosoft/dislinexport DISLIN
# 安装
./INSTALL
# 复制dislin_d.h 到dislin的文件下
cp ./example/dislin_d.h $DISLIN
# 删除安装文件(可选)
rm -rf dislin[color=#2aa198]-11.0[/color]

安装SHOREmap v3.x我这次安装的是当前最新的3.4版本,其他版本估计换汤不换药。
[Shell] 纯文本查看 复制代码
cd $HOME/biosoft
wget [url=http://bioinfo.mpipz.mpg.de/shoremap/SHOREmap_v3.4.tar.gz]http://bioinfo.mpipz.mpg.de/shoremap/SHOREmap_v3.4.tar.gz[/url]
# 替换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







总体流程OUTCROSS[td]
outcross的基本步骤描述
SHOREmap extract提取与SNP突变相关的重测序的一致的识别
SHOREmap create根据背景/亲本系的重测序质量创建SNP标记列表
SHOREmap outcross进行等位基因频率分析并定义mapping interval(也就是找到突变所在的大致区域)
SHOREmap annotate对mapping interval中的突变基因效应进行注释
BACKCROSS[td]
backcross的基本步骤描述
SHOREmap extract提取与SNP突变相关的重测序的一致的识别
SHOREmap backcross进行等位基因分析
SHOREmap annotate对mapping interval中的突变基因效应进行注
具体步骤下载数据
只安装软件,却没有数据,我们也只能干瞪眼。
oucross分析所需数据
[Shell] 纯文本查看 复制代码
# OCF2 
wget -4 -q [url=http://bioinfo.mpipz.mpg.de/shoremap/data/software/OC.fg.reads1.fq.gz]http://bioinfo.mpipz.mpg.de/shor ... /OC.fg.reads1.fq.gz[/url] &
wget -4 -q [url=http://bioinfo.mpipz.mpg.de/shoremap/data/software/OC.fg.reads2.fq.gz]http://bioinfo.mpipz.mpg.de/shor ... /OC.fg.reads2.fq.gz[/url] &
# Ler 
wget -4 -q [url=http://bioinfo.mpipz.mpg.de/shoremap/data/software/OC.bg.reads1.fq.gz]http://bioinfo.mpipz.mpg.de/shor ... /OC.bg.reads1.fq.gz[/url] &
wget -4 -qh ttp://bioinfo.mpipz.mpg.de/shoremap/data/software/OC.bg.reads2.fq.gz &

backcross分析所需数据
[Shell] 纯文本查看 复制代码
# BCF2
wget -4 -q [url=http://bioinfo.mpipz.mpg.de/shoremap/data/software/BC.fg.reads1.fq.gz]http://bioinfo.mpipz.mpg.de/shor ... /BC.fg.reads1.fq.gz[/url] &
wget -4 -q [url=http://bioinfo.mpipz.mpg.de/shoremap/data/software/BC.fg.reads2.fq.gz]http://bioinfo.mpipz.mpg.de/shor ... /BC.fg.reads2.fq.gz[/url] &
# mir159a (Col)
wget -4 -q [url=http://bioinfo.mpipz.mpg.de/shoremap/data/software/BC.bg.reads1.fq.gz]http://bioinfo.mpipz.mpg.de/shor ... /BC.bg.reads1.fq.gz[/url]
wget -4 -q [url=http://bioinfo.mpipz.mpg.de/shoremap/data/software/BC.bg.reads2.fq.gz]http://bioinfo.mpipz.mpg.de/shor ... /BC.bg.reads2.fq.gz[/url]

其他数据
除了最基本的测序数据外,我们可能还需要参考基因组,已有的注释数据等
[Shell] 纯文本查看 复制代码
# 参考基因组
wget -4 -q [url=http://bioinfo.mpipz.mpg.de/shoremap/data/software/TAIR10_chr_all.fas]http://bioinfo.mpipz.mpg.de/shor ... /TAIR10_chr_all.fas[/url] &
# 基因注释
wget -4 -q [url=http://bioinfo.mpipz.mpg.de/shoremap/data/software/TAIR10_GFF3_genes.gff]http://bioinfo.mpipz.mpg.de/shor ... IR10_GFF3_genes.gff[/url] &
# SHORE操作的结果数据
wget -4 -q [url=http://bioinfo.mpipz.mpg.de/shoremap/data/software/scoring_matrix_het.txt]http://bioinfo.mpipz.mpg.de/shor ... ring_matrix_het.txt[/url] &

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


[Shell] 纯文本查看 复制代码
# 混池测序结果,双端
BC.fg.reads1.fq.gzBC.fg.reads2.fq.gz
# 背景信息测序结果
BC.bg.reads1.fq.gzBC.bg.reads2.fa.gz
# 拟南芥参考基因组
TAIR10_chr_all.fas
# 拟南芥注释信息
TAIR10_GFF3_genes.gff


以下操作都是基于上述文件进行。
第一步: 序列比对,产生SAM文件
[Shell] 纯文本查看 复制代码
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预测突变位点
为了加快运算速度,可以先转换格式,并排序
[Shell] 纯文本查看 复制代码
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

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

[Shell] 纯文本查看 复制代码
# -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替代
[Shell] 纯文本查看 复制代码
# -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),这个没有找到合适的替代参数

这个没有找到合适的替代参数综上,推荐使用如下命令行
[AppleScript] 纯文本查看 复制代码
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

[AppleScript] 纯文本查看 复制代码
# 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
}

其实对于单个文件而言,可以直接用以下命令
[AppleScript] 纯文本查看 复制代码
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能认识的格式
[AppleScript] 纯文本查看 复制代码
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)
[AppleScript] 纯文本查看 复制代码
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
[AppleScript] 纯文本查看 复制代码
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

第四步:对结果进行注释
[AppleScript] 纯文本查看 复制代码
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: 如何调整阈值,如何辨认结果。






上一篇:Biostar:课程1、2-Linux中生信分析常用命令入门
下一篇:如何从cufflinks拼接的转录本中筛选出lncrna的转录本
回复

使用道具 举报

18

主题

56

帖子

403

积分

管理员

Rank: 9Rank: 9Rank: 9

积分
403
 楼主| 发表于 2017-3-17 20:55:41 | 显示全部楼层
必须要写一个专门的markdown转换程序了,不然发帖太累了
回复 支持 反对

使用道具 举报

您需要登录后才可以回帖 登录 | 立即注册

本版积分规则

QQ|手机版|小黑屋|生信技能树 ( 粤ICP备15016384号  

GMT+8, 2019-12-16 11:36 , Processed in 0.028127 second(s), 24 queries .

Powered by Discuz! X3.2

© 2001-2013 Comsenz Inc.