Population allelic richness by bootstrapping in R

[JackieBadger]

Hi,

I realize this may not be the best place for R help but thought I'd give it a shot.

I am trying to calculate allelic richness of MHC genes. To be more precise, because MHC alleles can not be assigned to loci I need to calculate richness on a population level. In addition because I have different n's for each population I need to standardize my allele counting using resampling.

So far in R:


Bootstrap.pop1<-replicate(1000,sample(POP1,10,replace=T))
table(Boot1)


This sub-samples all of the alleles found in my population with a sampling effort of 10, 1000 times.
Table then tells me how many times each allele was samples.

However, I need to count the number of unique alleles sampled at each of the 1000 iterations, and then find the mean of these values.

Any tips would be greatly appreciated....

Thanks,

J

[mastal]

Hi Jackie,

There seems to be a unique() function in R, although I've not used it before.

It removes duplicate elements from a vector, so then if the vector was made up of the alleles sampled in 1 iteration, length(unique.vector) would give you the number of unique elements sampled in the run.

But you might have to put the bootstrap procedure in a loop.

Maria

[JackieBadger]

Thanks for the suggestion to use a loop.

I only started using R last week, but what I came up works well so i thought I would share it incase others face the same problem.

R code:

###created an empty matrix of 1x 1000###
emptyvar<-numeric(1000)
emptyvar<-matrix(ncol=1000,nrow=1)

###Loop
#######Calculates the allelic richness in a population though bootstrapping. Needed when population n's are not equal. This is useful for MHC as most of the time alleles can not be assigned to loci.
Randomdly samples 10 alleles with replacement, 1000 times. Output is mean number of alleles within a population and standard deviation#######



for (i in 1:1000) {
sample1<-sample(MHC1,10,replace=T)


n.alleles1<-length(unique(sample1))
n.alleles1
emptyvar[i]<-n.alleles1
}
mean(emptyvar)
apply(emptyvar,1,sd)