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
标签:plink,--,PC1,test,recode,DOR,home,outcome 来源: https://www.cnblogs.com/liujiaxin2018/p/15501393.html