# for emacs: -*- mode: sh; -*-

#	Xenopus tropicalis
# ftp.ncbi.nlm.nih.gov/genbank/genomes/Eukaryotes/vertebrates_other/
#	Xenopus_tropicalis/v4.2/
#
# http://www.ncbi.nlm.nih.gov/Traces/wgs/?val=AAMC01
#	WGS: AAMC01000001:AAMC01190823

##########################################################################
# Download sequence (DONE - 2011-09-07 - Hiram)
    mkdir -p /hive/data/genomes/xenTro3/genbank
    cd /hive/data/genomes/xenTro3/genbank
    wget --timestamping -r --cut-dirs=6 --level=0 -nH -x \
	--no-remove-listing -np \
"ftp://ftp.ncbi.nlm.nih.gov/genbank/genomes/Eukaryotes/vertebrates_other/Xenopus_tropicalis/v4.2/*"

    faSize Primary_Assembly/unplaced_scaffolds/FASTA/unplaced.scaf.fa.gz 
1511717716 bases (153400444 N's 1358317272 real 1358317272 upper 0 lower) in 19549 sequences in 1 files

###########################################################################
# fixup contig names a little bit (DONE - 2011-09-07 - Hiram)
    cd /hive/data/genomes/xenTro3/genbank

    # removing the .1 from all the GL names
    
    zcat Primary_Assembly/unplaced_scaffolds/FASTA/unplaced.scaf.fa.gz \
	| sed -e "s/^>gi.[0-9][0-9]*.gb.GL/>GL/; s/\.1| .*//;" \
	| gzip -c > xenTro3.fa.gz
    #	real    6m24.993s

    zcat Primary_Assembly/unplaced_scaffolds/AGP/unplaced.scaf.agp.gz \
	| sed -e "s/^GL1\([0-9]*\)\.1/GL1\1/" | gzip -c > xenTro3.agp.gz

###########################################################################
#  Initial database build (DONE - 2011-09-08 - Hiram)

    cd /hive/data/genomes/xenTro3/genbank
    cat << '_EOF_' > xenTro3.config.ra
# Config parameters for makeGenomeDb.pl:
db xenTro3
clade vertebrate
genomeCladePriority 80
scientificName Xenopus (Silurana) tropicalis
commonName western clawed frog
assemblyDate Nov. 2009
assemblyLabel US DOE Joint Genome Institute (JGI-PGF) (GCA_000004195.1)
assemblyShortLabel X. tropicalis v4.2
orderKey 444
mitoAcc NC_006839
fastaFiles /cluster/data/xenTro3/genbank/xenTro3.fa.gz
agpFiles /cluster/data/xenTro3/genbank/xenTro3.agp.gz
dbDbSpeciesDir xenTro
taxId 8364
'_EOF_'
    # << happy emacs

    # verify sequence and agp are OK
    makeGenomeDb.pl -stop=agp xenTro3.config.ra > agp.log 2>&1

    time makeGenomeDb.pl -continue=db xenTro3.config.ra > db.log 2>&1
    #	real    12m7.144s

##########################################################################
# running repeat masker (DONE - 2011-09-08 - Hiram)
    mkdir /hive/data/genomes/xenTro3/bed/repeatMasker
    cd /hive/data/genomes/xenTro3/bed/repeatMasker
    time doRepeatMasker.pl -buildDir=`pwd` -noSplit \
	-species "Xenopus tropicalis" \
	-bigClusterHub=swarm -dbHost=hgwdev -workhorse=hgwdev \
	-smallClusterHub=memk xenTro3 > do.log 2>&1 &
    #	real    35m31.331s
    cat faSize.rmsk.txt
    #	1511735326 bases (153400444 N's 1358334882 real 1025832944 upper
    #	332501938 lower) in 19550 sequences in 1 files
    #	%21.99 masked total, %24.48 masked real

    grep -i versi do.log
# RepeatMasker version development-$Id: RepeatMasker,v 1.25 2010/09/08 21:32:26 angie Exp $
#    June 30 2010 (open-3-2-9) version of RepeatMasker
    featureBits -countGaps xenTro3 rmsk
    #	333219656 bases of 1511735326 (22.042%) in intersection
    # why is it different than the faSize above ?
    # because rmsk masks out some N's as well as bases, the count above
    #	separates out the N's from the bases, it doesn't show lower case N's

##########################################################################
# running simple repeat (DONE - 2011-09-08 - Hiram)
    mkdir /hive/data/genomes/xenTro3/bed/simpleRepeat
    cd /hive/data/genomes/xenTro3/bed/simpleRepeat
    time doSimpleRepeat.pl -buildDir=`pwd` -bigClusterHub=swarm \
	-dbHost=hgwdev -workhorse=hgwdev -smallClusterHub=memk \
	xenTro3 > do.log 2>&1 &
    # about 1 hour 3 minutes
    #	real    23m18.960s

    #	real    198m33.953s
    cat fb.simpleRepeat 
    #	103101528 bases of 1358346944 (7.590%) in intersection

    # add to rmsk after it is done:
    cd /hive/data/genomes/xenTro3
    twoBitMask xenTro3.rmsk.2bit \
	-add bed/simpleRepeat/trfMask.bed xenTro3.2bit
    #	you can safely ignore the warning about fields >= 13

    twoBitToFa xenTro3.2bit stdout | faSize stdin > faSize.xenTro3.2bit.txt
    cat faSize.xenTro3.2bit.txt
    #	1511735326 bases (153400444 N's 1358334882 real 1024824487 upper
    #	333510395 lower) in 19550 sequences in 1 files
    #	%22.06 masked total, %24.55 masked real

    rm /gbdb/xenTro3/xenTro3.2bit
    ln -s `pwd`/xenTro3.2bit /gbdb/xenTro3/xenTro3.2bit

#########################################################################
# Verify all gaps are marked, add any N's not in gap as type 'other'
#	(DONE - 2011-09-08 - Hiram)
    mkdir /hive/data/genomes/xenTro3/bed/gap
    cd /hive/data/genomes/xenTro3/bed/gap
    time nice -n +19 findMotif -motif=gattaca -verbose=4 \
	-strand=+ ../../xenTro3.unmasked.2bit > findMotif.txt 2>&1
    #	real    0m15.269s
    grep "^#GAP " findMotif.txt | sed -e "s/^#GAP //" > allGaps.bed
    featureBits xenTro3 -not gap -bed=notGap.bed
    #	1358346944 bases of 1358346944 (100.000%) in intersection
    featureBits xenTro3 allGaps.bed notGap.bed -bed=new.gaps.bed
    #	12062 bases of 1358346944 (0.001%) in intersection

    #	what is the highest index in the existing gap table:
    hgsql -N -e "select ix from gap;" xenTro3 | sort -n | tail -1
    #	1466
    cat << '_EOF_' > mkGap.pl
#!/bin/env perl

use strict;
use warnings;

my $ix=`hgsql -N -e "select ix from gap;" xenTro3 | sort -n | tail -1`;
chomp $ix;

open (FH,"<new.gaps.bed") or die "can not read new.gaps.bed";
while (my $line = <FH>) {
    my ($chrom, $chromStart, $chromEnd, $rest) = split('\s+', $line);
    ++$ix;
    printf "%s\t%d\t%d\t%d\tN\t%d\tother\tyes\n", $chrom, $chromStart,
        $chromEnd, $ix, $chromEnd-$chromStart;
}
close (FH);
'_EOF_'
    # << happy emacs
    chmod +x ./mkGap.pl
    ./mkGap.pl > other.bed
    featureBits -countGaps xenTro3 other.bed
    #	12062 bases of 1511735326 (0.001%) in intersection
    wc -l other.bed
    #	11242
    hgLoadBed -sqlTable=$HOME/kent/src/hg/lib/gap.sql \
	-noLoad xenTro3 otherGap other.bed
    #	starting with this many
    hgsql -e "select count(*) from gap;" xenTro3
    #	171274
    hgsql xenTro3 -e 'load data local infile "bed.tab" into table gap;'
    #	result count:
    hgsql -e "select count(*) from gap;" xenTro3
    #	182516
    # == 171274 + 11242
    # verify we aren't adding gaps where gaps already exist
    # this would output errors if that were true:
    gapToLift -minGap=1 xenTro3 nonBridged.lift -bedFile=nonBridged.bed
    # see example in danRer7.txt
    # there are no non-bridged gaps here:
    hgsql -N -e "select bridge from gap;" xenTro3 | sort | uniq -c
    #	182516 yes

##########################################################################
## WINDOWMASKER (DONE - 2011-09-08 - Hiram)
    mkdir /hive/data/genomes/xenTro3/bed/windowMasker
    cd /hive/data/genomes/xenTro3/bed/windowMasker
    time nice -n +19 doWindowMasker.pl -buildDir=`pwd` -workhorse=hgwdev \
	-dbHost=hgwdev xenTro3 > do.log 2>&1 &
    #	about 1h 52m

    # Masking statistics
    twoBitToFa xenTro3.wmsk.2bit stdout | faSize stdin
    #	1511735326 bases (153400444 N's 1358334882 real 845383044 upper
    #	512951838 lower) in 19550 sequences in 1 files
    #	%33.93 masked total, %37.76 masked real

    twoBitToFa xenTro3.wmsk.sdust.2bit stdout | faSize stdin
    #	1511735326 bases (153400444 N's 1358334882 real 835449906 upper
    #	522884976 lower) in 19550 sequences in 1 files
    #	%34.59 masked total, %38.49 masked real

    hgLoadBed xenTro3 windowmaskerSdust windowmasker.sdust.bed.gz
    #	Loaded 7010382 elements of size 3

    featureBits -countGaps xenTro3 windowmaskerSdust
    #	811865302 bases of 1799143587 (45.125%) in intersection

    #	eliminate the gaps from the masking
    featureBits xenTro3 -not gap -bed=notGap.bed
    #	1701353770 bases of 1701353770 (100.000%) in intersection
    time nice -n +19 featureBits xenTro3 windowmaskerSdust notGap.bed \
        -bed=stdout | gzip -c > cleanWMask.bed.gz
    #	714075485 bases of 1701353770 (41.971%) in intersection
    #	reload track to get it clean
    hgLoadBed xenTro3 windowmaskerSdust cleanWMask.bed.gz
    #	Loaded 8199653 elements of size 4
    featureBits -countGaps xenTro3 windowmaskerSdust
    #	714075485 bases of 1799143587 (39.690%) in intersection
    #	mask the sequence with this clean mask
    zcat cleanWMask.bed.gz \
	| twoBitMask ../../xenTro3.unmasked.2bit stdin \
	    -type=.bed xenTro3.cleanWMSdust.2bit
    twoBitToFa xenTro3.cleanWMSdust.2bit stdout | faSize stdin \
        > xenTro3.cleanWMSdust.faSize.txt
    cat xenTro3.cleanWMSdust.faSize.txt
    #	1799143587 bases (97789820 N's 1701353767 real 987278282 upper
    #	714075485 lower) in 6457 sequences in 1 files
    #	%39.69 masked total, %41.97 masked real

    # how much does this window masker and repeat masker overlap:
    featureBits -countGaps xenTro3 rmsk windowmaskerSdust
    #	280229203 bases of 1511735326 (18.537%) in intersection

#########################################################################
# MASK SEQUENCE WITH WM+TRF (DONE - 2011-09-09 - Hiram)
    cd /hive/data/genomes/xenTro3
    twoBitMask -add bed/windowMasker/xenTro3.cleanWMSdust.2bit \
	bed/simpleRepeat/trfMask.bed xenTro3.2bit
    #	safe to ignore the warnings about BED file with >=13 fields
    twoBitToFa xenTro3.2bit stdout | faSize stdin > faSize.xenTro3.txt
    cat faSize.xenTro3.txt
    #	1799143587 bases (97789820 N's 1701353767 real 987044085 upper
    #	714309682 lower) in 6457 sequences in 1 files
    #	%39.70 masked total, %41.98 masked real

    #	create symlink to gbdb
    ssh hgwdev
    rm /gbdb/xenTro3/xenTro3.2bit
    ln -s `pwd`/xenTro3.2bit /gbdb/xenTro3/xenTro3.2bit

########################################################################
# MAKE 11.OOC FILE FOR BLAT/GENBANK (DONE - 2011-09-09 - Hiram)
    # Use -repMatch=650, based on size -- for human we use 1024
    # use the "real" number from the faSize measurement,
    # hg19 is 2897316137, calculate the ratio factor for 1024:
    calc \( 1358317272 / 2897316137 \) \* 1024
    #	( 1358317272 / 2897316137 ) * 1024 = 480.070804
    # round up to 500

    cd /hive/data/genomes/xenTro3
    blat xenTro3.2bit /dev/null /dev/null -tileSize=11 \
      -makeOoc=jkStuff/xenTro3.11.ooc -repMatch=500
    #	Wrote 29991 overused 11-mers to jkStuff/xenTro3.11.ooc

    #	copy all of this stuff to the klusters:
    # there are no non-bridged gaps
    hgsql -N -e "select bridge from gap;" xenTro3 | sort | uniq -c
    #	182516 yes
    #	cd /hive/data/genomes/xenTro3/jkStuff
    #	gapToLift xenTro3 nonBridged.lift -bedFile=nonBridged.bed
    cd /hive/data/genomes/xenTro3
    mkdir /hive/data/staging/data/xenTro3
    cp -p jkStuff/xenTro3.11.ooc chrom.sizes \
	xenTro3.2bit /hive/data/staging/data/xenTro3
    # request rsync copy from cluster admin

#########################################################################
# AUTO UPDATE GENBANK (DONE - 2011-09-11 - Hiram)
    # examine the file:
    /cluster/data/genbank/data/organism.lst
    # for your species to see what counts it has for:
# organism       mrnaCnt estCnt  refSeqCnt
# Xenopus (Silurana) tropicalis   18786   1271379 8783
    # to decide which "native" mrna or ests you want to specify in genbank.conf
    # this appears that xenTro3 has plenty of native est's

    ssh hgwdev  
    cd $HOME/kent/src/hg/makeDb/genbank
    git pull
    # edit etc/genbank.conf to add xenTro3 before anoCar1 and commit to GIT
# xenTro3 (X. tropicalis) 19550 scaffolds
xenTro3.serverGenome = /hive/data/genomes/xenTro3/xenTro3.2bit
xenTro3.clusterGenome = /scratch/data/xenTro3/xenTro3.2bit
xenTro3.ooc = /scratch/data/xenTro3/xenTro3.11.ooc
xenTro3.lift = no
xenTro3.refseq.mrna.native.pslCDnaFilter  = ${lowCover.refseq.mrna.native.pslCDnaFilter}
xenTro3.refseq.mrna.xeno.pslCDnaFilter    = ${lowCover.refseq.mrna.xeno.pslCDnaFilter}
xenTro3.genbank.mrna.native.pslCDnaFilter = ${lowCover.genbank.mrna.native.pslCDnaFilter}
xenTro3.genbank.mrna.xeno.pslCDnaFilter   = ${lowCover.genbank.mrna.xeno.pslCDnaFilter}
xenTro3.genbank.est.native.pslCDnaFilter  = ${lowCover.genbank.est.native.pslCDnaFilter}
xenTro3.refseq.mrna.native.load = yes
xenTro3.genbank.est.native.load = yes
xenTro3.refseq.mrna.xeno.load = no
xenTro3.genbank.mrna.xeno.load = no
xenTro3.downloadDir = xenTro3
xenTro3.perChromTables = no
xenTro3.mgc = yes
# xenTro3.upstreamGeneTbl = mgcGenes
# xenTro3.upstreamMaf = multiz9way
# /hive/data/genomes/xenTro3/bed/multiz9way/species.lst

    # end of section added to etc/genbank.conf
    git commit -m "adding xenTro3 frog" genbank.conf
    git push
    make etc-update

    # ~/kent/src/hg/makeDb/genbank/src/lib/gbGenome.c already contains
    # anoCar genome information, if this is a new species, need to add stuff
    # there

    ssh hgwdev			# used to do this on "genbank" machine
    screen			# long running job managed in screen
    cd /cluster/data/genbank
    time nice -n +19 ./bin/gbAlignStep -initial xenTro3 &
    #	var/build/logs/2011.09.12-09:03:09.xenTro3.initalign.log
    #	real    360m20.657s

    # load database when finished
    ssh hgwdev
    cd /cluster/data/genbank
    time nice -n +19 ./bin/gbDbLoadStep -drop -initialLoad xenTro3 &
    #	logFile: var/dbload/hgwdev/logs/2011.09.12-20:47:05.dbload.log
    #	real    118m27.671s

    # enable daily alignment and update of hgwdev (DONE - 2011-09-13 - Hiram)
    cd ~/kent/src/hg/makeDb/genbank
    git pull
    # add xenTro3 to:
        etc/align.dbs
        etc/hgwdev.dbs
    git commit -m "Added xenTro3." etc/align.dbs etc/hgwdev.dbs
    git push
    make etc-update

#########################################################################
# construct ctgPos2 table (DONE - 2011-09-12 - Hiram)
    mkdir /hive/data/genomes/xenTro3/bed/ctgPos2
    cd /hive/data/genomes/xenTro3/bed/ctgPos2
    cat << '_EOF_' > chromSizesToCtgPos.pl
#!/usr/bin/env perl

use strict;
use warnings;

my %chromSizes;
open (FH, "<../../chrom.sizes") or die "can not read ../../chrom.sizes";
while (my $line = <FH>) {
    chomp $line;
    my ($contig, $size) = split('\s', $line);
    $chromSizes{$contig} = $size;
}
close (FH);

open (FH, "<../../genbank/Primary_Assembly/localID2acc") or
        die "can not read ../../genbank/Primary_Assembly/localID2acc";
while (my $line = <FH>) {
    chomp $line;
    my ($scaffold, $contig) = split('\s', $line);
    $scaffold =~ s/XENTR//;
    $contig =~ s/.1$//;
    my $size = $chromSizes{$contig};
    printf "%s\t%d\t%s\t0\t%d\tW\n", $scaffold, $size, $contig, $size;
}
close (FH);
'_EOF_'
    # << happy emacs
    chmod +x chromSizesToCtgPos.pl
    ./chromSizesToCtgPos.pl > xenTro3.ctgPos2.tab
    echo -e "AY789013\t17610\tchrM\t0\t17610\tF" >> xenTro3.ctgPos2.tab
    # check how big are the names:
    awk '{print length($1)}' xenTro3.ctgPos2.tab | sort -rn | head -1
    # 14
    awk '{print length($3)}' xenTro3.ctgPos2.tab | sort -rn | head
    # 8
    sed -e "s/20/14/; s/16/8/" $HOME/kent/src/hg/lib/ctgPos2.sql > ctgPos2.sql
    hgLoadSqlTab xenTro3 ctgPos2 ctgPos2.sql xenTro3.ctgPos2.tab

############################################################################
# running cpgIsland business (DONE - 2011-09-13 - Hiram)
    mkdir /hive/data/genomes/xenTro3/bed/cpgIsland
    cd /hive/data/genomes/xenTro3/bed/cpgIsland
    # use a previous binary for this program
    ln -s ../../../xenTro2/bed/cpgIsland/cpglh.exe .

    mkdir -p hardMaskedFa
    cut -f1 ../../chrom.sizes | while read C
do
    echo ${C}
    twoBitToFa ../../xenTro3.2bit:$C stdout \
	| maskOutFa stdin hard hardMaskedFa/${C}.fa
done

    ssh swarm
    cd /hive/data/genomes/xenTro3/bed/cpgIsland
    mkdir results
    cut -f1 ../../chrom.sizes > chr.list
    cat << '_EOF_' > template
#LOOP
./runOne $(root1) {check out exists results/$(root1).cpg}
#ENDLOOP
'_EOF_'
    # << happy emacs

    #	the faCount business is to make sure there is enough sequence to
    #	work with in the fasta.  cpglh.exe does not like files with too many
    #	N's - it gets stuck.
    cat << '_EOF_' > runOne
#!/bin/csh -fe
set C = `faCount hardMaskedFa/$1.fa | egrep -v "^#seq|^total" | awk '{print  $2 - $7 }'`
if ( $C > 200 ) then
    ./cpglh.exe hardMaskedFa/$1.fa > /scratch/tmp/$1.$$
    mv /scratch/tmp/$1.$$ $2
else
    touch $2
endif
'_EOF_'
    # << happy emacs
    chmod +x runOne

    gensub2 chr.list single template jobList
    para create jobList
    para try
    para check ... etc
    para time
# Completed: 1133 of 1133 jobs
# CPU time in finished jobs:        102s       1.70m     0.03h    0.00d  0.000 y
# IO & Wait Time:                  4469s      74.49m     1.24h    0.05d  0.000 y
# Average job time:                   4s       0.07m     0.00h    0.00d
# Longest finished job:              10s       0.17m     0.00h    0.00d
# Submission to last job:           485s       8.08m     0.13h    0.01d

    # Transform cpglh output to bed +
    catDir results | awk '{
$2 = $2 - 1;
width = $3 - $2;
printf("%s\t%d\t%s\t%s %s\t%s\t%s\t%0.0f\t%0.1f\t%s\t%s\n",
       $1, $2, $3, $5,$6, width,
       $6, width*$7*0.01, 100.0*2*$6/width, $7, $9);
}' > cpgIsland.bed

    # verify longest unique chrom name:
    cut -f1 cpgIsland.bed | awk '{print length($0)}' | sort -rn | head -1
    #	8
    # update the length 14 in the template to be 16:
    sed -e "s/14/8/" $HOME/kent/src/hg/lib/cpgIslandExt.sql > cpgIslandExt.sql

    cd /hive/data/genomes/xenTro3/bed/cpgIsland
    hgLoadBed xenTro3 cpgIslandExt -tab -sqlTable=cpgIslandExt.sql cpgIsland.bed
