shell高级技巧:提取vcf文件中一个contig

简介: 这是一个很小众的需求。大部分变异检测都是基于组装质量比较高的基因组,而不是那种初步拼接的contig。由于初步拼接的参考序列通常会有成千上万个contig序列,也就导致在VCF的头文件的##contig=部分会有成千上万个contig,将这个文件加载到IGV时, IGV会去解析VCF,这将会是非常缓慢的过程,最好的策略就是只提取其中一个contig查看。

这是一个很小众的需求。大部分变异检测都是基于组装质量比较高的基因组,而不是那种初步拼接的contig。

由于初步拼接的参考序列通常会有成千上万个contig序列,也就导致在VCF的头文件的##contig=<ID=xxx,length=xxx>部分会有成千上万个contig,将这个文件加载到IGV时, IGV会去解析VCF,这将会是非常缓慢的过程,最好的策略就是只提取其中一个contig查看。

一开始,我想到的是这个方法,直接用bcftools去提取目标contig

bcftools view some_variants.bcf contig_1

但是,你会发现HEADER部分还是会保留原来的信息。怎么办呢?

我们这里会用到一个Shell高级技巧--subshell。举个例子,比如说同时查看一个文件的头10行和后10行.

(head -n 10 ; tail -n10 ) < ~/referece/annotaiton/TAIR10/TAIR10_GFF3_genes.bed

这里的()就是一个subshell, 它的作用就是同时执行()中的所有命令。

用到此处,就是同时完成提取VCF的header部分并只保留目标contig和提取目标区间的所有记录这两个任务,然后保存为vcf.gz格式。

(bcftools view -h some_variants.bcf | grep -v "##contig" | sed "/^##reference/a $(bcftools view -h some_variants.bcf | grep -w "ID=contig_1")" ; bcftools view -H some_variants.bcf contig_1 ) | bcftools view -Oz -o contig_1.vcf.gz

上面是一个非常复杂的shell命令,希望你能看懂。

这里面重复出现了两个变量,一个是some_variants.bcf,一个是contig_1, 我们可以将其用$1$2代替,然后添加到你的~/.bashrc中,如下所示

# function
subvcf()
{
    (bcftools view -h $1 | grep -v "##contig" | sed "/^##reference/a $(bcftools view -h $1| grep -w "ID=$2")" ;  bcftools view -H $1 $2 ) | bcftools view -Oz -o $2.vcf.gz && bcftoos index $2.vcf.gz
}

最后用source ~/.bashrc添加刚才的函数到环境变量中,最后就是subvcf some_variants.bcf contig_1 代替之前的长串命令。

目录
相关文章
|
7月前
|
Shell Android开发
Android系统 adb shell push/pull 禁止特定文件
Android系统 adb shell push/pull 禁止特定文件
594 1
|
7月前
|
人工智能 机器人 Shell
【shell】文件读写及read用法
【shell】文件读写及read用法
|
7月前
|
Shell Linux API
【Shell 命令集合 备份压缩 】Linux 解压缩文件 unzip命令 使用指南
【Shell 命令集合 备份压缩 】Linux 解压缩文件 unzip命令 使用指南
262 0
|
7月前
|
缓存 Shell Linux
【Shell 命令集合 链接器(linker)工具】Linux ld命令 将目标文件与库链接为可执行文件或库文件
【Shell 命令集合 链接器(linker)工具】Linux ld命令 将目标文件与库链接为可执行文件或库文件
233 0
|
7月前
|
Shell Linux C语言
【Shell 命令集合 系统设置 】Linux 创建Kickstart文件mkkickstart命令 使用指南
【Shell 命令集合 系统设置 】Linux 创建Kickstart文件mkkickstart命令 使用指南
73 0
|
7月前
|
Shell Linux 编译器
【Shell 命令集合 备份压缩 】Linux 提取zip压缩文件的详细信息 zipinfo命令 使用指南
【Shell 命令集合 备份压缩 】Linux 提取zip压缩文件的详细信息 zipinfo命令 使用指南
119 0
|
7月前
|
存储 Shell Linux
【Shell 命令集合 备份压缩 】Linux 解码uuencode编码的文件 uudecode 命令 使用指南
【Shell 命令集合 备份压缩 】Linux 解码uuencode编码的文件 uudecode 命令 使用指南
84 0
|
7月前
|
监控 Shell Linux
【Shell 命令集合 系统管理 】Linux 自动轮转(log rotation)日志文件 logrotate命令 使用指南
【Shell 命令集合 系统管理 】Linux 自动轮转(log rotation)日志文件 logrotate命令 使用指南
166 0
|
2月前
|
Shell
Shell 文件包含
10月更文挑战第5天
35 4
|
7月前
|
Java 关系型数据库 MySQL
Elasticsearch【问题记录 01】启动服务&停止服务的2类方法【及 java.nio.file.AccessDeniedException: xx/pid 问题解决】(含shell脚本文件)
【4月更文挑战第12天】Elasticsearch【问题记录 01】启动服务&停止服务的2类方法【及 java.nio.file.AccessDeniedException: xx/pid 问题解决】(含shell脚本文件)
691 3