Jasieniuk:Implementation of a method to calculate microsatellite genotype distance in polyploids

Update

Note, April 2010: I'm working on developing an R package with these and other tools for polyploid microsatellite data. I have restructured the functions somewhat so that they will be more versatile (allow for other measures of distance, and import the genotypes as an R object so that other things can be done with them). I've written a function that will read the data more directly from the format produced by ABI GeneMapper and one that can process multiple loci at once and calculate the mean distances. Missing data is also allowed. There is also a function to export data to be read by Structure 2.3.3.

The functions are in a text file: Media:Polysatfun100409.R.txt

Documentation of the functions: Media:Polysatfun100409.Rd.txt

Overview

This is an implementation of the method of Bruvo et al. (2004). It calculates a pairwise measure of genetic distance between individuals at a single microsatellite locus. The allele copy number does not need to be known and the individuals do not even need to be the same ploidy level. Its primary advantage over band-sharing measures of genetic distance is that it allows for mutation, such that alleles aren't considered to be qualitatively different, but instead the difference in repeat number between two alleles is taken into account. The paper offers two options for dealing with virtual alleles (when one individual has more alleles than the other). My implementation uses the simpler option and sets virtual alleles as equal to infinity.

This is my first adventure into the R programming environment, so I apologize for any programming faux pas.

Arguments and Output

The function genotype.distance takes as arguments two single locus microsatellite genotypes. In the code these are genotype1 and genotype2. Each genotype should be a vector containing all of the alleles for a given individual at a given locus. The vectors do not need to be the same length, and should only be as long as needed to contain the alleles. The alleles should be expressed in repeat number rather than nucleotides (DNA fragment length). Because it is the difference in allele size that is measured, you can simply divide the nucleotide length by the repeat size (e.g. divide by two if you have dinucleotide repeats).

The output is a single number that is the distance between genotypes. In the paper, this is d_l.

Code for genotype.distance

