Difference between revisions of "Mixed models Association Analysis Exercise"

From Statistical Genetics Courses

Jump to: navigation, search
(Created page with "__NOTITLE__ ==Mixed models Association Analysis Exercise==")
 
Line 2: Line 2:
  
 
==Mixed models Association Analysis Exercise==
 
==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")

Revision as of 19:23, 23 January 2019

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")