分享

vcf 文件如何修改染色体修改样本名称提取样本

 育种数据分析 2023-10-10 发布于河南

大家好,我是邓飞。

对于vcf文件和plink文件是经常用的文件,对于基因型数据的处理,一般分为:

  • 数据质控
  • 数据提取
  • 染色体修改名称
  • 样本修改名称

今天介绍一下vcf文件的三个处理方法:

  • 1,染色体修改
  • 2,样本名称修改
  • 3,样本提取

用到的软件是bcftools,用到的系统是Linux

1. 数据描述

数据使用GWAS-Cookbook中的GWASdat1中的数据,将数据变为vcf格式。

相关基因型数据,可以在这里下载:快来领取 | 飞哥的GWAS分析教程更新啦

将plink的二进制文件,变为vcf的代码:

plink --bfile ../../gwas_cookbook/GWAS-dat1/1_QC_GWAS/HapMap_3_r3_1 --recode vcf-iid --out test1

vcf预览:

2. vcf文件修改染色体名称

整理染色体修改对应关系:

awk '{print $1}' HapMap_3_r3_1.bim |sort |uniq >chr.txt
awk '{print $1,"chr"$1}' chr.txt >chr_rename.txt

整理好的对应关系:

$ cat chr_rename.txt 
1 chr1
10 chr10
11 chr11
12 chr12
13 chr13
14 chr14
15 chr15
16 chr16
17 chr17
18 chr18
19 chr19
2 chr2
20 chr20
21 chr21
22 chr22
23 chr23
25 chr25
3 chr3
4 chr4
5 chr5
6 chr6
7 chr7
8 chr8
9 chr9

代码:

 bcftools annotate --rename-chrs chr_rename.txt -Oz -o test2.vcf.gz test1.vcf 
gzip -d test2.vcf.gz 
less -S test2.vcf 

修改后的结果:

可以看到,已经修改过了。

3. 修改样本的名称

样本对应关系txt文件整理:

awk '{print $2}' HapMap_3_r3_1.fam |head >sname.txt
awk '{print $1,"new_"$1}' sname.txt >sname_change.txt

对应关系文件内容:

$ head sname_change.txt 
NA06989 new_NA06989
NA11891 new_NA11891
NA11843 new_NA11843
NA12341 new_NA12341
NA12739 new_NA12739
NA10850 new_NA10850
NA06984 new_NA06984
NA12877 new_NA12877
NA12275 new_NA12275
NA06986 new_NA06986

代码:

bcftools reheader -s sname_change.txt test1.vcf -o test3.vcf.gz

修改后的vcf:

4. 提取vcf样本

代码:

bcftools view -S sname.txt test1.vcf -Ov > test4.vcf

提取后的文件:

搞定!

想要更好的学习和交流,快来加入飞哥的知识星球,这是一个生物统计+数量遗传学+GWAS+GS的社区,在这里你可以向飞哥提问、帮你定学习计划、跟着飞哥一起做实战项目,冲冲冲。点击这里加入吧:飞哥的学习圈子


分割线


大家好,我是邓飞,一个持续分享的农业数据分析师,这里我将自己公众号的干货内容挑重点罗列一下,方便大家阅读和使用。

1,GWAS学习教程(快来领取 | 飞哥的GWAS分析教程更新啦),这个pdf是我将公众号的内容进行了汇总,更方便从头学习GWAS分析,里面配套了数据、代码和讲解,属于干货推荐的Number 1。

2,农学人如何入门数据分析资料汇总(飞哥汇总 | 入门数据分析资源推荐),里面推荐了免费的教程,包括编程、统计和专业书籍。

3,数量遗传学电子书下载(数量遗传学,分享几本书的电子版

4,R语言电子书线上书籍推荐(学习R语言这几本电子书就够了!

    转藏 分享 献花(0

    0条评论

    发表

    请遵守用户 评论公约

    类似文章