linux 查看vcf文件,Linux生信练习4-vcf
原始数据准备
cd biosoft/bowtie2/bowtie2-2.4.1-linux-x86_64/example/reads/
../../bowtie2 -x ../index/lambda_virus -1 reads_1.fq -2 reads_2.fq | samtools sort -@ 5 -o tmp.bam -
bcftools mpileup -f ../reference/lambda_virus.fa tmp.bam |bcftools call -vm > tmp.vcf
Q1
把突变记录的vcf文件区分成 INDEL和SNP条目
grep -v "^#" tmp.vcf | grep INDEL > tmp.indel.vcf
grep -v "^#" tmp.vcf | grep -v INDEL > tmp.snp.vcf
Q2
统计INDEL和SNP条目的各自的平均测序深度
cut -f 8 tmp.indel.vcf | cut -d";" -f 4 | less -SN
number=$(grep -c gi tmp.indel.vcf)
sum=$(cut -f 8 tmp.indel.vcf | cut -d";" -f 4 | cut -d"=" -f 2 | paste -s -d +|bc)
echo "$sum/$number" | bc
cut -f 8 tmp.snp.vcf | cut -d";" -f 1 | less -SN
number=$(grep -c gi tmp.snp.vcf)
sum=$(cut -f 8 tmp.snp.vcf | cut -d";" -f 1 | cut -d"=" -f 2 | paste -s -d +|bc)
echo "$sum/$number" | bc
Q3
把INDEL条目再区分成insertion和deletion情况
head -1 tmp.indel.vcf | tr "\t" "\n" | cat -n
cat tmp.indel.vcf | awk 'length($4) > length($5){print}' > tmp.indel_del.vcf
cat tmp.indel.vcf | awk 'length($4) < length($5){print}' > tmp.indel_in.vcf
Q4
统计SNP条目的突变组合分布频率
cut -f 10 tmp.snp.vcf | cut -d":" -f 1 | sort | uniq -c
Q5
找到基因型不是 1/1 的条目,个数
grep -v "^#" tmp.vcf | grep -v 1/1 | wc
Q6
筛选测序深度大于20的条目
grep -E -o 'DP=[0-9]+' tmp.vcf | cut -d"=" -f 2 > DP_value
grep -v "^#" tmp.vcf | paste DP_value - | awk '$1>20{print}' | cut -f 2- | less -SN
#或者
grep 'DP=[2-9][0-9]' tmp.vcf #筛选大于20的巧妙思路
Q7
筛选变异位点质量值大于30的条目
cat tmp.vcf | awk '$6>30{print}' | less -SN
Q8
组合筛选变异位点质量值大于30并且深度大于20的条目
grep 'DP=[2-9][0-9]' tmp.vcf | awk '$6>30{print}' | less -SN
Q9
理解DP4=4,7,11,18 这样的字段,就是 Number of high-quality ref-forward , ref-reverse, alt-forward and alt-reverse bases 计算每个变异位点的 AF
grep -E -o 'DP[0-9]=[0-9]+,[0-9]+,[0-9]+,[0-9]+' tmp.vcf > DP4_value
cut -d"=" -f 2 DP4_value | tr "," "\t" | awk '{print ($3+$4)/($1+$2+$3+$4)}'| wc
#注意awk后花括号两边的引号要用单引号,而不是双引号
Q10
在前面步骤的bam文件里面找到这个vcf文件的某一个突变位点的测序深度表明的那些reads,并且在IGV里面可视化bam和vcf定位到该变异位点。
暂时由于电脑原因,igv可视化一段时间回校后学习
head -1 tmp.snp.vcf | tr "\t" "\n" | less -SN
#序列的第1104个碱基,由C变成了A,DP=43
samtools view tmp.sorted.bam | awk '$4<=1104 && length($10)>= 1104-$4{print}' | wc
# 返回50行,应该是vcf根据质量过滤的一小部分。
linux 查看vcf文件,Linux生信练习4-vcf相关推荐
- linux查看lib文件,linux下的 lib文件的学习思考
说到这个LIB文件,先从一个小故障说起. 某日开发说,一台测试用虚机可以PING通SSH不能连了.运维同学就赶紧去查,SSHD_CONFIG配置文件都正确啊,一点错误都没有,那为什么呢? 测试下,不管 ...
- linux查看usb文件,linux lsusb查看USB信息
linux lsusb查看USB信息 linux中lsusb用来显示系统中以及连接到系统的USB总线信息的工具,lsusb会显示驱动和内部连接到你系统的设备,包括PID和VID等,以及简单的设备描述. ...
- Linux查看时间段文件,Linux查看特定时间段内修改过的文件
一.Linux系统日志的一些信息,日志配置文件syslog.conf 系统日志一般都存在/var/log下 常用的系统日志如下: 核心启动日志:/var/log/dmesg 系统报错日志:/var/l ...
- linux 查看.img文件,linux img文件 分区挂载
首先是将制作的img文件比如hd5.img和loop设备建立联系. losetup /dev/loop0 hd5.img 然后用fdisk分区:fdisk /dev/loop0 mkfs.ext4 / ...
- linux查看usb文件,linux下查看usb个数
中午在图书馆闲逛,看到一本<linux指令语法辞典>,随手翻了翻,眼睛落在lsusb 列出所有usb设备那一行,顿时狂喊晕...于是回来试了下: [lbxwz@localhost ~]$ ...
- linux 查看nc文件,linux下nc的使用
今天在饮水思源上闲逛,看到了一个贴子关于Linux下nc命来实现文件传输,进行学习了解了一下.发送端:cattest.txt | nc -l -p 6666或者nc -l -p 6666 < ...
- linux 查看path文件,linux入门之环境变量与文件查找
环境变量 分类 当前 Shell 进程私有用户自定义变量,如上面我们创建的 temp 变量,只在当前 Shell 中有效. Shell 本身内建的变量. 从自定义变量导出的环境变量. declare ...
- linux查看pro文件,Linux下.pro文件的写法简介
1. 注释 从"#"开始,到这一行结束. 2. 指定源文件 SOURCES = *.cpp 对于多源文件,可用空格分开,如:SOURCES = 1.cpp 2.cpp 3.cpp ...
- linux 查看大文件,Linux 查看大文件内容的方法
查看文本文件内容的工具有很多,它们的实现方式和性能各有不同.当我们在大文件或者超大文本文件中查找内容时,考虑到执行效率,我们就要选择合适的方法和工具了. 一.文件大小介绍 英文的字母和标点占用一个字节 ...
- linux查看ld文件,Linux下库文件的设置 (/usr/bin/ld: cannot find -lxxx 的解决办法)
/usr/bin/ld: cannot find -lhdf5 这表示找不到库文件 libhdf5.so,若是其它库文件,则是 cannot find -lxxx 了,其中 xxx 是库文件的名字. ...
最新文章
- opencv判断 线夹角_opencv计算直线的斜率、截距,与水平线弧度值、角度值
- Windows域环境下部署ISA Server 2006防火墙(四)
- 简单介绍android studio中的Logcat
- Winmail邮件服务器
- 【ClickHouse 技术系列】- ClickHouse 中的嵌套数据结构
- redis学习-主从复制Master/slave
- 仿淘宝验证码 php,PHP中仿制 ecshop验证码实例
- python3 模块详细解释_详解Python3中的contextvars模块
- EnableQ在线问卷调查引擎在学校教学教评中的作用
- ARM嵌入式的位绑定原理
- Win11开始菜单怎么改成Win10模式?
- 韦东山第一二期衔接课程内容概要
- 洛谷试炼场 动态规划TG.lv(2)
- 统计笔记3:statistical inference
- 众外媒评析苹果2022秋季新品发布会:说服客户升级iPhone至关重要
- 2020-05-11
- 小红书流量高峰时间段是什么时候?早上发笔记好还是晚上好
- ROS2网络课程资料分享2019.10.26
- Linux中SELinux理解
- 22061周市场回顾