==GeneABEL==
plink --file GWAS_clean4 --pheno pheno.phen --pheno-name Aff --transpose --recode --out gwa_gabel
plink --file GWAS_clean4 --pheno pheno.phen --pheno-name systolic --transpose --recode --out gwa_gabel_qtl
plot(test.qt, col = "black")
add.plot(test.eg, col = "gray", pch = 3)
legend("topright", c("Original plot","After correction w/ EIGENSTRAT"), pch = c(1,3)) ==Imputation exercise== plink --file chr22_imputation_ex plink --file chr22_imputation_ex --maf 0.01 --mind 0.02 --geno 0.05 --hwe 0.001 --out qc_check plink --file chr22_imputation_ex --maf 0.01 --mind 0.02 --geno 0.05 --hwe 0.001 --recode --out chr22_clean1 plink --file chr22_clean1 --maf 0.01 --mind 0.02 --geno 0.05 --hwe 0.001 --out qc_check_2 plink --file chr22_clean1 --filter-cases --hwe 0.001 --recode --out chr22_cases_clean plink --file chr22_clean1 --filter-controls --recode --out chr22_controls_clean plink --file chr22_controls_clean --merge chr22_cases_clean.ped chr22_cases_clean.map -–hwe 0.001 --recode --out chr22_all_clean plink --file chr22_all_clean --logistic --out chr22_all_clean_geno R mydata = read.table(“chr22_all_clean_geno.assoc.logistic”, header=T) names(mydata) plot(mydata$BP, -log10(mydata$P)) smallp = mydata[which(mydata$P < 1E-4),] smallp smallp = smallp[order(smallp$BP),] smallp q() mach1 --hapmapFormat -d chr22_mach_merlin.map -p chr22_mach_merlin.ped --haps genotypes_chr22_CEU_r22_nr.b36_fwd.phase.gz –-snps genotypes_chr22_CEU_r22_nr.b36_fwd_legend.txt.gz --greedy --rounds 100 --mle --mldetails --autoflip -o chr22_HIHII plink –-dosage chr22_HIHII_dose_mach4plink.txt.gz Zin –-fam chr22_imputation_ex.fam –-map chr22_imputed_snps_positions.map --out chr22_HIHII_dosage R dosage = read.table("chr22_HIHII_dosage.assoc.dosage", header= T) names(dosage) plot(dosage$BP, -log10(dosage$P)) dosagep = dosage[which(dosage$P < 5E-8),] dosagep = dosagep[order(dosagep$BP),] dosagep interest = dosage[which(dosage$SNP=='rs715586'),] interest
)==Imputation exercise==
plink --file chr22_imputation_ex<br data-attributes="%20/" />plink --file chr22_imputation_ex --maf 0.01 --mind 0.02 --geno 0.05 --hwe 0.001 --out qc_check<br data-attributes="%20/" />plink --file chr22_imputation_ex --maf 0.01 --mind 0.02 --geno 0.05 --hwe 0.001 --recode --out chr22_clean1<br data-attributes="%20/" />plink --file chr22_clean1 --maf 0.01 --mind 0.02 --geno 0.05 --hwe 0.001 --out qc_check_2<br data-attributes="%20/" />plink --file chr22_clean1 --filter-cases --hwe 0.001 --recode --out chr22_cases_clean<br data-attributes="%20/" />plink --file chr22_clean1 --filter-controls --recode --out chr22_controls_clean<br data-attributes="%20/" />plink --file chr22_controls_clean --merge chr22_cases_clean.ped chr22_cases_clean.map -–hwe 0.001 --recode --out chr22_all_clean<br data-attributes="%20/" />plink --file chr22_all_clean --logistic --out chr22_all_clean_geno<br data-attributes="%20/" />R<br data-attributes="%20/" />mydata = read.table(“chr22_all_clean_geno.assoc.logistic”, header=T)<br data-attributes="%20/" />names(mydata)<br data-attributes="%20/" />plot(mydata$BP, -log10(mydata$P))<br data-attributes="%20/" />smallp = mydata[which(mydata$P < 1E-4),]<br data-attributes="%20/" />smallp<br data-attributes="%20/" />smallp = smallp[order(smallp$BP),]<br data-attributes="%20/" />smallp<br data-attributes="%20/" />q()<br data-attributes="%20/" />mach1 --hapmapFormat -d chr22_mach_merlin.map -p chr22_mach_merlin.ped --haps genotypes_chr22_CEU_r22_nr.b36_fwd.phase.gz –-snps genotypes_chr22_CEU_r22_nr.b36_fwd_legend.txt.gz --greedy --rounds 100 --mle --mldetails --autoflip -o chr22_HIHII<br data-attributes="%20/" />plink –-dosage chr22_HIHII_dose_mach4plink.txt.gz Zin –-fam chr22_imputation_ex.fam –-map chr22_imputed_snps_positions.map --out chr22_HIHII_dosage<br data-attributes="%20/" />R<br data-attributes="%20/" />dosage = read.table("chr22_HIHII_dosage.assoc.dosage", header= T)<br data-attributes="%20/" />names(dosage)<br data-attributes="%20/" />plot(dosage$BP, -log10(dosage$P))<br data-attributes="%20/" />dosagep = dosage[which(dosage$P < 5E-8),]<br data-attributes="%20/" />dosagep = dosagep[order(dosagep$BP),]<br data-attributes="%20/" />dosagep<br data-attributes="%20/" />interest = dosage[which(dosage$SNP=='rs715586'),]<br data-attributes="%20/" />interest
==PLINK_R==
Multifactorial
plink --ped dbp.cc.ped --map dbp.map --map3 --out logreg.add --logistic --noweb<br data-attributes="%20/" /> plink --ped dbp.cc.ped --map dbp.map --map3 --out logreg.add.ci --logistic --ci 0.95 --noweb<br data-attributes="%20/" /> plink --ped dbp.cc.ped --map dbp.map --map3 --out logreg.age.add --logistic --covar dbp.age.pheno --noweb<br data-attributes="%20/" /> plink --ped dbp.cc.ped --map dbp.map --map3 --out logreg.sex.add --logistic --sex --noweb<br data-attributes="%20/" /> plink --ped dbp.cc.ped --map dbp.map --map3 --out logreg.sexage.add --logistic --sex --covar dbp.age.pheno --noweb<br data-attributes="%20/" /> plink --ped dbp.cc.ped --map dbp.map --map3 --out logreg.snp1112.add --logistic --condition rs1112 --noweb<br data-attributes="%20/" /> plink --ped dbp.cc.ped --map dbp.map --map3 --out logreg.snp1117.add --logistic --condition rs1117<br data-attributes="%20/" /> plink --ped dbp.qt.ped --map dbp.map --map3 --out linreg.sex.add --linear --sex --noweb<br data-attributes="%20/" /> plink --ped dbp.cc.ped --map dbp.map --map3 --out logreg.sex.inter.add --logistic --sex --interaction --noweb<br data-attributes="%20/" /> plink --ped dbp.cc.ped --map dbp.map --map3 --out logreg.snp1112.inter.add --logistic --condition rs1112 --interaction --noweb<br data-attributes="%20/" /> R<br data-attributes="%20/" /> load("dbp.R")<br data-attributes="%20/" /> ls()<br data-attributes="%20/" /> dbp[1:5,]<br data-attributes="%20/" /> result.snp12 = glm (affection ~ rs1112, family=binomial("logit"), data=dpb)<br data-attributes="%20/" /> print (result.snp12)<br data-attributes="%20/" /> print ( class (result.snp12) )<br data-attributes="%20/" /> print ( summary(result.snp12) )<br data-attributes="%20/" /> dev.geno = anova (result.snp12, test="Chi")<br data-attributes="%20/" /> lrt.pvalue = pchisq(dev.geno[dim(dev.geno)[1],"Deviance"], df=2, ncp=0, FALSE)<br data-attributes="%20/" /> print ( lrt.pvalue )<br data-attributes="%20/" /> print ( summary(result.snp12)$coefficients )<br data-attributes="%20/" /> snp.beta = summary(result.snp12)$coefficients[2:3,1]<br data-attributes="%20/" /> print ( snp.beta )<br data-attributes="%20/" /> print ( exp(snp.beta) )<br data-attributes="%20/" /> ci = confint (result.snp12)<br data-attributes="%20/" /> print (ci)<br data-attributes="%20/" /> print ( exp(ci) )<br data-attributes="%20/" /> snp.data = dpb[,c("affection", "rs1112")]<br data-attributes="%20/" /> summary(snp.data)<br data-attributes="%20/" /> snp.data[,"rs1112"] summary(snp.data)<br data-attributes="%20/" /> result.all = glm (affection ~ rs1112, family=binomial("logit"), data=snp.data)<br data-attributes="%20/" /> dev.all = anova (result.all, test="Chi")<br data-attributes="%20/" /> summary(result.all)<br data-attributes="%20/" /> print(dev.all)<br data-attributes="%20/" /> snp.data = dpb[,c("affection", "trait","sex", "age", "rs1112", "rs1117")]<br data-attributes="%20/" /> summary(snp.data)<br data-attributes="%20/" /> snp.data[,"rs1112"] snp.data[,"rs1117"] result.adj = glm (affection ~ sex + rs1112 , family=binomial("logit"), data=snp.data)<br data-attributes="%20/" /> summary(result.adj)<br data-attributes="%20/" /> result.adj = glm (affection ~ age + rs1112 , family=binomial("logit"), data=snp.data)<br data-attributes="%20/" /> summary(result.adj)<br data-attributes="%20/" /> result.adj = glm (affection ~ sex + age + rs1112, family=binomial("logit"), data=snp.data)<br data-attributes="%20/" /> summary(result.adj)<br data-attributes="%20/" /> result.adj = glm (affection ~ rs1117 + rs1112, family=binomial("logit"), data=snp.data)<br data-attributes="%20/" /> summary(result.adj)<br data-attributes="%20/" /> anova (result.adj, test="Chi")<br data-attributes="%20/" /> result.adj = glm (affection ~ rs1112 + rs1117, family=binomial("logit"), data=snp.data)<br data-attributes="%20/" /> summary(result.adj)<br data-attributes="%20/" /> anova (result.adj, test="Chi")<br data-attributes="%20/" /> result.adj = lm (trait ~ rs1112, data=snp.data)<br data-attributes="%20/" /> summary(result.adj)<br data-attributes="%20/" /> result.adj = lm (trait ~ sex + rs1112, data=snp.data)<br data-attributes="%20/" /> summary(result.adj)<br data-attributes="%20/" /> result.inter = glm (affection ~ sex * rs1112, family=binomial("logit"), data=snp.data)<br data-attributes="%20/" /> summary(result.inter)<br data-attributes="%20/" /> result.inter = glm (affection ~ age * rs1112, family=binomial("logit"), data=snp.data)<br data-attributes="%20/" /> summary(result.inter)<br data-attributes="%20/" /> result.inter = glm (affection ~ rs1112 * rs1117, family=binomial("logit"), data=snp.data)<br data-attributes="%20/" /> summary(result.inter)<br data-attributes="%20/" /> result.reg = glm (affection ~ sex + age + rs1112 + rs1117, family=binomial("logit"), data=snp.data)<br data-attributes="%20/" /> summary(result.reg)<br data-attributes="%20/" /> modelchoice.result summary(modelchoice.result)<br data-attributes="%20/" /> q()
GWAS Control Substructure
Multiple Testing
plink --ped dbp.cc.ped --map dbp.map --map3 --out multtest --assoc --adjust --noweb<br data-attributes="%20/" /> plink --ped dbp.cc.ped --map dbp.map --map3 --out multperm5000 --assoc --mperm 5000 --noweb<br data-attributes="%20/" /> plink --ped dbp.cc.ped --map dbp.map --map3 --out multperm100000 --assoc --mperm 100000 --noweb<br data-attributes="%20/" /> R<br data-attributes="%20/" /> load("p.values.R")<br data-attributes="%20/" /> ls()<br data-attributes="%20/" /> p.values<br data-attributes="%20/" /> library (multtest)<br data-attributes="%20/" /> adj.p.values = mt.rawp2adjp(p.values,c("Bonferroni","Holm","SidakSS","BH"))<br data-attributes="%20/" /> adj.p.values<br data-attributes="%20/" /> rownames(adj.p.values$adjp) = names(p.values[adj.p.values$index])<br data-attributes="%20/" /> adj.p.values$adjp
==SEQPower==
spower -h
spower LOGIT -h
spower LOGIT Kryukov2009European1800.sfs --sample_size 1000 --OR_rare_detrimental $i --method "CFisher --name CMC$i" --title FixedOR$i -r 100 -j 4 -l 1 -o exercise2
done==Unphased==
unphased.sh
unphased mypeds.ped –marker 1 2 3 –missing –permutation 10
unphased all.ped -window 2 -LD >> results.txt
==VAT==
ulimit -s 8000
vtools -h