Changes

From Genome Analysis Wiki
Jump to navigationJump to search

Tutorial: GotCloud UW CMG (view source)

Revision as of 16:12, 6 August 2013

14,836 bytes added ,  16:12, 6 August 2013
→‎STEP 4 : EPACTS association analysis
Line 506: Line 506:  
  zcat /home/presenter02/test.v1/assoc/single.b.score.epacts.gz | awk '$9 != "NA" { print $0 }' | sort -g -k 9 | head -n 5000 > /home/presenter02/test.v1/assoc/single.b.score.epacts.top5000
 
  zcat /home/presenter02/test.v1/assoc/single.b.score.epacts.gz | awk '$9 != "NA" { print $0 }' | sort -g -k 9 | head -n 5000 > /home/presenter02/test.v1/assoc/single.b.score.epacts.top5000
 
  touch /home/presenter02/test.v1/assoc/single.b.score.epacts.OK
 
  touch /home/presenter02/test.v1/assoc/single.b.score.epacts.OK
 +
d-69-91-154-62:~ hmkang$ cat tmp.txt | perl -lane 's/test\.v1/test/g; print " $_"'
 +
uccessfully written phenotypes and 0 covariates across 99 individuals
 +
Processing chromosome 7...
 +
Finished generating EPACTS Makefile
 +
Running 1 parallel jobs of EPACTS
 +
make -f /home/presenter02/test/assoc/single.b.score.Makefile -j 1
 +
Loading required package: epactsR
 +
Sucessfully wrote ( 1388 * 10 ) matrix
 +
The following parameters are available.  Ones with "[]" are in effect:
 +
 +
Available Options
 +
    Required Parameters :
 +
                          -i [/home/presenter02/test/assoc/single.b.score.7.117000000.117500000.epacts]
 +
                          -o [/home/presenter02/test/assoc/single.b.score.7.117000000.117500000.epacts.tmp]
 +
          Gene Annotation
 +
              Parameters : -g [/home/presenter02/day2/session1/uwcmg_2013_08/epacts/bin/..//share/EPACTS/hg19_gencodeV14.txt.gz]
 +
                          -r [/home/presenter02/day2/session1/uwcmg_2013_08//examples/chr7Ref/hs37d5.chr7.fa]
 +
                          --inputFormat [epacts], --checkReference, -f [refGene]
 +
                          -p [/home/presenter02/day2/session1/uwcmg_2013_08/epacts/bin/..//share/EPACTS/priority.txt]
 +
                          -c [/home/presenter02/day2/session1/uwcmg_2013_08/epacts/bin/..//share/EPACTS/codon.txt]
 +
                          -u [], -d [], --se [], --si [], --outputFormat []
 +
  Other Annotation Tools : --genomeScore [], --bed [], --tabix []
 +
Load reference genome /home/presenter02/day2/session1/uwcmg_2013_08//examples/chr7Ref/hs37d5.chr7.fa...
 +
DONE: 1 chromosomes and 159138663 bases are loaded.
 +
Load codon file /home/presenter02/day2/session1/uwcmg_2013_08/epacts/bin/..//share/EPACTS/codon.txt...
 +
DONE: codon file loaded.
 +
Load priority file /home/presenter02/day2/session1/uwcmg_2013_08/epacts/bin/..//share/EPACTS/priority.txt...
 +
DONE: 24 priority annotation types loaded.
 +
Load gene file /home/presenter02/day2/session1/uwcmg_2013_08/epacts/bin/..//share/EPACTS/hg19_gencodeV14.txt.gz...
 +
DONE: 92627 gene loaded.
 +
DONE: Generated frequency of each annotype type in [ /home/presenter02/test/assoc/single.b.score.7.117000000.117500000.epacts.tmp.anno.frq ].
 +
DONE: Generated frequency of each highest priority annotation type in [ /home/presenter02/test/assoc/single.b.score.7.117000000.117500000.epacts.tmp.top.anno.frq ].
 +
Ts/Tv ratio: 2.02397
 +
Ts observed: 929  times; Tv observed: 459 times.
 +
DONE: Generated frequency of each base change in [ /home/presenter02/test/assoc/single.b.score.7.117000000.117500000.epacts.tmp.base.frq ].
 +
DONE: Generated frequency of each codon change in [ /home/presenter02/test/assoc/single.b.score.7.117000000.117500000.epacts.tmp.codon.frq ].
 +
DONE: Generated frequency of indel length in [ /home/presenter02/test/assoc/single.b.score.7.117000000.117500000.epacts.tmp.indel.frq ].
 +
nsnps = 866
 +
Rscript /home/presenter02/day2/session1/uwcmg_2013_08/epacts/bin/..//share/EPACTS/epactsSingle.R --vanilla /home/presenter02/day2/session1/uwcmg_2013_08/epacts/bin/../ /home/presenter02/test/assoc/single.b.score.phe NULL /home/presenter02/test/assoc/single.b.score.ind /home/presenter02/day2/session1/uwcmg_2013_08//out/snps/beagle/chr7/chr7.filtered.PASS.beagled.vcf.gz 7:117000000-117500000 /home/presenter02/test/assoc/single.b.score.7.117000000.117500000.epacts GT 1e-06 1 3 1000000000 0.5 0 FALSE single.b.score
 +
NULL
 +
NULL
 +
/home/presenter02/day2/session1/uwcmg_2013_08/epacts/bin/..//bin/epacts-anno --in /home/presenter02/test/assoc/single.b.score.7.117000000.117500000.epacts --ref /home/presenter02/day2/session1/uwcmg_2013_08//examples/chr7Ref/hs37d5.chr7.fa
 +
/home/presenter02/day2/session1/uwcmg_2013_08/epacts/bin/..//bin/anno -i /home/presenter02/test/assoc/single.b.score.7.117000000.117500000.epacts -r /home/presenter02/day2/session1/uwcmg_2013_08//examples/chr7Ref/hs37d5.chr7.fa -f refGene -g /home/presenter02/day2/session1/uwcmg_2013_08/epacts/bin/..//share/EPACTS/hg19_gencodeV14.txt.gz -c /home/presenter02/day2/session1/uwcmg_2013_08/epacts/bin/..//share/EPACTS/codon.txt -o /home/presenter02/test/assoc/single.b.score.7.117000000.117500000.epacts.tmp --inputFormat epacts -p /home/presenter02/day2/session1/uwcmg_2013_08/epacts/bin/..//share/EPACTS/priority.txt
 +
..............................................       
 +
  ...      Anno(tation)                      ...     
 +
  ...      Xiaowei Zhan, Goncalo Abecasis    ...     
 +
    ...      Speical Thanks:                    ...   
 +
    ...      Hyun Ming Kang, Yanming Li        ...   
 +
      ...      zhanxw@umich.edu                    ... 
 +
      ...      Sep 2011                            ...
 +
        ................................................
 +
                                                       
 +
DONE: 1388 varaints are annotated.
 +
DONE: Generated annotation output in [ /home/presenter02/test/assoc/single.b.score.7.117000000.117500000.epacts.tmp ].
 +
Annotation succeed!
 +
mv /home/presenter02/test/assoc/single.b.score.7.117000000.117500000.epacts.tmp /home/presenter02/test/assoc/single.b.score.7.117000000.117500000.epacts
 +
rm /home/presenter02/test/assoc/single.b.score.7.117000000.117500000.epacts.tmp.log /home/presenter02/test/assoc/single.b.score.7.117000000.117500000.epacts.tmp.top.anno.frq /home/presenter02/test/assoc/single.b.score.7.117000000.117500000.epacts.tmp.anno.frq /home/presenter02/test/assoc/single.b.score.7.117000000.117500000.epacts.tmp.base.frq /home/presenter02/test/assoc/single.b.score.7.117000000.117500000.epacts.tmp.codon.frq /home/presenter02/test/assoc/single.b.score.7.117000000.117500000.epacts.tmp.indel.frq
 +
/home/presenter02/day2/session1/uwcmg_2013_08/epacts/bin/..//bin/epacts-cat /home/presenter02/test/assoc/single.b.score.7.117000000.117500000.epacts | /home/presenter02/day2/session1/uwcmg_2013_08/epacts/bin/..//bin/bgzip -c > /home/presenter02/test/assoc/single.b.score.epacts.gz
 +
/home/presenter02/day2/session1/uwcmg_2013_08/epacts/bin/..//bin/tabix -f -pbed /home/presenter02/test/assoc/single.b.score.epacts.gz
 +
rm -f /home/presenter02/test/assoc/single.b.score.7.117000000.117500000.epacts
 +
/home/presenter02/day2/session1/uwcmg_2013_08/epacts/bin/..//bin/epacts-plot --in /home/presenter02/test/assoc/single.b.score.epacts.gz --region 7:0
 +
Rscript /home/presenter02/test/assoc/single.b.score.epacts.R --vanilla
 +
 +
export GDFONTPATH=/home/presenter02/day2/session1/uwcmg_2013_08/epacts/bin/../share/EPACTS; export GNUPLOT_FONTPATH=/home/presenter02/day2/session1/uwcmg_2013_08/epacts/bin/../share/EPACTS; export GNUPLOT_PS_DIR=/home/presenter02/day2/session1/uwcmg_2013_08/epacts/bin/../share/EPACTS; export PATH=$PATH:/home/presenter02/day2/session1/uwcmg_2013_08/epacts/bin/../bin/; export GNUPLOT_PFBTOPFA="pfbtops %s"; gnuplot /home/presenter02/test/assoc/single.b.score.epacts.cmd
 +
export GDFONTPATH=/home/presenter02/day2/session1/uwcmg_2013_08/epacts/bin/../share/EPACTS; export GNUPLOT_FONTPATH=/home/presenter02/day2/session1/uwcmg_2013_08/epacts/bin/../share/EPACTS; export GNUPLOT_PS_DIR=/home/presenter02/day2/session1/uwcmg_2013_08/epacts/bin/../share/EPACTS; export PATH=$PATH:/home/presenter02/day2/session1/uwcmg_2013_08/epacts/bin/../bin/; export GNUPLOT_PFBTOPFA="pfbtops %s"; /home/presenter02/day2/session1/uwcmg_2013_08/epacts/bin/../bin/epstopdf /home/presenter02/test/assoc/single.b.score.epacts.qq.eps
 +
export GDFONTPATH=/home/presenter02/day2/session1/uwcmg_2013_08/epacts/bin/../share/EPACTS; export GNUPLOT_FONTPATH=/home/presenter02/day2/session1/uwcmg_2013_08/epacts/bin/../share/EPACTS; export GNUPLOT_PS_DIR=/home/presenter02/day2/session1/uwcmg_2013_08/epacts/bin/../share/EPACTS; export PATH=$PATH:/home/presenter02/day2/session1/uwcmg_2013_08/epacts/bin/../bin/; export GNUPLOT_PFBTOPFA="pfbtops %s"; gnuplot /home/presenter02/test/assoc/single.b.score.epacts.cmd
 +
export GDFONTPATH=/home/presenter02/day2/session1/uwcmg_2013_08/epacts/bin/../share/EPACTS; export GNUPLOT_FONTPATH=/home/presenter02/day2/session1/uwcmg_2013_08/epacts/bin/../share/EPACTS; export GNUPLOT_PS_DIR=/home/presenter02/day2/session1/uwcmg_2013_08/epacts/bin/../share/EPACTS; export PATH=$PATH:/home/presenter02/day2/session1/uwcmg_2013_08/epacts/bin/../bin/; export GNUPLOT_PFBTOPFA="pfbtops %s"; /home/presenter02/day2/session1/uwcmg_2013_08/epacts/bin/../bin/epstopdf /home/presenter02/test/assoc/single.b.score.epacts.mh.eps
 +
rm /home/presenter02/test/assoc/single.b.score.epacts.cmd /home/presenter02/test/assoc/single.b.score.epacts.qq.eps /home/presenter02/test/assoc/single.b.score.epacts.mh.eps /home/presenter02/test/assoc/single.b.score.epacts.*.dat
 +
zcat /home/presenter02/test/assoc/single.b.score.epacts.gz | awk '$9 != "NA" { print $0 }' | sort -g -k 9 | head -n 5000 > /home/presenter02/test/assoc/single.b.score.epacts.top5000
 +
touch /home/presenter02/test/assoc/single.b.score.epacts.OK
 +
 +
To examine the top signals of association, try
 +
 +
% head assoc/single.b.score.epacts.top5000
 +
#CHROM BEGIN END MARKER_ID NS AC CALLRATE MAF PVALUE SCORE NS.CASE NS.CTRL AF.CASE AF.CTRL
 +
7 117130755 117130755 7:117130755_C/G_Intron:CFTR 99 9 1 0.045455 0.0018407 3.1148 50 49 0.09 0
 +
7 117149147 117149147 7:117149147_G/A_Nonsynonymous:CFTR 99 9 1 0.045455 0.0018407 3.1148 50 49 0.09 0
 +
7 117120141 117120141 7:117120141_G/C_Utr5:CFTR 99 14 1 0.070707 0.0034417 -2.9253 50 49 0.02 0.12245
 +
7 117356483 117356483 7:117356483_C/T_Intron:CTTNBP2 99 11 1 0.055556 0.0035695 -2.9139 50 49 0.01 0.10204
 +
7 117459023 117459023 7:117459023_A/G_Intron:CTTNBP2 99 121 1 0.38889 0.0047674 2.8223 5049 0.71 0.5102
 +
7 117430909 117430909 7:117430909_T/C_Intron:CTTNBP2 99 10 1 0.050505 0.0068926 -2.702 50 49 0.01 0.091837
 +
7 117456904 117456904 7:117456904_C/T_Intron:CTTNBP2 99 115 1 0.41919 0.0076959 2.6651 5049 0.68 0.47959
 +
7 117457141 117457141 7:117457141_G/C_Intron:CTTNBP2 99 115 1 0.41919 0.0076959 2.6651 5049 0.68 0.47959
 +
7 117467003 117467003 7:117467003_T/A_Intron:CTTNBP2 99 122 1 0.38384 0.0077344 2.6634 5049 0.71 0.52041
 +
 +
Also you can examine the PDF output file (although it is a tiny fraction of genome), if your X11 works.
 +
 +
% xpdf assoc/single.b.score.epacts.qq.pdf
 +
% xpdf assoc/single.b.score.epacts.mh.pdf
 +
 +
Running variant annotation using EPACTS is also possible
    
  % time ${GC}/epacts/bin/epacts anno --in ${GC}/out/snps/beagle/chr7/chr7.filtered.PASS.beagled.vcf.gz --out snps/chr7.filtered.PASS.beagled.anno.vcf.gz --ref ${GC}examples/chr7Ref/hs37d5.chr7.fa
 
  % time ${GC}/epacts/bin/epacts anno --in ${GC}/out/snps/beagle/chr7/chr7.filtered.PASS.beagled.vcf.gz --out snps/chr7.filtered.PASS.beagled.anno.vcf.gz --ref ${GC}examples/chr7Ref/hs37d5.chr7.fa
 +
/home/presenter02/day2/session1/uwcmg_2013_08/epacts/bin/..//bin/anno -i /home/presenter02/day2/session1/uwcmg_2013_08//out/snps/beagle/chr7/chr7.filtered.PASS.beagled.vcf.gz -r /home/presenter02/day2/session1/uwcmg_2013_08/examples/chr7Ref/hs37d5.chr7.fa -f refGene -g /home/presenter02/day2/session1/uwcmg_2013_08/epacts/bin/..//share/EPACTS/hg19_gencodeV14.txt.gz -c /home/presenter02/day2/session1/uwcmg_2013_08/epacts/bin/..//share/EPACTS/codon.txt -o snps/chr7.filtered.PASS.beagled.anno.vcf.gz --inputFormat vcf -p /home/presenter02/day2/session1/uwcmg_2013_08/epacts/bin/..//share/EPACTS/priority.txt
 +
The following parameters are available.  Ones with "[]" are in effect:
 +
 +
Available Options
 +
    Required Parameters :
 +
                          -i [/home/presenter02/day2/session1/uwcmg_2013_08//out/snps/beagle/chr7/chr7.filtered.PASS.beagled.vcf.gz]
 +
                          -o [snps/chr7.filtered.PASS.beagled.anno.vcf.gz]
 +
          Gene Annotation
 +
              Parameters : -g [/home/presenter02/day2/session1/uwcmg_2013_08/epacts/bin/..//share/EPACTS/hg19_gencodeV14.txt.gz]
 +
                          -r [/home/presenter02/day2/session1/uwcmg_2013_08/examples/chr7Ref/hs37d5.chr7.fa]
 +
                          --inputFormat [vcf], --checkReference, -f [refGene]
 +
                          -p [/home/presenter02/day2/session1/uwcmg_2013_08/epacts/bin/..//share/EPACTS/priority.txt]
 +
                          -c [/home/presenter02/day2/session1/uwcmg_2013_08/epacts/bin/..//share/EPACTS/codon.txt]
 +
                          -u [], -d [], --se [], --si [], --outputFormat []
 +
  Other Annotation Tools : --genomeScore [], --bed [], --tabix []
 +
Load reference genome /home/presenter02/day2/session1/uwcmg_2013_08/examples/chr7Ref/hs37d5.chr7.fa...
 +
DONE: 1 chromosomes and 159138663 bases are loaded.
 +
Load codon file /home/presenter02/day2/session1/uwcmg_2013_08/epacts/bin/..//share/EPACTS/codon.txt...
 +
DONE: codon file loaded.
 +
Load priority file /home/presenter02/day2/session1/uwcmg_2013_08/epacts/bin/..//share/EPACTS/priority.txt...
 +
DONE: 24 priority annotation types loaded.
 +
Load gene file /home/presenter02/day2/session1/uwcmg_2013_08/epacts/bin/..//share/EPACTS/hg19_gencodeV14.txt.gz...
 +
DONE: 92627 gene loaded.
 +
DONE: Generated frequency of each annotype type in [ snps/chr7.filtered.PASS.beagled.anno.vcf.gz.anno.frq ].
 +
DONE: Generated frequency of each highest priority annotation type in [ snps/chr7.filtered.PASS.beagled.anno.vcf.gz.top.anno.frq ].
 +
Ts/Tv ratio: 2.02397
 +
Ts observed: 929  times; Tv observed: 459 times.
 +
DONE: Generated frequency of each base change in [ snps/chr7.filtered.PASS.beagled.anno.vcf.gz.base.frq ].
 +
DONE: Generated frequency of each codon change in [ snps/chr7.filtered.PASS.beagled.anno.vcf.gz.codon.frq ].
 +
DONE: Generated frequency of indel length in [ snps/chr7.filtered.PASS.beagled.anno.vcf.gz.indel.frq ].
 +
..............................................       
 +
  ...      Anno(tation)                      ...     
 +
  ...      Xiaowei Zhan, Goncalo Abecasis    ...     
 +
    ...      Speical Thanks:                    ...   
 +
    ...      Hyun Ming Kang, Yanming Li        ...   
 +
      ...      zhanxw@umich.edu                    ... 
 +
      ...      Sep 2011                            ...
 +
        ................................................
 +
                                                       
 +
DONE: 1388 varaints are annotated.
 +
DONE: Generated annotation output in [ snps/chr7.filtered.PASS.beagled.anno.vcf.gz ].
 +
Annotation succeed!
 +
mv snps/chr7.filtered.PASS.beagled.anno.vcf.gz /home/presenter02/day2/session1/uwcmg_2013_08//out/snps/beagle/chr7/chr7.filtered.PASS.beagled.vcf.gz.tmp
 +
/home/presenter02/day2/session1/uwcmg_2013_08/epacts/bin/..//bin/bgzip -c /home/presenter02/day2/session1/uwcmg_2013_08//out/snps/beagle/chr7/chr7.filtered.PASS.beagled.vcf.gz.tmp > snps/chr7.filtered.PASS.beagled.anno.vcf.gz
 +
/home/presenter02/day2/session1/uwcmg_2013_08/epacts/bin/..//bin/tabix -pvcf -f snps/chr7.filtered.PASS.beagled.anno.vcf.gz
 +
rm /home/presenter02/day2/session1/uwcmg_2013_08//out/snps/beagle/chr7/chr7.filtered.PASS.beagled.vcf.gz.tmp
 +
rm snps/chr7.filtered.PASS.beagled.anno.vcf.gz.log snps/chr7.filtered.PASS.beagled.anno.vcf.gz.top.anno.frq snps/chr7.filtered.PASS.beagled.anno.vcf.gz.anno.frq snps/chr7.filtered.PASS.beagled.anno.vcf.gz.base.frq snps/chr7.filtered.PASS.beagled.anno.vcf.gz.codon.frq snps/chr7.filtered.PASS.beagled.anno.vcf.gz.indel.frq
 +
 +
You can look at example non-synonymous SNPs.
    
  % zcat snps/chr7.filtered.PASS.beagled.anno.vcf.gz | grep Nonsynonymous | grep CFTR | cut -f 1-8 | head -1
 
  % zcat snps/chr7.filtered.PASS.beagled.anno.vcf.gz | grep Nonsynonymous | grep CFTR | cut -f 1-8 | head -1
 +
 +
The following examples running gene-level burden test (skipped for the workshop),for your record.
    
  % ${GC}/epacts/bin/epacts make-group --vcf snps/chr7.filtered.PASS.beagled.anno.vcf.gz --out snps/chr7.filtered.PASS.beagled.anno.grp –nonsyn
 
  % ${GC}/epacts/bin/epacts make-group --vcf snps/chr7.filtered.PASS.beagled.anno.vcf.gz --out snps/chr7.filtered.PASS.beagled.anno.grp –nonsyn
Retrieved from "http://genome.sph.umich.edu/wiki/Special:MobileDiff/7758"

Navigation menu