genotype.distance <- function(genotype1, genotype2) {

require(combinat) # for permn function

if(identical(genotype1, genotype2)) {dist <- 0} else { if(length(genotype1) >= length(genotype2)) {genotypeL <- genotype1; genotypeS <- genotype2} else {genotypeL <- genotype2; genotypeS <- genotype1}

# if genotypes are identical, just return a distance of zero without doing the rest of the calculation

# whichever genotype has more alleles, make this genotypeL (long) and the other genotypeS (short)

kl <- length(genotypeL) # sets the ploidy level for this genotype comparison

ks <- length(genotypeS) # number of alleles in the shorter genotype

allele.distances <- array(0 , c(kl,ks)) # Create an empty matrix to contain the distances between alleles

for(n in 1:kl) { for(m in 1:ks) {allele.distances[n,m] <- genotypeL[n] - genotypeS[m]}} # fills the array with the differences in allele length

geometric.distances <- array(1 - 2^-abs(allele.distances) , c(kl,ks)) # geometric transformation based on mutation probabilities

#Next, find the minimum distance sum among all permutations

column <- 1:ks # an index of all columns (genotypeS alleles)

row <- 1:kl # an index of all rows (genotypeL alleles)

combinations <- combn(row, ks, FUN = NULL, simplify=FALSE) # all combinations of alleles in genotypeL that can be matched to non-virtual alleles in genotypeS

permutations <- permn(ks) # all possible orders that alleles within these combinations can go in

mindist <- Inf # this variable will store the minimum sum encountered so far.

for(i in 1:length(combinations)) { # the loop to go through every possible sum of compatible allele comparisons

rowcomb <- combinations[[i]] # choose one combination of rows for this round

for(l in 1:length(permutations)){ # go through all orders of this combinations of rows

sum <- 0 # this is si, the sum of allele comparisons

for(j in 1:ks) {sum <- sum + geometric.distances[rowcomb[permutations[[l]][j]],column[j]]} # the loop to calculate the sum for this permutation

if(sum < mindist) {mindist <- sum} # is this the minimum sum found so far?

}}

dist

}

Processing a file of genotypes

Explanation, arguments, and output

The function defined below utilizes the genotype.distance function and produces a symmetrical matrix of genotype distances, given a file of individual genotypes at a given locus. It takes three arguments:

• infile is the filepath of a tab-delimited text file containing the genotypes. The file should contain a header row. The first column should be labels for the individuals/samples, the second column should contain the number of alleles that the individual has at this locus, and columns 3-22 should be the alleles, expressed as the DNA fragment length. If, as is likely, a genotype contains less than 20 alleles, the allele fields on the right can be left blank or filled with zeros.
• outfile is the filepath where you want your output to be written. This will be a tab-delimited text file containing a labeled matrix of genotype distances. (This matrix is also returned by the function.)
• repeat.length is the length of the microsatellite repeat at this locus, for example 2 if this is a dinucleotide repeat.

entire.locus can be run individually on several different loci, then the matrices can be averaged.

I work with a dodecaploid organism, so the calculations can take awhile. For this reason, the function prints each pair of individuals as it works so that the user can monitor the progress.

Code for entire.locus

entire.locus <- function(infile,outfile,repeat.length){

usatdata<-read.table(infile, header=TRUE, sep = "\t", fill=TRUE)

labels <- row.names(usatdata)

genotypes <- list(0) #this will contain the genotype vectors

count <- dim(usatdata)[1] #the number of genotypes being read

for(i in 1:count){ #extract genotype vectors from the data table

genotypes[[i]] <- c(usatdata[i,2],usatdata[i,3],usatdata[i,4],usatdata[i,5],usatdata[i,6],usatdata[i,7],usatdata[i,8],usatdata[i,9],usatdata[i,10],usatdata[i,11],usatdata[i,12],usatdata[i,13],usatdata[i,14],usatdata[i,15],usatdata[i,16],usatdata[i,17],usatdata[i,18],usatdata[i,19],usatdata[i,20],usatdata[i,21])/repeat.length

length(genotypes[[i]]) <- usatdata[i,1]

}

distances<- array(0,c(count,count)) #array for containing genotype distances

for(m in 1:count){for(n in m:count){thisdistance <-genotype.distance(genotypes[[m]],genotypes[[n]]) #fill array with genotype distances

distances[m,n] <- thisdistance

distances[n,m] <- thisdistance

print(c(labels[m],labels[n]))}} #so the user can monitor the progress of the computation

rownames(distances) <- labels #add labels to the matrix of distances that was just produced

colnames(distances) <- labels

write.table(distances,file=outfile,sep="\t")

distances

}

Processing time

The amount of time that it takes genotype.distance to work is a function of the factorial of the number of alleles that the two genotypes have. If two individuals with four alleles each are being compared, there are only 24 sets of compatible allele comparisons. If the two individuals each have eight alleles, there are 8! or 40,320 sets of comparisons, and if the individuals each have twelve alleles there 479,001,600 sets of comparisons. I wrote the function to optimize for clonal reproduction and large genotypes, simply returning 0 without doing any calculation if the two genotypes are identical, and skipping all of the different ways that alleles can be matched to virtual alleles if one genotype has two or more virtual alleles (since the distance between any allele and a virtual allele is 1).

On my personal laptop (Windows Vista, 2 GHz processor, 4 Mb RAM), processing time became an issue when the allele distance grid was 10x10 (two genotypes of 10 alleles each) or larger. 9x9 grids took a few minutes to process, but anything larger could take over an hour. For my own use (not included in the above code), I added a clause in genotype.distance to return NA if both genotypes had more than nine alleles. I then calculated these allele distance grids in Excel and manually chose sets of compatible allele comparisons, since in these cases it was faster to visually pick out the short distances and work from there than it was to have the computer calculate every possible combination.

If your organism is tetraploid, processing time should not be an issue.

References

Bruvo R, Michiels NK, D'Souza TG, Schulenberg H (2004). "A simple method for the calculation of microsatellite genotype distances irrespective of ploidy level." Molecular Ecology 13:2101-2106.

Contact

User:Lindsay_V._Clark