如何将基因型数据转为 012 格式

生信菜鸟团

发布于 2020-06-04 10:12:05

10K4

发布于 2020-06-04 10:12:05

文章被收录于专栏：生信菜鸟团生信菜鸟团

最近碰到将基因型数据转为 012 格式的需求，就顺手总结了一些方法和大家分享，要是有更方便的法子欢迎大家多多补充~

012 格式一般就是：0:HOM_REF, 1:HET 2: HOM_VAR

plink1.9

直接用 plink 进行转换：

plink1.9 --bfile test.genotypes_no_missing_IDs --recodeA --out test

这一行命令会生成 .raw 文件，即为 012 矩阵，缺失的点标记为 NA：

head test.raw| cut -d" " -f 1-9

FID IID PAT MAT SEX PHENOTYPE rs74512038_T rs2977670_G rs4951859_C
19 19 0 0 1 1 0 0 0
26 26 0 0 1 2 0 0 0
39 39 0 0 1 2 0 0 1
44 44 0 0 2 1 0 0 0
46 46 0 0 1 2 0 0 0
47 47 0 0 1 1 1 1 1
49 49 0 0 2 1 1 1 1
52 52 0 0 2 1 0 0 0
55 55 0 0 1 2 1 1 1

下一步将 2-6 列删除，并格式化 RS ID：

cat test.raw | cut -d" " -f1,7- | sed 's/_[A-Z]//g' >genotype.txt

最后用 awk 进行矩阵转置：

awk '{i=1;while(i <= NF){col[i]=col[i] $i " ";i=i+1}} END {i=1;while(i<=NF){print col[i];i=i+1}}' genotype.txt | sed 's/[ \t]*$//g' > genotype_transpose.txt

head genotype_transpose.txt

FID 19 26 39 44 46 47 49 52 55 56 57 58 59 65 70 73 75 76 77 79 80 92 97 99 107 108 112 113 127 129 135 138 144 150 151 157 170 180 181 182 193 194 203 222 227 228 229 230 251 259 267 274 282 284 289 305 308 309 311 318 346 357 361 369 372 394 417 444 450 456 457 458 460 462 463 467 478 481 483 487 490 496 512
rs74512038 0 0 0 0 0 1 1 0 1 0 1 1 1 2 1 0 0 1 0 0 0 0 0 1 0 0 1 0 1 1 0 0 0 1 1 0 1 0 0 0 0 1 0 0 0 0 1 0 1 1 1 1 0 0 1 0 0 2 1 1 0 0 0 0 0 1 0 0 0 1 0 0 0 1 1 0 0 0 0 0 1 2 0
rs2977670 0 0 0 0 0 1 1 0 1 0 1 1 1 2 1 0 0 1 0 0 0 0 0 1 0 0 1 0 1 1 0 0 0 1 1 0 1 0 0 0 0 1 0 0 0 0 1 0 1 1 1 1 0 0 1 0 0 2 1 1 0 0 0 0 0 1 0 0 0 1 0 0 0 1 1 0 0 0 0 0 1 2 0
rs4951859 0 0 1 0 0 1 1 0 1 0 1 1 1 2 1 0 0 1 0 0 0 0 1 1 1 0 1 0 1 1 0 0 0 0 1 0 1 1 0 0 0 1 0 0 0 1 1 1 1 2 1 1 0 0 1 0 0 2 1 1 0 0 0 0 1 1 1 0 0 1 0 1 1 1 1 0 0 0 0 0 2 2 0
rs79010578 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 1 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 1 0 1 0 1 0 0 0 0 1 0 0 0 0 1 0 0 0 0 1 0 1 0 0 0 0 1 1 0 0 0 0 0 0 0 2 0 0
rs3094315 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 1 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 1 0 1 0 1 0 0 0 0 1 0 0 0 0 1 0 0 0 0 1 0 1 0 0 0 0 1 1 0 0 0 0 0 0 0 1 0 0
rs3131972 0 0 1 0 0 1 0 0 0 0 1 1 1 1 1 0 0 1 0 0 0 0 1 0 1 0 0 0 1 1 0 0 0 0 2 0 0 1 0 0 0 1 0 0 0 1 1 1 1 2 0 1 0 0 1 0 0 1 1 1 0 0 1 0 1 1 1 0 0 1 0 1 1 0 1 0 0 0 0 0 1 2 0
rs2073813 0 0 1 0 0 1 0 0 0 0 1 1 1 2 1 0 0 1 0 0 0 0 1 1 1 0 1 0 1 1 0 0 0 0 2 0 0 1 0 0 0 1 0 0 0 1 1 1 1 2 0 1 0 0 1 0 0 1 1 1 0 0 1 0 1 1 1 0 0 1 0 1 1 0 1 0 0 0 0 0 2 2 0
rs3131969 0 0 1 0 0 1 0 0 0 0 1 1 1 2 1 0 0 1 0 0 0 1 1 1 1 0 1 1 1 1 0 0 0 0 2 0 0 1 0 0 0 1 0 0 0 1 1 1 1 2 0 1 0 0 1 0 0 1 1 1 0 0 1 0 1 1 1 0 0 1 0 1 1 0 1 0 0 0 0 0 2 2 0
rs3131968 0 0 1 0 0 1 0 0 0 0 1 1 1 2 1 0 0 1 0 0 0 1 1 1 1 0 1 1 1 1 0 0 0 0 2 0 0 1 0 0 0 1 0 0 0 1 1 1 1 2 0 1 0 0 1 0 0 1 1 1 0 0 1 0 1 1 1 0 0 1 0 1 1 0 1 0 0 0 0 0 2 2 0

