2016-genetic-association-commands - Revision history

Serveradmin at 19:24, 19 May 2017

2017-05-19T19:24:40Z

Serveradmin at 19:15, 19 May 2017

2017-05-19T19:15:42Z

Serveradmin at 19:13, 19 May 2017

2017-05-19T19:13:17Z

Serveradmin at 19:12, 19 May 2017

2017-05-19T19:12:14Z

Serveradmin: /* PLINK_R */

2017-05-19T19:11:05Z

‎PLINK_R

Serveradmin: /* Unphased */

2017-01-04T20:46:56Z

‎Unphased

Serveradmin: /* GeneABEL */

2017-01-04T20:46:10Z

‎GeneABEL

Serveradmin: /* RV-TDT exercise */

2016-09-09T13:13:02Z

‎RV-TDT exercise

Serveradmin: /* RV-TDT exercise */

2016-09-09T13:12:13Z

‎RV-TDT exercise

Serveradmin at 12:25, 9 September 2016

2016-09-09T12:25:06Z

@@ Line 2: / Line 2: @@
 ==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
@@ Line 117: / Line 118: @@
   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))
+  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
@@ Line 143: / Line 143: @@
   names(dosage)
   plot(dosage$BP, -log10(dosage$P))
-  dosagep = dosage[which(dosage$P < 5E-8),]
+  dosagep = dosage[which(dosage$P < 5E-8),] dosagep = dosagep[order(dosagep$BP),]
- dosagep = dosagep[order(dosagep$BP),]
   dosagep
   interest = dosage[which(dosage$SNP=='rs715586'),]
   interest
 ==PLINK_R==
 Introduction
@@ Line 380: / Line 378: @@
   gws_unadj = aff_unadj[which(aff_unadj$P < 0.0000001),]
   gws_unadj
-  gws_adjusted = aff_C1C2[which(aff_C1C2$P < 0.0000001),] gws_adjusted q()
+  gws_adjusted = aff_C1C2[which(aff_C1C2$P < 0.0000001),]
  Multiple Testing