# Loaded 15341 elements of size 10

    featureBits xenTro3 cpgIslandExt
    #	9355659 bases of 1358334882 (0.689%) in intersection

    # there should be no output from checkTableCoords:
    checkTableCoords -verboseBlocks -table=cpgIslandExt xenTro3

    #	cleanup
    rm -fr hardMaskedFa

#########################################################################
# GENSCAN GENE PREDICTIONS (DONE - 2011-09-13 - Hiram)
    mkdir /hive/data/genomes/xenTro3/bed/genscan
    cd /hive/data/genomes/xenTro3/bed/genscan
    # use a previously existing genscan binary
    ln -s ../../../xenTro2/bed/genscan/hg3rdParty .

    # create hard masked .fa files
    mkdir -p hardMaskedFa
    cut -f1 ../../chrom.sizes | while read C
do
    echo ${C}
    twoBitToFa ../../xenTro3.2bit:$C stdout \
	| maskOutFa stdin hard hardMaskedFa/${C}.fa
done

    # Generate a list file, genome.list, of all the hard-masked contig chunks:
    find ./hardMaskedFa/ -type f | sed -e 's#^./##' > genome.list

    wc -l genome.list
    #	19550 genome.list

    # Run on small cluster (more mem than big cluster).
    ssh swarm
    cd /hive/data/genomes/xenTro3/bed/genscan
    # Make 3 subdirectories for genscan to put their output files in
    mkdir gtf pep subopt
    # Create template file, template, for gensub2.  For example (3-line file):
    cat << '_EOF_' > template
#LOOP
/cluster/bin/x86_64/gsBig {check in exists+ $(path1)} {check out exists gtf/$(root1).gtf} -trans={check out exists pep/$(root1).pep} -subopt={check out exists subopt/$(root1).bed} -exe=hg3rdParty/genscanlinux/genscan -par=hg3rdParty/genscanlinux/HumanIso.smat -tmp=/tmp -window=2400000
#ENDLOOP
'_EOF_'
    # << emacs
    gensub2 genome.list single template jobList
    para create jobList
    para try
    para check ... etc...
    para time
