R 脚本实现 xx.hmp.txt格式数据转换为plink格式

 

001、

dir()
dat <- fread("mdp_genotype_test.hmp.txt", header = F, data.table = F)   ## 读入文件

genotype <- t(dat[,12:ncol(dat)])[,-1]
genotype <- apply(genotype,2,function(x){gsub(pattern = "(^.)", replacement = "\\1 ", x)})
genotype <- apply(genotype,2,function(x){gsub(pattern = "N", replacement = "0", x)})

id <- t(dat[1,][,12:ncol(dat)])

ped <- data.frame(fid = id, id = id, pid = 0, mid = 0, gender = 0, phy = -9, genotype)
map <- data.frame(dat$V3, dat$V1, 0, dat$V4)[-1,]

write.table(ped, "result.ped", row.names = F, col.names = F, sep = "\t", quote = F)
write.table(map, "result.map", row.names = F, col.names = F, sep = "\t", quote = F)
dir()

 

 

002、plink验证

root@PC1:/home/test# ls
result.map  result.ped
root@PC1:/home/test# plink --file result --recode
PLINK v1.90b6.24 64-bit (6 Jun 2021)           www.cog-genomics.org/plink/1.9/
(C) 2005-2021 Shaun Purcell, Christopher Chang   GNU General Public License v3
Logging to plink.log.
Options in effect:
  --file result
  --recode

15969 MB RAM detected; reserving 7984 MB for main workspace.
.ped scan complete (for binary autoconversion).
Warning: Variant 424 quadallelic; setting rarest alleles missing.
Performing single-pass .bed write (3093 variants, 281 people).
--file: plink-temporary.bed + plink-temporary.bim + plink-temporary.fam
written.
3093 variants loaded from .bim file.
281 people (0 males, 0 females, 281 ambiguous) loaded from .fam.
Ambiguous sex IDs written to plink.nosex .
Using 1 thread (no multithreaded calculations invoked).
Before main variant filters, 281 founders and 0 nonfounders present.
Calculating allele frequencies... done.
Total genotyping rate is 0.963828.
3093 variants and 281 people pass filters and QC.
Note: No phenotypes present.
--recode ped to plink.ped + plink.map ... done.