大家好,我是邓飞。 对于vcf文件和plink文件是经常用的文件,对于基因型数据的处理,一般分为: 今天介绍一下vcf文件的三个处理方法: 用到的软件是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
提取后的文件:  搞定!
分割线
大家好,我是邓飞,一个持续分享的农业数据分析师,这里我将自己公众号的干货内容挑重点罗列一下,方便大家阅读和使用。
1,GWAS学习教程(快来领取 | 飞哥的GWAS分析教程更新啦),这个pdf是我将公众号的内容进行了汇总,更方便从头学习GWAS分析,里面配套了数据、代码和讲解,属于干货推荐的Number 1。 2,农学人如何入门数据分析资料汇总(飞哥汇总 | 入门数据分析资源推荐),里面推荐了免费的教程,包括编程、统计和专业书籍。
3,数量遗传学电子书下载(数量遗传学,分享几本书的电子版) 4,R语言电子书线上书籍推荐(学习R语言这几本电子书就够了!)
|