# Completed: 19550 of 19550 jobs
# CPU time in finished jobs:      25171s     419.52m     6.99h    0.29d  0.001 y
# IO & Wait Time:                 51770s     862.83m    14.38h    0.60d  0.002 y
# Average job time:                   4s       0.07m     0.00h    0.00d
# Longest finished job:             176s       2.93m     0.05h    0.00d
# Submission to last job:           336s       5.60m     0.09h    0.00d

    # this did not work, runs out of file handles ?
    find ./gtf -type f | xargs -n 256 endsInLf -zeroOk

    # Concatenate results:
    cd /hive/data/genomes/xenTro3/bed/genscan
    find ./gtf -type f | xargs cat > genscan.gtf
    find ./pep -type f | xargs cat > genscan.pep
    find ./subopt -type f | xargs cat > genscanSubopt.bed

    # Load into the database (without -genePredExt because no frame info):
    # Don't load the Pep anymore -- redundant since it's from genomic.
    ssh hgwdev
    cd /hive/data/genomes/xenTro3/bed/genscan
    # to construct a local file with the genePred business:
    gtfToGenePred genscan.gtf genscan.gp
    # this produces exactly the same thing and loads the table:
    ldHgGene -gtf xenTro3 genscan genscan.gtf
    #	Read 47412 transcripts in 293762 lines in 1 files
    #	47412 groups 9702 seqs 1 sources 1 feature types
    #	47412 gene predictions
    hgLoadBed xenTro3 genscanSubopt genscanSubopt.bed
    #	Loaded 369728 elements of size 6
    featureBits xenTro3 genscan
    #	62125158 bases of 1358334882 (4.574%) in intersection
    # previously:
    featureBits xenTro2 genscan
    #	62149060 bases of 1359412157 (4.572%) in intersection
    featureBits xenTro1 genscan
    #	63766483 bases of 1381238994 (4.617%) in intersection

#########################################################################
# CREATE MICROSAT TRACK (DONE - 2011-09-15 - Hiram
     ssh hgwdev
     mkdir /cluster/data/xenTro3/bed/microsat
     cd /cluster/data/xenTro3/bed/microsat
     awk '($5==2 || $5==3) && $6 >= 15 && $8 == 100 && $9 == 0 {printf("%s\t%s\t%s\t%dx%s\n", $1, $2, $3, $6, $16);}' \
	../simpleRepeat/simpleRepeat.bed > microsat.bed 
    hgLoadBed xenTro3 microsat microsat.bed
    #	Loaded 13365 elements of size 4

#########################################################################
# create ucscToEnsembl scaffold name mapping (DONE - 2011-09-15 - Hiram)
    # this allows the "ensembl" blue bar button to appear
    mkdir /hive/data/genomes/xenTro3/bed/ucscToEnsembl
    cd /hive/data/genomes/xenTro3/bed/ucscToEnsembl
    awk '{printf "%s\t%s.1\n", $1, $1}' ../../chrom.sizes > ucscToEnsembl.tab

    cat << '_EOF_' > ucscToEnsembl.sql
# UCSC to Ensembl chr name translation
CREATE TABLE ucscToEnsembl (
    ucsc varchar(255) not null,        # UCSC chromosome name
    ensembl varchar(255) not null,     # Ensembl chromosome name
              #Indices
    PRIMARY KEY(ucsc(8))
);
'_EOF_'

    hgsql xenTro3 < ucscToEnsembl.sql
    hgsql xenTro3 \
-e 'LOAD DATA LOCAL INFILE "ucscToEnsembl.tab" INTO TABLE ucscToEnsembl'

    # verify the blue bar "ensembl" link is now available
#########################################################################
# constructing downloads (DONE - 2011-09-16 - Hiram)
    cd /hive/data/genomes/xenTro3
    time makeDownloads.pl -workhorse=hgwdev xenTro3
    #	real    14m44.729s
    # edit the README files in goldenPath/*/README.txt

#########################################################################
# create pushQ entry (DONE - 2011-09-16 - Hiram)
    # first make sure all.joiner is up to date and has this new organism
    # a keys check should be clean:
    cd ~/kent/src/hg/makeDb/schema
    joinerCheck -database=xenTro3 -keys all.joiner

    mkdir /hive/data/genomes/xenTro3/pushQ
    cd /hive/data/genomes/xenTro3/pushQ
    makePushQSql.pl xenTro3 > xenTro3.sql 2> stderr.out
    # check stderr.out for no significant problems, it is common to see:
# WARNING: hgwdev does not have /gbdb/xenTro3/wib/gc5Base.wib
# WARNING: hgwdev does not have /gbdb/xenTro3/wib/quality.wib
# WARNING: hgwdev does not have /gbdb/xenTro3/bbi/quality.bw
# WARNING: xenTro3 does not have seq
# WARNING: xenTro3 does not have extFile
    # which are no real problem
    # if some tables are not identified:
# WARNING: Could not tell (from trackDb, all.joiner and hardcoded lists of
# supporting and genbank tables) which tracks to assign these tables to:
#	ensGtp
#	ensPep
#	ensemblSource
#	ensemblToGeneName
#	ucscToEnsembl
    # put them in manually after loading the pushQ entry
    scp -p xenTro3.sql hgwbeta:/tmp
    ssh hgwbeta
    cd /tmp
    hgsql qapushq < xenTro3.sql

#########################################################################
# HUMAN (hg18) PROTEINS TRACK (DONE 2011-09-17 braney )
    # bash  if not using bash shell already
    cd /cluster/data/xenTro3
    mkdir /cluster/data/xenTro3/blastDb

    awk '{if ($2 > 1000000) print $1}' chrom.sizes > 1meg.lst
    twoBitToFa -seqList=1meg.lst  xenTro3.2bit temp.fa
    faSplit gap temp.fa 1000000 blastDb/x -lift=blastDb.lft
# 1353 pieces of 1353 written
    rm temp.fa 1meg.lst

    awk '{if ($2 <= 1000000) print $1}' chrom.sizes > less1meg.lst
    twoBitToFa -seqList=less1meg.lst  xenTro3.2bit temp.fa
    faSplit about temp.fa 1000000 blastDb/y 
    rm temp.fa less1meg.lst

    cd blastDb
    for i in *.fa
    do
	/hive/data/outside/blast229/formatdb -i $i -p F
    done
    rm *.fa
    ls *.nsq | wc -l
