Skip to content

Instantly share code, notes, and snippets.

@monogenea
Last active December 2, 2019 14:21
Show Gist options
  • Star 0 You must be signed in to star a gist
  • Fork 0 You must be signed in to fork a gist
  • Save monogenea/fb2afba74367adbe63fa48d6dd9d31e4 to your computer and use it in GitHub Desktop.
Save monogenea/fb2afba74367adbe63fa48d6dd9d31e4 to your computer and use it in GitHub Desktop.
library(snpStats)
load("conversionTable.RData")
pathM <- paste("public/Genomics/108Malay_2527458snps", c(".bed", ".bim", ".fam"), sep = "")
SNP_M <- read.plink(pathM[1], pathM[2], pathM[3])
pathI <- paste("public/Genomics/105Indian_2527458snps", c(".bed", ".bim", ".fam"), sep = "")
SNP_I <- read.plink(pathI[1], pathI[2], pathI[3])
pathC <- paste("public/Genomics/110Chinese_2527458snps", c(".bed", ".bim", ".fam"), sep = "")
SNP_C <- read.plink(pathC[1], pathC[2], pathC[3])
# Ensure == number of markers across the three populations
if(ncol(SNP_C$genotypes) != ncol(SNP_I$genotypes)){
stop("Different number of columns in input files detected. This is not allowed.")
}
if(ncol(SNP_I$genotypes) != ncol(SNP_M$genotypes)){
stop("Different number of columns in input files detected. This is not allowed.")
}
# Merge the three SNP datasets
SNP <- SNP_M
SNP$genotypes <- rbind(SNP_M$genotypes, SNP_I$genotypes, SNP_C$genotypes)
colnames(SNP$map) <- c("chr", "SNP", "gen.dist", "position", "A1", "A2") # same for all three
SNP$fam<- rbind(SNP_M$fam, SNP_I$fam, SNP_C$fam)
# Rename SNPs present in the conversion table into rs IDs
mappedSNPs <- intersect(SNP$map$SNP, names(conversionTable))
newIDs <- conversionTable[match(SNP$map$SNP[SNP$map$SNP %in% mappedSNPs], names(conversionTable))]
SNP$map$SNP[rownames(SNP$map) %in% mappedSNPs] <- newIDs
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment