#No documentation of basic browser up til now...

# MULTIZ with P. abyssi, P. furi, T. kodakarensis
# This is redo of an earlier run, so some files are already in place
# DONE (10/8/05), kpollard

    mkdir -p /cluster/data/pyrHor1/bed/conservation
    cd /cluster/data/pyrHor1/bed/conservation
    cp /cluster/data/pyrFur2/bed/conservation/HoxD55.q .
    cp /cluster/data/pyrFur2/bed/conservation/*.chr .
    cp /cluster/data/pyrFur2/bed/conservation/*.nib .
    faSize -detailed *.chr > chrom.sizes
    #blastz
    blastz pyrHor1.chr pyrFur2.chr Q=HoxD55.q > pyrHor1-pyrFur2.lav
    blastz pyrHor1.chr pyrAby1.chr Q=HoxD55.q > pyrHor1-pyrAby1.lav
    blastz pyrHor1.chr therKoda1.chr Q=HoxD55.q > pyrHor1-therKoda1.lav
    #lavToAxt
    /cluster/bin/i386/lavToAxt pyrHor1-pyrFur2.lav . . pyrHor1-pyrFur2.axt
    /cluster/bin/i386/lavToAxt pyrHor1-pyrAby1.lav . . pyrHor1-pyrAby1.axt
    /cluster/bin/i386/lavToAxt pyrHor1-therKoda1.lav . . pyrHor1-therKoda1.axt
    #axtBest
    axtBest pyrHor1-pyrFur2.axt pyrHor1.chr -winSize=500 -minScore=5000 pyrHor1-pyrFur2-best.axt
    axtBest pyrHor1-pyrAby1.axt pyrHor1.chr -winSize=500 -minScore=5000 pyrHor1-pyrAby1-best.axt
    axtBest pyrHor1-therKoda1.axt pyrHor1.chr -winSize=500 -minScore=5000 pyrHor1-therKoda1-best.axt
    #axtToMaf
    axtToMaf pyrHor1-pyrFur2-best.axt chrom.sizes chrom.sizes pyrHor1-pyrFur2.maf
    axtToMaf pyrHor1-pyrAby1-best.axt chrom.sizes chrom.sizes pyrHor1-pyrAby1.maf
    axtToMaf pyrHor1-therKoda1-best.axt chrom.sizes chrom.sizes pyrHor1-therKoda1.maf
    #multiz
    #delete extra header lines from maf files
    multiz pyrHor1-pyrAby1.maf pyrHor1-pyrFur2.maf - > pyrHor1-pyrAby1-pyrFur2.maf
    multiz pyrHor1-therKoda1.maf pyrHor1-pyrAby1-pyrFur2.maf - > pyrHor1-pyrAby1-pyrFur2-therKoda1.maf
    #phyloHMM
    /cluster/bin/phast/msa_view -i MAF -M pyrHor1.chr -o SS pyrHor1-pyrAby1-pyrFur2-therKoda1.maf > pyrHor1.ss
    /cluster/bin/phast/phyloFit -i SS pyrHor1.ss -t "((pyrFur2,(pyrHor1,pyrAby1)),therKoda1)" 
    /cluster/bin/phast/msa_view -i SS pyrHor1.ss --summary-only
    /cluster/bin/phast/phastCons pyrHor1.ss phyloFit.mod --gc 0.4499 \
	--target-coverage 0.7 --estimate-trees pyr-tree \
	--expected-lengths 25 --no-post-probs --ignore-missing \
	--nrates 1,1
    /cluster/bin/phast/phastCons pyrHor1.ss \
	pyr-tree.cons.mod,pyr-tree.noncons.mod \
	--target-coverage 0.7 --expected-lengths 25 \
	--viterbi pyrHor1-elements.bed --score \
	--require-informative 0 --seqname chr > cons.dat
    wigEncode cons.dat phastCons.wig phastCons.wib
    draw_tree phyloFit.mod > pyr-tree.ps
    #check similar to pyrFur2 tree
    #move data
    mkdir wib
    mv phastCons.wib wib/phastCons.wib
    mv phastCons.wig wib/phastCons.wig
    ln -s /cluster/data/pyrHor1/bed/conservation/wib/phastCons.wib /gbdb/pyrHor1/wib
    mkdir -p /gbdb/pyrHor1/pwMaf
    mkdir -p otherSpp/pyrFur2 otherSpp/pyrAby1 otherSpp/therKoda1
    mv pyrHor1-pyrFur2.maf otherSpp/pyrFur2/chr.maf
    mv pyrHor1-pyrAby1.maf otherSpp/pyrAby1/chr.maf
    mv pyrHor1-therKoda1.maf otherSpp/therKoda1/chr.maf
    ln -s /cluster/data/pyrHor1/bed/conservation/otherSpp/pyrFur2 /gbdb/pyrHor1/pwMaf/pyrFur2_pwMaf
    ln -s /cluster/data/pyrHor1/bed/conservation/otherSpp/pyrAby1 /gbdb/pyrHor1/pwMaf/pyrAby1_pwMaf
    ln -s /cluster/data/pyrHor1/bed/conservation/otherSpp/therKoda1 /gbdb/pyrHor1/pwMaf/therKoda1_pwMaf
    mkdir multiz
    mv pyrHor1-pyrAby1-pyrFur2-therKoda1.maf multiz/chr.maf
    ln -s /cluster/data/pyrHor1/bed/conservation/multiz /gbdb/pyrHor1/multizPh1Pa1Pf2Tk1
    #load
    hgLoadWiggle pyrHor1 phastCons /cluster/data/pyrHor1/bed/conservation/wib/phastCons.wig
    hgLoadMaf -warn pyrHor1 multizPh1Pa1Pf2Tk1
    hgLoadMaf -warn pyrHor1 pyrFur2_pwMaf -pathPrefix=/gbdb/pyrHor1/pwMaf/pyrFur2_pwMaf
    hgLoadMaf -warn pyrHor1 pyrAby1_pwMaf -pathPrefix=/gbdb/pyrHor1/pwMaf/pyrAby1_pwMaf
    hgLoadMaf -warn pyrHor1 therKoda1_pwMaf -pathPrefix=/gbdb/pyrHor1/pwMaf/therKoda1_pwMaf
    hgLoadBed pyrHor1 phastConsElements pyrHor1-elements.bed 
    #trackDb
    cd ~/kent/src/hg/makeDb/trackDb/archae/pyrHor1
    #trackDb.ra entry
    # track multizPh1Pa1Pf2Tk1
    # shortLabel Conservation
    # longLabel Pyrococcus 4-way multiz alignments
    # group compGeno
    # priority 10.0
    # visibility pack
    # type wigMaf 0.0 1.0
    # maxHeightPixels 100:40:11
    # wiggle phastCons
    # yLineOnOff Off
    # autoScale Off
    # pairwise pwMaf
    # speciesOrder pyrAby1 pyrFur2 therKoda1
    cvs commit -m "New multiz track" trackDb.ra
    #html page for multizPh1Pa1Pf2Tk1
    cvs add multizPh1Pa1Pf2Tk1.html
    cvs commit -m "Details page for multiz track" multizPh1Pa1Pf2Tk1.html