R语言中实现plink --recode A指令

1、R实现

dat <- read.table("outcome.ped")

name1 <- c("FID", "IID", "PAT", "MAT","SEX", "PHENOTYPE" )
name2 <- dat[,1:6]
name3 <- rbind(name1, name2)

dat <- dat[,-(1:6)]
col2 <- vector()
for (i in 1:(ncol(dat)/2)) {
  temp1 <- vector()
  temp1 <- c(dat[,i * 2 - 1], dat[, i * 2])
  temp2 <- as.data.frame(table(temp1))
  temp2 <- temp2[order(temp2$Freq),]
  dim(temp2)
  if(nrow(temp2) == 1){
    col2 <- c(col2, "0")
  }else
  {
    col2 <- c(col2, as.character(temp2[1,1]))
  }
}

re1 <- data.frame()
for (i in 1:nrow(dat)) {
  vec = vector()
  for (j in 1:length(col2)) {
    count = 0
      if(dat[i,2 * j - 1] == col2[j] & dat[i, 2 * j - 1] != 0){
        count = count + 1}
      if(dat[i,2 * j] == col2[j] & dat[i, 2 * j] != 0){
        count = count + 1}
    vec <- c(vec, count)
  }
  re1 <- rbind(re1, vec)
}

map <- read.table("outcome.map")
col2 <- paste(map$V2,col2,sep = "_")

re1 <- rbind(col2,re1)
re1 <- cbind(name3, re1)
write.table(re1, "test.paw", col.names = F, row.names = F, quote = F)
re1

 

 

2、plink验证

root@PC1:/home/test# ls
outcome.map  outcome.ped  test
root@PC1:/home/test# plink --file outcome --recode A --out temp > /dev/null; rm *log *.nosex
root@PC1:/home/test# ls
outcome.map  outcome.ped  temp.raw  test
root@PC1:/home/test# cat temp.raw
FID IID PAT MAT SEX PHENOTYPE snp1_0 snp2_G snp3_0 snp4_0 snp5_A snp6_0 snp7_A snp8_G
DOR 1 0 0 0 -9 0 0 0 0 1 0 0 1
DOR 2 0 0 0 -9 0 1 0 0 0 0 1 0
DOR 3 0 0 0 -9 0 0 0 0 0 0 1 1
DOR 4 0 0 0 -9 0 0 0 0 0 0 0 2
DOR 5 0 0 0 -9 0 0 0 0 0 0 1 1
DOR 6 0 0 0 -9 0 0 0 0 0 0 2 0