AdvGeneMap2018Commands - Revision history

Serveradmin: /* GxG Interaction */

2018-01-23T19:46:57Z

‎GxG Interaction

Serveradmin at 19:45, 23 January 2018

2018-01-23T19:45:31Z

Serveradmin: /* Seqspark */

2018-01-23T16:36:50Z

‎Seqspark

Serveradmin: /* VAT */

2018-01-23T16:36:33Z

‎VAT

Serveradmin: /* Functional Annotation */

2018-01-23T15:55:52Z

‎Functional Annotation

Serveradmin: /* GenABEL */

2018-01-23T15:54:22Z

‎GenABEL

Show changes

Serveradmin: /* VAT */

2018-01-23T15:53:15Z

‎VAT

Serveradmin at 15:24, 23 January 2018

2018-01-23T15:24:13Z

Show changes

Serveradmin at 15:22, 23 January 2018

2018-01-23T15:22:58Z

Serveradmin at 15:21, 23 January 2018

2018-01-23T15:21:33Z

@@ Line 181: / Line 181: @@
 ===GxG Interaction===
   ./plink --noweb --ped simcasecon.ped --map simcasecon.map --assoc
   ./plink --noweb --ped simcasecon.ped --map simcasecon.map --fast-epistasis
-  ./plink --noweb --ped simcasecon.ped --map simcasecon.map --fast-epistasis --caseonly
+  ./plink --noweb --ped simcasecon.ped --map simcasecon.map --fast-epistasis --case-only
   ./plink --noweb --ped simcasecon.ped --map simcasecon.map --epistasis
   ./plink --noweb --ped simcasecon.ped --map simcasecon.map --recodeA --out recoded
@@ Line 189: / Line 190: @@
+  R
   # The following commands are in the R environment