SNP 位置信息可从 .map 文件获取：

plink1.9 --bfile test.genotypes_no_missing_IDs --recode --out test

head test.map

1    rs74512038    0    778597
1    rs2977670    0    788511
1    rs4951859    0    794299
1    rs79010578    0    800909
1    rs3094315    0    817186
1    rs3131972    0    817341
1    rs2073813    0    818161
1    rs3131969    0    818802
1    rs3131968    0    818812
1    rs3131967    0    818954

.map 文件共 4 列，每行为 SNP 位点，第一列为 SNP 所在的染色体，第二列为 SNP ID，第三列为 SNP 的遗传距离，没有实际值则用 0 填充，第四列为 SNP 的位置。可根据不同软件所需的格式提取数据：

echo "snpid chr pos" >snpsloc.txt
awk '{print $2,"chr"$1,$4}' test.map >>snpsloc.txt

vcftools

除了用 plink 转换外，还可以用 vcftools 得到 012 格式的矩阵。

用 plink 将二进制文件转为 VCF 格式：

plink1.9 --bfile test.genotypes_no_missing_IDs --recode vcf-iid --out test.genotypes_no_missing_IDs

用 vcftools 生成 012 矩阵：

vcftools --vcf test.genotypes_no_missing_IDs.vcf --012 --out snp_matrix

--012 参数会生成三个文件：

•后缀为 .012 的 012 基因型矩阵文件，每行为样本，列为基因型。缺失基因型用 -1 表示。

head snp_matrix.012 | cut -f 1-9

0    0    0    0    0    0    0    0    0
1    0    0    0    0    0    0    0    0
2    0    0    1    1    1    1    1    1
3    0    0    0    0    0    0    0    0
4    0    0    0    0    0    0    0    0
5    1    1    1    0    0    1    1    1
6    1    1    1    0    0    0    0    0
7    0    0    0    0    0    0    0    0
8    1    1    1    0    0    0    0    0
9    0    0    0    0    0    0    0    0

可以看到第一列为索引，并不是真实的样本名，在接下来的步骤中需要替换。

•后缀为 .012.indv 的样本信息文件。

head snp_matrix.012.indv

•后缀为 .012.pos 的 SNP 位点坐标文件。

head snp_matrix.012.pos

1    778597
1    788511
1    794299
1    800909
1    817186
1    817341
1    818161
1    818802
1    818812
1    818954

所以根据 vcftools 输出的 .012 矩阵，我们还需要给他加上 SNP ID 和样本 ID。首先将第一列替换为样本名：

awk '{$1=null;print $0}' snp_matrix.012 > tmp.txt
paste -d" " snp_matrix.012.indv tmp.txt | sed 's/  / /g' > snp_matrix_indv.txt

从 VCF 文件中提取 SNP ID：

echo "ID" > snpid.txt
grep -v "^#" test.genotypes_no_missing_IDs.vcf | cut -f3 >>snpid.txt

cat snpid.txt| tr "\n" " " >snpid2.txt
echo "" >> snpid2.txt

合并文件：

cat snpid2.txt snp_matrix_indv.txt> genotype.txt

最后同样进行矩阵转置即可：

awk '{i=1;while(i <= NF){col[i]=col[i] $i " ";i=i+1}} END {i=1;while(i<=NF){print col[i];i=i+1}}' genotype.txt | sed 's/[ \t]*$//g' > genotype_transpose.txt

组合 Linux 命令行工具

也可以直接根据 VCF 提取相应列，转换格式：

cat test.genotypes_no_missing_IDs.vcf | grep -v "^##" | cut -f3,10- | sed 's/0\/0/0/g' | sed 's/1\/1/2/g' | sed 's/0\/1/1/g' | sed 's/1\/0/1/g'| sed 's/\.\/\./-1/g' | tr "\t" " " >genotype_transpose.txt

生信技能树目前已经公开了三个生信知识库，记得来关注哦~

每周文献分享

https://www.yuque.com/biotrainee/weeklypaper

肿瘤外显子分析指南

https://www.yuque.com/biotrainee/wes

生物统计从理论到实践

https://www.yuque.com/biotrainee/biostat

本文参与腾讯云自媒体分享计划，分享自微信公众号。

原始发表：2020-05-31，如有侵权请联系 cloudcommunity@tencent.com 删除

https

网络安全

本文分享自生信菜鸟团微信公众号，前往查看

如有侵权，请联系 cloudcommunity@tencent.com 删除。

本文参与腾讯云自媒体分享计划，欢迎热爱写作的你一起参与！

https

网络安全

登录后参与评论

0 条评论

热度