使用LUMPY检测结构变异

简介: LUMPY是一款基于概率框架检测结构变异(structure variants)的软件, 它根据read-pair, split-read, read-depth和其他先验知识寻找基因组上可能的结构变异。

LUMPY是一款基于概率框架检测结构变异(structure variants)的软件, 它根据read-pair, split-read, read-depth和其他先验知识寻找基因组上可能的结构变异。

软件在编译的时候会先安装HTSLIB,根据Makefile, 需要预先安装好curl和zlib. 此外还推荐安装Python的Pysam和Numpy,Samtools(0.1.18+),SAMBLASTER(0.1.19+), sambamba。安装完成之后可以用测试数据进行测试, 数据下载地址为http://layerlab.org/lumpy/data.tar.gz, 要存放在/lumpy-sv/data

lumpy基于paired-end reads比对后得到的三类信息推断SV,局部异常的测序深度,不一致(discordant)的联配和断裂的联配(split-read alignment)。局部异常的测序深度比较容易理解,平均30X测序的地方,如果深度大于100X,意味着存在着拷贝数变异,如果深度程度非常低,可能意味着这里存在 大片段缺失。不一致的联配和断裂的联配能够提供的信息更多,如果基因组一个区域齐刷刷的截断(如下图),就意味着这个区域可能存在插入/缺失。当然也有其他可能,当两个read在不同链或者不同染色体时,可能是易位或倒置。

img_33c2dd603235cf965dc44d9f060957ff.jpe
read截断

因此,运行lumpy需要预先整理出不一致的短读以及断裂的联配, 如下是lumpy提供的数据预处理方式

# Align the data
bwa mem -R "@RG\tID:id\tSM:sample\tLB:lib" human_g1k_v37.fasta sample.1.fq sample.2.fq \
    | samblaster --excludeDups --addMateTags --maxSplitCount 2 --minNonOverlap 20 \
    | samtools view -S -b - \
    > sample.bam

# Extract the discordant paired-end alignments.
samtools view -b -F 1294 sample.bam > sample.discordants.unsorted.bam

# Extract the split-read alignments
samtools view -h sample.bam \
    | scripts/extractSplitReads_BwaMem -i stdin \
    | samtools view -Sb - \
    > sample.splitters.unsorted.bam

# Sort both alignments
samtools sort sample.discordants.unsorted.bam sample.discordants
samtools sort sample.splitters.unsorted.bam sample.splitters

让人感兴趣的是-F 1294用来提取不一致的联配,用samtools flags 1294可以发现1294表示"PROPER_PAIR,UNMAP,MUNMAP,SECONDARY,DUP",带上-F意味着以上这些标记在我们筛选的联配记录中都不会出现,也就意味着筛选的记录要符合下面要求

  • 不能是PROPER_PAIR: 就是比对工具认为都正确比对到基因组上,在同一条染色体,在同一条链的情况,常见的就是83,147和99,163
  • 不能是UNMAP和MUNMAP,也就是配对的短读至少有一个能够比对到参考基因组上
  • 也不能是SECONDARY, 也就是他必须是主要联配
  • 光学重复,DUP, 就更加不能要了

于是,经过上一步,那就得到了包含所有数据的sample.bam,不一致的联配sample.discordants.bam 和断裂联配sample.splitters.bam, 使用作者封装好的调用函数进行结构变异检测。

lumpyexpress \
    -B sample.bam \
    -S sample.splitters.bam \
    -D sample.discordants.bam \
    -o sample.vcf

在得到的结构变异基础上,作者推荐是用SVTyper进行基因型确定。

提高准确度的方法:根据先验剔除已知低复杂区域和高覆盖的区域,见参考资料的lumpy-sv教程。

参考资料

目录
相关文章
|
数据可视化 Linux 开发者
Processing有哪些常见的用途或者优势呢
Processing有哪些常见的用途或者优势呢
659 0
|
C++
如何使用MACS进行peak calling
MACS2是peak calling最常用的工具。 callpeak用法 这是MACS2的主要功能,因为MACS2的目的就是找peak,其他功能都是可有可无,唯独callpeak不可取代。
4304 0
如何解决 conda install 库时报错:The environment is inconsistent, please check the package plan carefully
如何解决 conda install 库时报错:The environment is inconsistent, please check the package plan carefully
如何解决 conda install 库时报错:The environment is inconsistent, please check the package plan carefully
|
监控 安全 Linux
在Linux中,如何查看和审计系统日志文件以检测异常活动?
在Linux中,如何查看和审计系统日志文件以检测异常活动?
|
搜索推荐 Docker 容器
生信分析代码之前还好好的,怎么就报错了 Error in Ops. data. frame(guide_loc, panel_loc) :'==' only defined for equally-sized data frames
执行 `DimPlot` 函数时遇到错误 `;Error in Ops. data. frame(g guides_loc, panel_loc) : '==' only defined for equally-sized data frames`。解决方案和办法
2298 0
生信分析代码之前还好好的,怎么就报错了 Error in Ops. data. frame(guide_loc, panel_loc) :'==' only defined for equally-sized data frames
|
Linux 数据处理
Linux中的pr命令:数据格式化与打印的艺术
`pr`命令是Linux下用于文本格式化的工具,擅长分页、设置页眉页脚及列宽,方便打印和阅读。它可以处理文件、管道输入,常用参数如 `-h` 设定页眉,`-t` 设置页脚,`-l` 控制每页行数,`-w` 设置列宽。例如,`pr -h "标题" -t "页码:%d" -l 2 file.txt` 可以将文本文件格式化并添加定制的页眉页脚。结合其他命令使用能增强文本处理能力。记得测试输出,了解详细帮助可查阅`man pr`。
|
Python
干货文:在 Mac 中卸载 Python 的方式
干货文:在 Mac 中卸载 Python 的方式
3024 1
|
机器学习/深度学习 数据可视化 数据挖掘
跟着Nature Genetics学数据分析:nucmer+lastz+svum流程全基因组比对鉴定CNV
跟着Nature Genetics学数据分析:nucmer+lastz+svum流程全基因组比对鉴定CNV
|
iOS开发 MacOS
ubuntu22/20 安装macos 主题桌面
ubuntu22/20 安装macos 主题桌面
651 0