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==
==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
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
==PLINK_R==
Introduction
snp.data = dbp[,c("affection", "rs1112")]
summary(snp.data)
snp.data[,"rs1112"] <- as.numeric(snp.data[,"rs1112"]) - 1
summary(snp.data)
result.all = glm (affection ~ rs1112, family=binomial("logit"), data=snp.data) dev.all = anova (result.all, test="Chi")
summary(result.all)
print(dev.all)
snp.data = dbp[,c("affection", "trait","sex", "age", "rs1112", "rs1117")]
summary(snp.data)
snp.data[,"rs1112"] <- as.numeric(snp.data[,"rs1112"]) - 1 snp.data[,"rs1117"] <- as.numeric(snp.data[,"rs1117"]) - 1 result.adj = glm (affection ~ sex + rs1112 , family=binomial("logit"), data=snp.data)
summary(result.adj)
result.adj = glm (affection ~ age + rs1112 , family=binomial("logit"), data=snp.data)
result.adj = lm (trait ~ sex + rs1112, data=snp.data)
summary(result.adj)
q()
Multifactorial Part 2
plink --ped dbp.cc.ped --map dbp.map --map3 --out logreg.sex.inter.add --logistic --sex --interaction --noweb
plink --ped dbp.cc.ped --map dbp.map --map3 --out logreg.snp1112.inter.add --logistic --condition rs1112 --interaction --noweb
R<br />load("dbp.R")<br />ls()<br />dbp[1:5,]<br />summary(dbp)
result.inter = glm (affection ~ sex * rs1112, family=binomial("logit"), data=snp.data)
summary(result.inter)
result.reg = glm (affection ~ sex + age + rs1112 + rs1117, family=binomial("logit"), data=snp.data)
summary(result.reg)
modelchoice.result<- step (result.reg) <br />summary(modelchoice.result) q()
GWAS Control Substructure
pseq myproj assoc --tests skat --phenotype BMI --covar SEX --mask include="DP>14" geno=DP:ge:11 file=YRI hwe=5.7e-7:1 "mac=1 maf=0.01" loc.group=refseq > SKAT_YRI.result
head -20 SKAT_YRI.result
cat SKAT_YRI.result | grep SKAT | grep -v "P=NA" | sort -k6 | head -15==RV-TDT exercise== vtools init rvtdt vtools import --format vcf data/data.vcf --build hg19 vtools phenotype --from_file data/phen.txt # variant selection vtools execute ANNOVAR geneanno vtools select variant "variant.region_type like '%splicing%'or variant.mut_type like 'nonsynonymous%' or variant.mut_type like 'frameshift%' or variant.mut_type like 'stop%'" -t func_variant # tped file vtools export func_variant --format tped --samples 'phenotype is not null' > vat_raw.tped sort -k4 -n vat_raw.tped | awk 'BEGIN{OFS="\t";prev="None";copy=1} {$2=$1"_"$4; $3=0; if($2==prev) {$2=$2"_"copy; copy=copy+1} else {prev=$2; copy=1}; print $0}' > vat_export.tped # tfam file vtools phenotype --out family sample_name pid mid sex phenotype > vat_export.tfam # anno file vtools use refGene-hg19_20130904 vtools update func_variant --set 'maf=0.001' vtools select func_variant -o chr pos refGene.name2 maf --header > vat_export.anno # Mendelian error and recode plink --noweb --tfile vat_export --recode12 --me 1 1 --set-me-missing --out "recode12_noME" sort -n -k1 -k6 -k2 recode12_noME.ped | sed 's/ /\t/g' | cut -f1,3,4,5 --complement > linkage.ped cut -f2 recode12_noME.map | awk 'BEGIN{OFS="\t";} {print "M",$0}' | sed '1i\I\tid\nA\tDisease' > linkage.dat java -Xmx10000m -jar java/linkage2beagle.jar linkage.dat linkage.ped > pre_beagle.bgl python script/pre_phase.py -i pre_beagle.bgl -a pre_beagle_withMissing.bgl java -Xmx10000m -jar java/beagle.jar missing=0 trios=pre_beagle.bgl out=bgl_phased verbose=false redundant=true gunzip bgl_phased.pre_beagle.bgl.phased.gz python script/post_phase.py -a vat_export.anno -b bgl_phased.pre_beagle.bgl.phased -o genes/ for g in `ls genes | grep tped | cut -d"." -f1 | head -20` do echo "runing rvTDT on gene "${g} rvTDT exercise_proj -G ./genes/${g}.tped -P ./data/rvtdt.phen -M ./genes/${g}.map --adapt 500 --alpha 0.00001 --permut 2000 --lower_cutoff 0 --upper_cutoff 100 --minVariants 3 --maxMissRatio 1 done # Answer vtools show tables ls genes/ | grep tped | wc cat exercise_proj_pval/*.pval | grep -v "^#" | sort -k2 cat exercise_proj_pval/*.pval | grep -v "^#" | sort -k3 # clean rm -r exercise_proj* genes/* bgl* linkage* recode12* pre_beagle* vat_export.*
==SEQPower==
spower -h
spower LOGIT Kryukov2009European1800.sfs --power 0.8 --OR_rare_detrimental 1.5 -j 1 -l 1 -o exercise
spower LNR Kryukov2009European1800.sfs --sample_size 1000 --meanshift_rare_detrimental 0.2 --method "CollapseQt --name CMC --alternative 2" -r 100 -j 4 -l 1 -o exercise
spower show exercise.csv sample* power
spower LNR Kryukov2009European1800.sfs --sample_size 1000 --meanshift_rare_detrimental 0.2 --meanshiftmax_rare_detrimental 0.5 --method "CollapseQt --alternative 2" -r 100 -j 4 -l 1 -o exercise
spower ELNR Kryukov2009European1800.sfs --sample_size 1000 --meanshift_rare_detrimental 0.2 --QT_thresholds 0.4 0.6 --method "CollapseQt --alternative 2" -r 100 -j 4 -l 1 -o exercise
spower show exercise.csv sample* power<br />spower ELNR Kryukov2009European1800.sfs --sample_size 1000 --p1 0.5 --meanshift_rare_detrimental 0.5 --QT_thresholds 0.4 0.6 --method "CollapseQt --alternative 2" -r 100 -j 4 -l 1 -o exercise
spower LOGIT Kryukov2009European1800.sfs --sample_size 1000 --OR_rare_detrimental 1.5 --method "GroupWrite ExerciseSimulation" -j 4 -o exercise -v1
spower show exercise.SEQPowerDB
for i in 1 1.5 2 2.5 3 3.5 4; do
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 <br />spower show exercise2.SEQPowerDB LOGIT method power title==Unphased==
unphased.sh
unphased mypeds.ped –marker 1 2 3 –missing –permutation 10
unphased all.ped -window 2 -LD
unphased all.ped -window 2 -LD >> results.txt
==VAT==
vtools -h
vtools init VATDemo