Rscripts - Revision history

Albrecht: /* Fst */

2015-03-12T13:25:29Z

Fst

2015-03-12T13:23:59Z

Fst

2015-03-12T13:23:24Z

Fst

← Older revision		Revision as of 15:23, 12 March 2015
Line 59:		Line 59:
	<pre>		<pre>
	WC84<-function(x,pop){		WC84<-function(x,pop){

			popNam<-sort(unique(pop))

	#number ind each population		#number ind each population
	n<-table(pop)		n<-table(pop)
Line 68:		Line 71:
	N<-length(pop)		N<-length(pop)
	#frequency in samples		#frequency in samples
	p<-apply(x,2,function(x,pop){tapply(x,pop,mean)/2},pop=pop)		p<-t(sapply(popNam,function(y) colMeans(x[pop==y,])/2))
			# p<-apply(x,2,function(x,pop){tapply(x,pop,mean)/2},pop=pop)
	#average frequency in all samples (apply(x,2,mean)/2)		#average frequency in all samples (apply(x,2,mean)/2)
	p_avg<-as.vector(n%*%p/N )		p_avg<-as.vector(n%*%p/N )

2015-02-13T08:35:17Z

PCA/Eigensoft/Eigenstrat

← Older revision		Revision as of 10:35, 13 February 2015
Line 11:		Line 11:
	nMis <- rowSums(is.na(geno))		nMis <- rowSums(is.na(geno))
	freq <- rowMeans(geno,na.rm=T)/2 # get allele frequency		freq <- rowMeans(geno,na.rm=T)/2 # get allele frequency
	keep <- freq>minMaf&freq<21-minMaf) & nMis<=maxMis # remove sites with non-polymorphic data		keep <- freq>minMaf&freq<(1-minMaf) & nMis<=maxMis # remove sites with non-polymorphic data
	freq<-freq[keep]		freq<-freq[keep]
	geno<-geno[keep,]		geno<-geno[keep,]

2015-02-13T08:34:54Z

PCA/Eigensoft/Eigenstrat

← Older revision		Revision as of 10:34, 13 February 2015
Line 9:		Line 9:
	##maxMis maximum allowed missing genotypes for a site		##maxMis maximum allowed missing genotypes for a site

	nMis<-rowSums(is.na(geno))		nMis <- rowSums(is.na(geno))
	~~avg~~<-rowMeans(geno,na.rm=T) # get allele frequency ~~times 2~~		freq <- rowMeans(geno,na.rm=T)/2 # get allele frequency
	keep<-~~avg~~>minMaf&~~avg~~<~~2*(1~~-minMaf)& nMis<=maxMis # remove sites with non-polymorphic data		keep <- freq>minMaf&freq<21-minMaf) & nMis<=maxMis # remove sites with non-polymorphic data
	~~avg~~<-~~avg~~[keep]		freq<-freq[keep]
	geno<-geno[keep,]		geno<-geno[keep,]
	snp<-nrow(geno) #number of snps used in analysis		snp<-nrow(geno) #number of snps used in analysis
	ind<-ncol(geno) #number of individuals used in analuysis		ind<-ncol(geno) #number of individuals used in analuysis
	~~freq<-avg/2 #frequency~~		M <- (geno-freq2)/sqrt(freq(1-freq)) #normalize the genotype matrix
	M <- (geno-~~avg~~)/sqrt(freq*(1-freq)) #normalize the genotype matrix
	M[is.na(M)]<-0		M[is.na(M)]<-0
	X<-t(M)%*%M #get the (almost) covariance matrix		X<-t(M)%*%M #get the (almost) covariance matrix

2015-02-13T08:30:15Z

PCA/Eigensoft/Eigenstrat

← Older revision		Revision as of 10:30, 13 February 2015
Line 45:		Line 45:
	Example		Example
	<pre>		<pre>
	ind<-c(100,100)		ind<-c(100,100,100)

	snp<-10000		snp<-10000
	freq=c(runif(snp),runif(snp))		freq<-cbind(runif(snp),runif(snp),runif(snp))
	~~geno~~<-~~c()~~		l<-lapply(1:length(ind),function(x) matrix(rbinom(snp*ind[x],2,freq[,x]),snp))
	~~for~~(~~pop in~~ 1:length(ind))		geno<-do.call(cbind,l)
	~~geno<-rbind(geno,~~matrix(rbinom(snp*ind[~~pop~~],2,freq[~~pop~~]),~~ind[pop]~~))
	geno<-t(~~geno~~)
	e<-eigenstrat(geno)		e<-eigenstrat(geno)

	plot(e,col=rep(1:length(ind),ind),xlab="PC1",ylab="PC2")		plot(e,col=rep(1:length(ind),ind),xlab="PC1",ylab="PC2")
	</pre>		</pre>

2015-02-13T08:22:13Z

2015-02-13T08:17:31Z

PCA/Eigensoft/Eigenstrat

← Older revision		Revision as of 10:17, 13 February 2015
Line 5:		Line 5:

	<pre>		<pre>
	eigenstrat<-function(geno){ #snp x ind matrix of genotypes \in 0,1,2		eigenstrat<-function(geno,maxMis=0,minMaf=0.01){
			## geno: snp x ind matrix of genotypes \in 0,1,2
			##maxMis maximum allowed missing genotypes for a site

	nMis<-rowSums(is.na(geno))		nMis<-rowSums(is.na(geno))
	~~geno<-geno[nMis==0,] #remove snps with missing data~~		avg<-rowMeans(geno,na.rm=T) # get allele frequency times 2
	avg<-~~rowSums(geno)/ncol~~(geno) # get allele frequency times 2		keep<-avg>minMaf&avg<2*(1-minMaf)& nMis<=maxMis # remove sites with non-polymorphic data
	keep<-avg~~!=0~~&avg!=2 # remove sites with non-polymorphic data
	avg<-avg[keep]		avg<-avg[keep]
	geno<-geno[keep,]		geno<-geno[keep,]
Line 16:		Line 18:
	freq<-avg/2 #frequency		freq<-avg/2 #frequency
	M <- (geno-avg)/sqrt(freq*(1-freq)) #normalize the genotype matrix		M <- (geno-avg)/sqrt(freq*(1-freq)) #normalize the genotype matrix
			M[is.na(M)]<-0
	X<-t(M)%*%M #get the (almost) covariance matrix		X<-t(M)%*%M #get the (almost) covariance matrix
	X<-X/(sum(diag(X))/(snp-1))		X<-X/(sum(diag(X))/(snp-1))
Line 25:		Line 28:
	E$mu<-mu		E$mu<-mu
	E$sigma<-sigma		E$sigma<-sigma
			E$nSNP <- nrow(geno)
			E$nInd <- ncol(geno)
	class(E)<-"eigenstrat"		class(E)<-"eigenstrat"
	E		E
	}		}



	plot.eigenstrat<-function(x,col=1,...)		plot.eigenstrat<-function(x,col=1,...)
	plot(x$vectors[,1:2],col=col,...)		plot(x$vectors[,1:2],col=col,...)

2013-08-15T10:33:29Z

@@ Line 1: / Line 1: @@
 = PCA/Eigensoft/Eigenstrat =
 <pre>
 </pre>
+Example
 <pre>
 </pre>
 = Fst =
 <pre>
 </pre>
@@ Line 20: / Line 99: @@
 <pre>
 </pre>
 = LD pruning in R =
+Install R package
 <pre>
+wget http://www.popgen.dk/albrecht/misc_Rpackages/Rpakker/pruning_0.51.tar.gz
 </pre>

2013-08-15T09:51:56Z

@@ Line 1: / Line 1: @@
 = PCA/Eigensoft/Eigenstrat =
 = Fst =
 = LD pruning in R =