# can't find earlier documentation

# MULTIZ with methKand1, methJann1, methTher1 
# DONE (10/11/05), kpollard

    cd /cluster/data/metMar1/bed/conservation
    cp metJan1.chr methJann1.chr
    cp metThe1.chr methTher1.chr
    cp /cluster/data/methJann1/bed/conservation/methKand1.chr .
    #fix names to be methTher1, methKand1 manually
    sed s/metJan/methJann/ methJann1.chr > temp
    mv temp methJann1.chr
    faToNib metMar1.chr metMar1.chr.nib
    faToNib methTher1.chr methTher1.chr.nib
    faToNib methKand1.chr methKand1.chr.nib
    faToTwoBit methJann1.chr methJann1.2bit

    #chrom sizes
    faSize -detailed metMar1.chr methKand1.chr methTher1.chr methJann1.chr > chrom.sizes

    #blastz 
    blastz metMar1.chr methJann1.chr Q=HoxD55.q > metMar1-methJann1.lav
    blastz metMar1.chr methTher1.chr Q=HoxD55.q > metMar1-methTher1.lav
    blastz metMar1.chr methKand1.chr Q=HoxD55.q > metMar1-methKand1.lav

    /cluster/bin/i386/lavToAxt metMar1-methJann1.lav . methJann1.2bit metMar1-methJann1.axt
    /cluster/bin/i386/lavToAxt metMar1-methTher1.lav . . metMar1-methTher1.axt
    /cluster/bin/i386/lavToAxt metMar1-methKand1.lav . . metMar1-methKand1.axt

    axtBest metMar1-methJann1.axt metMar1.chr -winSize=500 -minScore=5000 metMar1-methJann1-best.axt
    axtBest metMar1-methTher1.axt metMar1.chr -winSize=500 -minScore=5000 metMar1-methTher1-best.axt
    axtBest metMar1-methKand1.axt metMar1.chr -winSize=500 -minScore=5000 metMar1-methKand1-best.axt

    axtToMaf metMar1-methJann1-best.axt chrom.sizes chrom.sizes metMar1-methJann1.maf
    axtToMaf metMar1-methTher1-best.axt chrom.sizes chrom.sizes metMar1-methTher1.maf
    axtToMaf metMar1-methKand1-best.axt chrom.sizes chrom.sizes metMar1-methKand1.maf

    #multiz
    #remove extra header lines
    multiz metMar1-methJann1.maf metMar1-methTher1.maf - > metMar1-methJann1-methTher1.maf
    multiz metMar1-methKand1.maf metMar1-methJann1-methTher1.maf - > metMar1-methJann1-methTher1-methKand1.maf

    #phyloHMM
    /cluster/bin/phast/msa_view -i MAF -M metMar1.chr -o SS metMar1-methJann1-methTher1-methKand1.maf > metMar1.ss
    /cluster/bin/phast/phyloFit -i SS metMar1.ss -t "(methKand1,(methTher1,(methJann1,metMar1)))" -o MmMjMtMk
    /cluster/bin/phast/msa_view -i SS metMar1.ss --summary-only
    #add GC content to next call
    /cluster/bin/phast/phastCons metMar1.ss MmMjMtMk.mod --gc 0.3854 \
    --target-coverage 0.7 --estimate-trees met-tree \
    --expected-lengths 25 --no-post-probs --ignore-missing \
    --nrates 1,1
    /cluster/bin/phast/phastCons metMar1.ss \
    met-tree.cons.mod,met-tree.noncons.mod \
    --target-coverage 0.7 --expected-lengths 25 \
    --viterbi metMar1-elements.bed --score \
    --require-informative 0 --seqname chr > cons.dat
    wigEncode cons.dat phastCons.wig phastCons.wib
    /cluster/bin/phast/draw_tree MmMjMtMk.mod > met-tree.ps 
    #make ai and jpg files in Illustrator
    cp met-tree.jpg /usr/local/apache/htdocs/images/metMar1-tree.jpg

    #move data
    mkdir wib
    mv phastCons.wib wib/phastCons.wib
    mv phastCons.wig wib/phastCons.wig
    ln -s /cluster/data/metMar1/bed/conservation/wib/phastCons.wib /gbdb/metMar1/wib
    mkdir /gbdb/metMar1/pwMaf
    mkdir -p otherSpp/methJann1 otherSpp/methTher1 otherSpp/methKand1
    mv metMar1-methTher1.maf otherSpp/methTher1/chr.maf
    mv metMar1-methJann1.maf otherSpp/methJann1/chr.maf
    mv metMar1-methKand1.maf otherSpp/methKand1/chr.maf
    ln -s /cluster/data/metMar1/bed/conservation/otherSpp/methTher1 /gbdb/metMar1/pwMaf/methTher1_pwMaf
    ln -s /cluster/data/metMar1/bed/conservation/otherSpp/methKand1 /gbdb/metMar1/pwMaf/methKand1_pwMaf
    ln -s /cluster/data/metMar1/bed/conservation/otherSpp/methJann1 /gbdb/metMar1/pwMaf/methJann1_pwMaf
    mkdir multiz
    mv metMar1-methJann1-methTher1-methKand1.maf multiz/chr.maf
    ln -s /cluster/data/metMar1/bed/conservation/multiz /gbdb/metMar1/multizMmMjMtMk

    #load
    hgLoadWiggle metMar1 phastCons /cluster/data/metMar1/bed/conservation/wib/phastCons.wig
    hgLoadMaf -warn metMar1 multizMmMjMtMk
    hgLoadMaf -warn metMar1 methTher1_pwMaf -pathPrefix=/gbdb/metMar1/pwMaf/methTher1_pwMaf
    hgLoadMaf -warn metMar1 methKand1_pwMaf -pathPrefix=/gbdb/metMar1/pwMaf/methKand1_pwMaf
    hgLoadMaf -warn metMar1 methJann1_pwMaf -pathPrefix=/gbdb/metMar1/pwMaf/methJann1_pwMaf
    hgLoadBed metMar1 phastConsElements metMar1-elements.bed 

    #trackDb
    cd ~/kent/src/hg/makeDb/trackDb/archae/metMar1
    #trackDb.ra entry
    # track multizMmMjMtMk
    # shortLabel Conservation
    # longLabel Thermoplasma/Ferroplasma/Picrophilus 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 methJann1 methTher1 methKand1
    cvs add trackDb.ra
    cvs commit -m "New multiz track" trackDb.ra
    #html page
    cvs add multizMmMjMtMk.html
    cvs commit -m "Details page for multiz track" multizMmMjMtMk.html