Nothing
################################################################################
# Create a kinship matrix.
# Daniel Gatti
# Dan.Gatti@jax.org
# Oct. 27, 2012
################################################################################
# Create a kinship matrix based on allele sharing at the markers.
# Arguments: geno: character matrix containing allele calls for each sample and
# marker. num.snps x num.samples
kinship.alleles = function(geno) {
K = matrix(0, ncol(geno), ncol(geno), dimnames = list(colnames(geno),
colnames(geno)))
for(i in 1:ncol(geno)) {
K[,i] = colMeans(geno[,i] == geno)
} # for(i)
return(K)
} # kinship.alleles()
# This function takes the 8 state founder probabilities and calculates
# the cosine of the angle between each sample.
# Arguments: probs: 3D array containing founder haplotype contributions.
# num.samples x num.founders x num.snps.
# snps: data.frame containing the locations of the markers.
# SNP IDs, chromosomes, Mb & cM locations in columns
# 1:4. Optional.
# bychr: boolean that indicates if a separate kinship
# matrix should be made for each chromosome.
kinship.probs = function(probs, snps, bychr = FALSE) {
K = NULL
if(bychr) {
# Make a list of kinship matrices, one for each chromosome, in
# which the kinship on one chromosome uses the markers on the
# remaining chromosomes.
if(missing(snps)) {
stop(paste("snps cannot be missing if bychr is T. Please supply",
"the marker locations in the snps argument."))
} # if(missing(snps))
snps = snps[snps[,1] %in% dimnames(probs)[[3]],]
probs = probs[,,dimnames(probs)[[3]] %in% snps[,1]]
probs = probs[,,match(snps[,1], dimnames(probs)[[3]])]
snps[,2] = as.character(snps[,2])
chr = as.character(unique(snps[,2]))
Kbychr = as.list(chr)
names(Kbychr) = chr
keep = as.list(1:length(chr))
for(i in 1:length(chr)) {
print(paste("CHR", chr[i]))
Kbychr[[i]] = matrix(0, dim(probs)[1], dim(probs)[1], dimnames =
list(dimnames(probs)[[1]], dimnames(probs)[[1]]))
keep[[i]] = which(snps[,2] == chr[i])
res = .C(C_kinship,
dims = as.integer(dim(probs[,,keep[[i]]])),
probs = as.double(probs[,,keep[[i]]]),
K = as.double(Kbychr[[i]]))
# Restore the dimensions and dimnames. We have to multiply
# by the number of SNPs so that the means work out below.
Kbychr[[i]] = matrix(res$K * length(keep[[i]]), dim(probs)[1],
dim(probs)[1], dimnames =
list(dimnames(probs)[[1]], dimnames(probs)[[1]]))
} # for(i)
K = as.list(chr)
names(K) = chr
num.snps = sapply(keep, length)
for(i in 1:length(chr)) {
K[[i]] = matrix(0, dim(probs)[1], dim(probs)[1], dimnames =
list(dimnames(probs)[[1]], dimnames(probs)[[1]]))
for(j in chr[chr[i] != chr]) {
K[[i]] = K[[i]] + Kbychr[[j]]
} # for(j)
# Divide by the number of SNPs on the other (!= i) chromosomes.
K[[i]] = K[[i]] / sum(num.snps[chr != chr[i]])
} # for(i)
} else {
# Use all chromosomes to make a single kinship matrix.
K = matrix(0, dim(probs)[1], dim(probs)[1], dimnames = list(
dimnames(probs)[[1]], dimnames(probs)[[1]]))
probs[is.nan(probs) | is.na(probs) | is.infinite(probs)] = 0
res = .C(C_kinship,
dims = as.integer(dim(probs)),
probs = as.double(probs),
K = as.double(K))
# Restore the dimensions and dimnames.
K = matrix(res$K, dim(probs)[1], dim(probs)[1], dimnames = list(dimnames(probs)[[1]],
dimnames(probs)[[1]]))
} # else
return(K)
} # kinship.probs()
Any scripts or data that you put into this service are public.
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.