# 1808

    mkdir -p /cluster/data/xenTro3/bed/tblastn.hg18KG
    cd /cluster/data/xenTro3/bed/tblastn.hg18KG
    echo  ../../blastDb/*.nsq | xargs ls -S | sed "s/\.nsq//"  > query.lst
    wc -l query.lst
# 1808 query.lst

   # we want around 250000 jobs
   calc `wc /cluster/data/hg18/bed/blat.hg18KG/hg18KG.psl | awk '{print $1}'`/\(250000/`wc query.lst | awk '{print $1}'`\)
# 36727/(250000/1808) = 265.609664

   mkdir -p kgfa
   split -l 266 /cluster/data/hg18/bed/blat.hg18KG/hg18KG.psl  kgfa/kg
   cd kgfa
   for i in *; do 
     nice pslxToFa $i $i.fa; 
     rm $i; 
   done
   cd ..
   ls -1S kgfa/*.fa > kg.lst
   wc kg.lst
#  139  139 1807 kg.lst

   mkdir -p blastOut
   for i in `cat kg.lst`; do  mkdir blastOut/`basename $i .fa`; done
   tcsh
   cd /cluster/data/xenTro3/bed/tblastn.hg18KG
   cat << '_EOF_' > blastGsub
#LOOP
blastSome $(path1) {check in line $(path2)} {check out exists blastOut/$(root2)/q.$(root1).psl }
#ENDLOOP
'_EOF_'

   cat << '_EOF_' > blastSome
#!/bin/sh
BLASTMAT=/hive/data/outside/blast229/data
export BLASTMAT
g=`basename $2`
f=/tmp/`basename $3`.$g
for eVal in 0.01 0.001 0.0001 0.00001 0.000001 1E-09 1E-11
do
if /hive/data/outside/blast229/blastall -M BLOSUM80 -m 0 -F no -e $eVal -p tblastn -d $1 -i $2 -o $f.8
then
        mv $f.8 $f.1
        break;
fi
done
if test -f  $f.1
then
    if /cluster/bin/i386/blastToPsl $f.1 $f.2
    then
	liftUp -nosort -type=".psl" -nohead $f.3 /cluster/data/xenTro3/blastDb.lft carry $f.2
        liftUp -nosort -type=".psl" -pslQ -nohead $3.tmp /cluster/data/hg18/bed/blat.hg18KG/protein.lft warn $f.3
        if pslCheck -prot $3.tmp
        then                  
            mv $3.tmp $3     
            rm -f $f.1 $f.2 $f.3 $f.4
        fi
        exit 0               
    fi                      
fi                         
rm -f $f.1 $f.2 $3.tmp $f.8 $f.3 $f.4
exit 1
'_EOF_'
    # << happy emacs
    chmod +x blastSome
    exit 
    
    ssh swarm
    cd /cluster/data/xenTro3/bed/tblastn.hg18KG
    gensub2 query.lst kg.lst blastGsub blastSpec
    para create blastSpec
#    para try, check, push, check etc.

    para time

# Completed: 251312 of 251312 jobs
# CPU time in finished jobs:    8416814s  140280.23m  2338.00h   97.42d  0.267 y
# IO & Wait Time:               1310570s   21842.83m   364.05h   15.17d  0.042 y
# Average job time:                  39s       0.65m     0.01h    0.00d
# Longest finished job:             237s       3.95m     0.07h    0.00d
# Submission to last job:         19145s     319.08m     5.32h    0.22d

    cd /cluster/data/xenTro3/bed/tblastn.hg18KG
    mkdir chainRun
    cd chainRun
    tcsh
    cat << '_EOF_' > chainGsub
#LOOP
chainOne $(path1)
#ENDLOOP
'_EOF_'

    cat << '_EOF_' > chainOne
(cd $1; cat q.*.psl | simpleChain -prot -outPsl -maxGap=150000 stdin ../c.`basename $1`.psl)
'_EOF_'
    chmod +x chainOne
    ls -1dS ../blastOut/kg?? > chain.lst
    gensub2 chain.lst single chainGsub chainSpec
    # do the cluster run for chaining
    para create chainSpec
    para try, check, push, check etc.

# Completed: 139 of 139 jobs
# CPU time in finished jobs:        562s       9.37m     0.16h    0.01d  0.000 y
# IO & Wait Time:                 30974s     516.23m     8.60h    0.36d  0.001 y
# Average job time:                 227s       3.78m     0.06h    0.00d
# Longest finished job:             299s       4.98m     0.08h    0.00d
# Submission to last job:           310s       5.17m     0.09h    0.00d

    cd /cluster/data/xenTro3/bed/tblastn.hg18KG/blastOut
    for i in kg??
    do
       cat c.$i.psl | awk "(\$13 - \$12)/\$11 > 0.6 {print}" > c60.$i.psl
       sort -rn c60.$i.psl | pslUniq stdin u.$i.psl
       awk "((\$1 / \$11) ) > 0.60 { print   }" c60.$i.psl > m60.$i.psl
       echo $i
    done
    sort u.*.psl m60* | uniq | sort -T /tmp -k 14,14 -k 16,16n -k 17,17n > ../blastHg18KG.psl
    cd ..
    pslCheck blastHg18KG.psl
# checked: 41882 failed: 0 errors: 0

    # load table 
    ssh hgwdev
    cd /cluster/data/xenTro3/bed/tblastn.hg18KG
    hgLoadPsl xenTro3 blastHg18KG.psl

    # check coverage
    featureBits xenTro3 blastHg18KG 
# 20459423 bases of 1358334882 (1.506%) in intersection

    featureBits xenTro3 blastHg18KG ensGene  -enrichment
# blastHg18KG 1.506%, ensGene 2.717%, both 1.337%, cover 88.77%, enrich 32.67x

    rm -rf blastOut
#end tblastn

##########################################################################
#  BLATSERVERS ENTRY (DONE - 2011-09-20 - Hiram)
#	After getting a blat server assigned by the Blat Server Gods,
    ssh hgwdev

    hgsql -e 'INSERT INTO blatServers (db, host, port, isTrans, canPcr) \
	VALUES ("xenTro3", "blat12", "17828", "1", "0"); \
	INSERT INTO blatServers (db, host, port, isTrans, canPcr) \
	VALUES ("xenTro3", "blat12", "17829", "0", "1");' \
	    hgcentraltest
    #	test it with some sequence

############################################################################
# set default position the same as was xenTro2 via blat
#	(DONE - 2011-09-20 - Hiram)
    hgsql -e \
'update dbDb set defaultPos="GL172663:565941-578374" where name="xenTro3";' \
	hgcentraltest

############################################################################
# Swap lastz galGal3 Chicken (DONE - 2011-09-20 - Hiram)
    # original alignment
    cd /hive/data/genomes/galGal3/bed/lastzXenTro3.2011-09-20
    cat fb.galGal3.chainXenTro3Link.txt 
    #	54978665 bases of 1042591351 (5.273%) in intersection

    # and this swap
    mkdir /hive/data/genomes/xenTro3/bed/blastz.galGal3.swap
    cd /hive/data/genomes/xenTro3/bed/blastz.galGal3.swap
    time nice -n +19 doBlastzChainNet.pl -verbose=2 \
	/hive/data/genomes/galGal3/bed/lastzXenTro3.2011-09-20/DEF \
	-noLoadChainSplit -chainMinScore=5000 -chainLinearGap=loose \
	-swap -workhorse=hgwdev -smallClusterHub=memk -bigClusterHub=swarm \
	> swap.log 2>&1 &
    #	real    13m36.708s

    cat fb.xenTro3.chainGalGal3Link.txt
    #	64250673 bases of 1358334882 (4.730%) in intersection

#########################################################################
# Swap lastz anoCar2 Lizard (DONE - 2011-09-21 - Hiram)
    # original alignment
    cd /cluster/data/anoCar2/bed/blastz.xenTro3.2011-09-19
    cat fb.anoCar2.chainXenTro3Link.txt
    #	86932478 bases of 1701353770 (5.110%) in intersection

    # and this swap
    mkdir /hive/data/genomes/xenTro3/bed/blastz.anoCar2.swap
    cd /hive/data/genomes/xenTro3/bed/blastz.anoCar2.swap

    time doBlastzChainNet.pl -verbose=2 \
	/cluster/data/anoCar2/bed/blastz.xenTro3.2011-09-19/DEF \
	-bigClusterHub=swarm -workhorse=hgwdev -smallClusterHub=memk \
	-swap -chainMinScore=5000 \
	-chainLinearGap=loose > swap.log 2>&1 &
    #	Elapsed time: 88m21s
    cat fb.xenTro3.chainAnoCar2Link.txt
    #	92033566 bases of 1358334882 (6.775%) in intersection

    cd /hive/data/genomes/xenTro3/bed
    ln -s blastz.anoCar2.swap lastz.anoCar2

##########################################################################
# Swap lastz hg19 Human (DONE - 2011-09-21 - Hiram)
    # original alignment
    cd /hive/data/genomes/hg19/bed/lastzXenTro3.2011-09-20
    cat fb.hg19.chainXenTro3Link.txt
    #	87928753 bases of 2897316137 (3.035%) in intersection

    #	running the swap - DONE - 2011-09-21
    mkdir /hive/data/genomes/xenTro3/bed/blastz.hg19.swap
    cd /hive/data/genomes/xenTro3/bed/blastz.hg19.swap
    time nice -n +19 doBlastzChainNet.pl -verbose=2 \
	/hive/data/genomes/hg19/bed/lastzXenTro3.2011-09-20/DEF \
	-chainMinScore=5000 -chainLinearGap=loose \
	-workhorse=hgwdev -smallClusterHub=memk -bigClusterHub=swarm \
	-swap > swap.log 2>&1 &
    #	real    37m7.001s
    cat fb.xenTro3.chainHg19Link.txt 
    #	90929066 bases of 1358334882 (6.694%) in intersection

##############################################################################
# Swap lastz danRer7 Zebrafish (DONE - 2011-09-21 - Hiram)
    # original alignment
    cd /hive/data/genomes/danRer7/bed/lastzXenTro3.2011-09-20
    cat fb.danRer7.chainXenTro3Link.txt 
    #	91771881 bases of 1409770109 (6.510%) in intersection
    cd /hive/data/genomes/danRer7/bed
    ln -s lastzXenTro3.2011-09-20 lastz.xenTro3

    #	running the swap
    mkdir /hive/data/genomes/xenTro3/bed/blastz.danRer7.swap
    cd /hive/data/genomes/xenTro3/bed/blastz.danRer7.swap
    time nice -n +19 doBlastzChainNet.pl -verbose=2 \
	/hive/data/genomes/danRer7/bed/lastzXenTro3.2011-09-20/DEF \
	-chainMinScore=5000 -chainLinearGap=loose \
	-workhorse=hgwdev -smallClusterHub=memk -bigClusterHub=swarm \
	-swap > swap.log 2>&1 &
    #	real    67m46.506s
    cat fb.xenTro3.chainDanRer7Link.txt 
    #	96539932 bases of 1358334882 (7.107%) in intersection

#########################################################################
# Swap lastz mm9 Mouse (DONE - 2011-09-21 - Hiram)
    # original alignment
    cd /hive/data/genomes/mm9/bed/lastzXenTro3.2011-09-20
    cat fb.mm9.chainXenTro3Link.txt 
    #	81920795 bases of 2620346127 (3.126%) in intersection

    #	running the swap - DONE - 2011-09-21
    mkdir /hive/data/genomes/xenTro3/bed/blastz.mm9.swap
    cd /hive/data/genomes/xenTro3/bed/blastz.mm9.swap
    time nice -n +19 doBlastzChainNet.pl -verbose=2 \
	/hive/data/genomes/mm9/bed/lastzXenTro3.2011-09-20/DEF \
	-chainMinScore=5000 -chainLinearGap=loose \
	-workhorse=hgwdev -smallClusterHub=memk -bigClusterHub=swarm \
	-swap > swap.log 2>&1 &
    #	real    57m36.529s
    cat fb.xenTro3.chainMm9Link.txt 
    #	89770014 bases of 1358334882 (6.609%) in intersection

##############################################################################
# Swap lastz melGal1 Turkey (DONE - 2011-09-23 - Hiram)
    # repeated a second time with better lastz parameters, see below
    # original alignment
    cd /hive/data/genomes/melGal1/bed/lastzXenTro3.2011-09-22
    cat fb.melGal1.chainXenTro3Link.txt 
    #	35993896 bases of 935922386 (3.846%) in intersection

    #   running the swap 
    mkdir /hive/data/genomes/xenTro3/bed/blastz.melGal1.swap
    cd /hive/data/genomes/xenTro3/bed/blastz.melGal1.swap
    time nice -n +19 doBlastzChainNet.pl -verbose=2 \
        /hive/data/genomes/melGal1/bed/lastzXenTro3.2011-09-22/DEF \
        -swap \
        -workhorse=hgwdev -smallClusterHub=memk -bigClusterHub=swarm \
        -chainMinScore=5000 -chainLinearGap=loose > swap.log 2>&1
    #	real    3m29.407s

    cat fb.xenTro3.chainMelGal1Link.txt 
    #	42384698 bases of 1358334882 (3.120%) in intersection

    cd /hive/data/genomes/xenTro3/bed
    ln -s blastz.melGal1.swap lastz.melGal1

#######################################################################
# Swap lastz monDom5 Opossum (DONE - 2011-09-23 - Hiram)
    # the original alignment
    cd /hive/data/genomes/monDom5/bed/lastzXenTro3.2011-09-22
    cat fb.monDom5.chainXenTro3Link.txt
    #	75079446 bases of 3501660299 (2.144%) in intersection

    # and the swap
    mkdir /hive/data/genomes/xenTro3/bed/blastz.monDom5.swap
    cd /hive/data/genomes/xenTro3/bed/blastz.monDom5.swap
    time nice -n +19 doBlastzChainNet.pl -verbose=2 \
        /hive/data/genomes/monDom5/bed/lastzXenTro3.2011-09-22/DEF \
        -swap -syntenicNet \
        -workhorse=hgwdev -smallClusterHub=memk -bigClusterHub=swarm \
        -chainMinScore=5000 -chainLinearGap=loose > swap.log 2>&1 &
    #   real    52m10.229s

    cat fb.xenTro3.chainMonDom5Link.txt
    #	74593612 bases of 1358334882 (5.492%) in intersection
    cd /hive/data/genomes/xenTro3/bed
    ln -s blastz.monDom5.swap lastz.monDom5

#########################################################################
# Swap lastz rn4 Rat (DONE - 2011-09-23 - Hiram)
    # original alignment
    cd /hive/data/genomes/rn4/bed/lastzXenTro3.2011-09-22
    cat fb.rn4.chainXenTro3Link.txt 
    #	66202042 bases of 2571531505 (2.574%) in intersection

    #	and the swap
    mkdir /hive/data/genomes/xenTro3/bed/blastz.rn4.swap
    cd /hive/data/genomes/xenTro3/bed/blastz.rn4.swap

    time nice -n +19 doBlastzChainNet.pl -verbose=2 \
        /hive/data/genomes/rn4/bed/lastzXenTro3.2011-09-22/DEF \
        -swap -syntenicNet \
        -workhorse=hgwdev -smallClusterHub=memk -bigClusterHub=swarm \
        -chainMinScore=5000 -chainLinearGap=loose > swap.log 2>&1 &
    #	Elapsed time: 20m33s

    cat fb.xenTro3.chainRn4Link.txt 
    #	74593612 bases of 1358334882 (5.492%) in intersection

    cd /hive/data/genomes/xenTro3/bed
    ln -s blastz.rn4.swap lastz.rn4

#############################################################################
# lastz Turkey melGal1 again (DONE - 2011-09-28 - Hiram)
    # original alignment
    cat fb.melGal1.chainXenTro3Link.txt
    #	49081729 bases of 935922386 (5.244%) in intersection

    #   running the swap
    mkdir /hive/data/genomes/xenTro3/bed/blastz.melGal1.swap
    cd /hive/data/genomes/xenTro3/bed/blastz.melGal1.swap
    time nice -n +19 doBlastzChainNet.pl -verbose=2 \
        /hive/data/genomes/melGal1/bed/lastzXenTro3.2011-09-27/DEF \
        -swap \
        -workhorse=hgwdev -smallClusterHub=memk -bigClusterHub=swarm \
        -chainMinScore=5000 -chainLinearGap=loose > swap.log 2>&1
    #	real    4m9.023s

    cat fb.xenTro3.chainMelGal1Link.txt
    #	57087653 bases of 1358334882 (4.203%) in intersection

    cd /hive/data/genomes/xenTro3/bed
    ln -s blastz.melGal1.swap lastz.melGal1

#######################################################################
## 9-Way Multiz (DONE - 2011-09-28 - Hiram)
    ssh hgwdev
    mkdir /hive/data/genomes/xenTro3/bed/multiz9way
    cd /hive/data/genomes/xenTro3/bed/multiz9way

    # working on: xenTro3 anoCar2 danRer7 galGal3 hg19 melGal1 mm9 monDom5 rn4

kent/src/hg/utils/phyloTrees/50way.nh
    xenTro3,xenTro2,galGal3,mm9,hg19,gasAcu1,ornAna1

    #	All distances remain as specified in the 46way.nh
    /cluster/bin/phast/tree_doctor --prune-all-but \
	anoCar1,xenTro2,danRer6,galGal3,melGal1,mm9,hg19,monDom5,rn4 \
	$HOME/kent/src/hg/utils/phyloTrees/50way.nh \
	| sed -e "s/anoCar1/anoCar2/; s/xenTro2/xenTro3/; s/danRer6/danRer7/" \
	> 9way.nh
    #	what that looks like:
    cat 9way.nh
# (((((hg19:0.144018,(mm9:0.084509,rn4:0.091589):0.271974):0.278985,
# monDom5:0.340786):0.181168,((melGal1:0.100000,galGal3:0.165536):0.299223,
# anoCar2:0.489241):0.105143):0.172371,xenTro3:0.855573):0.311354,
# danRer7:0.878701);

    #	rearrange to get xenTro3 on top:
    cat << '_EOF_' > xenTro3.9way.nh
((xenTro3:0.855573,(((hg19:0.144018,(mm9:0.084509,
rn4:0.091589):0.271974):0.278985,monDom5:0.340786):0.181168,
((melGal1:0.100000,galGal3:0.165536):0.299223,
anoCar2:0.489241):0.105143):0.172371):0.311354,danRer7:0.878701);
'_EOF_'
    # << happy emacs

    # convert to species names
    /cluster/bin/phast/tree_doctor --rename \
"xenTro3->X_tropicalis; hg19->Human; mm9->Mouse; rn4->Rat; monDom5->Opossum; melGal1->Turkey; galGal3->Chicken; anoCar2->Lizard; danRer7->Zebrafish" \
xenTro3.9way.nh > xenTro3.commonNames.9way.nh

# ((X_tropicalis:0.855573,(((Human:0.144018,
# (Mouse:0.084509,Rat:0.091589):0.271974):0.278985,
# Opossum:0.340786):0.181168,((Turkey:0.100000,
# Chicken:0.165536):0.299223,Lizard:0.489241):0.105143):0.172371):0.311354,
# Zebrafish:0.878701);

    #	Use this specification in the phyloGif tool:
    #	http://genome.ucsc.edu/cgi-bin/phyloGif
    #	to obtain a png image for src/hg/htdocs/images/phylo/xenTro3_9way.png

    /cluster/bin/phast/all_dists xenTro3.9way.nh > 9way.distances.txt
    #	Use this output to create the table below
    grep -i xenTro3 9way.distances.txt | sort -k3,3n
# xenTro3 melGal1 1.532310
# xenTro3 monDom5 1.549898
# xenTro3 galGal3 1.597846
# xenTro3 anoCar2 1.622328
# xenTro3 hg19    1.632115
# xenTro3 mm9     1.844580
# xenTro3 rn4     1.851660
# xenTro3 danRer7 2.045628

#	If you can fill in all the numbers in this table, you are ready for
#	the multiple alignment procedure

#                         featureBits chainLink measures
#                                   chainAnoCar2Link         chain linearGap
#    distance                    on xenTro3  on other     minScore
#  1  1.532 - turkey melGal1     (%  4.20) (%  5.24)       5000     loose
#  2  1.550 - opossum monDom5    (%  5.49) (%  2.14)       5000     loose
#  3  1.598 - chicken galGal3    (%  4.73) (%  5.27)       5000     loose
#  4  1.622 - lizard anoCar2     (%  6.78) (%  5.11)       5000     loose
#  5  1.632 - human hg19         (%  6.69) (%  3.04)       5000     loose
#  6  1.845 - mouse mm9          (%  6.61) (%  3.13)       5000     loose
#  7  1.852 - rat rn4            (%  5.49) (%  2.57)       5000     loose
#  8  2.046 - zebrafish danRer7  (%  7.11) (%  6.51)       5000     loose

# None of this concern for distances matters in building the first step, the
# maf files.

    # create species list and stripped down tree for autoMZ
    sed 's/[a-z][a-z]*_//g; s/:[0-9\.][0-9\.]*//g; s/;//; /^ *$/d' \
	xenTro3.9way.nh > tmp.nh
    echo `cat tmp.nh` > tree-commas.nh
    echo `cat tree-commas.nh` | sed 's/ //g; s/,/ /g' > tree.nh
    sed 's/[()]//g; s/,/ /g' tree.nh > species.list

    # split the maf files into a set of hashed named files
    # this hash named split keeps the same chr/contig names in the same
    # named hash file.
    mkdir /hive/data/genomes/xenTro3/bed/multiz9way/mafSplit
    cd /hive/data/genomes/xenTro3/bed/multiz9way/mafSplit
    for D in `sed -e "s/xenTro3 //" ../species.list`
