使用SHOREmap做mapping-by-sequencing

简介: 简介SHOREmap可以用来分析传统作图群体(自然系natural strains和分化系,diverged accession杂交,或outcrossing)或近等作图群体(isogenic mapping population, 诱变后代与未诱变亲本进行杂交,即会交,backcrossing)所产生的重测序数据。

简介

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验证

img_a2d24ee799a6c799e76e1818cf24bb7c.png

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

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格式描述了

目录
相关文章
|
2月前
|
算法 数据挖掘
文献解读-Consistency and reproducibility of large panel next-generation sequencing: Multi-laboratory assessment of somatic mutation detection on reference materials with mismatch repair and proofreading deficiency
Consistency and reproducibility of large panel next-generation sequencing: Multi-laboratory assessment of somatic mutation detection on reference materials with mismatch repair and proofreading deficiency,大panel二代测序的一致性和重复性:对具有错配修复和校对缺陷的参考物质进行体细胞突变检测的多实验室评估
26 6
文献解读-Consistency and reproducibility of large panel next-generation sequencing: Multi-laboratory assessment of somatic mutation detection on reference materials with mismatch repair and proofreading deficiency
|
3月前
|
机器学习/深度学习 算法
【文献学习】Channel Estimation Method Based on Transformer in High Dynamic Environment
一种基于CNN和Transformer的信道估计方法,用于在高度动态环境中跟踪信道变化特征,并通过实验结果展示了其相比传统方法的性能提升。
56 0
|
算法 计算机视觉 知识图谱
ACL2022:A Simple yet Effective Relation Information Guided Approach for Few-Shot Relation Extraction
少样本关系提取旨在通过在每个关系中使用几个标记的例子进行训练来预测句子中一对实体的关系。最近的一些工作引入了关系信息
123 0
|
算法 数据可视化 机器人
Object SLAM: An Object SLAM Framework for Association, Mapping, and High-Level Tasks 论文解读
Object SLAM: An Object SLAM Framework for Association, Mapping, and High-Level Tasks 论文解读
84 0
|
机器学习/深度学习 人工智能 自然语言处理
One SPRING to Rule Them Both Symmetric AMR Semantic Parsing and Generation without Complex Pipeline
在文本到AMR解析中,当前最先进的语义解析器集成了几个不同模块或组件的繁琐管道,并利用图重新分类,即在训练集的基础上开发的一组特定内容的启发式方法。
127 0
|
自然语言处理
Re24:读论文 IOT-Match Explainable Legal Case Matching via Inverse Optimal Transport-based Rationale Ext
Re24:读论文 IOT-Match Explainable Legal Case Matching via Inverse Optimal Transport-based Rationale Ext
Re24:读论文 IOT-Match Explainable Legal Case Matching via Inverse Optimal Transport-based Rationale Ext
|
机器学习/深度学习 人工智能 开发框架
IJCAI_2020_Channel Pruning via Automatic Structure Search
1 摘要 通道修剪是压缩深层神经网络的主要方法之一。
183 0
IJCAI_2020_Channel Pruning via Automatic Structure Search
annotationProcessor‘ dependencies won‘t be recognized as kapt annotation processors. Please change
annotationProcessor‘ dependencies won‘t be recognized as kapt annotation processors. Please change
355 0
MapStruct - Ambiguous constructors found
MapStruct - Ambiguous constructors found
455 0
一起谈.NET技术,NHibernate3剖析:Mapping篇之ConfORM实战(4):ManyToMany语义
  ConfORM概述   如果你不熟悉ConfORM请查看前几篇文章,你可以到http://code.google.com/p/codeconform/获取ConfORM最新版本。   在Domain设计中经常使用集合,在.
1202 0