Changes

From Statistical Genetics Courses

Jump to: navigation, search

AdvGeneMap2017Commands

1,677 bytes removed, 20:48, 4 January 2017
__NOTITLE__
===GeneABEL=== 
plink --file GWAS_clean4 --pheno pheno.phen --pheno-name Aff --transpose --recode --out gwa_gabel --noweb
plink --file GWAS_clean4 --pheno pheno.phen --pheno-name systolic --transpose --recode --out gwa_gabel_qtl --noweb
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 --noweb plink --file chr22_imputation_ex --maf 0.01 --mind 0.02 --geno 0.05 --hwe 0.001 --out qc_check --noweb plink --file chr22_imputation_ex --maf 0.01 --mind 0.02 --geno 0.05 --hwe 0.001 --recode --out chr22_clean1 --noweb plink --file chr22_clean1 --maf 0.01 --mind 0.02 --geno 0.05 --hwe 0.001 --out qc_check_2 --noweb plink --file chr22_clean1 --filter-cases --hwe 0.001 --recode --out chr22_cases_clean --noweb plink --file chr22_clean1 --filter-controls --recode --out chr22_controls_clean --noweb plink --file chr22_controls_clean --merge chr22_cases_clean.ped chr22_cases_clean.map --hwe 0.001 --recode --out chr22_all_clean --noweb plink --file chr22_all_clean --logistic --out chr22_all_clean_geno --noweb R mydata = read.table("chr22_all_clean_geno.assoc.logistic", header=T) names(mydata) plot(mydata$BP, -log10(mydata$P)) smallp GWAS Data QC= 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 --noweb 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
===GWAS Data QC===
plink --file GWAS --noweb
plink --file GWAS --mind 0.10 --recode --out GWAS_clean_mind --noweb
plink --file GWAS_clean3 --exclude HWE_out.txt --recode --out GWAS_clean4 --noweb
===GWAS Control Substructure=== 
plink --file GWAS_clean4 --genome --mds-plot 10 --noweb
R