do
    echo "${D}"
    mkdir $D
    cd $D
    echo "mafSplit -byTarget -useHashedName=10 /dev/null . ../../../lastz.${D}/axtChain/xenTro3.${D}.maf.gz"
    mafSplit -byTarget -useHashedName=8 /dev/null . \
	../../../lastz.${D}/mafNet/xenTro3.${D}.net.maf.gz
    cd ..
done

    # construct a list of all possible maf file names.
    # they do not all exist in each of the species directories
    find . -type f | wc -l
    # 2048
    find . -type f | grep ".maf$" | xargs -L 1 basename | sort -u > maf.list
    wc -l maf.list
    # 256 maf.list

    mkdir /hive/data/genomes/xenTro3/bed/multiz9way/splitRun
    cd /hive/data/genomes/xenTro3/bed/multiz9way/splitRun
    mkdir maf run
    cd run
    mkdir penn
    cp -p /cluster/bin/penn/multiz.2009-01-21/multiz penn
    cp -p /cluster/bin/penn/multiz.2009-01-21/maf_project penn
    cp -p /cluster/bin/penn/multiz.2009-01-21/autoMZ penn

    #	set the db and pairs directories here
    cat > autoMultiz.csh << '_EOF_'
#!/bin/csh -ef
set db = xenTro3
set c = $1
set result = $2
set run = `/bin/pwd`
set tmp = /scratch/tmp/$db/multiz.$c
set pairs = /hive/data/genomes/xenTro3/bed/multiz9way/mafSplit
/bin/rm -fr $tmp
/bin/mkdir -p $tmp
/bin/cp -p ../../tree.nh ../../species.list $tmp
pushd $tmp > /dev/null
foreach s (`/bin/sed -e "s/$db //" species.list`)
    set in = $pairs/$s/$c.maf
    set out = $db.$s.sing.maf
    if (-e $in.gz) then
        /bin/zcat $in.gz > $out
	if (! -s $out) then
	    echo "##maf version=1 scoring=autoMZ" > $out
	endif
    else if (-e $in) then
        /bin/ln -s $in $out
    else
        echo "##maf version=1 scoring=autoMZ" > $out
    endif
end
set path = ($run/penn $path); rehash
$run/penn/autoMZ + T=$tmp E=$db "`cat tree.nh`" $db.*.sing.maf $c.maf \
	> /dev/null
popd > /dev/null
/bin/rm -f $result
/bin/cp -p $tmp/$c.maf $result
/bin/rm -fr $tmp
/bin/rmdir --ignore-fail-on-non-empty /scratch/tmp/$db
'_EOF_'
# << happy emacs
    chmod +x autoMultiz.csh

    cat  << '_EOF_' > template
#LOOP
./autoMultiz.csh $(root1) {check out line+ /hive/data/genomes/xenTro3/bed/multiz9way/splitRun/maf/$(root1).maf}
#ENDLOOP
'_EOF_'
# << happy emacs

    ln -s ../../mafSplit/maf.list maf.list
    ssh swarm
    cd /hive/data/genomes/xenTro3/bed/multiz9way/splitRun/run
    # the tac reverses the list to get the small jobs first
    gensub2 maf.list single template stdout | tac > jobList
    para -ram=8g create jobsList
