# Load files
library(GenABEL)
convert.snp.tped(tped &eq#61; "gwa_gabel_qtl.tped", tfam &eq#61; "gwa_gabel_qtl.tfam", out &eq#61; "gwa_gabel_qtl.raw", strand &eq#61; "u") g.dat <- load.gwaa.data(phen &eq#61; "gwa_gabel_qtl.praw", gen &eq#61; "gwa_gabel_qtl.raw", force &eq#61; T)
slotNames(g.dat)
slotNames(g.dat@gtdata)
# Trait
summary(g.dat@phdata$disease)
hist(g.dat@phdata$disease, main&eq#61;"Quantitative Phenotype data summary", xlab &eq#61; "Systolic pressure measure", freq &eq#61; F,breaks&eq#61;20, col&eq#61;"gray")
rug(g.dat@phdata$disease)
###
###
# GLM test
test.snp <- scan.glm('disease ~ CRSNP', family &eq#61; gaussian(), data &eq#61; g.dat)
names(test.snp)
alpha <- 5e-8
test.snp$P1df[test.snp$P1df < alpha]
# Score test
test.qt <- qtscore(disease, data &eq#61; g.dat, trait &eq#61; "gaussian")
slotNames(test.qt)
names(test.qt@results)
obs <- sort(results(test.qt)$P1df)
ept <- c(1:length(obs)) / (length(obs) + 1)
plot(-log10(ept), -log10(obs), main &eq#61; "GWAS QQ plot, qtl", xlab&eq#61;"Expected -log10(pvalue)", ylab&eq#61;"Observed -log10(pvalue)") abline(0, 1, col &eq#61; "red") abline(h &eq#61; 8, lty &eq#61; 2)
# Manhattan plot
plot(test.qt, col &eq#61; "black")
# Adding confounders
test.qt.sex <- qtscore(disease ~ sex, data &eq#61; g.dat, trait &eq#61; "gaussian")
rownames(results(test.qt.sex))[results(test.qt)$P1df < alpha]
summary(lm(disease ~ sex, data &eq#61; g.dat))
###
# MDS
###
gkin <- ibs(g.dat, weight &eq#61; "freq")
gkin[1:10,1:10]
cps.full <- cmdscale(as.dist(.5 - gkin), eig &eq#61; T, k &eq#61; 10)
names(cps.full)
cps <- cps.full$points
plot(cps[,1], cps[,2], pch &eq#61; g.dat@phdata$popn) legend(-0.16, 0.06, c("TSI","MEX", "CEU"), pch &eq#61; c(1,2,3))
###
# Corrected test
gpc.dat <- g.dat
gpc.dat@phdata<-cbind(g.dat@phdata, cps)
test.pc.a <- scan.glm('disease ~ CRSNP + C1 + C2 + C3 + C4 + C5', family&eq#61;gaussian(), data &eq#61; gpc.dat)
test.pc.a$snpnames[test.pc.a$P1df < alpha]
test.pc.a$P1df[test.pc.a$P1df < alpha]
test.pc.b <- qtscore(disease ~ C1 + C2 + C3 + C4 + C5, data &eq#61; gpc.dat, trait &eq#61; "gaussian")
test.pc.b@lambda
# scree plot
plot(cps.full$eig[1:10]/sum(cps.full$eig), axes &eq#61; F, type &eq#61; "b", xlab &eq#61; "Components", ylim &eq#61; c(0,0.05), ylab &eq#61; "Proportion of Variations", main &eq#61; "MDS analysis scree plot")
axis(1, 1:10)
axis(2)
# cumulative plot
plot(cumsum(cps.full$eig[1:10])/sum(cps.full$eig), axes &eq#61; F, type &eq#61; "b", ylim &eq#61; c(0,0.2), xlab &eq#61; "Components", ylab &eq#61; "Proportion of Variations", main &eq#61; "MDS analysis cumulative plot")
axis(1, 1:10)
axis(2)
# Check for inflation of statistic
obs <- sort(results(test.qt)$chi2.1df)
ept <- sort(qchisq(1:length(obs) / (length(obs) + 1), df &eq#61; 1)) plot(ept, obs, main &eq#61; "Genomic control (slope is the inflation factor)", xlab&eq#61;"Expected chisq, 1df", ylab&eq#61;"Observed chisq, 1df") abline(0, 1, col &eq#61; "red") abline(0, test.qt@lambda[1], lty &eq#61; 2)
# Definition of GIF
# Conventional definition
obs <- sort(results(test.qt)$Pc1df)
ept <- c(1:length(obs)) / (length(obs) + 1)
plot(-log10(ept), -log10(obs), main &eq#61; "GWAS QQ plot adj. via Genomic Control", xlab&eq#61;"Expected -log10(pvalue)", ylab&eq#61;"Observed -log10(pvalue)") abline(0, 1, col &eq#61; "red") abline(h &eq#61; 8, lty &eq#61; 2)
# EIGENSTRAT
adj.gkin &eq#61; gkin diag(adj.gkin) &eq#61; hom(g.dat)$Var # naxes &eq#61; 3 is default value test.eg <- egscore(disease, data &eq#61; g.dat, kin &eq#61; adj.gkin, naxes &eq#61; 2)
descriptives.scan(test.eg)
snp.eg <- row.names(results(test.eg))[results(test.eg)$P1df < alpha]
# Change #PCs
for (k in 1:10){
test.tmp <- egscore(disease, data &eq#61; g.dat, kin &eq#61; adj.gkin, naxes &eq#61; k)
print(test.tmp@lambda$estimate)
}
obs <- sort(results(test.eg)$Pc1df)
ept <- c(1:length(obs)) / (length(obs) + 1)
qqplot(-log10(ept), -log10(obs), main &eq#61; "GWAS QQ plot adj. w/ EIGENSTRAT", xlab&eq#61;"Expected -log10(pvalue)", ylab&eq#61;"Observed -log10(pvalue)") abline(0, 1, col &eq#61; "red") abline(h &eq#61; 8, lty &eq#61; 2)
# Manhattan plot comparison
plot(test.qt, col &eq#61; "black") add.plot(test.eg, col &eq#61; "gray", pch &eq#61; 3) legend("topright", c("Original plot","After correction w/ EIGENSTRAT"), pch &eq#61; c(1,3))
###
# Basic test, binary trait
###
# load files to GenABEL
convert.snp.tped(tped &eq#61; "gwa_gabel.tped", tfam &eq#61; "gwa_gabel.tfam", out &eq#61; "gwa_gabel.raw", strand &eq#61; "u") b.dat <- load.gwaa.data(phen &eq#61; "gwa_gabel.praw", gen &eq#61; "gwa_gabel.raw", force &eq#61; T)
slotNames(b.dat)
slotNames(b.dat@gtdata)
b.dat@gtdata@nids
# number of cases and controls
case.size <- length(which(b.dat@phdata$disease &eq#61;&eq#61; 1)) control.size <- length(which(b.dat@phdata$disease &eq#61;&eq#61; 0))
case.size
control.size
snpsb.total <- b.dat@gtdata@nsnps
# GLM test
testb.snp <- scan.glm('disease ~ CRSNP', family &eq#61; binomial(), data &eq#61; b.dat)
names(testb.snp)
alpha <- 5e-8
testb.snp$P1df[testb.snp$P1df < alpha]
# Score test
testb.qt <- qtscore(disease, data &eq#61; b.dat, trait &eq#61; "binomial")
slotNames(testb.qt)
descriptives.scan(testb.qt)
#### in R - open R by simply typing R
setwd("to_your_working_directory/")
sexcheck &eq#61; read.table("GWAS_sex_checking.sexcheck", header&eq#61;T)
names(sexcheck)
sex_problem &eq#61; sexcheck[which(sexcheck$STATUS&eq#61;&eq#61;"PROBLEM"),]
sex_problem
q()
#### in R
setwd("to_your_working_directory/")
dups &eq#61; read.table("duplicates.genome", header &eq#61; T) problem_pairs &eq#61; dups[which(dups$PI_HAT > 0.4),]
problem_pairs
problem_pairs &eq#61; dups[which(dups$PI_HAT > 0.05),] myvars &eq#61; c("FID1", "IID1", "FID2", "IID2", "PI_HAT")
problem_pairs[myvars]
q()
plink --file GWAS_clean3 --het
###### in R
Dataset <- read.table("plink.het", header&eq#61;TRUE, sep&eq#61;"", na.strings&eq#61;"NA", dec&eq#61;".", strip.white&eq#61;TRUE)
mean(Dataset$F)
sd(Dataset$F)
jpeg("hist.jpeg", height&eq#61;1000, width&eq#61;1000) hist(scale(Dataset$F), xlim&eq#61;c(-4,4))
dev.off()
q()
plink --file GWAS_clean3 --pheno pheno.txt --pheno-name Aff --hardy
##### in R
hardy &eq#61; read.table("plink.hwe", header &eq#61; T)
names(hardy)
hwe_prob &eq#61; hardy[which(hardy$P < 0.0000009),]
hwe_prob
q()
plink --file GWAS_clean4 --genome --cluster --mds-plot 10
#### in R
mydata &eq#61; read.table("mds_components.txt", header&eq#61;T) mydata$pch[mydata$Group&eq#61;&eq#61;1 ] <-15 mydata$pch[mydata$Group&eq#61;&eq#61;2 ] <-16 mydata$pch[mydata$Group&eq#61;&eq#61;3 ] <-2 jpeg("mds.jpeg", height&eq#61;500, width&eq#61;500) plot(mydata$C1, mydata$C2 ,pch&eq#61;mydata$pch)
dev.off()
q()
expected <- c(1:length(observed))
lexp <- -(log10(expected / (length(expected)+1)))
plot(c(0,7), c(0,7), col&eq#61;"red", lwd&eq#61;3, type&eq#61;"l", xlab&eq#61;"Expected (-logP)", ylab&eq#61;"Observed (-logP)", xlim&eq#61;c(0,max(lobs)), ylim&eq#61;c(0,max(lobs)), las&eq#61;1, xaxs&eq#61;"i", yaxs&eq#61;"i", bty&eq#61;"l", main &eq#61; title) points(lexp, lobs, pch&eq#61;23, cex&eq#61;.4, bg&eq#61;"black") } jpeg("qqplot_compare.jpeg", height&eq#61;1000, width&eq#61;500) par(mfrow&eq#61;c(2,1)) aff_unadj<-read.table("unadj.assoc.logistic", header&eq#61;TRUE) aff_unadj.add.p<-aff_unadj[aff_unadj$TEST&eq#61;&eq#61;c("ADD"),]$P
broadqq(aff_unadj.add.p,"Some Trait Unadjusted")
aff_C1C2<-read.table("PC1-PC2.assoc.logistic", header&eq#61;TRUE) aff_C1C2.add.p<-aff_C1C2[aff_C1C2$TEST&eq#61;&eq#61;c("ADD"),]$P
broadqq(aff_C1C2.add.p, "Some Trait Adjusted for PC1 and PC2")
dev.off()
gws_unadj &eq#61; aff_unadj[which(aff_unadj$P < 0.0000001),]
gws_unadj
gws_adjusted &eq#61; aff_C1C2[which(aff_C1C2$P < 0.0000001),]
gws_adjusted
head GenotypeSummary.txt
vtools output variant "max(DP)" "min(DP)" "avg(DP)" "stdev(DP)" "lower_quartile(DP)" "upper_quartile(DP)" --header
vtools select variant "filter&eq#61;’PASS’" --count vtools select variant "filter&eq#61;’PASS’" -o "max(DP)" "min(DP)" "avg(DP)" "stdev(DP)" "lower_quartile(DP)" "upper_quartile(DP)" --header vtools update variant --from_stat ’total&eq#61;#(GT)’ ’num&eq#61;#(alt)’ ’het&eq#61;#(het)’ ’hom&eq#61;#(hom)’ ’other&eq#61;#(other)’ ’minDP&eq#61;min(DP_geno)’ ’maxDP&eq#61;max(DP_geno)’ ’meanDP&eq#61;avg(DP_geno)’ ’maf&eq#61;maf()’
vtools show fields
vtools show table variant
vtools update variant --from_stat ’totalGD10&eq#61;#(GT)’ ’numGD10&eq#61;#(alt)’ ’hetGD10&eq#61;#(het)’ ’homGD10&eq#61;#(hom)’ ’otherGD10&eq#61;#(other)’ ’mafGD10&eq#61;maf()’ --genotypes "DP_geno > 10"
vtools show fields
vtools show table variant
vtools output variant chr pos maf mafGD10 --header --limit 20
vtools phenotype --set "RACE&eq#61;0" --samples "filename like ’YRI%’" vtools phenotype --set "RACE&eq#61;1" --samples "filename like ’CEU%’"
vtools show samples --limit 10
vtools update variant --from_stat ’CEU_mafGD10&eq#61;maf()’ --genotypes ’DP_geno>10’ --samples "RACE&eq#61;1" vtools update variant --from_stat ’YRI_mafGD10&eq#61;maf()’ --genotypes ’DP_geno>10’ --samples "RACE&eq#61;0"
vtools output variant chr pos mafGD10 CEU_mafGD10 YRI_mafGD10 --header --limit 10
vtools phenotype --from_stat ’CEU_totalGD10&eq#61;#(GT)’ ’CEU_numGD10&eq#61;#(alt)’ --genotypes ’DP_geno>10’ --samples "RACE&eq#61;1" vtools phenotype --from_stat ’YRI_totalGD10&eq#61;#(GT)’ ’YRI_numGD10&eq#61;#(alt)’ --genotypes ’DP_geno>10’ --samples "RACE&eq#61;0"
vtools phenotype --output sample_nameCEU_totalGD10CEU_numGD10YRI_totalGD10YRI_numGD10 --header
vtools execute ANNOVAR geneanno
vtools show tables
vtools remove genotypes "DP_geno<10" -v0
vtools select variant "mut_type like ’non%’ or mut_type like ’stop%’ or region_type&eq#61;’splicing’" -t v_funct
vtools show tables
vtools show samples --limit 5
vtools select variant --samples "RACE&eq#61;1" -t CEU
mkdir -p ceu
cd ceu
vtools init ceu --parent ../ --variants CEU --samples "RACE&eq#61;1" --build hg19 vtools show project vtools select variant "CEU_mafGD10>&eq#61;0.05" -t common_ceu
vtools select v_funct "CEU_mafGD10<0.01" -t rare_ceu
vtools use refGene
less EA_RV_VT.asso.res
sort -g -k6 EA_RV_VT.asso.res | head
vtools select rare_ceu "refGene.name2&eq#61;’ABCC1’" -o chr pos ref alt CEU_mafGD10 numGD10 mut_type --header
cd ..
vtools select variant --samples "RACE&eq#61;0" -t YRI
mkdir -p yri
cd yri
vtools init yri --parent ../ --variants YRI --samples "RACE&eq#61;0" --build hg19 vtools select variant "YRI_mafGD10>&eq#61;0.05" -t common_yri vtools select v_funct "YRI_mafGD10<0.01" -t rare_yri
vtools use refGene
vtools associate common_yri BMI --covariate SEX -m "LinRegBurden --alternative 2" -j1 --to_db YA_CV > YA_CV.asso.res