-  je &lt;-read.table("cassi.out", header&#61;T)
+  je &lt;-read.table("cassi.out", header=T)
   je
   library(ORMDR)
-  recoded&lt;-read.table("recoded.raw", header&#61;T)
+  recoded&lt;-read.table("recoded.raw", header=T)
   head(recoded)
-  newdata&lt;-recoded&#91;7:106&#93;
+  newdata&lt;-recoded[7:106]
   ormdrdata&lt;-cbind(newdata,recoded$PHENOTYPE-1)
-  names(ormdrdata)&#91;101&#93;&lt;-"casestatus"
+  names(ormdrdata)[101]&lt;-"casestatus"
   head(ormdrdata)
-  mdr1&lt;-mdr.c(ormdrdata, colresp&#61;101, cs&#61;1, combi&#61;1, cv.fold &#61; 10)
+  mdr1&lt;-mdr.c(ormdrdata, colresp=101, cs=1, combi=1, cv.fold = 10)
   mdr1$min.comb
-  mdr2&lt;-mdr.c(ormdrdata, colresp&#61;101, cs&#61;1, combi&#61;2, cv.fold &#61; 10)
+  mdr2&lt;-mdr.c(ormdrdata, colresp=101, cs=1, combi=2, cv.fold = 10)
   mdr2$min.comb
-  mdr3&lt;-mdr.c(ormdrdata, colresp&#61;101, cs&#61;1, combi&#61;3, cv.fold &#61; 10)
+  mdr3&lt;-mdr.c(ormdrdata, colresp=101, cs=1, combi=3, cv.fold = 10)
   mdr3$min.comb
   mdr1$test.erate
@@ Line 217: / Line 218: @@
   mdr3$best.combi
   mdr3$min.comb
-  logreg12&lt;-glm(casestatus ~ factor(snp1_2)*factor(snp2_1), family&#61;binomial,
+  logreg12&lt;-glm(casestatus ~ factor(snp1_2)*factor(snp2_1), family=binomial,
-  data&#61;ormdrdata)
+  data=ormdrdata)
   summary(logreg12)
   anova(logreg12)
-  pchisq(701.68,4,lower.tail&#61;F)
+  pchisq(701.68,4,lower.tail=F)
-  pchisq(703.82,8,lower.tail&#61;F)
+  pchisq(703.82,8,lower.tail=F)
   logreg345&lt;-glm(casestatus ~ factor(snp3_2)*factor(snp4_2)*factor(snp5_2),
-  family&#61;binomial, data&#61;ormdrdata)
+  family=binomial, data=ormdrdata)
   summary(logreg345)
   anova(logreg345)
-  pchisq(45.6,8,lower.tail&#61;F)
+  pchisq(45.6,8,lower.tail=F)
   q()
   ### The following commands are in the linux shell
   ./BEAM3 beam3data.txt -o beam3results
   ./BEAM3 beam3data.txt -o beam3results -T 10
 ===Plink - Part 1 - Data QC===

@@ Line 178: / Line 178: @@
   results(testb.qt)$P1df[results(testb.qt)$P1df < alpha]
   results(testb.qt)$Pc1df[results(testb.qt)$Pc1df < alpha]
 ===Plink - Part 1 - Data QC===

@@ Line 308: / Line 308: @@
 ===Seqspark===
   hdfs dfs -put demo.vcf.bz2
   hdfs dfs -put demo.tsv
@@ Line 314: / Line 313: @@
   seqspark qc.conf
   seqspark demo.conf
 ===VAT===

@@ Line 305: / Line 305: @@
   done
 ===VAT===
   vtools -h
   vtools init VATDemo

@@ Line 5: / Line 5: @@
 ===Functional Annotation===
+ table_annovar.pl
 ===GenABEL===

@@ Line 259: / Line 259: @@
-===VAT===
+===RV-TDT===
   vtools -h
   vtools init VATDemo
@@ Line 277: / Line 309: @@
   head GenotypeSummary.txt
   vtools output variant "max(DP)" "min(DP)" "avg(DP)" "stdev(DP)" "lower_quartile(DP)" "upper_quartile(DP)" --header
-  vtools select variant "filter&#61;’PASS’" --count
+  vtools select variant "filter=’PASS’" --count
-  vtools select variant "filter&#61;’PASS’" -o "max(DP)" "min(DP)" "avg(DP)" "stdev(DP)" "lower_quartile(DP)" "upper_quartile(DP)" --header
+  vtools select variant "filter=’PASS’" -o "max(DP)" "min(DP)" "avg(DP)" "stdev(DP)" "lower_quartile(DP)" "upper_quartile(DP)" --header
-  vtools update variant --from_stat ’total&#61;#(GT)’ ’num&#61;#(alt)’ ’het&#61;#(het)’ ’hom&#61;#(hom)’ ’other&#61;#(other)’ ’minDP&#61;min(DP_geno)’ ’maxDP&#61;max(DP_geno)’ ’meanDP&#61;avg(DP_geno)’ ’maf&#61;maf()’
+  vtools update variant --from_stat ’total=#(GT)’ ’num=#(alt)’ ’het=#(het)’ ’hom=#(hom)’ ’other=#(other)’ ’minDP=min(DP_geno)’ ’maxDP=max(DP_geno)’ ’meanDP=avg(DP_geno)’ ’maf=maf()’
   vtools show fields
   vtools show table variant
-  vtools update variant --from_stat ’totalGD10&#61;#(GT)’ ’numGD10&#61;#(alt)’ ’hetGD10&#61;#(het)’ ’homGD10&#61;#(hom)’ ’otherGD10&#61;#(other)’ ’mafGD10&#61;maf()’ --genotypes "DP_geno &gt; 10"
+  vtools update variant --from_stat ’totalGD10=#(GT)’ ’numGD10=#(alt)’ ’hetGD10=#(het)’ ’homGD10=#(hom)’ ’otherGD10=#(other)’ ’mafGD10=maf()’ --genotypes "DP_geno &gt; 10"
   vtools show fields
   vtools show table variant
   vtools output variant chr pos maf mafGD10 --header --limit 20
-  vtools phenotype --set "RACE&#61;0" --samples "filename like ’YRI%’"
+  vtools phenotype --set "RACE=0" --samples "filename like ’YRI%’"
-  vtools phenotype --set "RACE&#61;1" --samples "filename like ’CEU%’"
+  vtools phenotype --set "RACE=1" --samples "filename like ’CEU%’"
   vtools show samples --limit 10
-  vtools update variant --from_stat ’CEU_mafGD10&#61;maf()’ --genotypes ’DP_geno&gt;10’ --samples "RACE&#61;1"
+  vtools update variant --from_stat ’CEU_mafGD10=maf()’ --genotypes ’DP_geno&gt;10’ --samples "RACE=1"
-  vtools update variant --from_stat ’YRI_mafGD10&#61;maf()’ --genotypes ’DP_geno&gt;10’ --samples "RACE&#61;0"
+  vtools update variant --from_stat ’YRI_mafGD10=maf()’ --genotypes ’DP_geno&gt;10’ --samples "RACE=0"
   vtools output variant chr pos mafGD10 CEU_mafGD10 YRI_mafGD10 --header --limit 10
-  vtools phenotype --from_stat ’CEU_totalGD10&#61;#(GT)’ ’CEU_numGD10&#61;#(alt)’ --genotypes ’DP_geno&gt;10’ --samples "RACE&#61;1"
+  vtools phenotype --from_stat ’CEU_totalGD10=#(GT)’ ’CEU_numGD10=#(alt)’ --genotypes ’DP_geno&gt;10’ --samples "RACE=1"
-  vtools phenotype --from_stat ’YRI_totalGD10&#61;#(GT)’ ’YRI_numGD10&#61;#(alt)’ --genotypes ’DP_geno&gt;10’ --samples "RACE&#61;0"
+  vtools phenotype --from_stat ’YRI_totalGD10=#(GT)’ ’YRI_numGD10=#(alt)’ --genotypes ’DP_geno&gt;10’ --samples "RACE=0"
   vtools phenotype --output sample_nameCEU_totalGD10CEU_numGD10YRI_totalGD10YRI_numGD10 --header
   vtools execute ANNOVAR geneanno
@@ Line 304: / Line 336: @@
   vtools show tables
   vtools remove genotypes "DP_geno&lt;10" -v0
-  vtools select variant "mut_type like ’non%’ or mut_type like ’stop%’ or region_type&#61;’splicing’" -t v_funct
+  vtools select variant "mut_type like ’non%’ or mut_type like ’stop%’ or region_type=’splicing’" -t v_funct
   vtools show tables
   vtools show samples --limit 5
-  vtools select variant --samples "RACE&#61;1" -t CEU
+  vtools select variant --samples "RACE=1" -t CEU
   mkdir -p ceu
   cd ceu
-  vtools init ceu --parent ../ --variants CEU --samples "RACE&#61;1" --build hg19 vtools show project
+  vtools init ceu --parent ../ --variants CEU --samples "RACE=1" --build hg19 vtools show project
-  vtools select variant "CEU_mafGD10&gt;&#61;0.05" -t common_ceu
+  vtools select variant "CEU_mafGD10&gt;=0.05" -t common_ceu
   vtools select v_funct "CEU_mafGD10&lt;0.01" -t rare_ceu
   vtools use refGene
@@ Line 331: / Line 363: @@
   less EA_RV_VT.asso.res
   sort -g -k6 EA_RV_VT.asso.res | head
-  vtools select rare_ceu "refGene.name2&#61;’ABCC1’" -o chr pos ref alt CEU_mafGD10 numGD10 mut_type --header
+  vtools select rare_ceu "refGene.name2=’ABCC1’" -o chr pos ref alt CEU_mafGD10 numGD10 mut_type --header
   cd ..
-  vtools select variant --samples "RACE&#61;0" -t YRI
+  vtools select variant --samples "RACE=0" -t YRI
   mkdir -p yri
   cd yri
-  vtools init yri --parent ../ --variants YRI --samples "RACE&#61;0" --build hg19 vtools select variant "YRI_mafGD10&gt;&#61;0.05" -t common_yri vtools select v_funct "YRI_mafGD10&lt;0.01" -t rare_yri
+  vtools init yri --parent ../ --variants YRI --samples "RACE=0" --build hg19 vtools select variant "YRI_mafGD10&gt;=0.05" -t common_yri vtools select v_funct "YRI_mafGD10&lt;0.01" -t rare_yri
   vtools use refGene
   vtools associate common_yri BMI --covariate SEX -m "LinRegBurden --alternative 2" -j1 --to_db YA_CV &gt; YA_CV.asso.res

@@ Line 236: / Line 236: @@
   broadqq &lt;-function(pvals, title)
+  {
-observed &lt;- sort(pvals)
+     observed &lt;- sort(pvals)
-lobs &lt;- -(log10(observed))
+     lobs &lt;- -(log10(observed))
-expected &lt;- c(1:length(observed))
+     expected &lt;- c(1:length(observed))
-lexp &lt;- -(log10(expected / (length(expected)+1)))
+     lexp &lt;- -(log10(expected / (length(expected)+1)))
-plot(c(0,7), c(0,7), col&eq;"red", lwd&eq;3, type&eq;"l", xlab&eq;"Expected (-logP)", ylab&eq;"Observed (-logP)", xlim&eq;c(0,max(lobs)), ylim&eq;c(0,max(lobs)), las&eq;1, xaxs&eq;"i", yaxs&eq;"i", bty&eq;"l", main &eq; title)
+     plot(c(0,7), c(0,7), col&eq;"red", lwd&eq;3, type&eq;"l", xlab&eq;"Expected (-logP)", ylab&eq;"Observed (-logP)", xlim&eq;c(0,max(lobs)), ylim&eq;c(0,max(lobs)), las&eq;1, xaxs&eq;"i", yaxs&eq;"i", bty&eq;"l", main &eq; title)
-points(lexp, lobs, pch&eq;23, cex&eq;.4, bg&eq;"black") }
+     points(lexp, lobs, pch&eq;23, cex&eq;.4, bg&eq;"black") }
-−
   jpeg("qqplot_compare.jpeg", height&eq;1000, width&eq;500)
   par(mfrow&eq;c(2,1))

← Older revision		Revision as of 15:21, 23 January 2018
Line 216:		Line 216:


		+	===Plink - Part 2 - Controlling for Substructure===

−	~~===Plink - Part 2 - Controlling for Substructure===~~
−
	plink --file GWAS_clean4 --genome --cluster --mds-plot 10		plink --file GWAS_clean4 --genome --cluster --mds-plot 10
	#### in R		#### in R