Changes

AdvGeneMap2018Commands

615 bytes added, 15:02, 23 January 2018
===Plink - Part 1 - Data QC===
  plink --file GWAS
plink --file GWAS --mind 0.10 --recode --out GWAS_clean_mind
plink --file GWAS_clean_mind --maf 0.05 --recode --out MAF_greater_5
sexcheck = read.table("GWAS_sex_checking.sexcheck", header=T)
names(sexcheck)
sex_problem = sexcheck[[which(sexcheck$STATUS=="PROBLEM"),]]
sex_problem
q()
setwd("to_your_working_directory/")
dups = read.table("duplicates.genome", header = T)
problem_pairs = dups[[which(dups$PI_HAT > 0.4),]]
problem_pairs
problem_pairs = dups[[which(dups$PI_HAT > 0.05),]]
myvars = c("FID1", "IID1", "FID2", "IID2", "PI_HAT")
problem_pairs[[myvars]]
q()
######
plink --file GWAS_clean3 --het
###### in R
Dataset<- read.table("plink.het", header=TRUE, sep="", na.strings="NA", dec=".", strip.white=TRUE) mean(Dataset$F) sd(Dataset$F) jpeg("hist.jpeg", height=1000, width=1000) hist(scale(Dataset$F), xlim=c(-4,4)) dev.off() q() ###### plink --file GWAS_clean3 --pheno pheno.txt --pheno-name Aff --hardy ##### in R hardy = read.table("plink.hwe", header = T) names(hardy) hwe_prob = hardy&#91;which(hardy$P < 0.0000009),&#93; hwe_prob q() ########## plink --file GWAS_clean3 --exclude HWE_out.txt --recode --out GWAS_clean4 
===Plink - Part 2 - Controlling for Substructure===
<nowiki>
plink --file GWAS_clean4 --genome --cluster --mds-plot 10
#### in R
mydata = read.table("mds_components.txt", header=T)
mydata$pch[&#91;mydata$Group==1 ] &#93; <-15 mydata$pch[&#91;mydata$Group==2 ] &#93; <-16 mydata$pch[&#91;mydata$Group==3 ] &#93; <-2
jpeg("mds.jpeg", height=500, width=500)
plot(mydata$C1, mydata$C2 ,pch=mydata$pch)
par(mfrow=c(2,1))
aff_unadj<-read.table("unadj.assoc.logistic", header=TRUE)
aff_unadj.add.p<-aff_unadj[&#91;aff_unadj$TEST==c("ADD"),]&#93;$P
broadqq(aff_unadj.add.p,"Some Trait Unadjusted")
aff_C1C2<-read.table("PC1-PC2.assoc.logistic", header=TRUE)
aff_C1C2.add.p<-aff_C1C2[&#91;aff_C1C2$TEST==c("ADD"),]&#93;$P
broadqq(aff_C1C2.add.p, "Some Trait Adjusted for PC1 and PC2")
dev.off()
gws_unadj = aff_unadj[&#91;which(aff_unadj$P < 0.0000001),]&#93;
gws_unadj
gws_adjusted = aff_C1C2[&#91;which(aff_C1C2$P < 0.0000001),]&#93;
gws_adjusted
<nowiki>
Bureaucrat, administrator
1,252
edits