@@ Line 502: / Line 502: @@
   rm -r exercise_proj* genes/* bgl* linkage* recode12* pre_beagle* vat_export.*
 ==SEQPower==
   spower -h
   spower LOGIT -h

@@ Line 2: / Line 2: @@
 ==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
@@ Line 121: / Line 120: @@
 ==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
@@ Line 382: / Line 380: @@
   gws_unadj = aff_unadj[which(aff_unadj$P < 0.0000001),]
   gws_unadj
-  gws_adjusted = aff_C1C2[which(aff_C1C2$P < 0.0000001),] gws_adjusted
+  gws_adjusted = aff_C1C2[which(aff_C1C2$P < 0.0000001),] gws_adjusted q()
- q()
  Multiple Testing
@@ Line 495: / Line 492: @@
   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
+      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
@@ Line 505: / Line 502: @@
   rm -r exercise_proj* genes/* bgl* linkage* recode12* pre_beagle* vat_export.*
 ==SEQPower==
   spower -h
   spower LOGIT -h
@@ Line 535: / Line 533: @@
   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
@@ Line 543: / Line 540: @@
   unphased all.ped -window 2 -LD
   unphased all.ped -window 2 -LD &gt;&gt; results.txt
   vtools -h
   vtools init VATDemo

@@ Line 2: / Line 2: @@
 ==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
@@ Line 118: / Line 119: @@
   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
@@ Line 144: / Line 145: @@
   names(dosage)
   plot(dosage$BP, -log10(dosage$P))
-  dosagep = dosage[which(dosage$P < 5E-8),] dosagep = dosagep[order(dosagep$BP),] dosagep
+  dosagep = dosage[which(dosage$P < 5E-8),]
   interest = dosage[which(dosage$SNP=='rs715586'),]
   interest
@@ Line 379: / Line 382: @@
   gws_unadj = aff_unadj[which(aff_unadj$P < 0.0000001),]
   gws_unadj
-  gws_adjusted = aff_C1C2[which(aff_C1C2$P < 0.0000001),]
+  gws_adjusted = aff_C1C2[which(aff_C1C2$P < 0.0000001),] gws_adjusted
- gws_adjusted
   q()
@@ Line 466: / Line 468: @@
   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
@@ Line 504: / Line 505: @@
   rm -r exercise_proj* genes/* bgl* linkage* recode12* pre_beagle* vat_export.*
 ==SEQPower==
   spower -h
   spower LOGIT -h
@@ Line 534: / Line 534: @@
   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
+  done<br />spower show exercise2.SEQPowerDB LOGIT method power title==Unphased==
   unphased.sh
   unphased mypeds.ped –marker 1 2 3 –missing –permutation 10
@@ Line 545: / Line 544: @@
   unphased all.ped -window 2 -LD &gt;&gt; results.txt
   vtools -h
   vtools init VATDemo

@@ Line 2: / Line 2: @@
 ==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
@@ Line 122: / Line 121: @@
 ==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
@@ Line 146: / Line 144: @@
   names(dosage)
   plot(dosage$BP, -log10(dosage$P))
-  dosagep = dosage[which(dosage$P < 5E-8),] dosagep = dosagep[order(dosagep$BP),]
+  dosagep = dosage[which(dosage$P < 5E-8),] dosagep = dosagep[order(dosagep$BP),] dosagep
- dosagep
   interest = dosage[which(dosage$SNP=='rs715586'),]
   interest
@@ Line 382: / Line 379: @@
   gws_unadj = aff_unadj[which(aff_unadj$P < 0.0000001),]
   gws_unadj
-  gws_adjusted = aff_C1C2[which(aff_C1C2$P < 0.0000001),] gws_adjusted
+  gws_adjusted = aff_C1C2[which(aff_C1C2$P < 0.0000001),]
   q()
@@ Line 467: / Line 465: @@
   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
+  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
@@ Line 506: / Line 503: @@
   # clean
   rm -r exercise_proj* genes/* bgl* linkage* recode12* pre_beagle* vat_export.*
 ==SEQPower==
   spower -h
   spower LOGIT -h
@@ Line 540: / Line 537: @@
 ==Unphased==
   unphased.sh
   unphased mypeds.ped –marker 1 2 3 –missing –permutation 10
@@ Line 551: / Line 547: @@
 ==VAT==
   vtools -h
   vtools init VATDemo

@@ Line 305: / Line 305: @@
   snp.data = dbp[,c("affection", "rs1112")]
   summary(snp.data)
-  snp.data[,"rs1112"]<-as.numeric(snp.data[,"rs1112"])-1<br />summary(snp.data)<br />result.all=glm(affection ~ rs1112, family=binomial("logit"), data=snp.data)<br />dev.all = anova(result.all, test="Chi")<br />print(dev.all)<br />snp.data = dbp[,c("affection","trait","sex","age","rs1112","rs1117")]<br />summary(snp.data)<br />snp.data[,"rs1112"]<-as.numeric(snp.data[,"rs1112"])-1<br />snp.data[,"rs1117"]<-as.numeric(snp.data[,"rs1117"])-1<br />result.adj = glm(affection ~ sex + rs1112, family=binomial("logit"), data=snp.data)<br />summary(result.adj)
+  snp.data[,"rs1112"]<-as.numeric(snp.data[,"rs1112"])-1
   result.adj = glm (affection ~ age + rs1112 , family=binomial("logit"), data=snp.data)
   summary(result.adj)

@@ Line 454: / Line 454: @@
 ==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
@@ Line 495: / Line 483: @@
       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
@@ Line 501: / Line 488: @@
   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.*

← Older revision		Revision as of 13:12, 9 September 2016
Line 454:		Line 454:

	==RV-TDT exercise==		==RV-TDT exercise==
−	~~for g in `ls rvtdt_exercise_data \| grep tped \| cut -d"." -f1`~~
−	do
−	~~echo "runing rvTDT on gene "${g}~~
−	~~./rvTDT exercise_proj -G ./rvtdt_exercise_data/${g}.tped \~~
−	~~-P ./rvtdt_exercise_data/rvtdt_exercise.phen \~~
−	~~-M ./rvtdt_exercise_data/${g}.map \~~
−	~~--adapt 500 --alpha 0.00001 --permut 2000 \~~
−	~~--lower_cutoff 0 --upper_cutoff 100 \~~
−	~~--minVariants 3 \~~
−	~~--maxMissRatio 1~~
−	~~done~~

		+	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==		==SEQPower==