View source for Mixed models Association Analysis Exercise

__NOTITLE__

==Mixed models Association Analysis Exercise==

 head FLMMresults
 more topsnps.txt
 plink --bfile quantfamdata --extract topsnps.txt --recodeA --out topsnpsFLMM
 R
 fam<-read.table("quantfamdata.fam", header=F)
 pheno=data.frame(fam[,1:2],fam[,6])
 write.table(pheno,file="phenos.txt",col.names=F,row.names=F,quote=F)
 gcta64 --mlma --bfile quantfamdata --pheno phenos.txt --out GCTAresults
 R
 source("qqmanHJCupdated.R")
 res3<-read.table("GCTAresults.mlma", header=T)
 head(res3)
 chi<-(qchisq(1-res3$p,1))
 lambda=median(chi)/0.456
 lambda
 new3<-data.frame(res3$SNP, res3$Chr, res3$bp, res3$p)
 names(new3)<-c("SNP", "CHR", "BP", "P")
 head(new3)
 qq(new3$P)
 manhattan(new3, pch=20, suggestiveline=F, genomewideline=F, ymin=2, cex.x.axis=0.65, colors=c("black","dodgerblue"), cex=0.5)
 res2<-read.table("FLMMresults", header=T)
 new2<-data.frame(res2$SNP, res2$Chromosome, res2$Position, res2$Pvalue)
 names(new2)<-c("SNP", "CHR", "BP", "P")
 merged=merge(new3,new2, by="SNP", sort=F)
 head(res2)
 head(new2)
 head(new3)
 head(merged)
 plot(-log10(merged$P.x),-log10(merged$P.y))
 abline(0,1)
 q()
 gcta64 --bfile quantfamdata --autosome --make-grm-bin --out GCTAgrm
 gcta64 --reml --grm-bin GCTAgrm --pheno phenos.txt --out GCTAherit
 gcta64 --bfile quantfamdata --chr 1 --make-grm-bin --out GCTAgrmchr1
 gcta64 --reml --grm-bin GCTAgrmchr1 --pheno phenos.txt --out GCTAheritchr1
 gcta64 --bfile quantfamdata --chr 2 --make-grm-bin --out GCTAgrmchr2
 gcta64 --reml --grm-bin GCTAgrmchr2 --pheno phenos.txt --out GCTAheritchr2
 gcta64 --bfile quantfamdata --chr 6 --make-grm-bin --out GCTAgrmchr6
 gcta64 --reml --grm-bin GCTAgrmchr6 --pheno phenos.txt --out GCTAheritchr6
 gcta64 --bfile quantfamdata --chr 12 --make-grm-bin --out GCTAgrmchr12
 gcta64 --reml --grm-bin GCTAgrmchr12 --pheno phenos.txt --out GCTAheritchr12
 gcta64 --bfile quantfamdata --chr 1 --make-grm-bin --out GCTAgrmchr1
 gcta64 --bfile quantfamdata --chr 2 --make-grm-bin --out GCTAgrmchr2
 gcta64 --bfile quantfamdata --chr 3 --make-grm-bin --out GCTAgrmchr3
 gcta64 --bfile quantfamdata --chr 4 --make-grm-bin --out GCTAgrmchr4
 gcta64 --bfile quantfamdata --chr 5 --make-grm-bin --out GCTAgrmchr5
 gcta64 --bfile quantfamdata --chr 6 --make-grm-bin --out GCTAgrmchr6
 gcta64 --bfile quantfamdata --chr 7 --make-grm-bin --out GCTAgrmchr7
 gcta64 --bfile quantfamdata --chr 8 --make-grm-bin --out GCTAgrmchr8
 gcta64 --bfile quantfamdata --chr 9 --make-grm-bin --out GCTAgrmchr9
 gcta64 --bfile quantfamdata --chr 10 --make-grm-bin --out GCTAgrmchr10
 gcta64 --bfile quantfamdata --chr 11 --make-grm-bin --out GCTAgrmchr11
 gcta64 --bfile quantfamdata --chr 12 --make-grm-bin --out GCTAgrmchr12
 gcta64 --bfile quantfamdata --chr 13 --make-grm-bin --out GCTAgrmchr13
 gcta64 --bfile quantfamdata --chr 14 --make-grm-bin --out GCTAgrmchr14
 gcta64 --bfile quantfamdata --chr 15 --make-grm-bin --out GCTAgrmchr15
 gcta64 --bfile quantfamdata --chr 16 --make-grm-bin --out GCTAgrmchr16
 gcta64 --bfile quantfamdata --chr 17 --make-grm-bin --out GCTAgrmchr17
 gcta64 --bfile quantfamdata --chr 18 --make-grm-bin --out GCTAgrmchr18
 gcta64 --bfile quantfamdata --chr 19 --make-grm-bin --out GCTAgrmchr19
 gcta64 --bfile quantfamdata --chr 20 --make-grm-bin --out GCTAgrmchr20
 gcta64 --bfile quantfamdata --chr 21 --make-grm-bin --out GCTAgrmchr21
 gcta64 --bfile quantfamdata --chr 22 --make-grm-bin --out GCTAgrmchr22
 gcta64 --reml --mgrm-bin multipleGRMs.txt --pheno phenos.txt --out GCTAherit22GRMs
 gcta64 --mgrm-bin multipleGRMsnot6.txt --make-grm --out GCTAgrmnot6
 gcta64 --reml --mgrm-bin multipleGRMs6andnot6.txt --pheno phenos.txt --out GCTAherit6andnot6
 gcta64 --reml --reml-pred-rand --grm-bin GCTAgrm --pheno phenos.txt --out GCTAblups
 R
 blups<-read.table("GCTAblups.indi.blp", header=F)
 names(blups)=c("FID","IID","intV1","BLUP", "intV2", "Residual")
 head(blups)
 hist(blups$BLUP)
 topgenos<-read.table("topsnpsFLMM.raw",header=T)
 blupsgenos=data.frame(blups, topgenos)
 blupsgenos$combo=10*blupsgenos$rs9271209_G+blupsgenos$rs233722_G
 head(blupsgenos)
 par(mfrow=c(3,3))
 hist(blupsgenos$BLUP[blupsgenos$combo==0], main="genotype=0", xlab="BLUP")
 hist(blupsgenos$BLUP[blupsgenos$combo==1], main="genotype=1", xlab="BLUP") 
 hist(blupsgenos$BLUP[blupsgenos$combo==2], main="genotype=2", xlab="BLUP")
 hist(blupsgenos$BLUP[blupsgenos$combo==10], main="genotype=10", xlab="BLUP")
 hist(blupsgenos$BLUP[blupsgenos$combo==11], main="genotype=11", xlab="BLUP")
 hist(blupsgenos$BLUP[blupsgenos$combo==12], main="genotype=12", xlab="BLUP")
 hist(blupsgenos$BLUP[blupsgenos$combo==20], main="genotype=20", xlab="BLUP")
 hist(blupsgenos$BLUP[blupsgenos$combo==21], main="genotype=21", xlab="BLUP")
 hist(blupsgenos$BLUP[blupsgenos$combo==22], main="genotype=22", xlab="BLUP")
 q()
 gcta64 --bfile quantfamdata --blup-snp GCTAblups.indi.blp --out GCTAsnpblups
 R
 snpblups<-read.table("GCTAsnpblups.snp.blp", header=F)
 names(snpblups)=c("SNP","Allele","snpBLUP","snpResidual")
 head(snpblups)
 res2<-read.table("FLMMresults", header=T)
 new2<-data.frame(res2$SNP, res2$Chromosome, res2$Position, res2$Pvalue)
 names(new2)<-c("SNP", "CHR", "BP", "P")
 head(new2)
 mergeblup=merge(snpblups,new2, by="SNP",sort=F)
 head(mergeblup)
 par(mfrow=c(1,1))
 plot(mergeblup$snpBLUP,-log10(mergeblup$P))