1. 程式人生 > 其它 >R語言中實現plink --recode A指令

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