# Completed: 256 of 256 jobs
# CPU time in finished jobs:      14427s     240.45m     4.01h    0.17d  0.000 y
# IO & Wait Time:                  1927s      32.12m     0.54h    0.02d  0.000 y
# Average job time:                  64s       1.06m     0.02h    0.00d
# Longest finished job:             280s       4.67m     0.08h    0.00d
# Submission to last job:          1536s      25.60m     0.43h    0.02d

    # assemble into a single maf file
    cd /hive/data/genomes/xenTro3/bed/multiz9way
    head -1 splitRun/maf/001.maf > multiz9way.maf
    for F in splitRun/maf/*.maf
do
    egrep -v "^#" ${F}
done >> multiz9way.maf
    tail -1 splitRun/maf/001.maf >> multiz9way.maf

# -rw-rw-r--  1 1524521286 Sep 28 13:28 multiz9way.maf

    # Load into database
    ssh hgwdev
    cd /hive/data/genomes/xenTro3/bed/multiz9way
    mkdir /gbdb/xenTro3/multiz9way
    ln -s `pwd`/multiz9way.maf /gbdb/xenTro3/multiz9way
    cd /scratch/tmp
    time nice -n +19 hgLoadMaf xenTro3 multiz9way
    #	Indexing and tabulating /gbdb/xenTro3/multiz9way/multiz9way.maf
    #	Loaded 2591342 mafs in 1 files from /gbdb/xenTro3/multiz9way
    #	real    1m4.383s

    time nice -n +19 hgLoadMafSummary -verbose=2 -minSize=30000 \
	-mergeGap=1500 -maxSize=200000 xenTro3 multiz9waySummary \
	/gbdb/xenTro3/multiz9way/multiz9way.maf
    #	Created 516059 summary blocks from 5782695 components and 2591342 mafs
    #	from /gbdb/xenTro3/multiz9way/multiz9way.maf
    #	real    1m2.197s

    wc -l multiz9way*.tab
    #	2591342 multiz9way.tab
    #	516059 multiz9waySummary.tab
    #	3107401 total

    rm multiz9way*.tab

#######################################################################
# GAP ANNOTATE MULTIZ9WAY MAF AND LOAD TABLES (DONE - 2011-09-28 - Hiram)
    # mafAddIRows has to be run on single chromosome maf files, it does not
    #	function correctly when more than one reference sequence
    #	are in a single file.
    mkdir -p /hive/data/genomes/xenTro3/bed/multiz9way/anno/mafSplit
    cd /hive/data/genomes/xenTro3/bed/multiz9way/anno/mafSplit
    time mafSplit -byTarget -useFullSequenceName \
	/dev/null . ../../multiz9way.maf
    #	real    0m53.756s

    ls | wc
    #	9202    9202  119622

    cd /hive/data/genomes/xenTro3/bed/multiz9way/anno
    # a couple of these do not yet have an N.bed file
    for DB in xenTro3 melGal1 anoCar2
do
    cd /hive/data/genomes/${DB}
    twoBitInfo -nBed ${DB}.2bit ${DB}.N.bed
done

for DB in `cat ../species.list`
do
    echo "${DB} "
    ln -s  /hive/data/genomes/${DB}/${DB}.N.bed ${DB}.bed
    echo ${DB}.bed  >> nBeds
    ln -s  /hive/data/genomes/${DB}/chrom.sizes ${DB}.len
    echo ${DB}.len  >> sizes
done
    # make sure they all are successful symLinks:
    ls -ogrtL

    ssh swarm
    cd /hive/data/genomes/xenTro3/bed/multiz9way/anno
    mkdir result
    # NEXT TIME: this template should have a check out exists+ statement
    cat << '_EOF_' > template
#LOOP
mafAddIRows -nBeds=nBeds $(path1) /hive/data/genomes/xenTro3/xenTro3.2bit result/$(file1)
#ENDLOOP
'_EOF_'
    # << happy emacs

    ls mafSplit/*.maf > maf.list
    gensub2 maf.list single template jobList
    # limit jobs to one per node with the ram=8g requirement
    para -ram=8g create jobList
    para try ... check ... push ...
# Completed: 9202 of 9202 jobs
# CPU time in finished jobs:       4291s      71.52m     1.19h    0.05d  0.000 y
# IO & Wait Time:                 24581s     409.68m     6.83h    0.28d  0.001 y
# Average job time:                   3s       0.05m     0.00h    0.00d
# Longest finished job:              10s       0.17m     0.00h    0.00d
# Submission to last job:           146s       2.43m     0.04h    0.00d

    # verify all result files have some content, look for 0 size files:
    find . -type f -size 0
    # should see none

    # combine into one file
    head -q -n 1 result/GL172637.maf > xenTro3.9way.maf
    for F in result/*.maf
do
    grep -h -v "^#" ${F}
done >> xenTro3.9way.maf
    #	these maf files do not have the end marker, this does nothing:
    #	tail -q -n 1 result/GL172637.maf >> xenTro3.9way.maf
    # How about an official end marker:
    echo "##eof maf" >> xenTro3.9way.maf

    # Load into database
    rm /gbdb/xenTro3/multiz9way/multiz9way.maf   # remove old symlink
    ln -s `pwd`/xenTro3.9way.maf /gbdb/xenTro3/multiz9way/multiz9way.maf
    cd /scratch/tmp
    time nice -n +19 hgLoadMaf xenTro3 multiz9way
    #	Loaded 2929343 mafs in 1 files from /gbdb/xenTro3/multiz9way
    #	real    1m21.568s

    time nice -n +19 hgLoadMafSummary -verbose=2 -minSize=30000 \
	-mergeGap=1500 -maxSize=200000 xenTro3 multiz9waySummary \
	/gbdb/xenTro3/multiz9way/multiz9way.maf
    #	Created 516059 summary blocks from 5782695 components and 2929343 mafs
    #	from /gbdb/xenTro3/multiz9way/multiz9way.maf
    #	real    1m6.705s

    wc -l multiz9way*.tab
    #	2929343 multiz9way.tab
    #	516059 multiz9waySummary.tab
    #	3445402 total

    rm multiz9way*.tab

#########################################################################
# MULTIZ9WAY MAF FRAMES (DONE - 2011-05-16 - Hiram)
    ssh hgwdev
    mkdir /hive/data/genomes/xenTro3/bed/multiz9way/frames
    cd /hive/data/genomes/xenTro3/bed/multiz9way/frames
    mkdir genes

    #	we get more annotations using the ensGene on mm9 and hg19
    #	instead of knownGene, so, using ensGene for hg19 and mm9

xenTro3	ensGene
hg19	knownGene
mm9	knownGene
rn4	rgdGene2
monDom5	ensGene
melGal1	ensGene
galGal3	ensGene
anoCar2	ensGene
danRer7	ensGene

    #------------------------------------------------------------------------
    # get the genes for all genomes
    # mRNAs with CDS.  single select to get cds+psl, then split that up and
    # create genePred
    # using ensGene for:
for qDB in xenTro3 monDom5 melGal1 galGal3 anoCar2 danRer7
do
  echo hgsql -N -e \"'select * from 'ensGene\;\" ${qDB}
      hgsql -N -e "select * from ensGene" ${qDB} | cut -f 2-11 \
      | genePredSingleCover stdin stdout | gzip -2c > genes/$qDB.gp.gz
done

    # if you were using knownGene for mm9 hg19
    # genePreds; (must keep only the first 10 columns for knownGene)
    for qDB in mm9 hg19
    do
      echo hgsql -N -e \"'select * from 'knownGene\;\" ${qDB}
      hgsql -N -e "select * from knownGene" ${qDB} | cut -f 1-10 \
      | genePredSingleCover stdin stdout | gzip -2c > genes/$qDB.gp.gz
    done

      hgsql -N -e "select * from rgdGene2" rn4 | cut -f 1-10 \
      | genePredSingleCover stdin stdout | gzip -2c > genes/rn4.gp.gz

    # verify counts for genes are reasonable:
    for T in genes/*.gz
do
    echo -n "# $T: "
    zcat $T | cut -f1 | sort | uniq -c | wc -l
done

# genes/anoCar2.gp.gz: 17757
# genes/danRer7.gp.gz: 25874
# genes/galGal3.gp.gz: 16491
# genes/hg19.gp.gz: 20709
# genes/melGal1.gp.gz: 14040
# genes/mm9.gp.gz: 20905
# genes/monDom5.gp.gz: 19188
# genes/rn4.gp.gz: 18093
# genes/xenTro3.gp.gz: 18416

    ssh hgwdev
    cd /hive/data/genomes/xenTro3/bed/multiz9way/frames
    time (cat ../anno/xenTro3.9way.maf \
	| nice -n +19 genePredToMafFrames xenTro3 stdin stdout \
	    `sed -e "s#\([a-zA-Z0-9]*\)#\1 genes/\1.gp.gz#g" ../species.list` \
		| gzip > multiz9way.mafFrames.gz)
    #	real    2m22.546s

    # verify there are frames on everything:
    zcat multiz9way.mafFrames.gz | awk '{print $4}' | sort | uniq -c
    #	279231 anoCar2
    #	372399 danRer7
    #	258019 galGal3
    #	341339 hg19
    #	249007 melGal1
    #	322663 mm9
    #	257623 monDom5
    #	214301 rn4
    #	191614 xenTro3

    ssh hgwdev
    cd /hive/data/genomes/xenTro3/bed/multiz9way/frames
    time hgLoadMafFrames xenTro3 multiz9wayFrames multiz9way.mafFrames.gz
    #	real    0m21.425s

############################################################################
# creating upstream mafs (DONE - 2011-10-04 - Hiram)
    ssh hgwdev
    mkdir /hive/data/genomes/xenTro3/goldenPath/multiz9way
    cd /hive/data/genomes/xenTro3/goldenPath/multiz9way
    # run this bash script
    cat << '_EOF_' > mkUpstream.sh
DB=xenTro3
GENE=ensGene
NWAY=multiz9way
export DB GENE

for S in 1000 2000 5000
do
    echo "making upstream${S}.${GENE}.maf"
    echo "featureBits ${DB} ${GENE}:upstream:${S} -fa=/dev/null -bed=stdout"
    featureBits ${DB} ${GENE}:upstreamAll:${S} -fa=/dev/null -bed=stdout \
        | perl -wpe 's/_up[^\t]+/\t0/' | sort -k1,1 -k2,2n \
        | mafFrags ${DB} ${NWAY} stdin stdout \
                -orgs=/hive/data/genomes/${DB}/bed/${NWAY}/species.list \
        | gzip -c > upstream${S}.${GENE}.maf.gz
    echo "done upstream${S}.${GENE}.maf.gz"
done
'_EOF_'
    # << happy emacs
    chmod +x mkUpstream.sh
    time ./mkUpstream.sh
    #	real    17m4.037s

    # FIXUP/VERIFY the genbank.conf file to indicate this gene choice for the
    #	upstream automatic generation

#########################################################################
# Phylogenetic tree from 9-way (DONE - 2011-09-28 - Hiram)
    # this is a chancy operation with this business.  These alignments are
    #	too distant and not much of the gene sequence is matching.
    mkdir /hive/data/genomes/xenTro3/bed/multiz9way/4d
    cd /hive/data/genomes/xenTro3/bed/multiz9way/4d

    # the split annotated maf's are in:
    ../anno/result/*.maf

    cd /hive/data/genomes/xenTro3/bed/multiz9way/4d

    # using ensGene for xenTro3, only transcribed genes and nothing
    #	from the randoms and other misc.
    hgsql xenTro3 -Ne \
    "select * from ensGene WHERE cdsEnd > cdsStart;" | cut -f 2-20 \
	> ensGene.gp
    #	22705 ensGene.gp

    genePredSingleCover ensGene.gp stdout | sort > ensGeneNR.gp
    wc -l ensGeneNR.gp
    #	18416 ensGeneNR.gp

    ssh encodek
    mkdir /hive/data/genomes/xenTro3/bed/multiz9way/4d/run
    cd /hive/data/genomes/xenTro3/bed/multiz9way/4d/run
    mkdir ../mfa

    # newer versions of msa_view have a slightly different operation
    # the sed of the gp file inserts the reference species in the chr name
    cat << '_EOF_' > 4d.csh
#!/bin/csh -fe
set PHASTBIN = /cluster/bin/phast.build/cornellCVS/phast.2010-12-30/bin
set r = "/hive/data/genomes/xenTro3/bed/multiz9way"
set c = $1
set infile = $r/anno/result/$2
set outfile = $3
cd /scratch/tmp
# 'clean' maf
perl -wpe 's/^s ([^.]+)\.\S+/s $1/' $infile > $c.maf    
awk -v C=$c '$2 == C {print}' $r/4d/ensGeneNR.gp | sed -e "s/\t$c\t/\txenTro3.$c\t/" > $c.gp
set NL=`wc -l $c.gp| gawk '{print $1}'`
if ("$NL" != "0") then
    $PHASTBIN/msa_view --4d --features $c.gp -i MAF $c.maf -o SS > $c.ss
    $PHASTBIN/msa_view -i SS --tuple-size 1 $c.ss > $r/4d/run/$outfile
else
    echo "" > $r/4d/run/$outfile
endif
rm -f $c.gp $c.maf $c.ss
'_EOF_'
    # << happy emacs
    chmod +x 4d.csh

    ls -1S /hive/data/genomes/xenTro3/bed/multiz9way/anno/result/*.maf \
	| sed -e "s#.*multiz9way/anno/result/##" \
	> maf.list

    cat << '_EOF_' > template
#LOOP
4d.csh $(root1) $(path1) {check out line+ ../mfa/$(root1).mfa}
#ENDLOOP
'_EOF_'
    # << happy emacs

    # the tac puts the quick jobs at the front
    gensub2 maf.list single template stdout | tac > jobList
    para create jobList
    para try ... check
    para -maxJob=5 push
    para time
    # there are numerous problems in some of the result.  Not a whole lot
    #	of sequence is matching in these distant organisms
# Completed: 9143 of 9202 jobs
# Crashed: 30 jobs
# Other count: 29 jobs
# CPU time in finished jobs:        711s      11.85m     0.20h    0.01d  0.000 y
# IO & Wait Time:                 23273s     387.89m     6.46h    0.27d  0.001 y
# Average job time:                   3s       0.04m     0.00h    0.00d
# Longest finished job:               7s       0.12m     0.00h    0.00d
# Submission to last job:          2270s      37.83m     0.63h    0.03d

    # combine mfa files
    ssh hgwdev
    cd /hive/data/genomes/xenTro3/bed/multiz9way/4d
    # remove the broken empty files, size 0 and size 1:
    find ./mfa -type f -size 0 | xargs rm -f
    # most interesting, this did not identify files of size 1:
#    find ./mfa -type f -size 1
    ls -og mfa | awk '$3 == 1' > empty.list
    sed -e "s#^#mfa/##" empty.list | xargs rm -f
    #want comma-less species.list
    /cluster/bin/phast.build/cornellCVS/phast.2010-12-30/bin/msa_view \
	--aggregate "`cat ../species.list`" mfa/*.mfa | sed s/"> "/">"/ \
	    > 4d.all.mfa
    # check they are all in there:
    grep "^>" 4d.all.mfa
    #	>xenTro3
    #	>hg19
    #	>mm9
    #	>rn4
    #	>monDom5
    #	>melGal1
    #	>galGal3
    #	>anoCar2
    #	>danRer7

    # tree-commas.nh:
#((xenTro3,(((hg19,(mm9, rn4)),monDom5), ((melGal1,galGal3), anoCar2))),danRer7)

    # use phyloFit to create tree model (output is phyloFit.mod)
    time nice -n +19 \
	/cluster/bin/phast.build/cornellCVS/phast.2010-12-30/bin/phyloFit \
	    --EM --precision MED --msa-format FASTA --subst-mod REV \
		--tree ../tree-commas.nh 4d.all.mfa
    #	real    2m10.469s
    mv phyloFit.mod all.mod

    grep TREE all.mod
    # TREE: ((xenTro3:0.678022,(((hg19:0.160837,(mm9:0.0938355,rn4:0.110316):0.205024):0.213345,monDom5:0.346247):0.158444,((melGal1:0.0616202,galGal3:0.0595626):0.332209,anoCar2:0.500881):0.098368):0.212664):0.492003,danRer7:0.492003);

    # it actually looks pretty good as a tree estimate ...

#########################################################################
# phastCons 9-way (DONE - 2011-09-29 - Hiram)

    # split 9way mafs into 10M chunks and generate sufficient statistics 
    # files for # phastCons
    ssh encodek
    mkdir -p /hive/data/genomes/xenTro3/bed/multiz9way/cons/SS
    cd /hive/data/genomes/xenTro3/bed/multiz9way/cons/SS
    mkdir result done

    cat << '_EOF_' > mkSS.csh
#!/bin/csh -ef
set c = $1
set MAF = /hive/data/genomes/xenTro3/bed/multiz9way/anno/result/$c.maf
set WINDOWS = /hive/data/genomes/xenTro3/bed/multiz9way/cons/SS/result/$c
set WC = `cat $MAF | wc -l`
set NL = `grep "^#" $MAF | wc -l`
if ( -s $2 ) then
    exit 0
endif
if ( -s $2.running ) then
    exit 0
endif

date >> $2.running

rm -fr $WINDOWS
mkdir $WINDOWS
pushd $WINDOWS > /dev/null
if ( $WC != $NL ) then
/cluster/bin/phast.build/cornellCVS/phast.2010-12-30/bin/msa_split \
    $MAF -i MAF -o SS -r $WINDOWS/$c -w 10000000,0 -I 1000 -B 5000
endif
popd > /dev/null
date >> $2
rm -f $2.running
'_EOF_'
    # << happy emacs
    chmod +x mkSS.csh

    cat << '_EOF_' > template
#LOOP
mkSS.csh $(root1) {check out line+ done/$(root1)}
#ENDLOOP
'_EOF_'
    # << happy emacs

    #	do the easy ones first to see some immediate results
    ls -1S -r ../../anno/result | sed -e "s/.maf//" > maf.list

    gensub2 maf.list single template jobList
    para create jobList
    para try ... check ... etc
# Completed: 9202 of 9202 jobs
# CPU time in finished jobs:        816s      13.60m     0.23h    0.01d  0.000 y
# IO & Wait Time:                 24389s     406.49m     6.77h    0.28d  0.001 y
# Average job time:                   3s       0.05m     0.00h    0.00d
# Longest finished job:              10s       0.17m     0.00h    0.00d
# Submission to last job:          1276s      21.27m     0.35h    0.01d

    # not all of them produce results
    #	 there isn't much aligned in these MAF files
    find ./result -type f | wc
    #	 5690    5690  214806

    # Run phastCons
    #	This job is I/O intensive in its output files, beware where this
    #	takes place or do not run too many at once.
    ssh encodek
    mkdir -p /hive/data/genomes/xenTro3/bed/multiz9way/cons/run.cons
    cd /hive/data/genomes/xenTro3/bed/multiz9way/cons/run.cons

    #	there are going to be only one phastCons run using
    #	this same script.  It triggers off of the current working directory
    #	$cwd:t which is the "grp" in this script.  It is:
    #	all 

    cat << '_EOF_' > doPhast.csh
#!/bin/csh -fe
set PHASTBIN = /cluster/bin/phast.build/cornellCVS/phast.2010-12-30/bin
set c = $1
set cX = $1:r
set f = $2
set len = $3
set cov = $4
set rho = $5
set grp = $cwd:t
set cons = /hive/data/genomes/xenTro3/bed/multiz9way/cons
set tmp = $cons/tmp/$f
mkdir -p $tmp
set ssSrc = $cons
set useGrp = "$grp.mod"
ln -s $ssSrc/SS/result/$c/$f.ss $tmp
ln -s $cons/$grp/$grp.mod $tmp
pushd $tmp > /dev/null
$PHASTBIN/phastCons $f.ss $useGrp \
    --rho $rho --expected-length $len --target-coverage $cov --quiet \
    --seqname $c --idpref $c --most-conserved $f.bed --score > $f.pp
popd > /dev/null
mkdir -p pp/$c bed/$c
sleep 4
touch pp/$c bed/$c
rm -f pp/$c/$f.pp
rm -f bed/$c/$f.bed
mv $tmp/$f.pp pp/$c
mv $tmp/$f.bed bed/$c
rm -fr $tmp
'_EOF_'
    # << happy emacs
    chmod +x doPhast.csh

    #	this template will serve for all runs
    #	root1 == chrom name, file1 == ss file name without .ss suffix
    cat << '_EOF_' > template
#LOOP
../run.cons/doPhast.csh $(root1) $(file1) 45 0.3 0.3 {check out line+ pp/$(root1)/$(file1).pp}
#ENDLOOP
'_EOF_'
    # << happy emacs

    find ../SS -type f | grep ".ss$" | sed -e 's/.ss$//' > ss.list
    wc -l ss.list
    #	3291 ss.list

    # Create parasol batch and run it
    # run for all species
    cd /hive/data/genomes/xenTro3/bed/multiz9way/cons
    mkdir -p all
    cd all
    #	Using the .mod tree
    cp -p ../../4d/all.mod ./all.mod

    gensub2 ../run.cons/ss.list single ../run.cons/template jobList
    para create jobList
    para try ... check ...
    para push
# Completed: 5690 of 5690 jobs
# CPU time in finished jobs:       2411s      40.18m     0.67h    0.03d  0.000 y
# IO & Wait Time:                 37365s     622.76m    10.38h    0.43d  0.001 y
# Average job time:                   7s       0.12m     0.00h    0.00d
# Longest finished job:              20s       0.33m     0.01h    0.00d
# Submission to last job:          2030s      33.83m     0.56h    0.02d

    cd /hive/data/genomes/xenTro3/bed/multiz9way/cons/all

    find ./bed -type f | grep ".bed$" | xargs cat | sort -k1,1 -k2,2n \
	| awk '{printf "%s\t%d\t%d\tlod=%d\t%s\n", $1, $2, $3, $5, $5;}' \
	    > tmpMostConserved.bed

    /cluster/bin/scripts/lodToBedScore tmpMostConserved.bed > mostConserved.bed

    # load into database
    time nice -n +19 hgLoadBed xenTro3 phastConsElements9way mostConserved.bed
    #	Loaded 827762 elements of size 5
    #	real    0m4.880s

    # on human we often try for 5% overall cov, and 70% CDS cov 
    # most bets are off here for that goal, these alignments are too few
    #	and too far between
    featureBits xenTro3 -enrichment ensGene:cds phastConsElements9way
    #	--rho 0.3 --expected-length 45 --target-coverage 0.3
# ensGene:cds 2.244%, phastConsElements9way 10.055%, both 1.669%, cover 74.39%, enrich 7.40x
    featureBits xenTro3 -enrichment mgcGenes:cds phastConsElements9way
# mgcGenes:cds 0.741%, phastConsElements9way 10.055%, both 0.619%, cover 83.53%, enrich 8.31x

    # Create merged posterier probability file and wiggle track data files
    cd /hive/data/genomes/xenTro3/bed/multiz9way/cons/all
    mkdir downloads

    find ./pp -type f | sed -e "s#^./##; s#\.# d #g; s#-# m #;" \
	| sort -k1,1 -k3,3n | sed -e "s# d #.#g; s# m #-#g;" | xargs cat \
	| gzip -c > downloads/phastCons9way.wigFix.gz

    # check integrity of data with wigToBigWig
    time (zcat downloads/phastCons9way.wigFix.gz \
	| wigToBigWig -verbose=2 stdin /hive/data/genomes/xenTro3/chrom.sizes \
	    phastCons9way.bw) > bigWig.log 2>&1 &
    tail bigWig.log
    #	pid=32277: VmPeak:     2219472 kB
    #	real    2m45.147s

    bigWigInfo phastCons9way.bw
# version: 4
# isCompressed: yes
# isSwapped: 0
# primaryDataSize: 489,270,288
# primaryIndexSize: 14,137,652
# zoomLevels: 10
# chromCount: 5690
# basesCovered: 202,032,632
# mean: 0.630799
# min: 0.000000
# max: 1.000000
# std: 0.378048

    #	encode those files into wiggle data
    time (zcat downloads/phastCons9way.wigFix.gz \
	| wigEncode stdin phastCons9way.wig phastCons9way.wib)
    #	Converted stdin, upper limit 1.00, lower limit 0.00
    #	real    1m2.832s

    du -hsc *.wi?
    #	193M    phastCons9way.wib
    #	43M     phastCons9way.wig
    #	236M    total

    # Load gbdb and database with wiggle.
    ln -s `pwd`/phastCons9way.wib /gbdb/xenTro3/multiz9way/phastCons9way.wib
    time nice -n +19 hgLoadWiggle -pathPrefix=/gbdb/xenTro3/multiz9way \
	xenTro3 phastCons9way phastCons9way.wig
    #	real    0m4.207s

    # use to set trackDb.ra entries for wiggle min and max
    # and verify table is loaded correctly

    wigTableStats.sh xenTro3 phastCons9way
# db.table      min max mean count sumData stdDev viewLimits
#xenTro3.phastCons9way    0 1 0.630799 202032632 1.27442e+08 0.378048 viewLimits=0:1

    #  Create histogram to get an overview of all the data
    time nice -n +19 hgWiggle -doHistogram -db=xenTro3 \
	-hBinSize=0.001 -hBinCount=1000 -hMinVal=0.0 -verbose=2 \
	    phastCons9way > histogram.data 2>&1
    #	real    0m42.544s

    #	create plot of histogram:

    cat << '_EOF_' | gnuplot > histo.png
set terminal png small x000000 xffffff xc000ff x66ff66 xffff00 x00ffff
set size 1.4, 0.8
set key left box
set grid noxtics
set grid ytics
set title " X. tropicalis xenTro3 Histogram phastCons9way track"
set xlabel " phastCons9way score"
set ylabel " Relative Frequency"
set y2label " Cumulative Relative Frequency (CRF)"
set y2range [0:1]
set y2tics
set yrange [0:0.02]

plot "histogram.data" using 2:5 title " RelFreq" with impulses, \
        "histogram.data" using 2:7 axes x1y2 title " CRF" with lines
'_EOF_'
    #	<< happy emacs

    display histo.png &

#############################################################################
# phyloP for 7-way (DONE - 2011-09-29 - Hiram)
    # run phyloP with score=LRT 
    ssh encodek
    mkdir /cluster/data/xenTro3/bed/multiz9way/consPhyloP
    cd /cluster/data/xenTro3/bed/multiz9way/consPhyloP

    mkdir run.phyloP
    cd run.phyloP
    # Adjust model file base composition background and rate matrix to be
    # representative of the chromosomes in play
    grep BACKGROUND ../../cons/all/all.mod | awk '{printf "%0.3f\n", $3 + $4}'
    #	0.495
    /cluster/bin/phast.build/cornellCVS/phast.2010-12-30/bin/modFreqs \
	../../cons/all/all.mod 0.495 > all.mod
    # verify, the BACKGROUND should now be paired up:
    grep BACK all.mod 
    #	BACKGROUND: 0.252500 0.247500 0.247500 0.252500 

    #	following the pattern from panTro3 with grp: "all" (no other grp)
    cat << '_EOF_' > doPhyloP.csh
#!/bin/csh -fe
set PHASTBIN = /cluster/bin/phast.build/cornellCVS/phast.2010-12-30/bin
set f = $1
set file1 = $f:t
set out = $2
set cName = $f:t:r
set n = $f:r:e
set grp = $cwd:t
set cons = /hive/data/genomes/xenTro3/bed/multiz9way/consPhyloP
set tmp = $cons/tmp/$grp/$f
rm -fr $tmp
mkdir -p $tmp
set ssSrc = "/hive/data/genomes/xenTro3/bed/multiz9way/cons/SS/result/$f"
set useGrp = "$grp.mod"
ln -s $cons/run.phyloP/$grp.mod $tmp
pushd $tmp > /dev/null
$PHASTBIN/phyloP --method LRT --mode CONACC --wig-scores --chrom $cName \
    -i SS $useGrp $ssSrc.ss > $file1.wigFix
popd > /dev/null
mkdir -p $out:h
sleep 4
mv $tmp/$file1.wigFix $out
rm -fr $tmp
rmdir --ignore-fail-on-non-empty $cons/tmp/$grp/$cName
rmdir --ignore-fail-on-non-empty $cons/tmp/$grp
rmdir --ignore-fail-on-non-empty $cons/tmp
'_EOF_'
    # << happy emacs

    # Create list of chunks
    find ../../cons/SS/result -type f | grep ".ss$" \
	| sed -e "s/.ss$//; s#^../../cons/SS/result/##" > ss.list
    # make sure the list looks good
    wc -l ss.list
    #	5690 ss.list

    # Create template file
    #	file1 == $chr/$chunk/file name without .ss suffix
    cat << '_EOF_' > template
#LOOP
../run.phyloP/doPhyloP.csh $(path1) {check out line+ wigFix/$(dir1)/$(file1).wigFix}
#ENDLOOP
'_EOF_'
    # << happy emacs

    ######################   Running all species  #######################
    # setup run for all species
    mkdir /hive/data/genomes/xenTro3/bed/multiz9way/consPhyloP/all
    cd /hive/data/genomes/xenTro3/bed/multiz9way/consPhyloP/all
    rm -fr wigFix
    mkdir wigFix

    gensub2 ../run.phyloP/ss.list single ../run.phyloP/template jobList
    # the -ram=8g will allow only one job per node to slow this down since
    #	it would run too fast otherwise.  Either run on one of the small
    #	klusters or use the -ram=8g on the para create
    para -ram=8g create jobList
    para try ... check ... etc ...
    para -maxJob=64 push
    para time > run.time
# Completed: 5690 of 5690 jobs
# CPU time in finished jobs:       9670s     161.17m     2.69h    0.11d  0.000 y
# IO & Wait Time:                 38146s     635.76m    10.60h    0.44d  0.001 y
# Average job time:                   8s       0.14m     0.00h    0.00d
# Longest finished job:              63s       1.05m     0.02h    0.00d
# Submission to last job:          2448s      40.80m     0.68h    0.03d

    # make downloads
    mkdir downloads

    time (find ./wigFix -type f | sed -e "s#^./##; s#\.# d #g; s#-# m #;" \
	| sort -k1,1 -k3,3n | sed -e "s# d #.#g; s# m #-#g;" | xargs cat \
	| gzip -c > downloads/phyloP9way.wigFix.gz) &
    #	real    5m28.973s

    # check integrity of data with wigToBigWig
    time (zcat downloads/phyloP9way.wigFix.gz \
	| wigToBigWig -verbose=2 stdin /hive/data/genomes/xenTro3/chrom.sizes \
	phyloP9way.bw) > bigWig.log 2>&1 &
    egrep "real|VmPeak" bigWig.log
    # pid=17458: VmPeak:     2219472 kB
    #	real    2m53.180s

    bigWigInfo phyloP9way.bw 
# version: 4
# isCompressed: yes
# isSwapped: 0
# primaryDataSize: 366,163,491
# primaryIndexSize: 14,137,652
# zoomLevels: 10
# chromCount: 5690
# basesCovered: 202,032,632
# mean: 0.447745
# min: -3.316000
# max: 2.298000
# std: 0.772169

    #	encode those files into wiggle data
    time (zcat downloads/phyloP9way.wigFix.gz \
	| wigEncode stdin phyloP9way.wig phyloP9way.wib) &
    #	Converted stdin, upper limit 2.30, lower limit -3.32
    #	real    1m3.694s

    du -hsc *.wi?
    #	193M    phyloP9way.wib
    #	43M     phyloP9way.wig
    #	236M    total

    # Load gbdb and database with wiggle.
    ln -s `pwd`/phyloP9way.wib /gbdb/xenTro3/multiz9way/phyloP9way.wib
    nice hgLoadWiggle -pathPrefix=/gbdb/xenTro3/multiz9way xenTro3 \
	phyloP9way phyloP9way.wig

    # use to set trackDb.ra entries for wiggle min and max
    # and verify table is loaded correctly

    wigTableStats.sh xenTro3 phyloP9way
# db.table      min max mean count sumData stdDev viewLimits
# xenTro3.phyloP9way      -3.316 2.298 0.447745 202032632 9.04591e+07 0.772169
# viewLimits=-3.316:2.298

    #  Create histogram to get an overview of all the data
    time nice -n +19 hgWiggle -doHistogram -db=xenTro3 \
	-hBinSize=0.001 -hBinCount=1000 -hMinVal=0.0 -verbose=2 \
	    phyloP9way > histogram.data 2>&1
    #	real    0m19.396s

    # find out the range for the 2:5 graph
    grep -v chrom histogram.data | grep "^[0-9]" | ave -col=5 stdin
# Q1 0.000240
# median 0.000383
# Q3 0.000726
# average 0.001000
# min 0.000107
# max 0.041706
# count 1000
# total 1.000011
# standard deviation 0.002765

    #	create plot of histogram:

    cat << '_EOF_' | gnuplot > histo.png
set terminal png small x000000 xffffff xc000ff x66ff66 xffff00 x00ffff
set size 1.4, 0.8
set key left box
set grid noxtics
set grid ytics
set title " X. tropicalis xenTro3 Histogram phyloP9way track"
set xlabel " phyloP9way score"
set ylabel " Relative Frequency"
set y2label " Cumulative Relative Frequency (CRF)"
set y2range [0:1]
set y2tics
set yrange [0:0.007]

plot "histogram.data" using 2:5 title " RelFreq" with impulses, \
        "histogram.data" using 2:7 axes x1y2 title " CRF" with lines
'_EOF_'
    #	<< happy emacs

    display histo.png &

#############################################################################
# download data for 7-way (DONE - 2011-09-29 - Hiram)
    mkdir /usr/local/apache/htdocs-hgdownload/goldenPath/xenTro3/phastCons9way
    cd /usr/local/apache/htdocs-hgdownload/goldenPath/xenTro3/phastCons9way
    ln -s \
/hive/data/genomes/xenTro3/bed/multiz9way/cons/all/downloads/phastCons9way.wigFix.gz \
	./phastCons9way.wigFix.gz
    ln -s /hive/data/genomes/xenTro3/bed/multiz9way/cons/all/phastCons9way.bw \
	./phastCons9way.bw

    ln -s /hive/data/genomes/xenTro3/bed/multiz9way/cons/all/all.mod \
	./vertebrate.mod
    # use a README from a recent multiz like this, this one is from anoCar2
    cd /hive/data/genomes/xenTro3/bed/multiz9way/cons
    cp -p /hive/data/genomes/anoCar2/bed/multiz7way/cons/README.txt .
    # need to edit this to finish it off

    ln -s /hive/data/genomes/xenTro3/bed/multiz9way/cons/README.txt .
    md5sum * > /hive/data/genomes/xenTro3/bed/multiz9way/cons/md5sum.txt
    ln -s /hive/data/genomes/xenTro3/bed/multiz9way/cons/md5sum.txt .

    mkdir /usr/local/apache/htdocs-hgdownload/goldenPath/xenTro3/phyloP9way
    cd /usr/local/apache/htdocs-hgdownload/goldenPath/xenTro3/phyloP9way
    ln -s \
/hive/data/genomes/xenTro3/bed/multiz9way/consPhyloP/all/downloads/phyloP9way.wigFix.gz \
	./phyloP9way.wigFix.gz
    ln -s \
/hive/data/genomes/xenTro3/bed/multiz9way/consPhyloP/all/phyloP9way.bw \
	./phyloP9way.bw
    ln -s \
/hive/data/genomes/xenTro3/bed/multiz9way/consPhyloP/run.phyloP/all.mod \
	./vertebrate.mod
    # use a README from a recent multiz like this, this one is from anoCar2
    /hive/data/genomes/xenTro3/bed/multiz9way/consPhyloP
    cp /hive/data/genomes/anoCar2/bed/multiz7way/consPhyloP/README.txt .
    # need to edit this to finish it off

    ln -s /hive/data/genomes/xenTro3/bed/multiz9way/consPhyloP/README.txt .
    md5sum * > /hive/data/genomes/xenTro3/bed/multiz9way/consPhyloP/md5sum.txt
    ln -s /hive/data/genomes/xenTro3/bed/multiz9way/consPhyloP/md5sum.txt .

    mkdir /hive/data/genomes/xenTro3/bed/multiz9way/downloads
    cd /hive/data/genomes/xenTro3/bed/multiz9way/downloads
    # this was already done elsewhere:
    grep TREE ../cons/all/all.mod | sed -e "s/TREE: //" > 9way.nh
    /cluster/bin/phast/tree_doctor --rename \
"xenTro3->X_tropicalis; hg19->Human; mm9->Mouse; rn4->Rat; monDom5->Opossum; melGal1->Turkey; galGal3->Chicken; anoCar2->Lizard; danRer7->Zebrafish" \
9way.nh > xenTro3.commonNames.9way.nh

    time cp -p ../anno/xenTro3.9way.maf ./multiz9way.maf
    #	real    1m13.654s
    time gzip multiz9way.maf
    #	real    7m12.108s

    mkdir /usr/local/apache/htdocs-hgdownload/goldenPath/xenTro3/multiz9way
    cd /usr/local/apache/htdocs-hgdownload/goldenPath/xenTro3/multiz9way
    ln -s /hive/data/genomes/xenTro3/bed/multiz9way/downloads/9way.nh .
    ln -s \
/hive/data/genomes/xenTro3/bed/multiz9way/downloads/xenTro3.commonNames.9way.nh .
    ln -s \
/hive/data/genomes/xenTro3/bed/multiz9way/downloads/multiz9way.maf.gz .

    # use a README from a recent multiz like this, this one is from danRer7
    cd /hive/data/genomes/xenTro3/bed/multiz9way/downloads
    cp -p /hive/data/genomes/anoCar2/bed/multiz7way/downloads/README.txt .
    # need to edit this to finish it off
    cd /usr/local/apache/htdocs-hgdownload/goldenPath/xenTro3/multiz9way
    ln -s /hive/data/genomes/xenTro3/bed/multiz9way/downloads/README.txt .

    md5sum *.nh *.gz *.txt \
	> /hive/data/genomes/xenTro3/bed/multiz9way/downloads/md5sum.txt
    ln -s /hive/data/genomes/xenTro3/bed/multiz9way/downloads/md5sum.txt .

#############################################################################
# hgPal downloads (DONE braney 2011-09-30)
#   FASTA from 9way for ensGene

    ssh hgwdev
    screen
    bash
    rm -rf /cluster/data/xenTro3/bed/multiz9way/pal
    mkdir /cluster/data/xenTro3/bed/multiz9way/pal
    cd /cluster/data/xenTro3/bed/multiz9way/pal
    for i in `cat ../species.list`; do echo $i; done > order.lst

    mz=multiz9way
    gp=ensGene
    db=xenTro3
    mkdir exonAA exonNuc ppredAA ppredNuc
    for j in `sort -nk 2 /cluster/data/$db/chrom.sizes | awk '{print $1}'`
    do
	echo "date"
	echo "mafGene -chrom=$j  $db $mz $gp order.lst stdout | \
	    gzip -c > ppredAA/$j.ppredAA.fa.gz"
	echo "mafGene -chrom=$j -noTrans $db $mz $gp order.lst stdout | \
	    gzip -c > ppredNuc/$j.ppredNuc.fa.gz"
	echo "mafGene -chrom=$j -exons -noTrans $db $mz $gp order.lst stdout | \
	    gzip -c > exonNuc/$j.exonNuc.fa.gz"
	echo "mafGene -chrom=$j -exons $db $mz $gp order.lst stdout | \
	    gzip -c > exonAA/$j.exonAA.fa.gz"
    done > $gp.jobs

    time sh -x $gp.jobs > $gp.jobs.log 2>&1 &
    sleep 1
    tail -f $gp.jobs.log

# real    29m58.498s
# user    8m8.312s
# sys     6m7.591s

    mz=multiz9way
    gp=ensGene
    db=xenTro3
    zcat exonAA/*.gz | gzip -c > $gp.$mz.exonAA.fa.gz
    zcat exonNuc/*.gz | gzip -c > $gp.$mz.exonNuc.fa.gz
    zcat ppredAA/*.gz | gzip -c > $gp.$mz.ppredAA.fa.gz
    zcat ppredNuc/*.gz | gzip -c > $gp.$mz.ppredNuc.fa.gz

    rm -rf exonAA exonNuc ppredAA ppredNuc

    # we're only distributing exons at the moment
    mz=multiz9way
    gp=ensGene
    db=xenTro3
    pd=/usr/local/apache/htdocs-hgdownload/goldenPath/$db/$mz/alignments
    mkdir -p $pd
    ln -s `pwd`/$gp.$mz.exonAA.fa.gz $pd/$gp.exonAA.fa.gz
    ln -s `pwd`/$gp.$mz.exonNuc.fa.gz $pd/$gp.exonNuc.fa.gz

