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

# MULTIZ with P. furi, P. hori, 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/pyrAby1/bed/conservation
    cd /cluster/data/pyrAby1/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 pyrAby1.chr pyrFur2.chr Q=HoxD55.q > pyrAby1-pyrFur2.lav
    blastz pyrAby1.chr pyrHor1.chr Q=HoxD55.q > pyrAby1-pyrHor1.lav
    blastz pyrAby1.chr therKoda1.chr Q=HoxD55.q > pyrAby1-therKoda1.lav
    #lavToAxt
    /cluster/bin/i386/lavToAxt pyrAby1-pyrFur2.lav . . pyrAby1-pyrFur2.axt
    /cluster/bin/i386/lavToAxt pyrAby1-pyrHor1.lav . . pyrAby1-pyrHor1.axt
    /cluster/bin/i386/lavToAxt pyrAby1-therKoda1.lav . . pyrAby1-therKoda1.axt
    #axtBest
    axtBest pyrAby1-pyrFur2.axt pyrAby1.chr -winSize=500 -minScore=5000 pyrAby1-pyrFur2-best.axt
    axtBest pyrAby1-pyrHor1.axt pyrAby1.chr -winSize=500 -minScore=5000 pyrAby1-pyrHor1-best.axt
    axtBest pyrAby1-therKoda1.axt pyrAby1.chr -winSize=500 -minScore=5000 pyrAby1-therKoda1-best.axt
    #axtToMaf
    axtToMaf pyrAby1-pyrFur2-best.axt chrom.sizes chrom.sizes pyrAby1-pyrFur2.maf
     axtToMaf pyrAby1-pyrHor1-best.axt chrom.sizes chrom.sizes pyrAby1-pyrHor1.maf
     axtToMaf pyrAby1-therKoda1-best.axt chrom.sizes chrom.sizes pyrAby1-therKoda1.maf
    #multiz
    #delete extra header lines from maf files
    multiz pyrAby1-pyrHor1.maf pyrAby1-pyrFur2.maf - > pyrAby1-pyrHor1-pyrFur2.maf
    multiz pyrAby1-therKoda1.maf pyrAby1-pyrHor1-pyrFur2.maf - > pyrAby1-pyrHor1-pyrFur2-therKoda1.maf
    #phyloHMM
    /cluster/bin/phast/msa_view -i MAF -M pyrAby1.chr -o SS pyrAby1-pyrHor1-pyrFur2-therKoda1.maf > pyrAby1.ss
    /cluster/bin/phast/phyloFit -i SS pyrAby1.ss -t "((pyrFur2,(pyrAby1,pyrHor1)),therKoda1)" 
    /cluster/bin/phast/msa_view -i SS pyrAby1.ss --summary-only
    /cluster/bin/phast/phastCons pyrAby1.ss phyloFit.mod --gc 0.4522 \
	--target-coverage 0.7 --estimate-trees pyr-tree \
	--expected-lengths 25 --no-post-probs --ignore-missing \
	--nrates 1,1
    /cluster/bin/phast/phastCons pyrAby1.ss \
	pyr-tree.cons.mod,pyr-tree.noncons.mod \
	--target-coverage 0.7 --expected-lengths 25 \
	--viterbi pyrAby1-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/pyrAby1/bed/conservation/wib/phastCons.wib /gbdb/pyrAby1/wib
    mkdir -p /gbdb/pyrAby1/pwMaf
    mkdir -p otherSpp/pyrFur2 otherSpp/pyrHor1 otherSpp/therKoda1
    mv pyrAby1-pyrFur2.maf otherSpp/pyrFur2/chr.maf
    mv pyrAby1-pyrHor1.maf otherSpp/pyrHor1/chr.maf
    mv pyrAby1-therKoda1.maf otherSpp/therKoda1/chr.maf
    ln -s /cluster/data/pyrAby1/bed/conservation/otherSpp/pyrFur2 /gbdb/pyrAby1/pwMaf/pyrFur2_pwMaf
    ln -s /cluster/data/pyrAby1/bed/conservation/otherSpp/pyrHor1 /gbdb/pyrAby1/pwMaf/pyrHor1_pwMaf
    ln -s /cluster/data/pyrAby1/bed/conservation/otherSpp/therKoda1 /gbdb/pyrAby1/pwMaf/therKoda1_pwMaf
    mkdir multiz
    mv pyrAby1-pyrHor1-pyrFur2-therKoda1.maf multiz/chr.maf
    ln -s /cluster/data/pyrAby1/bed/conservation/multiz /gbdb/pyrAby1/multizPa1Ph1Pf2Tk1
    #load
    hgLoadWiggle pyrAby1 phastCons /cluster/data/pyrAby1/bed/conservation/wib/phastCons.wig
    hgLoadMaf -warn pyrAby1 multizPa1Ph1Pf2Tk1
    hgLoadMaf -warn pyrAby1 pyrFur2_pwMaf -pathPrefix=/gbdb/pyrAby1/pwMaf/pyrFur2_pwMaf
    hgLoadMaf -warn pyrAby1 pyrHor1_pwMaf -pathPrefix=/gbdb/pyrAby1/pwMaf/pyrHor1_pwMaf
    hgLoadMaf -warn pyrAby1 therKoda1_pwMaf -pathPrefix=/gbdb/pyrAby1/pwMaf/therKoda1_pwMaf
    hgLoadBed pyrAby1 phastConsElements pyrAby1-elements.bed 
    #trackDb
    cd ~/kent/src/hg/makeDb/trackDb/archae/pyrAby1
    #trackDb.ra entry
    # track multizPa1Ph1Pf2Tk1
    # 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 pyrHor1 pyrFur2 therKoda1
    cvs commit -m "New multiz track" trackDb.ra
    #html page for multizPa1Ph1Pf2Tk1
    cvs add multizPa1Ph1Pf2Tk1.html
    cvs commit -m "Details page for multiz track" multizPa1Ph1Pf2Tk1.html