src/hg/makeDb/doc/panTro2.txt 1.16

1.16 2009/08/21 06:25:21 markd
re-annotate/reload unrelease maf to fix bug
Index: src/hg/makeDb/doc/panTro2.txt
===================================================================
RCS file: /projects/compbio/cvsroot/kent/src/hg/makeDb/doc/panTro2.txt,v
retrieving revision 1.15
retrieving revision 1.16
diff -b -B -U 1000000 -r1.15 -r1.16
--- src/hg/makeDb/doc/panTro2.txt	21 Jul 2009 21:01:45 -0000	1.15
+++ src/hg/makeDb/doc/panTro2.txt	21 Aug 2009 06:25:21 -0000	1.16
@@ -1,3226 +1,3251 @@
 # for emacs: -*- mode: sh; -*-
 
 # This file describes browser build for the panTro2 chimp genome: Feb 2006
 #
 #	"$Id$"
 #
 
 #  NOTE:  this doc may have genePred loads that fail to include
 #  the bin column.  Please correct that for the next build by adding
 #  a bin column when you make any of these tables:
 
 #  mysql> SELECT tableName, type FROM trackDb WHERE type LIKE "%Pred%";
 #  +-------------+-------------------------+
 #  | tableName   | type                    |
 #  +-------------+-------------------------+
 #  | refGene     | genePred refPep refMrna |
 #  | xenoRefGene | genePred                |
 #  | nscanGene   | genePred nscanPep       |
 #  | genscan     | genePred genscanPep     |
 #  +-------------+-------------------------+
 
 
 #######################################################################
 # DOWNLOAD 6X CHIMP SEQUENCE FROM LaDeana HIllier, WUSTL (2006-02-27 - kate)
 #
 # panTro2
 # 6X Chimp from LaDeana Hillier, WUSLT 
 
 # NCBI accession for the project is:  AADA01000000  (02-SEP-2005)
 # NOTE: chrom numbering follows hunan orthology, a change from panTro1.
 # Human chrom 2 is orthologus to chimp chroms 2a and 2b.
 
 # SETUP BUILD AREA AND DOWNLOAD ASSEMBLY (DONE 2006-02-27 kate)
 
     ssh kkstore01
     mkdir /cluster/store8/panTro2
     ln -s /cluster/store8/panTro2 /cluster/data
     cd /cluster/data/panTro2
 
     # annotations dir
     mkdir bed
 
     # temp dir for lift files and 
     # scripts copy-pasted from this doc
     mkdir jkStuff
 
     # make download dir
     wget ftp://genome.wustl.edu/pub/user/lhillier/panTro2.tar.gz
 
     tar xvfz *.gz
     mv panTro2 wustl
 
     # for finished chr21, there is chr21.agp, chr21.clones.agp,
     #   and chr21.fa (complete chrom fasta file)
     # for partially finished chrY, there is:
     #   .agp, .clones.agp, and .fa for chrY and chrY_random
     # for chr7, that includes a 5M finished region, there
     # are:  chr7.agp, chr7.contigs.fa, and chr7.clones.fa
     # (also chr7_random.{agp,contigs.fa}
     # For all other chroms there are AGP files plus contig fasta:
     #   chr*.agp, chr*.contigs.fa files
 
     grep '>' wustl/*.contigs.fa | wc -l
     # 246371
     # NOTE: Release notes list 265882 contigs 
     # (due to replacement of contigs by clone sequence in finished chroms?)
     grep '>' wustl/*.clones.fa | wc -l
     # 505
 
 
 # BUILD CHROM FASTA FILES (2006-02-28 kate)
 #  Most chroms are built from contigs, a few have
 #       clones in the AGP, chr7 has both!
 
     cut -f1 wustl/*.agp | uniq | grep -v '^$' | grep -v random | \
         sed 's/chr//' > chrom.lst
     wc -l chrom.lst
     # 28
     # 1,2a,2b,3-22,X,Y,M,Un, and chr6_hla_hap1
 
 cat > jkStuff/makeChroms.csh << 'EOF'
     date
     foreach d (`cat chrom.lst`)
         set c = chr$d
         echo $c
         mkdir $d
         cp wustl/$c.{agp,*.fa} $d
         if (-e $d/$c.clones.fa) then
             set fa = clones.fa
             # chr7 has both contigs and clones and clone file
             # doesn't have genbank sequence header -- merge them
             if ($c == "chr7") then
                 cat $d/$c.contigs.fa  >> $d/$c.clones.fa
                 rm $d/$c.contigs.fa
             else
                 # extract clone id from sequence name line to please agpToFa
                 awk -F\| '{if (/^>/) printf(">%s\n", $4); else print}' \
                 wustl/$c.clones.fa > $d/$c.clones.fa
             endif
         else
             set fa = contigs.fa
         endif
         agpToFa -simpleMultiMixed $d/$c.agp $c $d/$c.fa $d/$c.$fa
         if (-e wustl/${c}_random.agp) then
             set r = ${c}_random
             echo $r
             cp wustl/$r.{agp,*.fa} $d
             if (-e $d/$r.clones.fa) then
                 set fa = clones.fa
             awk -F\| '{if (/^>/) printf(">%s\n", $4); else print}' \
                 wustl/$r.clones.fa > $d/$r.clones.fa
             else
                 set fa = contigs.fa
             endif
             agpToFa -simpleMultiMixed $d/$r.agp $r $d/$r.fa $d/$r.$fa
         endif
     end
     date
 'EOF'
     # << happy emacs
     csh jkStuff/makeChroms.csh >&! jkStuff/makeChroms.log &
         # 5 minutes
 
 cat > jkStuff/checkChroms.csh << 'EOF'
     date
     foreach d (`cat chrom.lst`)
         set c = chr$d
         echo $c
         faSize $d/$c.fa
         checkAgpAndFa $d/$c.agp $d/$c.fa | tail -1
         if (-e $d/${c}_random.fa) then
             set r = ${c}_random
             echo $r
             faSize $d/$r.fa
             checkAgpAndFa $d/$r.agp $d/$r.fa | tail -1
         endif
     end
     date
 'EOF'
     # << happy emacs
     csh jkStuff/checkChroms.csh >&! jkStuff/checkChroms.log &
         # 3 minutes
 
     # Replace chrM and chrM_random from this assembly with
     # chrM from NCBI.
     # so retrieving it from NCBI, via method described in canFam2 make doc.
     # (search "pan troglodytes mitochondrion genome" finds NC_001643
     rm -fr M
     mkdir M
     cd M
     wget -O chrM.fa 'http://www.ncbi.nlm.nih.gov/entrez/viewer.fcgi?db=nucleotide&val=5835121&dopt=FASTA'
     # Edit chrM.fa: make sure the long fancy header line says it's the 
     # Canis familiaris mitochondrion complete genome, and then replace the 
     # header line with just ">chrM".
     faSize chrM.fa
         # 16554 bases
 
     # cleanup chrom dirs
     cd ..
     foreach d (`cat chrom.lst`)
         set c = chr$d
         mkdir $d/parts
         mv $d/*.clones.fa $d/*.contigs.fa $d/parts
     end
 
     # retrieve fixed chr9.agp from LaDeana (2006-03-02)
     mkdir wustl.chr9
     cd wustl.chr9
     wget ftp://genome.wustl.edu/pub/user/lhillier/pub/chr9.agp.gz
     gunzip *.gz
     cp -p chr9.agp ../9 
     mv wustl/chr9.agp wustl/chr9.agp.old
     cp -p chr9.agp wustl/chr9.agp
     cd ../9
     agpToFa -simpleMultiMixed chr9.agp chr9 chr9.fa parts/chr9.contigs.fa
     checkAgpAndFa chr9.agp chr9.fa | tail -1
 
 
 ########################
 # QUALITY SCORES  (2006-03-15 kate)
 # From Joanne Nelson, WUSTL
 # Redo chr7 with finished clone scores newly available (2006-07-13 kate)
 # Finished chr21 scores are being submitted to Genbank by Riken,
 # but are not yet available
 
     ssh kkstore04
     cd /cluster/data/panTro2
     mkdir bed/quality
     cd bed/quality
     wget -r ftp://genome.wustl.edu/private/963082470306293/upload_for_kate
     mkdir wustl
     mv genome.wustl.edu/private/*/upload*/*.qvl wustl
     rm -fr genome.wustl.edu
     grep '>' wustl/chr*.qvl | wc -l
         # 246371
         # same as #contigs in sequence files, above
     gzip wustl/*.qvl
 
 cat > makeQuals.csh << 'EOF'
     date
     mkdir -p qac
     set b = /cluster/data/panTro2
     foreach f (`ls wustl/*.qvl.gz | grep -v chrM | grep -v ContigF`)
         set c = $f:t:r:r:r
         echo $c
         set d = `echo $c | sed 's/chr//;s/_random//'`
         # contig names in quality files are slightly different from AGP's
         sed 's/Contig/bld2_Cont/' $b/$d/$c.agp > agp.tmp
         nice zcat wustl/$c.contigs.qvl | \
             nice qaToQac stdin stdout | \
             nice qacAgpLift agp.tmp stdin qac/$c.qac
     end
     date
 'EOF'
     # << happy emacs
     csh makeQuals.csh >&! makeQuals.log &
 #chr7
 #AC161123.3 not found
     # need special handling for chr7, that has both contigs
     # (with quality scores), and clones (w/o)
 
     # For chimpHiQualDiffs, construct a chr7 in one piece, with
     # missing data filled in with score=98 (manually assigned "low quality"),
     # although these are clones, and should have high quality.
     # By request of Daryl Thomas
 
     sed 's/Contig/bld2_Cont/' /cluster/data/panTro2/7/chr7.agp > agp.tmp
     nice zcat wustl/chr7.contigs.qvl.gz | \
         nice qaToQac stdin stdout | \
         nice qacAgpLift -mScore=98 agp.tmp stdin qac/chr7.qac
 
     # Redo chr7 with finished clone scores newly available (2006-07-13)
     wget -nd -r ftp://genome.wustl.edu/private/723378195268200/chimpqual
     # oops -- no longer at this site -- inquiring with Joanne Nelson
 
     # Also, create dummy chrY, chrY_random,  chr21, all of which
     # are finished sequence from clones
 cat > makeFakeQuals.csh << 'EOF'
     date
     foreach c (chr21 chrY chrY_random)
         echo $c
         set d = `echo $c | sed 's/chr//;s/_random//'`
         # pick a small set of quality scores, needed so utils don't complain
         nice zcat wustl/chr22_random.contigs.qvl.gz | \
             nice qaToQac stdin stdout | \
             nice qacAgpLift -mScore=98 /cluster/data/panTro2/$d/$c.agp \
                 stdin qac/$c.qac
     end
     date
 'EOF'
 
     # And a dummy chrM
     echo 'chrM	1	16554	1       W       Unknown  1	16554    +' \
         > agp.tmp
     nice zcat wustl/chr7.contigs.qvl.gz | \
         nice qaToQac stdin stdout | \
         nice qacAgpLift -mScore=98 agp.tmp stdin qac/chrM.qac
     
     # Load track table, excluding "dummy" chroms
 cat > toWig.csh << 'EOF'
     foreach f (`ls qac/chr*.qac | egrep -v 'chrY|chr21|chrM'`)
         qacToWig -fixed $f stdout
     end 
 'EOF'
     # << happy emacs
     date
     csh toWig.csh | nice wigEncode stdin quality.wig quality.wib
     date
 
     ssh hgwdev
     cd /cluster/data/panTro2/bed/quality
     ln -s /cluster/data/panTro2/bed/quality/quality.wib /gbdb/panTro2/wib
     hgLoadWiggle panTro2 quality quality.wig
 
     # create downloads for quality scores
     ssh kkstore04
     cd /cluster/data/panTro2/bed/quality
     mkdir qa
 cat > makeDownload.csh >> 'EOF'
     foreach f (qac/chr*.qac)
         set c = $f:t:r
         echo $c
         qacToQa $f stdout | gzip > qa/$c.qa.gz
     end
     cd qa
     tar cvfz chromQuals.tar.gz chr*.qa.gz
 'EOF'
     rm -fr qa
     ssh hgwdev
     ln -s /cluster/data/panTro2/bed/quality/chromQuals.tar.gz \
         /usr/local/apache/htdocs/goldenPath/panTro2/bigZips
     # edit README.txt and update md5sum.txt
 
 
 ########################
 # CREATE DATABASE AND GRP TABLE (DONE 2006-02-28 kate)
     ssh hgwdev
     hgsql '' -e 'create database panTro2'
     # Use df to make sure there is at least 75G free on hgwdev:/var/lib/mysql
     df -h /var/lib/mysql
         # /dev/sdc1             1.8T  1.6T   71G  96% /var/lib/mysql
     # oops -- not enough ?
 
     hgsql panTro2 -e \
       "create table grp (PRIMARY KEY(NAME)) select * from mm8.grp"
     # remove encode track groups
     hgsql panTro2 -e "delete from grp where name like 'encode%'"
     
 
 ########################
 # MAKE CHROMINFO TABLE WITH (TEMPORARILY UNMASKED) 2BIT (DONE 2006-02-28 kate)
     # Redo to fix chr9 (2006-03-02 kate)
     # Make .2bit, unmasked until RepeatMasker and TRF steps are done.
     # Do this now so we can load up RepeatMasker and run featureBits; 
     # can also load up other tables that don't depend on masking.  
     ssh kkstore01
     cd /cluster/data/panTro2
     nice faToTwoBit ?{,?}/chr*.fa 6_hla_hap1/chr*.fa panTro2.2bit
 
     mkdir bed/chromInfo
     twoBitInfo panTro2.2bit stdout \
     | awk '{print $1 "\t" $2 "\t/gbdb/panTro2/panTro2.2bit";}' \
       > bed/chromInfo/chromInfo.tab
     
     # Make symbolic links from /gbdb/panTro2/ to the real .2bit.
     ssh hgwdev
     mkdir /gbdb/panTro2
     ln -s /cluster/data/panTro2/panTro2.2bit /gbdb/panTro2/
     # Load /gbdb/panTro2/panTro2.2bit paths into database and save size info.
     cd /cluster/data/panTro2
     hgsql panTro2  < $HOME/kent/src/hg/lib/chromInfo.sql
     hgsql panTro2 -e 'load data local infile \
       "/cluster/data/panTro2/bed/chromInfo/chromInfo.tab" \
       into table chromInfo;'
     echo "select chrom,size from chromInfo" | hgsql -N panTro2 > chrom.sizes
     wc chrom.sizes
        # 52
 
 # GOLD AND GAP TRACKS (DONE 2006-02-28 kate)
     # Redo to fix chr9 (2006-03-02 kate)
     # Change chrY centromere to unbridged (2006-07-11 kate)
     ssh hgwdev
     cd /cluster/data/panTro2
     cp -p wustl/chrY.agp wustl/chrY.agp.old
     sed '/centromere/s/yes/no/' wustl/chrY.agp.old > wustl/chrY.agp
     cp wustl/chrY.agp Y
     cat wustl/*.agp | grep -v chrM | hgGoldGapGl -noGl panTro2 stdin
 
 [kate@hgwdev panTro2]$ nice featureBits panTro2 -countGaps -noRandom gap
 423043673 bases of 3175632892 (13.322%) in intersection
 [kate@hgwdev panTro2]$ nice featureBits panTro1 -countGaps -noRandom gap
 671583599 bases of 3083993401 (21.776%) in intersection
 [kate@hgwdev panTro2]$ nice featureBits hg15 -countGaps -noRandom gap
 238329157 bases of 3070537687 (7.762%) in intersection
 [kate@hgwdev panTro2]$ nice featureBits hg18 -countGaps -noRandom gap
 222401287 bases of 3091592211 (7.194%) in intersection
 
 chr1   12771775 bases of 229974691 (5.554%) in intersection
 chr10   9298149 bases of 135001995 (6.887%) in intersection
 chr11   10601065 bases of 134204764 (7.899%) in intersection
 chr12   5495847 bases of 135371336 (4.060%) in intersection
 chr13   28069601 bases of 115868456 (24.225%) in intersection
 chr14   21093350 bases of 107349158 (19.649%) in intersection
 chr15   23087195 bases of 100063422 (23.073%) in intersection
 chr16   16171896 bases of 90682376 (17.834%) in intersection
 chr17   9948345 bases of 83384210 (11.931%) in intersection
 chr18   3076876 bases of 77261746 (3.982%) in intersection
 chr19   12470245 bases of 64473437 (19.342%) in intersection
 chr20   4187763 bases of 62293572 (6.723%) in intersection
 chr21   13764311 bases of 46489110 (29.608%) in intersection
 chr22   17821597 bases of 50165558 (35.526%) in intersection
 chr2a   8581091 bases of 114460064 (7.497%) in intersection
 chr2b   120728260 bases of 248603653 (48.563%) in intersection
 chr3   8990241 bases of 203962478 (4.408%) in intersection
 chr4   7932593 bases of 194897272 (4.070%) in intersection
 chr5   8761164 bases of 183994906 (4.762%) in intersection
 chr6   9202606 bases of 173908612 (5.292%) in intersection
 chr6_hla_hap1   0 bases of 34169 (0.000%) in intersection
 chr7   9183564 bases of 160261443 (5.730%) in intersection
 chr8   6927299 bases of 145085868 (4.775%) in intersection
 #chr9   115285361 bases of 224587821 (51.332%) in intersection
 #  Seems too big -- notified LaDeana of huge gap (80MB)
 #  She fixed the AGP, and here's the result
 chr9    29207531 bases of 138509991 (21.087%) in intersection
 chrM   0 bases of 16554 (0.000%) in intersection
 chrUn   8188633 bases of 58616431 (13.970%) in intersection
 chrX   24409881 bases of 155361357 (15.712%) in intersection
 chrY   1261428 bases of 23952694 (5.266%) in intersection
 
 
 # HGCENTRALTEST ENTRY AND TRACKDB TABLE (DONE 2006-03-03 kate)
     ssh hgwdev
     cd $HOME/kent/src/hg/makeDb/trackDb
     cvs up -d -P
     # Edit that makefile to add in all the right places and do
     make update
     cvs commit makefile
     mkdir -p chimp/panTro2
     cvs add chimp/panTro2
     cvs ci -m "trackDb dir for second chimp genome(s)" chimp/panTro2
     # Do this in a clean (up-to-date, no edits) tree:
     make alpha DBS=panTro2
 
     # Add dbDb entry (not a new organism so defaultDb and genomeClade already 
     # have entries):
     hgsql -h genome-testdb hgcentraltest \
       -e 'insert into dbDb (name, description, nibPath, organism,  \
           defaultPos, active, orderKey, genome, scientificName,  \
           htmlPath, hgNearOk, hgPbOk, sourceName)  \
           values("panTro2", "Jan. 2006", \
           "/gbdb/panTro2", "Chimp", "chr7:115705331-115981791", 1, \
           15, "Chimp", "Pan troglodytes", \
           "/gbdb/panTro2/html/description.html", 0, 0, \
           "Washington University Build 2 Version 1");'
 
 
 # SPLIT SEQUENCE FOR REPEATMASKER (2006-03-01 kate)
 # Split into 500K chunks, at gaps if possible
 
     ssh kkstore01
     cd /cluster/data/panTro2
     mkdir split500K
     foreach d (`cat chrom.lst`)
         set c = chr$d
         mkdir split500K/$c
         faSplit gap -minGapSize=100 -verbose=2 $d/$c.fa 500000 \
                 split500K/$c/${c}_ -lift=split500K/$c.lft
         set r = ${c}_random
         if (-e $d/$r.fa) then
             mkdir split500K/$r
             faSplit gap -minGapSize=100 -verbose=2 $d/$r.fa 500000 \
                 split500K/$r/${r}_ -lift=split500K/$r.lft
         endif
     end
     mkdir lifts
     cat split500K/*.lft > lifts/split500K.lft
 
     # Redo chr9 (2006-03-02 kate)
     rm -fr split500K/chr9
     mkdir split500K/chr9
     faSplit gap -minGapSize=100 -verbose=2 9/chr9.fa 500000 \
                 split500K/chr9/chr9_ -lift=split500K/chr9.lft
     # forgot to redo this -- doing it now (2006-03-21 kate)
     cat split500K/*.lft > lifts/split500K.lft
 
 
 # REPEATMASKER RUN (2006-03-01 kate)
 # Using -species pan option to get chimp repeats
 # Redone 3/10 with additional chimp repeats from
 # Evan Eichler (left out of -pan lib)
 
     ssh pk
     cd /cluster/data/panTro2
     mkdir RMRun
     cd RMRun
 
     # Record RM version used
     ls -l /cluster/bluearc/RepeatMasker
         #lrwxrwxrwx    1 hiram    protein        18 Jan 20 13:13 /cluster/bluearc/RepeatMasker -> RepeatMasker060120/
     cat /cluster/bluearc/RepeatMasker/Libraries/version > RM.version
         #RM database version 20060120
     # Run RepeatMasker on a dummy input, just to make it initialize its chimp 
     # libraries once before the cluster run
     /cluster/bluearc/RepeatMasker/RepeatMasker -spec pan /dev/null
         # Building species libraries in: /cluster/bluearc/RepeatMasker060120/Libraries/20060120/pan
 
     pushd  /cluster/data/panTro2/split500K
     ls -1S */*.fa | sort > ../RMRun/split.lst
     popd
 
     cat << 'EOF' > template
 #LOOP
 ./RMChimp.csh $(dir1) $(root1) {check out line out/$(dir1)/$(root1).out}
 #ENDLOOP
 'EOF'
     # << for emacs
     
 cat > RMChimp.csh << 'EOF'
 #!/bin/csh  -ef
 set d = /cluster/data/panTro2
 set tmp = /scratch/tmp/panTro2/$2
 mkdir -p $tmp
 mkdir -p out/$1
 cp $d/split500K/$1/$2.fa $tmp
 pushd $tmp
 /cluster/bluearc/RepeatMasker060120/RepeatMasker -ali -s -species pan $2.fa
 popd
 cp -p $tmp/$2.fa.out $3
 if (-e $tmp/$2.fa.align) cp $tmp/$2.fa.align out/$1
 if (-e $tmp/$2.fa.tbl) cp $tmp/$2.fa.tbl out/$1
 if (-e $tmp/$2.fa.cat) cp $tmp/$2.fa.cat out/$1
 rm -fr $tmp/*
 rmdir --ignore-fail-on-non-empty $tmp
 rmdir --ignore-fail-on-non-empty /scratch/tmp/panTro2
 'EOF'
     # << for emacs
     chmod +x RMChimp.csh
 
     gensub2 split.lst single template jobList
     para create jobList
         # 6538 jobs
     para try
 
 # CPU time in finished jobs:   32493557s  541559.29m  9025.99h  376.08d  1.030 y
 # IO & Wait Time:                 40999s     683.31m    11.39h    0.47d  0.001 y
 # Average job time:                4976s      82.94m     1.38h    0.06d
 # Longest running job:                0s       0.00m     0.00h    0.00d
 # Longest finished job:            7607s     126.78m     2.11h    0.09d
 # Submission to last job:        121890s    2031.50m    33.86h    1.41d
 
 #  34 hours
 
     # Lift up the 500KB .out's to chrom level
     ssh kkstore01
     cd /cluster/data/panTro2/RMRun
 cat > lift.csh << 'EOF'
     foreach d (`cat ../chrom.lst`)
         set c = chr$d
         echo $c
         liftUp ../$d/$c.fa.out ../lifts/split500K.lft error out/$c/*.out
         set r = ${c}_random
         if (-e out/$r) then
             echo $r
             liftUp ../$d/$r.fa.out ../lifts/split500K.lft error out/$r/*.out
         endif
     end
 'EOF'
     # << for emacs
     csh lift.csh >&! lift.log &
 
     # Redo chr9 after remaking split500K.lft (2006-03-24 kate)
     liftUp ../9/chr9.fa.out ../lifts/split500K.lft error out/chr9/*.out
 
     # Load .outs into database
     # Redone after relift (2006-03-24 kate)
     ssh hgwdev
     cd /cluster/data/panTro2
     ls */*.fa.out | wc -l
         # 52 (one per chrom)... load them up
     hgLoadOut panTro2 */*.fa.out
 # Several messages like: Strange perc. field -3.1 line 327855 of 1/chr1.fa.out
 # And at end: 953 records dropped due to repStart > repEnd
 # These are noted in rn4 and other make docs.
 # TODO: run it with -verbose=2 and send output to Robert Hubley/Arian Smit
 
 
 # VERIFY REPEATMASKER RESULTS (2006-03-03 kate)
     # TODO: Eyeball some repeat annotations in the browser, compare to lib seqs.
     # Run featureBits on previous genome build, and compare:
     #nice featureBits panTro2 rmsk
         # w/o 2 Eichler repeats:
         # 1396630330 bases of 2909512873 (48.002%) in intersection
         # with 2 Eichler repeats:
         # 1401134418 bases of 2909512873 (48.157%) in intersection
     nice featureBits panTro1 rmsk
         # 1311281819 bases of 2733948177 (47.963%) in intersection
         # after chr9 .out reload...
         # 1311281819 bases of 2733948177 (47.963%) in intersection
 
     # looks good -- similar but just a bit higher
 
     ssh hgwdev
     cd /cluster/data/panTro2/RMRun
     # take a look at 
      hgsql panTro2 -e "select distinct(repFamily) from chr1_rmsk" > repFamilies.panTro2
      hgsql panTro2 -e "select distinct(repName) from chr1_rmsk where repFamily <> 'Simple_repeat' and repFamily <> 'Low_complexity' order by repName" > repNames.panTro2
     wc -l repNames.panTro2
      hgsql panTro1 -e "select distinct(repFamily) from chr1_rmsk" > repFamilies.panTro1
      hgsql panTro1 -e "select distinct(repName) from chr1_rmsk" > repNames.panTro1
      hgsql panTro1 -e "select distinct(repName) from chr1_rmsk where repFamily <> 'Simple_repeat' and repFamily <> 'Low_complexity' order by repName" > repNames.panTro1
     wc -l *Names*
         # 864 repNames.panTro1
         # 906 repNames.panTro2
             # 40 new repeats (+ 2 added from Eichler)
     wc -l *Families*
         # 38 repFamilies.panTro1
         # 38 repFamilies.panTro2
             # no new repeat families
 
 
 # RUN TRF TO IDENTIFY SIMPLE REPEATS (2006-03-01 kate)
 
     ssh kkstore01
     cd /cluster/data/panTro2
 
     # First break up sequence into 5MB chunks at contigs/gaps
     # NOTE: Probably unnecessary -- Hiram has run on whole chroms
     # using small cluster
     mkdir split5M
     foreach d (`cat chrom.lst`)
         set c = chr$d
         mkdir split5M/$c
         faSplit gap -minGapSize=100 -verbose=2 $d/$c.fa 5000000 \
                 split5M/$c/${c}_ -lift=split5M/$c.lft
         set r = ${c}_random
         if (-e $d/$r.fa) then
             mkdir split5M/$r
             faSplit gap -minGapSize=100 -verbose=2 $d/$r.fa 5000000 \
                 split5M/$r/${r}_ -lift=split5M/$r.lft
         endif
     end
     mkdir lifts
     cat split5M/*.lft > lifts/split5M.lft
 
     # Run TRF
     mkdir -p bed/simpleRepeat/trf
     cd bed/simpleRepeat
 cat > runTrf.csh << 'EOF'
     date
     foreach d (`find  /cluster/data/panTro2/split5M/* -type d`)
         cd $d
         foreach f (*.fa)
             set b = $f:r
             echo $f
             trfBig -trf=/cluster/bin/i386/trf $f /dev/null -tempDir=/tmp \
                     -bedAt=/cluster/data/panTro2/bed/simpleRepeat/trf/$b.bed
         end
     end
     date
 'EOF'
     csh runTrf.csh >&! runTrf.log &
 
     # ~6 hours run-time
 
     # Redo chr9 (2006-03-02 kate)
     cd /cluster/data/panTro2
     rm -fr split5M/chr9
     mkdir split5M/chr9
     faSplit gap -minGapSize=100 -verbose=2 9/chr9.fa 5000000 \
                 split5M/chr9/chr9_ -lift=split5M/chr9.lft
     cd bed/simpleRepeat
     cp runTrf.csh runTrf.chr9.csh
      # edit for just chr9
     csh runTrf.chr9.csh >&! runTrf.chr9.log &
     # check for EOF in output files to give some assurance
     endsInLf trf/*
         # trf/chrM_0.bed zero length
 
     # lift to chromosome coordinates
     liftUp simpleRepeat.bed /cluster/data/panTro2/lifts/split5M.lft \
          warn trf/*.bed
 
     # forgot to include the new lift in the all split (2006-03-21 kate)
     # so do it now and reload table
     cd /cluster/data/panTro2
     cat split5M/*.lft > lifts/split5M.lft
     cd bed/simpleRepeat
     liftUp simpleRepeat.bed /cluster/data/panTro2/lifts/split5M.lft \
          warn trf/*.bed
 
     # load into database
     ssh hgwdev
     cd /cluster/data/panTro2/bed/simpleRepeat
     hgLoadBed panTro2 simpleRepeat simpleRepeat.bed \
       -sqlTable=$HOME/kent/src/hg/lib/simpleRepeat.sql
   
     nice featureBits panTro2 simpleRepeat
         # 82897828 bases of 2909512873 (2.849%) in intersection
     nice featureBits panTro1 simpleRepeat
         # 53732632 bases of 2733948177 (1.965%) in intersection
     # seems quite a bit higher -- compare it to hg17 chrom1
 
         # nice featureBits panTro2 simpleRepeat -chrom=chr1
             # 3068104 bases of 217202916 (1.413%) in intersection
         nice featureBits hg17 simpleRepeat -chrom=chr1
             # 3438627 bases of 222827847 (1.543%) in intersection
                     # darn similar
 
 
 # PROCESS SIMPLE REPEATS INTO MASK
     # After the simpleRepeats track has been built, make a filtered version 
     # of the trf output: keep trf's with period <= 12:
 
     ssh kkstore01
     cd /cluster/data/panTro2/bed/simpleRepeat
     mkdir -p trfMask trfMaskChrom
     foreach f (trf/chr*.bed)
         awk '{if ($5 <= 12) print;}' $f > trfMask/$f:t
     end
 
     # lift filtered simple repeats to chrom coordinates, by chrom
 cat > liftChroms.csh << 'EOF'
     foreach d (`cat /cluster/data/panTro2/chrom.lst`)
         set c = chr$d
         echo $c
         liftUp trfMaskChrom/$c.bed /cluster/data/panTro2/split5M/$c.lft \
                 warn trfMask/${c}_?{,?}.bed 
         set r = ${c}_random
         if (-e /cluster/data/panTro2/$d/$r.fa) then
             liftUp trfMaskChrom/$r.bed /cluster/data/panTro2/split5M/$r.lft \
                 warn trfMask/${r}_?{,?}.bed 
         endif
     end
 'EOF'
     # << for emacs
     csh liftChroms.csh >&! liftChroms.log &
 
     # check coverage of filtered TRF
     hgwdev
     cd /cluster/data/panTro2/bed/simpleRepeat
     cat trfMaskChrom/*.bed > /tmp/filtTrf.bed
     featureBits panTro2 /tmp/filtTrf.bed
         # 28800863 bases of 2909512873 (0.990%) in intersection
 
 
 # MASK SEQUENCE WITH REPEATMASKER AND SIMPLE REPEAT/TRF  (2006-03-03)
     
     # soft-mask (lower-case) the chrom fasta files
     # then make hard-masked versions from the soft-masked
     # note: next time simplify this (as below for chr9)
     ssh kkstore01
     cd /cluster/data/panTro2
     set trfChr=bed/simpleRepeat/trfMaskChrom
     foreach f (?{,?}/chr*.fa 6_hla_hap1/chr*.fa)
       cp $f $f.unmasked
       echo "repeat- and trf-masking $f"
       maskOutFa -soft $f $f.out $f
       cp $f $f.rmsk
       set c = $f:t:r
       maskOutFa -softAdd $f $trfChr/$c.bed $f
       echo "hard-masking $f"
       maskOutFa $f hard $f.masked
     end
 
     # remask chr9 after relifting .out's (2006-03-24)
     set f = 9/chr9.fa
     set c = chr9
     maskOutFa -soft $f.unmasked $f.out $f.rmsk
     maskOutFa -softAdd $f.rmsk $trfChr/$c.bed $f
     maskOutFa $f hard $f.masked
 
     # a few messages (e.g. 10 per chrom) like:
     # WARNING: negative rEnd: -531 chr1:14090705-14090763 L1MEc
 
     # after checking, remove .unmasked and .rmsk
 
     # make 2bit for browser and blat
     # redone for chr9 correction (2006-03-24 kate)
     faToTwoBit ?{,?}/chr*.fa 6_hla_hap1/chr*.fa panTro2.2bit
     mkdir /cluster/bluearc/scratch/hg/panTro2
     ssh kkr1u00
     cp -p /cluster/data/panTro2/panTro2.2bit /iscratch/i/panTro2/panTro2.2bit
     ~kent/bin/iSync
 
 
     # after chr9 remask:
     cp -p /cluster/data/panTro2/panTro2.2bit /iscratch/i/panTro2
     cp -p /cluster/data/panTro2/nib/chr9.nib /iscratch/i/panTro2/nib
     ~kent/bin/iSync
 
     # make nibs for blastz w/linSpecRep
     mkdir nib
     foreach f (?{,?}/chr*.fa 6_hla_hap1/chr*.fa)
       echo $f:t:r
       faToNib -softMask $f nib/$f:t:r.nib
     end
     # remake chr9 (2006-03-24 kate)
     faToNib -softMask 9/chr9.fa nib/chr9.nib
 
 
 # MAKE 11.OOC 
 #   #Remade after chr9 remask (2006-03-24 kate)
 
     ssh kolossus
     cd /cluster/data/panTro2
     mkdir /cluster/bluearc/panTro2
     blat panTro2.2bit /dev/null /dev/null \
         -tileSize=11 -makeOoc=/cluster/bluearc/panTro2/11.ooc -repMatch=1024
     cp -p /cluster/bluearc/panTro2/11.ooc .
 
 
 # COPY SEQUENCE TO CLUSTER NODES FOR CLUSTER RUNS
 
     set scratch = /cluster/bluearc/scratch/hg/panTro2
     mkdir $scratch
     cp -p panTro2.2bit $scratch
     cp -rp nib $scratch
     cp -rp chrom.sizes $scratch
     cp -p 11.ooc $scratch
     # request cluster-admin sync /scratch/hg/panTro2 to cluster nodes
 
     # update 2bit and chr9 nib after updating chr9 masking
     cp -p panTro2.2bit $scratch
     cp -p nib/chr9.nib $scratch
     cp -p 11.ooc $scratch
     # request cluster-admin sync /scratch/hg/panTro2 to cluster nodes
 
 
 
 # GENBANK ALIGNMENTS (DONE 2006-03-25 kate/markd)
     # restarted 3/13 4pm
 
     # Make a lift file that identifies gap locations for Genbank
     # This is a pseudo "liftAll" made from the AGP files, however negative
     #  strand contigs are renamed to indicate they need to be
     #  reversed. 
     
     ssh kkstore01
     cd /cluster/data/panTro2
     cat ?{,?}/*.agp 6_hla_hap1/*.agp | agpToLift -revStrand > \
         lifts/genbank.lft
     
     ssh hgwdev
     cd ~/kent/src/hg/makeDb/genbank
     cvsup
     # edit etc/genbank.conf to add panTro2
     # the hack in the genbank code for panTro1 as changed for
     # panTro2.  panTro2 treats both chimp and human cDNAs as
     # native, and all else as xeno.  Don't do xeno ESTs any more.
 
 panTro2.serverGenome = /cluster/data/panTro2/panTro2.2bit
 panTro2.clusterGenome = /scratch/hg/panTro2/panTro2.2bit
 panTro2.ooc = /cluster/bluearc/panTro2/11.ooc
 panTro2.align.unplacedChroms = chrUn,chr*_random
 panTro2.lift = /cluster/data/panTro2/lifts/genbank.lft
 panTro2.refseq.mrna.native.pslCDnaFilter  = ${ordered.refseq.mrna.native.pslCDnaFilter}
 panTro2.refseq.mrna.xeno.pslCDnaFilter    = ${ordered.refseq.mrna.xeno.pslCDnaFilter}
 panTro2.genbank.mrna.native.pslCDnaFilter = ${ordered.genbank.mrna.native.pslCDnaFilter}
 panTro2.genbank.mrna.xeno.pslCDnaFilter   = ${ordered.genbank.mrna.xeno.pslCDnaFilter}
 panTro2.genbank.est.native.pslCDnaFilter  = ${ordered.genbank.est.native.pslCDnaFilter}
 panTro2.genbank.est.xeno.pslCDnaFilter    = ${ordered.genbank.est.xeno.pslCDnaFilter}
 panTro2.downloadDir = panTro2
 panTro2.genbank.est.xeno.load = no
 panTro2.refseq.mrna.native.load  = yes
 panTro2.refseq.mrna.xeno.load = yes
 panTro2.refseq.mrna.xeno.loadDesc  = yes
 
     cvs ci etc/genbank.conf
     make etc-update
 
     ssh kkstore02
     cd /cluster/data/genbank
     nice bin/gbAlignStep -initial panTro2 &
     # various problems with ssh, bluearc, etc....
     nice bin/gbAlignStep -initial -continue=copy panTro2 &
 
     ssh hgwdev
     nice bin/gbDbLoadStop -initialLoad -drop panTro2 &
   
 
 
 # MAKE LINEAGE-SPECIFIC REPEATS FOR BLASTZ (2006-03-15 kate)
 
     ssh kkstore01
     cd /cluster/data/panTro2
     mkdir linSpecRep
     cd linSpecRep
     mkdir dates
     cd dates
 
     # run Arian's script to annotate .out files with species-specificity
     # include all currently supported comparison species
 cat > outRepeats.csh << 'EOF'
     date
     set d = /cluster/data/panTro2
     foreach f ($d/?{,?}/chr*.fa.out $d/6_hla_hap1/chr*.fa.out)
         set c = $f:t:r:r
         echo $c
         cp $f .
         /cluster/bluearc/RepeatMasker060120/DateRepeats $c.fa.out -query human \
             -comp mouse -comp rat -comp dog -comp cow -comp rabbit
         mv $c.fa.*cuniculus $c.fa.dates.out
         end
     date
 'EOF'
     # << for emacs
     csh outRepeats.csh >&! outRepeats.log &
         # 25 minutes
 
     # Redo for fixed chr9 (2006-03-24 kate)
     cp /cluster/data/panTro2/9/chr9.fa.out .
     set c = chr9
     /cluster/bluearc/RepeatMasker060120/DateRepeats $c.fa.out -query human \
             -comp mouse -comp rat -comp dog -comp cow -comp rabbit
         mv $c.fa.*cuniculus $c.fa.dates.out
     
     # run our script to extract lineage-specific by column 
     #   from the annotated .outs
     # just run on chr1 to assess which variants produce any
     #   differences.
     mkdir testChr1
     cd testChr1
 cat > testChr1.csh << 'EOF'
     set script = /cluster/bin/scripts/extractRepeats
     set f = ../chr1.fa.dates.out
     echo mouse
     $script 1  $f > notInMouse.out.spec
     echo rat
     $script 2 $f > notInRat.out.spec
     echo dog
     $script 3 $f > notInDog.out.spec
     echo cow
     $script 4 $f > notInCow.out.spec
     echo rabbit
     $script 5 $f > notInRabbit.out.spec
 'EOF'
     csh testChr1.csh >&! testChr1.log &
     ls -l *.spec
 -rw-rw-r--  1 kate protein 17678799 Mar 15 09:51 notInCow.out.spec
 -rw-rw-r--  1 kate protein 17678799 Mar 15 09:51 notInDog.out.spec
 -rw-rw-r--  1 kate protein 17773939 Mar 15 09:48 notInHuman.out.spec
 -rw-rw-r--  1 kate protein 17773939 Mar 15 09:50 notInMouse.out.spec
 -rw-rw-r--  1 kate protein 17773939 Mar 15 09:51 notInRabbit.out.spec
 -rw-rw-r--  1 kate protein 17773939 Mar 15 09:50 notInRat.out.spec
 
     # this indicates there are 2 unique LSR files:
     # 2) notInRodent (includes Mouse, Rat, Rabbit)
     # 3) notInOther  (includes Cow & Dog)
 
     # Extract mouse (column 1) and dog (column 3)
     cd /cluster/data/panTro2/linSpecRep
     mkdir notInRodent notInOthers
 cat > makeLSR.csh << 'EOF'
     foreach f (dates/*.fa.dates.out)
         set c = $f:t:r:r:r
         echo $c
         set out = $c.out.spec
         /cluster/bin/scripts/extractRepeats 1 $f > notInRodent/$out
         /cluster/bin/scripts/extractRepeats 3 $f > notInOthers/$out
     end
 'EOF'
     # << happy emacs
     csh makeLSR.csh >&! makeLSR.log &
 
     # copy to bluearc for blastz runs
     set d = /cluster/bluearc/panTro2/linSpecRep
     mkdir -p $d
     cp -rp notInRodent $d
     cp -rp notInOthers $d
 
     # update for chr9 (2006-03-24 kate)
     set c = chr9
     extractRepeats 1 dates/$c.fa.dates.out > \
                 notInRodent/$c.out.spec
     extractRepeats 3 dates/$c.fa.dates.out > \
                 notInOthers/$c.out.spec
     cp -rp notInRodent $d
     cp -rp notInOthers $d
 
 
 # CREATE LINEAGE-SPECIFIC REPEATS FOR BLASTZ WITH ZEBRAFISH (2006-04-07 kate)
 # Treat all repeats as lineage-specific
     ssh kkstore04
     cd /cluster/data/panTro2
     set d = /cluster/bluearc/panTro2/linSpecRep/notInZebrafish
     mkdir $d
     foreach f ({{?{,?}},6_hla_hap1}/chr*.fa.out)
     set c = $f:t:r:r
     echo $c
     cp -p $f $d/$c.out.spec
     end
 
 
 # SET UP BLAT SERVER (2006-03-06 kate)
 # Request blat server from cluster-admin, then enter ports in database
 
     ssh hgwdev
     hgsql -e 'INSERT INTO blatServers (db, host, port, isTrans, canPcr) \
         VALUES ("panTro2", "blat12", "17790", "1", "0"); \
         INSERT INTO blatServers (db, host, port, isTrans, canPcr) \
         VALUES ("panTro2", "blat12", "17791", "0", "1");' \
             hgcentraltest
     #   test it with some sequence
 
 
 # GC 5 BASE TRACK (2006-03-06 kate)
 
     ssh kkstore01
     cd /cluster/data/panTro2/bed
     mkdir gc5base
     cd gc5base
     time hgGcPercent -wigOut -doGaps -file=stdout -win=5 panTro2 \
               /cluster/data/panTro2 | wigEncode stdin gc5Base.wig gc5Base.wib
 
     sh hgwdev
     cd /cluster/data/panTro2/bed/gc5base
     mkdir -p /gbdb/panTro2/wib
     ln -s `pwd`/*.wib /gbdb/panTro2/wib
     hgLoadWiggle panTro2 gc5Base gc5Base.wig
 
 
 # CPGISLANDS (2006-03-06 - kate)
 #       Redo 3/13/06 after second Repeatmasker run
 #       Redone for chr9 remasking (2006-03-24 kate)
     ssh hgwdev
     mkdir /cluster/data/panTro2/bed/cpgIsland
     cd /cluster/data/panTro2/bed/cpgIsland
 
     # Build software from Asif Chinwalla (achinwal@watson.wustl.edu)
     cvs co hg3rdParty/cpgIslands
     cd hg3rdParty/cpgIslands
     make
     #   gcc readseq.c cpg_lh.c -o cpglh.exe
     cd ../..
     ln -s hg3rdParty/cpgIslands/cpglh.exe .
 
     # cpglh.exe requires hard-masked (N) .fa's.  
     # There may be warnings about "bad character" for IUPAC ambiguous 
     # characters like R, S, etc.  Ignore the warnings.  
     ssh kkstore01
     cd /cluster/data/panTro2/bed/cpgIsland
 cat > doCpg.csh << 'EOF'
     foreach f (../../?{,?}/chr*.fa.masked ../../6_hla_hap1/chr*.fa.masked)
         set c = $f:t:r:r
         echo $c
         ./cpglh.exe $f > $c.cpg
     end
 'EOF'
     # << happy emacs
     csh doCpg.csh >&! doCpg.log &
 
     #   Several chroms have 0 results:
     #   -rw-rw-r--  1     0 Feb 16 15:19 chr10_random.cpg
     #   -rw-rw-r--  1     0 Feb 16 15:20 chr15_random.cpg
     #   -rw-rw-r--  1     0 Feb 16 15:22 chr8_random.cpg
     #   -rw-rw-r--  1     0 Feb 16 15:22 chr9_random.cpg
     #   -rw-rw-r--  1     0 Feb 16 15:22 chrM.cpg
     #   -rw-rw-r--  1     0 Feb 16 15:22 chrX_random.cpg
     #   -rw-rw-r--  1     0 Feb 16 15:22 chrY.cpg
 
     # Transform cpglh output to bed +
     cat << '_EOF_' > filter.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);
 }
 '_EOF_'
     #   happy emacs
     awk -f filter.awk chr*.cpg | sort -k1,1 -k2,2n > cpgIsland.bed
 
     ssh hgwdev
     cd /cluster/data/panTro2/bed/cpgIsland
     hgLoadBed -strict panTro2 cpgIslandExt -tab -noBin \
       -sqlTable=$HOME/kent/src/hg/lib/cpgIslandExt.sql cpgIsland.bed
     featureBits panTro2 cpgIslandExt
         # 19412837 bases of 2909512873 (0.667%) in intersection
 
 
 # BAC Ends (2006-03-08 IN PROGRESS kate)
 
     # download BAC ends from NCBI.  
     ssh kkstore
     cd /cluster/data/ncbi/bacends/chimp
     mkdir bacends.chimp.2
     cd bacends.chimp.2
     # wget ftp.ncbi.nih.gov/genomes/CLONEEND/pan_troglodytes/*
     # trouble with wget on files in this dir -- they're hanging,
     # so I used ftp mget w/o prompt
     # creates files:  9598_clone_end???.mfa and 95998_clone_info???.txt
     # 10 files of each type
     # NOTE: the *.info files contain a superset of the data in the
     # previous cl_acc_gi_len file -- extract fields 1-5 and 8 to
     # use as input to Terry's script that generaes the .pair file
     # for compatibility with downstream tools
     ftp ftp.ncbi.nih.gov
     cd /genomes/CLONEEND/pan_troglodytes
     prompt
     mget *
     gunzip *.gz
 
     # setup for track build area
     cd /cluster/data/panTro2
     mkdir -p bed/bacends
     cd bed/bacends
     ln -s  /cluster/data/ncbi/bacends/chimp/bacends.chimp.2 ncbi
 
     # generate pairs info file by/for processing by Terry Furey's scripts
     grep -h -v '^#' ncbi/*.txt | \
         awk '{printf ("%s\t%s\t%s\t%s\t%s\t%s\n", $1, $2, $3, $4, $5, $8)}' \
                 > cl_acc_gi_len
     wc -l cl_acc_gi_len
         # 175757 cl_acc_gi_len
         # previous load (2004): 158588 entries
     faSize ncbi/*.mfa
         # 175757 sequences in 10 files
         # Info and sequence catch match -- good!
      /cluster/bin/scripts/convertBacEndPairInfo cl_acc_gi_len
         # 85911 pairs and 3935 singles
         # previous: 78846 pairs and 896 singles
         # Note large increase in singles -- either lower quality in
         # new data set, or script could use some updating.
         # Just note this for now.
     wc -l bacEndPairs.txt bacEndSingles.txt
         #  85911 bacEndPairs.txt
         #   3935 bacEndSingles.txt
 
     # create sequence file
     cp /cluster/data/ncbi/bacends/chimp/bacends.chimp.1/convert.pl .
     cat ncbi/*.mfa | ./convert.pl > BACends.fa
     faSize BACends.fa
         # 175757 sequences
 
     # make accessible to track
     ssh hgwdev
     mkdir -p /gbdb/panTro2/bacends
     ln -s /cluster/data/panTro2/bed/bacends/BACends.fa /gbdb/panTro1/bacends
 
     # split for aligning on cluster
     ssh kkstore01
     cd /cluster/data/panTro2/bed/bacends
     set tmp = /san/sanvol1/scratch/panTro2/bacends/
     mkdir -p $tmp
     faSplit BACends.fa sequence 10 $tmp
 
     # blat vs. unmasked sequence in 5M chunks
     # 652 chunks * 10 bacends files = 6520 jobs
     ssh pk
     cd /cluster/data/panTro2/bed/bacends
     mkdir run
     cd run
     mkdir ../out
 
     # list chrom chunks and bacends chunks
     set dir = /cluster/data/panTro2/split5M
     ls $dir/*/*.fa | sed "s^$dir/^^" > chromSplit.lst
     set dir = /san/sanvol1/scratch/panTro2/bacends
     ls $dir/*.fa | sed "s^$dir/^^" > bacends.lst
 
 cat > align.csh << 'EOF'
 #!/bin/csh -ef 
     set d = $1:h
     set f = $1:t
     set e = $2
     mkdir -p ../out/$d
     set tmp = /scratch/tmp/panTro2/align.$$
     mkdir -p $tmp
     cp /cluster/data/panTro2/split5M/$d/$f $tmp
     cp /san/sanvol1/scratch/panTro2/bacends/$e $tmp
     blat $tmp/$f $tmp/$e -ooc=/cluster/bluearc/panTro2/11.ooc ../out/$d/$f:r.$e:r.psl
     rm -fr $tmp
 'EOF'
     # << happy emacs
     chmod +x align.csh
 
 cat > gsub << 'EOF'
 #LOOP
 ./align.csh $(path1) $(path2) {check out line+ ../out/$(dir1)/$(root1).$(root2).psl}
 #ENDLOOP
 'EOF'
     # << happy emacs
 
     gensub2 chromSplit.lst bacends.lst gsub jobList
     para create jobList
         # 6520 jobs
     para try
     # NOTE: these run quickly
 
 
 ########################
 # MAKE DOWNLOADABLE SEQUENCE FILES (2006-04-02 kate)
 
     ssh kkstore04
     cd /cluster/data/panTro2
     mkdir -p downloads/bigZips downloads/chromosomes
     cd downloads
 cat > makeChroms.csh << 'EOF'
     date
     foreach f (../?{,?}/*.fa ../6_hla_hap1/*.fa)
         echo $f:t:r
         nice gzip -c $f > chromosomes/$f:t.gz
     end
     date
 'EOF'
     csh makeChroms.csh >&! makeChroms.log &
         # 14 minutes
 
     cat > bigZips.csh << 'EOF'
 cd /cluster/data/panTro2
 tar cvzf downloads/bigZips/chromAgp.tar.gz \
         ?{,?}/chr*.agp 6_hla_hap1/*.agp
 tar cvzf downloads/bigZips/chromFa.tar.gz \
         ?{,?}/chr*.fa 6_hla_hap1/chr*.fa
 tar cvzf downloads/bigZips/chromFaMasked.tar.gz \
         ?{,?}/chr*.fa.masked 6_hla_hap1/chr*.fa.masked 
 tar cvzf downloads/bigZips/chromOut.tar.gz \
         ?{,?}/chr*.fa.out 6_hla_hap1/chr*.fa.out
 cd /cluster/data/panTro2/bed/simpleRepeat
 tar cvzf ../../downloads/bigZips/chromTrf.tar.gz ./trfMaskChrom
 'EOF'
     csh bigZips.csh >&! bigZips.log &
 
     # get GenBank native mRNAs and refGene 
     ssh hgwdev
     cd /cluster/data/genbank
     time ./bin/i386/gbGetSeqs -db=panTro2 -native GenBank mrna \
         /cluster/data/panTro2/downloads/bigZips/mrna.fa
     cd /cluster/data/panTro2/downloads/bigZips
     cd /cluster/data/panTro2/downloads/bigZips
     foreach I (1000 2000 5000) 
         echo "upstream${I} working ... "
         featureBits panTro2 refGene:upstream:${I} -fa=stdout \
                 | gzip -c > upstream${I}.fa.gz
         echo "upstream${I} done"
     end
     #   real    11m25.493s
 
     ssh kkstore04
     cd /cluster/data/panTro2/downloads/chromosomes
     #   copy previous release README.txt
     scp hgwdev:/usr/local/apache/htdocs/goldenPath/panTro1/chromosomes/README.txt .
     #   edit it to bring it up to date
     md5sum *.gz > md5sum.txt
 
     cd /cluster/data/panTro2/downloads/bigZips
     gzip mrna.fa
     cp -p ../../panTro2.2bit .
     #   copy previous release README.txt
     scp hgwdev:/usr/local/apache/htdocs/goldenPath/panTro1/bigZips/README.txt .
     # edit README.txt to indicate proper version of sequence and
     #   RepeatMasker
     md5sum *.gz *.2bit README.txt > md5sum.txt
 
     ssh hgwdev
     set d = /usr/local/apache/htdocs/goldenPath/panTro2
     mkdir $d
     ln -s /cluster/data/panTro2/downloads/bigZips $d
     ln -s /cluster/data/panTro2/downloads/chromosomes $d
 
 
 ########################
 # BLASTZ HUMAN hg18 (2006-05-23 kate)
 # Do not use lineage-specific repeats, on recommendation of Arian Smit
 # Redo with human/chimp scoring matrix
 
     ssh pk
     cd /cluster/data/panTro2/bed
     mkdir blastz.hg18.2006-05-23
     rm blastz.hg18
     ln -s blastz.hg18.2006-05-23 blastz.hg18
     cd blastz.hg18
 
     cat << '_EOF_' > DEF
 # chimp vs human
 export PATH=/usr/bin:/bin:/usr/local/bin:/cluster/bin/penn:/cluster/bin/x86_64:/cluster/home/angie/schwartzbin:/parasol/bin
 
 BLASTZ=blastz.v7.x86_64
 BLASTZ_ABRIDGE_REPEATS=0
 
 # TARGET: Chimmp panTro2
 SEQ1_DIR=/scratch/hg/panTro2/panTro2.2bit
 SEQ1_CHUNK=10000000
 SEQ1_LAP=10000
 SEQ1_LEN=/cluster/bluearc/panTro2/chrom.sizes
 
 # QUERY: Human hg18 - single chunk big enough to run entire genome
 SEQ2_DIR=/scratch/hg/hg18/hg18.2bit
 SEQ2_LEN=/scratch/hg/hg18/chrom.sizes
 SEQ2_CHUNK=3000000000
 SEQ2_LAP=0
 
 BASE=/cluster/data/panTro2/bed/blastz.hg18.2006-05-23
 TMPDIR=/scratch/tmp
 '_EOF_'
 
     # << happy emacs
     /cluster/bin/scripts/doBlastzChainNet.pl -verbose=2 \
         -chainMinScore=5000 -chainLinearGap=medium -bigClusterHub=pk \
         -blastzOutRoot /san/sanvol1/scratch/panTro2/blastz.hg18 \
         `pwd`/DEF >&! blastz.out &
 
 ########################
 # BLASTZ HUMAN hg18 (2006-03-13 kate)
 # Do not use lineage-specific repeats, on recommendation of Arian Smit
 
     ssh pk
     cd /cluster/data/panTro2/bed
     mkdir blastz.hg18.2006-03-13
     ln -s blastz.hg18.2006-03-13 blastz.hg18
     cd blastz.hg18
 
     cat << '_EOF_' > DEF
 # chimp vs human
 export PATH=/usr/bin:/bin:/usr/local/bin:/cluster/bin/penn:/cluster/bin/x86_64:/cluster/home/angie/schwartzbin:/parasol/bin
 
 BLASTZ=blastz.v7.x86_64
 BLASTZ_ABRIDGE_REPEATS=0
 
 # TARGET: Chimmp panTro2
 SEQ1_DIR=/scratch/hg/panTro2/panTro2.2bit
 SEQ1_CHUNK=10000000
 SEQ1_LAP=10000
 SEQ1_LEN=/cluster/bluearc/panTro2/chrom.sizes
 
 # QUERY: Human hg18 - single chunk big enough to run entire genome
 SEQ2_DIR=/scratch/hg/hg18/hg18.2bit
 SEQ2_LEN=/scratch/hg/hg18/chrom.sizes
 SEQ2_CHUNK=3000000000
 SEQ2_LAP=0
 
 BASE=/cluster/data/panTro2/bed/blastz.hg18.2006-03-13
 TMPDIR=/scratch/tmp
 '_EOF_'
 
     # << happy emacs
     /cluster/bin/scripts/doBlastzChainNet.pl -verbose=2 \
         -chainMinScore=5000 -chainLinearGap=medium -bigClusterHub=pk \
         -blastzOutRoot /san/sanvol1/scratch/panTro2/blastz.hg18 \
         `pwd`/DEF >&! blastz.out &
 
     # failed during chaining step, because /scratch/hg/panTro2
     #   doesn't exist on minicluster nodes
     # for now, just symlink /iscratch/i/panTro2 there 
     #   so that the 2bit is available to this pipeline -- later
     # we'll make sure the nodes are completely synced up.
     /cluster/bin/scripts/doBlastzChainNet.pl -verbose=2 \
         -chainMinScore=5000 -chainLinearGap=medium -bigClusterHub=pk \
         -blastzOutRoot /san/sanvol1/scratch/panTro2/blastz.hg18 \
         -continue=chainRun \
         `pwd`/DEF >&! blastz.2.out &
 
     # chr19 failed with out-of-mem, so try it on kolossus
     # moved to hgwdev by request....
     ssh kolossus
     cd /cluster/data/panTro2/bed/blastz.hg18/axtChain/run
     csh chain.csh panTro2.2bit:chr19: chain/panTro2.2bit:chr19:.chain >&! \
         chainChr19.log &
     # Create run.time file and remove liftedChain so we can continue
 
     /cluster/bin/scripts/doBlastzChainNet.pl -verbose=2 \
         -chainMinScore=5000 -chainLinearGap=medium -bigClusterHub=pk \
         -continue=chainMerge \
         `pwd`/DEF >&! blastz.3.out &
 
     # SWAP ALIGNMENTS TO PANTRO2 (2006-03-19 kate)
     /cluster/bin/scripts/doBlastzChainNet.pl -swap \
         -workhorse=hgwdev64 \
         /cluster/data/panTro2/bed/blastz.hg18/DEF >&! swap.log &
 
     # Redo of chr9 (2006-03-28 kate)
     # Create blastz run dir for just chr9, and run manually on pk.
     ssh pk
     cd /cluster/data/panTro2/bed/blastz.hg18
     mkdir run.blastz.chr9
     cd run.blastz.chr9
     cp ../run.blastz/*.csh .
     grep chr9 ../run.blastz/jobList | grep -v random > jobList
     ln -s ../run.blastz/qParts .
     para create jobList
         # 28 jobs
     para push
     para time > run.time
 
     # merge chr9 psl's on small cluster ("cat" step in automated script)
     ssh kki
     cd /cluster/data/panTro2/bed/blastz.hg18
     mkdir run.cat.chr9
     cd run.cat.chr9
     cp ../run.cat/cat.csh .
     grep chr9 ../run.cat/jobList > jobList
     para create jobList
         # 14 jobs
     para push
     para time > run.time
 
     cd ../axtChain
     mkdir run.chr9
     cd run.chr9
     cp ../run/chain.csh .
     grep chr9 ../run/jobList > jobList
     mkdir chain
     para create jobList
         # 1 job
     para push
     para time > run.time
 
     # replace chr9 chain in all.chain
     ssh hgwdev64 # kolossus busy
     cd /cluster/data/panTro2/bed/blastz.hg18/axtChain/run
     mkdir chain
     set all = panTro2.hg18.all.chain.gz
     chainSplit chain/ ../$all
     rm chain/chr9.chain
     cp ../run.chr9/chain/*.chain chain
     chainMergeSort chain/*.chain | nice gzip -c > ../$all
     rm -fr chain/*
     cd ..
     mkdir chain
     chainSplit chain/ $all
 
     # automate the rest
     ssh pk
     cd /cluster/data/panTro2/bed/blastz.hg18
     rm -fr axtNet mafNet
     ssh hgwdev "rm -fr /usr/local/apache/htdocs/goldenPath/panTro2/vsHg18"
     /cluster/bin/scripts/doBlastzChainNet.pl -verbose=2 \
         -bigClusterHub=pk -workhorse=hgwdev64 \
         -continue=net \
         `pwd`/DEF >&! blastz.net.chr9.out &
     # done
 
     # swap after chr9 fix (2006-03-30 kate)
     rm -fr /cluster/data/hg18/bed/blastz.panTro2.swap
     /cluster/bin/scripts/doBlastzChainNet.pl -swap \
         -workhorse=kkr4u00 \
         /cluster/data/panTro2/bed/blastz.hg18/DEF >&! swap.log &
 
     # downloads failed -- need to remove the old dir
      ssh hgwdev "rm -fr /usr/local/apache/htdocs/goldenPath/hg18/vsPanTro2"
     /cluster/bin/scripts/doBlastzChainNet.pl -swap -continue=download\
         /cluster/data/panTro2/bed/blastz.hg18/DEF >&! swap.download.log &
 
     # Reciprocal best (DONE 2006-04-02 kate)
     # These are best alignments in both directions
     # Start with the best chains from the net (these are in the liftOver chain).
     # Swap to human-reference, and net them
     # Finally extract the best chains from this net, and make axt's
     # REDO with improved blastz params (2006-07-27 kate)
     ssh kolossus
     cd /cluster/data/hg18/bed
     ln -s blastz.panTro2.swap blastz.panTro2
     cd blastz.panTro2/axtChain
     mkdir rBestNet
 cat > rbest.csh << 'EOF'
     set swapChain = panTro2ToHg18.swap.chain
     echo "merge and swap $swapChain"
     zcat /cluster/data/panTro2/bed/liftOver/panTro2ToHg18.over.chain.gz | \ 
         chainSwap stdin stdout | \
         chainMergeSort stdin > $swapChain
     echo "split chain"
     chainSplit swapChain/ $swapChain
     echo "run netter"
     chainPreNet $swapChain /cluster/bluearc/hg18/chrom.sizes \
                           /cluster/bluearc/panTro2/chrom.sizes stdout | \
         chainNet stdin -minSpace=1 /cluster/bluearc/hg18/chrom.sizes \
                                    /cluster/bluearc/panTro2/chrom.sizes \
                         stdout /dev/null | \
         netSyntenic stdin rbest.noClass.net
     echo "split net"
     netSplit rbest.noClass.net rBestNet
     echo "extract chains from net"
     netChainSubset -verbose=0 rbest.noClass.net $swapChain stdout | \
         chainMergeSort stdin | \
         gzip -c > hg18.panTro2.rbest.chain.gz
     echo "split extracted chain"
     chainSplit rBestChain/ hg18.panTro2.rbest.chain.gz 
     cd ..
     mkdir axtRBestNet
     echo "extract axts"
     foreach f (axtChain/rBestNet/*.net)
         echo $f:t:r
         netToAxt $f axtChain/swapChain/$f:t:r.chain \
                 /san/sanvol1/scratch/hg18/hg18.2bit  \
                 /cluster/bluearc/panTro2/panTro2.2bit stdout | \
         axtSort stdin stdout | \
         gzip -c > axtRBestNet/$f:t:r.hg18.panTro2.net.axt.gz
     end
 'EOF'
     csh rbest.csh >&! rbest.log &
 
     # load chains for viewing on genome-test
     ssh hgwdev
     cd /cluster/data/hg18/bed/blastz.panTro2/axtChain
     set chain = hg18.panTro2.rbest.chain
     # skip for now -- load by request
     cd rBestChain
 cat > loadRBestChain.csh << 'EOF'
     date
     foreach c (`awk '{print $1;}' /cluster/bluearc/hg18/chrom.sizes`)
         echo $c
         set f = $c.chain
         if (! -e $f) then
           echo no chains for $c
           set f = /dev/null
         endif
         hgLoadChain hg18 ${c}_chainRBestPanTro2 $f
     end
     date
 'EOF'
     cd ..
     #end skip for now
 
     # compare to previous set (other blastz run) and also full set
     featureBits hg18 chainRBestOldPanTro2Link -chrom=chr7
         # 146029338 bases of 154952424 (94.241%) in intersection
     featureBits hg18 chainRBestPanTro2Link -chrom=chr7
         # 145403551 bases of 154952424 (93.838%) in intersection
     featureBits hg18 chainPanTro2Link -chrom=chr7
         # 148298785 bases of 154952424 (95.706%) in intersection
 
     # some additional processing of nets
     # add classification info
     netClass -noAr rbest.noClass.net hg18 panTro2 hg18.panTro2.rbest.net
     gzip hg18.panTro2.rbest.net
 
     # link to downloads
     md5sum *.gz > md5sum.txt
     cd ..
     md5sum axtNet/*.gz >> axtChain/md5sum.txt
     md5sum axtRBestNet/*.gz >> axtChain/md5sum.txt
     cd  /usr/local/apache/htdocs/goldenPath/hg18/vsPanTro2
     set d = /cluster/data/hg18/bed/blastz.panTro2
     ln -s $d/axtChain/*.gz .
     mkdir axtRBestNet
     ln -s /cluster/data/hg18/bed/blastz.panTro2/axtRBestNet/*.axt.gz axtRBestNet
     # edit README.txt
 
     # cleanup: TODO after release
     cd /cluster/data/hg18/bed/blastz.panTro2/axtChain
     rm -f rbest.noClass.net panTro2ToHg18.swap.chain
     rm -fr swapChain rBestChain rBestNet 
 
 
 ########################
 # BLASTZ MACAQUE rheMac2 (2006-03-14 kate)
 
     ssh pk
     cd /cluster/data/panTro2/bed
     mkdir blastz.rheMac2.2006-03-14
     ln -s blastz.rheMac2.2006-03-14 blastz.rheMac2
     cd blastz.rheMac2
 
     cat << '_EOF_' > DEF
 # chimp vs monkey
 export PATH=/usr/bin:/bin:/usr/local/bin:/cluster/bin/penn:/cluster/bin/x86_64:/cluster/home/angie/schwartzbin:/parasol/bin
 
 BLASTZ=blastz.v7.x86_64
 BLASTZ_ABRIDGE_REPEATS=0
 
 # TARGET: Chimp panTro2
 SEQ1_DIR=/scratch/hg/panTro2/panTro2.2bit
 SEQ1_CHUNK=10000000
 SEQ1_LAP=10000
 SEQ1_LEN=/cluster/bluearc/panTro2/chrom.sizes
 
 # QUERY: Macaque rheMac2 - single chunk big enough to run entire genome
 SEQ2_DIR=/scratch/hg/rheMac2/rheMac2.2bit
 SEQ2_LEN=/scratch/hg/rheMac2/chrom.sizes
 SEQ2_CHUNK=3000000000
 SEQ2_LAP=0
 
 BASE=/cluster/data/panTro2/bed/blastz.rheMac2.2006-03-14
 TMPDIR=/scratch/tmp
 '_EOF_'
 
     # NOTE: probably should have used chainMinScore=3000
     # << happy emacs
     /cluster/bin/scripts/doBlastzChainNet.pl -verbose=2 \
         -chainMinScore=5000 -chainLinearGap=medium -bigClusterHub=pk \
         -blastzOutRoot /san/sanvol1/scratch/panTro2/blastz.rheMac2 \
         `pwd`/DEF >&! blastz.out &
 
     # Started Mar 14 18:40
     # Done Mar 15 18:50
 
     # Copy MAF files to san for use in multiple alignment
     set dir = /san/sanvol1/scratch/panTro2/mafNet
     mkdir -p $dir
     cp -rp mafNet $dir/rheMac2
 
     ssh hgwdev featureBits panTro2 chainRheMac2Link
         # crashes out-of-mem on hgwdev (try on kolossus)
 
     # SWAP ALIGNMENTS TO PANTRO2 (2006-03-15 kate)
 
     /cluster/bin/scripts/doBlastzChainNet.pl -swap \
         /cluster/data/panTro2/bed/blastz.rheMac2/DEF >&! swap.log &
 
         # failed due to missing file hgwdev:/scratch/hg/rheMac2/chrom.sizes
         # so I copied the file and ran the "loadUp" script
     ssh hgwdev
     cd /cluster/data/rheMac2/bed/blastz.panTro2.swap/axtChain
     nice loadUp.csh >&! loadUp.log &
 
     # Redo of chr9 (2006-03-28 kate)
     # Create blastz run dir for just chr9, and run manually on pk.
     ssh pk
     cd /cluster/data/panTro2/bed/blastz.rheMac2
     mkdir run.blastz.chr9
     cd run.blastz.chr9
     cp ../run.blastz/*.csh .
     grep chr9 ../run.blastz/jobList | grep -v random > jobList
     ln -s ../run.blastz/qParts .
     para create jobList
         # 28 jobs
     para try
     para push
     para time > run.time
 
     # merge chr9 psl's on small cluster ("cat" step in automated script)
     ssh kki
     cd /cluster/data/panTro2/bed/blastz.rheMac2
     mkdir run.cat.chr9
     cd run.cat.chr9
     cp ../run.cat/cat.csh .
     grep chr9 ../run.cat/jobList > jobList
     para create jobList
         # 14 jobs
     para push
     para time > run.time
 
     cd ../axtChain
     mkdir run.chr9
     cd run.chr9
     cp ../run/chain.csh .
     grep chr9 ../run/jobList > jobList
     mkdir chain
     para create jobList
         # 1 job
     para push
     para time > run.time
 
     # replace chr9 chain in all.chain
     ssh hgwdev64 # kolossus busy
     cd /cluster/data/panTro2/bed/blastz.rheMac2/axtChain/run
     mkdir chain
     set all = panTro2.rheMac2.all.chain.gz
     chainSplit chain/ ../$all
     rm chain/chr9.chain
     cp ../run.chr9/chain/*.chain chain
     chainMergeSort chain/*.chain | nice gzip -c > ../$all
     rm -fr chain/*
     cd ..
     mkdir chain
     chainSplit chain/ $all
         # GOT HERE
 
     # automate the rest
     ssh pk
     cd /cluster/data/panTro2/bed/blastz.rheMac2
     rm -fr axtNet mafNet
     ssh hgwdev "rm -fr /usr/local/apache/htdocs/goldenPath/panTro2/vsRheMac2"
     /cluster/bin/scripts/doBlastzChainNet.pl -verbose=2 \
         -bigClusterHub=pk -workhorse=hgwdev64 \
         -continue=net \
         `pwd`/DEF >&! blastz.net.chr9.out &
 
     # gotta remove another file for the automation
     rm axtChain/panTro2.rheMac2.net.gz
     /cluster/bin/scripts/doBlastzChainNet.pl -verbose=2 \
         -bigClusterHub=pk -workhorse=kkr1u00 \
         -continue=load \
         `pwd`/DEF >&! blastz.load.chr9.out &
 
     # missing noClass.net, work around it
     ssh hgwdev
     cd /cluster/data/panTro2/bed/blastz.rheMac2/axtChain
     cat net/*.net > noClass.net
     # create a script from loadUp.csh, containing only netClass and hgLoadNet
     csh loadUp2.csh >&! ../blastz.load2.chr9.out &
 
     /cluster/bin/scripts/doBlastzChainNet.pl -verbose=2 \
         -bigClusterHub=pk -workhorse=kkr1u00 \
         -continue=download \
         `pwd`/DEF >&! blastz.download.chr9.out &
 
     featureBits panTro2 chainRheMac2Link >&! chainRheMac2Link.fb
         # 2452377221 bases of 2909512873 (84.288%) in intersection
 
     # NOTE: had to run this on kolossus as hgwdev ran out of mem.
     # (needed to push chromInfo and chr*_gap tables over first).
 
     # swap alignments after chr9 fix
     ssh pk
     rm -fr /cluster/data/rheMac2/bed/blastz.panTro2.swap
     cd /cluster/data/panTro2/bed/blastz.rheMac2
     /cluster/bin/scripts/doBlastzChainNet.pl -swap \
         -workhorse=kkr1u00 \
        `pwd`/DEF >&! swap.log &
 
 
 ########################
 # BLASTZ DOG canFam2 (2006-03-15 kate)
 
     ssh pk
     cd /cluster/data/panTro2/bed
     mkdir blastz.canFam2.2006-03-15
     ln -s blastz.canFam2.2006-03-15 blastz.canFam2
     cd blastz.canFam2
 
     cat << '_EOF_' > DEF
 # chimp vs dog
 export PATH=/usr/bin:/bin:/usr/local/bin:/cluster/bin/penn:/cluster/bin/x86_64:/cluster/home/angie/schwartzbin:/parasol/bin
 
 BLASTZ=blastz.v7.x86_64
 
 # Specific settings for dog (per Webb email to Brian Raney)
 BLASTZ_ABRIDGE_REPEATS=1
 # NOTE: must use nibs with repeat abridging
 
 # TARGET: Chimp panTro2
 SEQ1_DIR=/scratch/hg/panTro2/nib
 SEQ1_SMSK=/cluster/bluearc/panTro2/linSpecRep/notInOthers
 SEQ1_LEN=/scratch/hg/panTro2/chrom.sizes
 SEQ1_IN_CONTIGS=0 
 SEQ1_CHUNK=10000000
 SEQ1_LAP=10000
 
 # QUERY: Dog CanFam2 - single chunk big enough to run entire genome
 SEQ2_DIR=/scratch/hg/canFam2/nib
 SEQ2_LEN=/cluster/bluearc/canFam2/chrom.sizes
 SEQ2_SMSK=/san/sanvol1/scratch/canFam2/linSpecRep.notInHuman
 SEQ2_IN_CONTIGS=0
 SEQ2_CHUNK=200000000
 SEQ2_LAP=0
 
 BASE=/cluster/data/panTro2/bed/blastz.canFam2.2006-03-15
 TMPDIR=/scratch/tmp
 '_EOF_'
     # << happy emacs
 
     /cluster/bin/scripts/doBlastzChainNet.pl -verbose=2 \
 	-bigClusterHub=pk -chainMinScore=3000 -chainLinearGap=medium \
         -blastzOutRoot /san/sanvol1/scratch/panTro2/blastz.canFam2 \
 	`pwd`/DEF >&! blastz.out &
     #	Started 2006-03-15 15:40
 
     # Failed due to missing chimp nibs on small cluster /scratch/hg
     # copied them over manually, then continued, running from
     # fileserver, as pk is down
 
     /cluster/bin/scripts/doBlastzChainNet.pl -verbose=2 \
 	-bigClusterHub=kk -chainMinScore=3000 -chainLinearGap=medium \
         -continue=chainRun \
         -blastzOutRoot /san/sanvol1/scratch/panTro2/blastz.canFam2 \
 	`pwd`/DEF >&! blastz.2.out &
     #	Started 2006-03-16 17:50
 
     # Failed due to missing chrom.sizes on hgwdev:/scratch/hg/panTro2
     # Copied manually then continued
 
     /cluster/bin/scripts/doBlastzChainNet.pl -verbose=2 \
 	-bigClusterHub=pk -continue=load \
 	`pwd`/DEF >&! blastz.3.out &
 
     # Redo of chr9 (2006-03-28 kate)
     # Create blastz run dir for just chr9, and run manually on pk.
     ssh pk
     cd /cluster/data/panTro2/bed/blastz.canFam2
     mkdir run.blastz.chr9
     cd run.blastz.chr9
     cp ../run.blastz/*.csh .
     grep panTro2/nib/chr9.nib ../run.blastz/jobList > jobList
     para create jobList
         # 14 jobs
     para try
     para push
     para time > run.time
 
     # merge chr9 psl's on small cluster ("cat" step in automated script)
     ssh kki
     cd /cluster/data/panTro2/bed/blastz.canFam2
     mkdir run.cat.chr9
     cd run.cat.chr9
     cp ../run.cat/cat.csh .
     grep chr9 ../run.cat/jobList > jobList
     para create jobList
         # 15 jobs
     para push
     para time > run.time
 
     cd ../axtChain
     mkdir run.chr9
     cd run.chr9
     cp ../run/chain.csh .
     grep chr9 ../run/jobList | grep -v random > jobList
     mkdir chain
     para create jobList
         # 1 job
     para push
     para time > run.time
 
     # replace chr9 chain in all.chain
     ssh hgwdev64 # kolossus busy
     cd /cluster/data/panTro2/bed/blastz.canFam2/axtChain/run
     mkdir chain
     set all = panTro2.canFam2.all.chain.gz
     chainSplit chain/ ../$all
     rm chain/chr9.chain
     cp ../run.chr9/chain/*.chain chain
     chainMergeSort chain/*.chain | nice gzip -c > ../$all
     rm chain/*
     chainSplit chain/ ../$all
     mv chain ..
 
     # automate the rest
     ssh pk
     cd /cluster/data/panTro2/bed/blastz.canFam2
     rm -fr axtNet mafNet
     rm axtChain/{noClass.net,panTro2.canFam2.net.gz}
     ssh hgwdev "rm -fr /usr/local/apache/htdocs/goldenPath/panTro2/vsCanFam2"
     /cluster/bin/scripts/doBlastzChainNet.pl -verbose=2 \
         -bigClusterHub=pk -workhorse=kkr8u00 \
         -continue=net \
         `pwd`/DEF >&! blastz.net.chr9.out &
     # done
 
     featureBits panTro2 chainCanFam2Link > chainCanFam2Link.fb
         # 1516576565 bases of 2909512873 (52.125%) in intersection
 
 
 
 #########################################################################
 # BLASTZ galGal2 (2006-03-18 kate)
 
     ssh pk
     cd /cluster/data/panTro2/bed/
     mkdir blastz.galGal2.2006-03-18
     ln -s blastz.galGal2.2006-03-18 blastz.galGal2
     cd blastz.galGal2
 
     cat << '_EOF_' > DEF
 # chimp vs chicken
 export PATH=/usr/bin:/bin:/usr/local/bin:/cluster/bin/penn:/cluster/bin/i386:/cluster/home/angie/schwartzbin:/parasol/bin
     
 BLASTZ=blastz.v7
 BLASTZ_ABRIDGE_REPEATS=1 
     
 # TARGET: Chimp panTro2
 SEQ1_DIR=/scratch/hg/panTro2/nib
 SEQ1_SMSK=/cluster/bluearc/panTro2/linSpecRep/notInOthers
 SEQ1_LEN=/scratch/hg/panTro2/chrom.sizes
 SEQ1_CHUNK=50000000
 SEQ1_LAP=10000
     
 # QUERY: Chicken galGal2 - single chunk big enough for whole chroms at once
 SEQ2_DIR=/scratch/hg/galGal2/nib
 SEQ2_LEN=/scratch/hg/galGal2/chrom.sizes
 SEQ2_SMSK=/scratch/hg/galGal2/linSpecRep
 SEQ2_IN_CONTIGS=0
 SEQ2_CHUNK=200000000
 SEQ2_LAP=0
 
 BASE=/cluster/data/panTro2/bed/blastz.galGal2.2006-03-18
 TMPDIR=/scratch/tmp
 '_EOF_'
     #   << happy emacs
 
     /cluster/bin/scripts/doBlastzChainNet.pl -verbose=2 \
         -bigClusterHub=pk -chainMinScore=5000 -chainLinearGap=loose \
         -blastzOutRoot /san/sanvol1/scratch/panTro2/blastz.galGal2 \
         `pwd`/DEF >&! blastz.out &
 
     # error in net step -- can't find all.chain for some reason
     # try restarting
     /cluster/bin/scripts/doBlastzChainNet.pl -verbose=2 \
         -bigClusterHub=pk -chainMinScore=5000 -chainLinearGap=loose \
         -continue=net \
         `pwd`/DEF >&! blastz.2.out &
 
     # Redo of chr9 (2006-03-28 kate)
     # Create blastz run dir for just chr9, and run manually on pk.
     ssh pk
     cd /cluster/data/panTro2/bed/blastz.galGal2
     mkdir run.blastz.chr9
     cd run.blastz.chr9
     cp ../run.blastz/*.csh .
     grep panTro2/nib/chr9.nib ../run.blastz/jobList > jobList
         # 162 jobs
     para create jobList
     para push
     para time > run.time
 
     # merge chr9 psl's on small cluster ("cat" step in automated script)
     ssh kki
     cd /cluster/data/panTro2/bed/blastz.galGal2
     mkdir run.cat.chr9
     cd run.cat.chr9
     cp ../run.cat/cat.csh .
     grep chr9 ../run.cat/jobList > jobList
     para create jobList
         # 15 jobs
     para push
     para time > run.time
 
     cd ../axtChain
     mkdir run.chr9
     cd run.chr9
     cp ../run/chain.csh .
     grep chr9 ../run/jobList | grep -v random > jobList
     mkdir chain
     para create jobList
         # 1 job
     para push
         # GOT HERE
     para time > run.time
 
     # replace chr9 chain in all.chain
     ssh kkstore01 # kolossus busy
     cd /cluster/data/panTro2/bed/blastz.galGal2/axtChain/run
     mkdir chain
     set all = panTro2.galGal2.all.chain.gz
     nice chainSplit chain/ ../$all
     rm chain/chr9.chain
     cp ../run.chr9/chain/*.chain chain
     nice chainMergeSort chain/*.chain | nice gzip -c > ../$all
     rm chain/*
     cd ..
     mkdir chain
     nice chainSplit chain/ $all
 
     # automate the rest
     ssh pk
     cd /cluster/data/panTro2/bed/blastz.galGal2
     rm -fr axtNet mafNet
     rm axtChain/{noClass.net,panTro2.galGal2.net.gz}
     ssh hgwdev "rm -fr /usr/local/apache/htdocs/goldenPath/panTro2/vsGalGal2"
     /cluster/bin/scripts/doBlastzChainNet.pl -verbose=2 \
         -bigClusterHub=pk -workhorse=kkr7u00 \
         -continue=net \
         `pwd`/DEF >&! blastz.net.chr9.out &
     
     featureBits panTro2 chainGalGal2Link > chainGalGal2Link.fb
         # 79103883 bases of 2909512873 (2.719%) in intersection
 
     # swap alignments to other browser
     /cluster/bin/scripts/doBlastzChainNet.pl -verbose=2 \
         -swap `pwd`/DEF >&! swap.out &
 
 
 ############################################################################
 # GENSCAN PREDICTIONS (2006-03-18 kate)
 
     ssh kkstore04
     cd /cluster/data/panTro2
     mkdir bed/genscan
 
     # setup to break hardmasked chrom sequence into 5MB chunks at 
     # unbridged contigs/gaps
     set sdir = /san/sanvol1/scratch/panTro2/splitHard
     foreach d (`cat chrom.lst | grep -v M`)
         mkdir -p $sdir/$d
         foreach agp ($d/chr*.agp)
             set c = $agp:t:r
             echo $c
             cp $agp $sdir/$d
             cp $d/$c.fa.masked $sdir/$d
         end
     end
 
     ssh pk
     cd /cluster/data/panTro2
     cd bed/genscan
     # Make 3 subdirectories for genscan output files
     mkdir gtf pep subopt
 
     # finish split
 cat > bed/genscan/doSplit.csh << 'EOF'
     set sdir = /san/sanvol1/scratch/panTro2/splitHard
     cd /san/sanvol1/scratch/panTro2/splitHard
     foreach agp (*/chr*.agp)
         set d = $agp:h
         set c = $agp:t:r
         set fa = $d/$c.fa.masked
         echo splitting $agp $fa
         nice splitFaIntoContigs $agp $fa . -nSize=5000000
     end
 'EOF'
     csh bed/genscan/doSplit.csh >&! bed/genscan/doSplit.log &
 
     # add chrM
         
     # create list of chunks for cluster run
     cd bed/genscan
     rm -f genome.list
     set sdir = /san/sanvol1/scratch/panTro2/splitHard
 
     # slip in chrM
     mkdir -p $sdir/M/chrM_1
     cp ../../M/chrM.fa.masked $sdir/M/chrM_1/chrM_1.fa
     mkdir $sdir/M/lift
     awk '/chrM/ {printf("%d\t%s\t%d\t%s\t%d\n", 0, "M/chrM", $2, "chrM", $2)}' \
         /cluster/data/panTro2/chrom.sizes > $sdir/M/lift/ordered.lft
     
     cd $sdir
     foreach f (*/chr*_*/*.fa) 
       egrep '[ACGT]' $f > /dev/null
       if ($status == 0) echo $sdir/$f >> \
                 /cluster/data/panTro2/bed/genscan/genome.list
     end
 
     # retrieve executable
     ssh hgwdev
     cd /cluster/data/panTro2/bed/genscan
     cvs co hg3rdParty/genscanlinux
 
     ssh kki
     cd /cluster/data/panTro2/bed/genscan
 
     # Create template file, gsub, for gensub2.  For example (3-line file):
     cat << '_EOF_' > gsub
 #LOOP
 /cluster/bin/x86_64/gsBig {check in line+ $(path1)} {check out line gtf/$(root1).gtf} -trans={check out line pep/$(root1).pep} -subopt={check out line subopt/$(root1).bed} -exe=hg3rdParty/genscanlinux/genscan -par=hg3rdParty/genscanlinux/HumanIso.smat -tmp=/tmp -window=2400000
 #ENDLOOP
 '_EOF_'
     # << this line makes emacs coloring happy
     gensub2 genome.list single gsub jobList
     para create jobList
         # 646 jobs written to batch
     para try
     para check
     para push
 # CPU time in finished jobs:     444979s    7416.32m   123.61h    5.15d  0.014 y
 # IO & Wait Time:                  1898s      31.63m     0.53h    0.02d  0.000 y
 # Time in running jobs:          150791s    2513.18m    41.89h    1.75d  0.005 y
 # Average job time:                 693s      11.55m     0.19h    0.01d
 # Longest running job:           150791s    2513.18m    41.89h    1.75d
 # Longest finished job:           86693s    1444.88m    24.08h    1.00d
 # Submission to last job:         87402s    1456.70m    24.28h    1.01d
 
    # Excessively long on chunks with centromere not at start or end of chunk
    # Next time, try to break these.
 
     # Convert these to chromosome level files
     ssh kkstore04
     cd /cluster/data/panTro2/bed/genscan
     set sdir = /san/sanvol1/scratch/panTro2/splitHard
     cat $sdir/*/lift/*.lft > $sdir/liftAll.lft
     liftUp genscan.gtf $sdir/liftAll.lft warn gtf/*.gtf
     liftUp genscanSubopt.bed $sdir/liftAll.lft warn subopt/*.bed
     cat pep/*.pep > genscan.pep
 
     # Load into the database
     ssh hgwdev
     cd /cluster/data/panTro2/bed/genscan
     ldHgGene panTro2 genscan genscan.gtf
         # 42844 gene predictions
     hgsql panTro1 -Ne "select count(*) from genscan"
         # 41471
     hgPepPred panTro2 generic genscanPep genscan.pep
     hgLoadBed -strict panTro2 genscanSubopt genscanSubopt.bed
 
     featureBits panTro2 genscan
         # 53758249 bases of 2909485072 (1.848%) in intersection
     featureBits panTro1 genscan
         # 53757565 bases of 2909485072 (1.848%) in intersection
 
 QA NOTE:
 
 Ran mytouch for genscanPep table and changed the Update_time to 2006-07-19
 18:00:00, because runJoiner was complaining. 
 
 #######################################################################
 #  OPOSSUM BLASTZ - (WORKING - 2006-03-20 - Hiram)
     ssh pk
     mkdir /cluster/data/panTro2/bed/blastzMonDom4.2006-03-20
     cd /cluster/data/panTro2/bed
     ln -s blastzMonDom4.2006-03-20 blastz.monDom4
     cd /cluster/data/panTro2/bed/blastzMonDom4.2006-03-20
 
     cat << '_EOF_' > DEF
 # chimp vs. opossum
 export PATH=/usr/bin:/bin:/usr/local/bin:/cluster/bin/penn:/cluster/bin/x86_64:/parasol/bin
 
 BLASTZ=blastz.v7.x86_64
 
 # settings for more distant organism alignments
 BLASTZ_H=2000
 BLASTZ_Y=3400
 BLASTZ_L=10000
 BLASTZ_K=2200
 BLASTZ_M=50
 BLASTZ_Q=/cluster/data/blastz/HoxD55.q
 BLASTZ_ABRIDGE_REPEATS=0
 
 # TARGET: Chimp (panTro2)
 SEQ1_DIR=/scratch/hg/panTro2/nib
 SEQ1_LEN=/scratch/hg/panTro2/chrom.sizes
 SEQ1_CHUNK=10000000
 SEQ1_LAP=10000
 
 # QUERY: Opossum monDom4
 SEQ2_DIR=/san/sanvol1/scratch/monDom4/monDom4.2bit
 SEQ2_LEN=/san/sanvol1/scratch/monDom4/chrom.sizes
 SEQ2_IN_CONTIGS=0
 SEQ2_CHUNK=30000000
 SEQ2_LAP=0
 
 BASE=/cluster/data/panTro2/bed/blastzMonDom4.2006-03-20
 TMPDIR=/scratch/tmp
 '_EOF_'
     #	<< happy emacs
 
     time /cluster/bin/scripts/doBlastzChainNet.pl -verbose=2 \
 	-bigClusterHub=pk -chainMinScore=5000 -chainLinearGap=loose \
 	`pwd`/DEF > blastz.out 2>&1 &
     #	running 2006-03-20 11:30
 
 
 # SWAP CHAINS/NET RN4 (DONE 3/27/06 angie)
     # This run used the updated chr9.nib.
     ssh kkstore01
     mkdir /cluster/data/panTro2/bed/blastz.rn4.swap
     cd /cluster/data/panTro2/bed/blastz.rn4.swap
     doBlastzChainNet.pl -swap /cluster/data/rn4/bed/blastz.panTro2/DEF \
       -workhorse kkr8u00 >& do.log & tail -f do.log
     ln -s blastz.rn4.swap /cluster/data/panTro2/bed/blastz.rn4
 
 
 ############################################################################
 ### swap of Mm8 alignments - (DONE - 2006-03-30 - Hiram)
     ssh pk
     mv /cluster/data/panTro2/bed/blastz.mm8.swap \
 	/cluster/data/panTro2/bed/blastz.mm8.swap.2006-03-21
     mkdir /cluster/data/panTro2/bed/blastz.mm8.swap
     cd /cluster/data/panTro2/bed/blastz.mm8.swap
 
     time /cluster/bin/scripts/doBlastzChainNet.pl -verbose=2 \
 	-swap -bigClusterHub=pk -chainMinScore=3000 -chainLinearGap=medium \
 	/cluster/data/mm8/bed/blastzPanTro2.2006-03-28/DEF \
 	> blastz.out 2>&1 &
     #  XXX started 2006-03-30
 
     time nice -n +19 featureBits panTro2 chainMm8Link \
 	> fb.panTro2.chainMm8Link
     # XXXX first time was:
     #	986978326 bases of 2909512873 (33.922%) in intersection
 
 
 ############################################################################
 # SPLIT SEQUENCE FOR LIFTOVER CHAINS FROM PANTRO1  (2006-03-13 kate)
 
     ssh kkr1u00
     cd /cluster/data/panTro2/bed
     mkdir liftOver
     cd liftOver
     makeLoChain-split panTro2 /cluster/data/panTro2/nib >&! split.log &
 
 ########################
 # BLASTZ HUMAN HG17 (2006-06-17 kate)
 # Do not use lineage-specific repeats, on recommendation of Arian Smit
 # Do use latest human-chimp scoring matrix and gap penalties from Webb Miller
 # via rico (Richard Burhans) at PSU
 
 # NOTE: these alignments have been extensively reviewed by Kateryna 
 # Makerova, Richard Burnhams and Webb Miller at Penn state to 
 # verify parameters are suitable.  The same parameters were
 # used for the HG18 alignments, below.
 
     cat /cluster/data/penn/human_chimp.v2.q
 #    A    C    G    T
 #    90 -330 -236 -356
 #  -330  100 -318 -236
 #  -236 -318  100 -330
 #  -356 -236 -330   90
 # O = 600 E = 150
 
     ssh pk
     cd /cluster/data/panTro2/bed
     mkdir blastz.hg17.2006-06-17-05
     ln -s blastz.hg17.2006-06-17-05 blastz.hg17
     cd blastz.hg17
 
     ssh kkr1u00
     cp -p /cluster/data/hg17/hg17.2bit /scratch/hg/hg17
     ~kent/bin/iSync
        # failed for some reason, so I copied the files manually to
        # kkr*u00
 
     cat << '_EOF_' > DEF
 # chimp vs human
 # Specially tuned blastz parameters from Webb Miller
 
 export PATH=/usr/bin:/bin:/usr/local/bin:/cluster/bin/penn:/cluster/bin/x86_64:/cluster/home/angie/schwartzbin:/parasol/bin
 
 BLASTZ=blastz.v7.x86_64
 BLASTZ_ABRIDGE_REPEATS=0
 BLASTZ_O=600
 BLASTZ_E=150
 BLASTZ_Y=15000
 BLASTZ_T=2
 BLASTZ_K=4500
 BLASTZ_Q=/cluster/data/blastz/human_chimp.v2.q
 
 # TARGET: Chimmp panTro2
 SEQ1_DIR=/scratch/hg/panTro2/panTro2.2bit
 SEQ1_CHUNK=10000000
 SEQ1_LAP=10000
 SEQ1_LEN=/cluster/bluearc/panTro2/chrom.sizes
 
 # QUERY: Human hg17 - single chunk big enough to run entire genome
 SEQ2_DIR=/scratch/hg/hg17/hg17.2bit
 SEQ2_LEN=/cluster/bluearc/hg17/chrom.sizes
 SEQ2_CHUNK=3000000000
 SEQ2_LAP=0
 
 BASE=/cluster/data/panTro2/bed/blastz.hg17.2006-06-17-05
 TMPDIR=/scratch/tmp
 '_EOF_'
 
 # << happy emacs
 /cluster/bin/scripts/doBlastzChainNet.pl -verbose=2 \
     -chainMinScore=5000 -chainLinearGap=medium -bigClusterHub=pk \
     -blastzOutRoot /san/sanvol1/scratch/panTro2/blastz.hg17 \
     `pwd`/DEF >&! blastz.out &
 
 # failed due to badly formatted scoring matrix file, which
 # fortuitously pointed out the need for a trailing line
 # containing gap penalties in the scoring matrix file.
 # If it's missing (as in all of our current .q files),
 # axtChain uses defaults, not the penalties used for the
 # blastz run (although it could now get these from the axt files).
 # I'm adding this to the matrix file, also will change library
 # (axt.c) to add the gap values used by axtChain to the header comment
 
 /cluster/bin/scripts/doBlastzChainNet.pl -verbose=2 \
     -chainMinScore=5000 -chainLinearGap=medium -bigClusterHub=pk \
     -blastzOutRoot /san/sanvol1/scratch/panTro2/blastz.hg17 \
     -continue=chainRun \
     `pwd`/DEF >&! blastz.2.out &
 
 # failed due to ssh key failure on hgwdev
 
 /cluster/bin/scripts/doBlastzChainNet.pl -verbose=2 \
     -chainMinScore=5000 -chainLinearGap=medium -bigClusterHub=pk \
     -blastzOutRoot /san/sanvol1/scratch/panTro2/blastz.hg17 \
     -continue=load \
     `pwd`/DEF >&! blastz.3.out &
 
     nice featureBits panTro2 chainHg17Link
         # 2759308860 bases of 2909512873 (94.837%) in intersection
     nice featureBits panTro2 chainHg17OldLink
         # from alignments with default matrix
         #2782705676 bases of 2909512873 (95.642%) in intersection
 
     hgsql panTro2 -e "select count(*) from chr7_chainHg17"
         # 75376
     hgsql panTro2 -e "select count(*) from chr7_chainHg17Old"
         # 241200
 
     # swap alignments to hg17 reference
     /cluster/bin/scripts/doBlastzChainNet.pl -verbose=2 \
         -swap -bigClusterHub=pk \
         `pwd`/DEF >&! swap.out &
 
     nice featureBits hg17 chainPanTro2Link
         # 2719566113 bases of 2866216770 (94.883%) in intersection
     nice featureBits hg17 chainPanTro2LinkOld
         # 2739413164 bases of 2866216770 (95.576%) in intersection
 
     hgsql panTro2 -e "select count(*) from chr7_chainPanTro2"
         # 105398
     hgsql panTro2 -e "select count(*) from chr7_chainPanTro2Old"
         # 333069
 
     nice featureBits panTro2 chainHg17Link | tee fb.panTro2.chainHg17Link
         # 2759308860 bases of 2823407242 (97.730%) in intersection
 
 
 ########################
 # BLASTZ HUMAN HG18 (2006-07-9 kate)
 # Do not use lineage-specific repeats, on recommendation of Arian Smit
 # Do use latest human-chimp scoring matrix and gap penalties from Webb Miller
 # via rico (Richard Burhans) at PSU
 
     cat /cluster/data/penn/human_chimp.v2.q
 #    A    C    G    T
 #    90 -330 -236 -356
 #  -330  100 -318 -236
 #  -236 -318  100 -330
 #  -356 -236 -330   90
 # O = 600 E = 150
 
     ssh pk
     cd /cluster/data/panTro2/bed
     mkdir blastz.hg18.2006-07-09
     ln -s blastz.hg18.2006-07-09 blastz.hg18
     cd blastz.hg18
 
     cat << '_EOF_' > DEF
 # chimp vs human
 # Specially tuned blastz parameters from Webb Miller
 
 export PATH=/usr/bin:/bin:/usr/local/bin:/cluster/bin/penn:/cluster/bin/x86_64:/cluster/home/angie/schwartzbin:/parasol/bin
 
 BLASTZ=blastz.v7.x86_64
 BLASTZ_ABRIDGE_REPEATS=0
 BLASTZ_O=600
 BLASTZ_E=150
 BLASTZ_Y=15000
 BLASTZ_T=2
 BLASTZ_K=4500
 BLASTZ_Q=/cluster/data/blastz/human_chimp.v2.q
 
 # TARGET: Chimmp panTro2
 SEQ1_DIR=/scratch/hg/panTro2/panTro2.2bit
 SEQ1_CHUNK=10000000
 SEQ1_LAP=10000
 SEQ1_LEN=/cluster/bluearc/panTro2/chrom.sizes
 
 # QUERY: Human hg18 - single chunk big enough to run entire genome
 SEQ2_DIR=/scratch/hg/hg18/hg18.2bit
 SEQ2_LEN=/cluster/bluearc/hg18/chrom.sizes
 SEQ2_CHUNK=3000000000
 SEQ2_LAP=0
 
 BASE=/cluster/data/panTro2/bed/blastz.hg18.2006-07-09
 TMPDIR=/scratch/tmp
 '_EOF_'
 
 # << happy emacs
 /cluster/bin/scripts/doBlastzChainNet.pl -verbose=2 \
     -chainMinScore=5000 -chainLinearGap=medium -bigClusterHub=pk \
     -blastzOutRoot /san/sanvol1/scratch/panTro2/blastz.hg18 \
     `pwd`/DEF >&! blastz.out &
 
 # Mini-cluster nodes were dead, so restart required after admins fixed.
 
 /cluster/bin/scripts/doBlastzChainNet.pl -verbose=2 \
     -chainMinScore=5000 -chainLinearGap=medium -bigClusterHub=pk \
     -blastzOutRoot /san/sanvol1/scratch/panTro2/blastz.hg18 \
     -continue=cat \
     `pwd`/DEF >&! blastz.2.out &
 
     hgsql panTro2 -Nse "select count(*) from chr7_chainHg18"
         # 76122
     hgsql panTro2 -Nse "select count(*) from chr7_chainHg17"
         # 75376
 
     featureBits panTro2 -chrom=chr7 chr7_chainHg18Link
         # 146031865 bases of 151077879 (96.660%) in intersection
     featureBits panTro2 -chrom=chr7 chr7_chainHg17Link
         # 145990849 bases of 151077879 (96.633%) in intersection
 
     nice featureBits panTro2 chainHg18Link >& fb.panTro2.chainHg18Link &
     cat  fb.panTro2.chainHg18Link
         # 2760856313 bases of 2823407242 (97.785%) in intersection
     cat ../blastz.hg17/fb.panTro2.chainHg17Link
         # 2759308860 bases of 2823407242 (97.730%) in intersection
 
     # SWAP ALIGNMENTS (2006-07-11 kate)
     ssh kkstore02
     cd /cluster/data/hg18/bed
     mkdir blastz.panTro2.swap
     cd blastz.panTro2
     /cluster/bin/scripts/doBlastzChainNet.pl -swap \
         /cluster/data/panTro2/bed/blastz.hg18/DEF >&! swap.log &
 
     # failed during liftover chain copy (netChains.csh)
 
     /cluster/bin/scripts/doBlastzChainNet.pl -swap -continue=download\
         /cluster/data/panTro2/bed/blastz.hg18/DEF >&! swap.download.log &
 
 
 #########################################################################
 # BLASTZ danRer3 (2006-04-07 kate)
 # Includes both randoms
 
     ssh pk
     mkdir /cluster/data/panTro2/bed/blastz.danRer3.2006-04-07
     cd /cluster/data/panTro2/bed
     ln -s blastz.danRer3.2006-04-07 blastz.danRer3
     cd blastz.danRer3
 
 cat << 'EOF' > DEF
 # human target, zebrafish query
 export PATH=/usr/bin:/bin:/usr/local/bin:/cluster/bin/penn:/cluster/bin/x86_64:/cluster/home/angie/schwartzbin:/parasol/bin
 
 BLASTZ=blastz.v7.x86_64
 BLASTZ_ABRIDGE_REPEATS=1
 
 # use parameters suggested for human-fish evolutionary distance
 # recommended in doBlastzChainNet.pl help
 # (previously used for  hg16-fr1, danrer1-mm5)
 BLASTZ_H=2000
 BLASTZ_Y=3400
 BLASTZ_L=6000
 BLASTZ_K=2200
 BLASTZ_Q=/san/sanvol1/scratch/blastz/HoxD55.q
 
 # TARGET: Human panTro2
 SEQ1_DIR=/scratch/hg/panTro2/nib
 SEQ1_SMSK=/cluster/bluearc/panTro2/linSpecRep/notInZebrafish
 SEQ1_CHUNK=10000000
 SEQ1_LAP=10000
 SEQ1_LEN=/scratch/hg/panTro2/chrom.sizes
 
 # QUERY: zebrafish danRer3
 # Use all chroms, including both randoms (chrUn and chrNA)
 SEQ2_DIR=/san/sanvol1/scratch/danRer3/nib
 SEQ2_SMSK=/san/sanvol1/scratch/danRer3/linSpecRep.notInOthers
 SEQ2_LEN=/cluster/bluearc/danRer3/chrom.sizes
 SEQ2_CHUNK=30000000
 SEQ2_LAP=1000
 
 BASE=/cluster/data/panTro2/bed/blastz.danRer3.2006-04-07
 TMPDIR=/scratch/tmp
 'EOF'
 # << happy emacs
 
 /cluster/bin/scripts/doBlastzChainNet.pl -verbose=2 \
     -bigClusterHub=pk -chainMinScore=5000 -chainLinearGap=loose  \
     -blastzOutRoot /san/sanvol1/scratch/panTro2/blastz.danRer3 \
     `pwd`/DEF >& blastz.out &
 
 # failed on 'cat.run' step, for no good reason.
 # continuing manually with 'para shove -retries=20'...
 
     ssh kki
     cd /cluster/data/panTro2/bed/blastz.danRer3/run.cat
     para time > run.time
 
     ssh pk
     cd /cluster/data/panTro2/bed/blastz.danRer3
     /cluster/bin/scripts/doBlastzChainNet.pl -verbose=2 \
     -bigClusterHub=pk -chainMinScore=5000 -chainLinearGap=loose  \
     -continue=chainMerge \
     `pwd`/DEF >& blastz.chainMerge.out &
 
 # failed due to bad NFS mount on one of the mini-cluster nodes.
 # finished manually, and created run,time
 
 #########################################################################
 # BLASTZ danRer4 (2006-07-10 kate)
 # Includes both randoms
 
     ssh pk
     mkdir /cluster/data/panTro2/bed/blastz.danRer4.2006-07-10
     cd /cluster/data/panTro2/bed
     ln -s blastz.danRer4.2006-07-10 blastz.danRer4
     cd blastz.danRer4
 
     cat << 'EOF' > DEF
 # human target, zebrafish query
 export PATH=/usr/bin:/bin:/usr/local/bin:/cluster/bin/penn:/cluster/bin/x86_64:/cluster/home/angie/schwartzbin:/parasol/bin
 
 BLASTZ=blastz.v7.x86_64
 BLASTZ_ABRIDGE_REPEATS=1
 
 # use parameters suggested for human-fish evolutionary distance
 # recommended in doBlastzChainNet.pl help
 # (previously used for  hg16-fr1, danrer1-mm5)
 BLASTZ_H=2000
 BLASTZ_Y=3400
 BLASTZ_L=6000
 BLASTZ_K=2200
 BLASTZ_Q=/san/sanvol1/scratch/blastz/HoxD55.q
 
 # TARGET: Human panTro2
 SEQ1_DIR=/scratch/hg/panTro2/nib
 SEQ1_SMSK=/cluster/bluearc/panTro2/linSpecRep/notInZebrafish
 SEQ1_CHUNK=10000000
 SEQ1_LAP=10000
 SEQ1_LEN=/scratch/hg/panTro2/chrom.sizes
 
 # QUERY: zebrafish danRer4
 # Use all chroms, including both randoms (chrUn and chrNA)
 SEQ2_DIR=/san/sanvol1/scratch/danRer4/nib
 SEQ2_SMSK=/san/sanvol1/scratch/danRer4/linSpecRep.notInOthers
 SEQ2_LEN=/san/sanvol1/scratch/danRer4/chrom.sizes
 SEQ2_CHUNK=30000000
 SEQ2_LAP=1000
 
 BASE=/cluster/data/panTro2/bed/blastz.danRer4.2006-07-10
 TMPDIR=/scratch/tmp
 'EOF'
 # << happy emacs
 
     /cluster/bin/scripts/doBlastzChainNet.pl -verbose=2 \
         -bigClusterHub=pk -chainMinScore=5000 -chainLinearGap=loose  \
         -blastzOutRoot /san/sanvol1/scratch/panTro2/blastz.danRer4 \
         `pwd`/DEF >& blastz.out &
 
 # 2 jobs crashed on kki, pk was rebooted
 # finish jobs with para push on kki, then continue
 
     /cluster/bin/scripts/doBlastzChainNet.pl -verbose=2 \
         -bigClusterHub=pk -chainMinScore=5000 -chainLinearGap=loose  \
         -blastzOutRoot /san/sanvol1/scratch/panTro2/blastz.danRer4 \
         -continue=chainRun \
         `pwd`/DEF >& blastz.2.out &
 
     # swap alignemnts to other assembly
     /cluster/bin/scripts/doBlastzChainNet.pl -verbose=2 \
         -swap `pwd`/DEF >& swap.out &
 
 #########################################################################
 # LifUunder to panTro1 (2005-07-11 kate)
 
     ssh pk
     cd /cluster/data/panTro2
     doSameSpeciesLiftOver.pl panTro2 panTro1 >& liftUnder.out &
     ln -s panTro1.2005-07-14 blat.panTro1
 
     # stalled out looking for workhorse, restarting
     doSameSpeciesLiftOver.pl panTro2 panTro1 \
         -buildDir=/cluster/data/panTro2/bed/blat.panTro1 \
         -continue=net -workhorse=kolossus >& liftOver.2.out &
 
 #############################################################################
 # 8-WAY MULTIZ ALIGNMENTS (2006-07-15 kate)
 
 # hg18
 # rheMac2
 # mm8
 # canFam2
 # monDom4
 # galGal2
 # danRer4
 
     # copy net mafs to cluster-friendly storage for multiz run
     ssh kkstore04
     cd /cluster/data/panTro2/bed
     set d = /san/sanvol1/scratch/panTro2/mafNet
     mkdir -p $d
     foreach s (hg18 rheMac2 mm8 canFam2 monDom4 galGal2 danRer4)
         echo $s
         cp -Rp blastz.$s/mafNet $d/$s
     end
 
     # make output dir and run dir
     ssh pk
     cd /cluster/data/panTro2/bed
     mkdir -p multiz8way.2006-07-15
     ln -s multiz8way.2006-07-15 multiz8way
     cd multiz8way
 
     cp /san/sanvol1/scratch/mm7/cons/elliotsEncode.mod .
     cat elliotsEncode.mod
 ALPHABET: A C G T 
 ORDER: 0
 SUBST_MOD: REV
 TRAINING_LNL: -988246.132962
 BACKGROUND: 0.295 0.205 0.205 0.295
 RATE_MAT:
   -1.165221    0.315494    0.589884    0.259843 
    0.189778   -0.878194    0.208718    0.479698 
    0.444622    0.261535   -0.885604    0.179447 
    0.234867    0.720815    0.215191   -1.170872 
 TREE: (((((((((((((hg17:0.006690,panTro1:0.007571):0.024272,(colobus_monkey:0.015404,(baboon:0.008258,rheMac2:0.028617):0.008519):0.022120):0.023960,(dusky_titi:0.025662,(owl_monkey:0.012151,marmoset:0.029549):0.008236):0.027158):0.066101,(mouse_lemur:0.059024,galago:0.121375):0.032386):0.017073,((rn3:0.081728,mm7:0.077017):0.229273,oryCun1:0.206767):0.023340):0.023026,(((bosTau2:0.159182,canFam2:0.147731):0.004946,rfbat:0.138877):0.010150,(hedgehog:0.193396,shrew:0.261724):0.054246):0.024354):0.028505,dasNov1:0.149862):0.015994,(loxAfr1:0.104891,echTel1:0.259797):0.040371):0.218400,monDom2:0.371073):0.065268,platypus:0.468116):0.123856,galGal2:0.454691):0.123297,xenTro1:0.782453):0.156067,((tetNig1:0.199381,fr1:0.239894):0.492961,danRer3:0.782561):0.156067);
 
     # create species list and stripped down tree for autoMZ
     sed -n '/TREE: /s///p' elliotsEncode.mod > elliotsEncode.nh
 
     /cluster/bin/phast/tree_doctor elliotsEncode.nh \
      --prune-all-but=hg17,panTro1,rheMac1,mm7,canFam2,monDom2,galGal2,danRer3 \
      --rename="hg17->hg18;panTro1->panTro2;rheMac1->rheMac2;mm7->mm8;monDom2->monDom4;danRer3->danRer4" > 8way.nh
 
     sed 's/[a-z][a-z]*_//g; s/:[0-9\.][0-9\.]*//g; s/;//' 8way.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.lst
 
     /cluster/bin/phast/all_dists 8way.nh > 8way.distances.txt
     grep panTro2 8way.distances.txt | sort -k3,3n | \
         awk '{printf ("%.4f\t%s\t%s\n", $3, $2, $1)}' > distances.txt
     cat distances.txt
 # 0.0143  chimp_panTro2   human_hg18
 # 0.0910  macaque_rheMac2 chimp_panTro2
 # 0.2660  dog_canFam2     chimp_panTro2
 # 0.4686  mouse_mm8       chimp_panTro2
 # 0.7128  monodelphis_monDom4     chimp_panTro2
 # 0.9855  chicken_galGal2 chimp_panTro2
 # 1.7488  zebrafish_danRer4       chimp_panTro2
 
     mkdir -p maf run
     cd run
 
     # stash binaries 
     mkdir penn
     cp -p /cluster/bin/penn/multiz.v11.x86_64/multiz-tba/multiz penn
     cp -p /cluster/bin/penn/multiz.v11.x86_64/multiz-tba/maf_project penn
     cp -p /cluster/bin/penn/multiz.v11.x86_64/multiz-tba/autoMZ penn
 
 cat > autoMultiz.csh << 'EOF'
 #!/bin/csh -ef
     set db = panTro2
     set c = $1
     set maf = $2
     set run = `pwd`
     set tmp = /scratch/tmp/$db/multiz.$c
     set pairs = /san/sanvol1/scratch/$db/mafNet
     rm -fr $tmp
     mkdir -p $tmp
     cp ../{tree.nh,species.lst} $tmp
     pushd $tmp
     foreach s (`cat species.lst`)
         set in = $pairs/$s/$c.maf
         set out = $db.$s.sing.maf
         if ($s == panTro2) then
             continue
         endif
         if (-e $in.gz) then
             zcat $in.gz > $out
         else if (-e $in) then
             cp $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
     popd
     cp $tmp/$c.maf $maf
     rm -fr $tmp
 # << happy emacs
     chmod +x autoMultiz.csh
 
 cat  << 'EOF' > spec
 #LOOP
 ./autoMultiz.csh $(root1) {check out line+ /cluster/data/panTro2/bed/multiz8way.2006-07-15/maf/$(root1).maf}
 #ENDLOOP
 'EOF'
 # << happy emacs
 
     awk '{print $1}' /cluster/data/panTro2/chrom.sizes > chrom.lst
 
     gensub2 chrom.lst single spec jobList
     para create jobList
         # 52 files
     para try
     para check
     para push
     para time > run.time
 # CPU time inn finished jobs:      99662s    1661.03m    27.68h    1.15d  0.003 y
 # IO & Wait Time:                  2009s      33.49m     0.56h    0.02d  0.000 y
 # Average job time:                1955s      32.59m     0.54h    0.02d
 # Longest running job:                0s       0.00m     0.00h    0.00d
 # Longest finished job:            8172s     136.20m     2.27h    0.09d
 # Submission to last job:          8172s     136.20m     2.27h    0.09d
 
     # annotate mafs
     ssh kolossus
     cd /cluster/data/panTro2/bed/multiz8way
     mkdir anno 
     cd anno
     mkdir maf run
     cd run
     rm -f sizes nBeds
 
     foreach i (`cat /cluster/data/panTro2/bed/multiz8way/species.lst`)
         ln -s  /cluster/data/$i/chrom.sizes $i.len
         ln -s  /cluster/data/$i/$i.N.bed $i.bed
         echo $i.bed  >> nBeds
         echo $i.len  >> sizes
     end
 
     echo date > jobs.csh
     foreach i (../../maf/*.maf)
         echo nice mafAddIRows -nBeds=nBeds -sizes=sizes $i /cluster/data/panTro2/panTro2.2bit ../maf/`basename $i` >> jobs.csh
         echo "echo $i" >> jobs.csh
     end 
     echo date >> jobs.csh
 
     # do smaller jobs first
     tac jobs.csh > jobsRev.csh
     mv jobsRev.csh jobs.csh
     
     csh jobs.csh > jobs.log &
 
     # Load into database
     ssh hgwdev
     cd /cluster/data/panTro2/bed/multiz8way/anno/maf
     mkdir -p /gbdb/panTro2/multiz8way/anno/maf
     ln -s /cluster/data/panTro2/bed/multiz8way/anno/maf/*.maf \
         /gbdb/panTro2/multiz8way/anno/maf
 cat > loadMaf.csh << 'EOF'
     nice hgLoadMaf -pathPrefix=/gbdb/panTro2/multiz8way/anno/maf panTro2 multiz8way
     cat *.maf | \
         nice hgLoadMafSummary panTro2 -minSize=30000 -mergeGap=1500 -maxSize=200000  multiz8waySummary stdin
 'EOF'
     # 16414516 mafs, 1281405 summary blocks
 
 #<< happy emacs
     csh loadMaf.csh  >&! loadMaf.log &
 
     # generate .gif of tree
     phyloGif
 
+----------------------------------------------------------------------------
+2009-08-20 markd with guidance from braney
+    ###
+    # mafs not correctly displayed due to incorrectly using  -sizes=sizes
+    # with mafAddIRows
+    ###
+
+    # annotate mafs
+    ssh kolossus
+    cd /cluster/data/panTro2/bed/multiz8way/anno
+    rm -f maf/*.maf run/*.maf
+    cd run
+    # edit jobs.csh to remove -sizes=sizes
+
+    csh jobs.csh >& jobs.log &
+
+>>>>>>>>>>>>>>>>>
+
+    # Load into database
+    ssh hgwdev
+    cd /cluster/data/panTro2/bed/multiz8way/anno/maf
+    mkdir -p /gbdb/panTro2/multiz8way/anno/maf
+    csh loadMaf.csh  >& loadMaf.log &
+
+
 #############################################################################
 # PHASTCONS CONSERVATION (2006-07-18 kate)
 
     # create cluster-friendly copies of chrom fasta files
     ssh kkstore04
     cd /cluster/data/panTro2
     set sdir = /san/sanvol1/scratch/panTro2/chrom
     mkdir -p $sdir
     foreach d (`cat chrom.lst`)
         echo $d
         cp -p $d/chr*.fa $sdir
     end
 
     # split mafs
     ssh pk
     cd /cluster/data/panTro2/bed/multiz8way
     mkdir cons
     cd cons
     mkdir run.split
     cd run.split
     set WINDOWS = /san/sanvol1/scratch/panTro2/multiz8way/cons/ss
     rm -fr $WINDOWS
     mkdir -p $WINDOWS
 
     cat << 'EOF' > doSplit.csh
 #!/bin/csh -ef
     set MAFS = /cluster/data/panTro2/bed/multiz8way/maf
     set WINDOWS = /san/sanvol1/scratch/panTro2/multiz8way/cons/ss
     set sdir = /san/sanvol1/scratch/panTro2/chrom
     cd $WINDOWS
     set c = $1
     echo $c
     rm -fr $c
     mkdir $c
     /cluster/bin/phast/$MACHTYPE/msa_split $MAFS/$c.maf -i MAF \
         -M $sdir/$c.fa -o SS -r $c/$c -w 10000000,0 -I 1000 -B 5000
     echo "Done" >> $c.done
 'EOF'
 # << happy emacs
     chmod +x doSplit.csh
 
     rm -f jobList
     foreach f (../../maf/*.maf) 
         set c = $f:t:r
 	echo "doSplit.csh $c {check out line+ $WINDOWS/$c.done}" >> jobList
     end
     
     para create jobList
         # 52 jobs
     para try
     para check
     para push
 #CPU time in finished jobs:       2426s      40.44m     0.67h    0.03d  0.000 y
 #IO & Wait Time:                   988s      16.46m     0.27h    0.01d  0.000 y
 #Average job time:                  66s       1.09m     0.02h    0.00d
 #Longest running job:                0s       0.00m     0.00h    0.00d
 #Longest finished job:             234s       3.90m     0.07h    0.00d
 #Submission to last job:           346s       5.77m     0.10h    0.00d
 
     # run phastCons
     #	This job is I/O intensive in its output files, thus it is all
     #	working over in /scratch/tmp/
     cd ..
     # create tree model file from ENCODE model file, with tree
     # pruned and renamed for this set of alignments
 
     sed '/TREE: /d' ../elliotsEncode.mod > 8way.mod
     echo "TREE: `cat ../8way.nh`" >> 8way.mod
     cat 8way.mod
 #ALPHABET: A C G T 
 #ORDER: 0
 #SUBST_MOD: REV
 #TRAINING_LNL: -988246.132962
 #BACKGROUND: 0.295 0.205 0.205 0.295
 #RATE_MAT:
 #   -1.165221    0.315494    0.589884    0.259843 
 #    0.189778   -0.878194    0.208718    0.479698 
 #    0.444622    0.261535   -0.885604    0.179447 
 #    0.234867    0.720815    0.215191   -1.170872 
 #TREE:
 #(((((((hg18:0.006690,panTro2:0.007571):0.024272,rheMac2:0.059256):0.107134,mm8:0.329630):0.023026,canFam2:0.187181):0.262899,monDom4:0.371073):0.189124,galGal2:0.454691):0.279364,danRer4:0.938628);
 
     # calculate GC in model file
     awk '$1 == "BACKGROUND:" {print $3 + $4;}' 8way.mod
         #0.41
 
     mkdir run.cons
     cd run.cons
     cat > doPhast.csh << 'EOF'
 #!/bin/csh -fe
 set c = $1
 set f = $2
 set len = $3
 set cov = $4
 set rho = $5
 set tmp = /scratch/tmp/$f
 mkdir -p $tmp
 set san = /san/sanvol1/scratch/panTro2/multiz8way/cons
 cp -p $san/ss/$c/$f.ss ../8way.mod $tmp
 pushd $tmp > /dev/null
 /cluster/bin/phast/$MACHTYPE/phastCons $f.ss 8way.mod \
    --rho $rho --expected-length $len --target-coverage $cov --quiet \
     --not-informative hg18,rheMac2 \
 	--seqname $c --idpref $c --viterbi $f.bed --score > $f.pp
 popd > /dev/null
 mkdir -p $san/pp/$c $san/bed/$c
 sleep 1
 mv $tmp/$f.pp $san/pp/$c
 mv $tmp/$f.bed $san/bed/$c
 rm -fr $tmp
 'EOF'
     # emacs happy
     chmod a+x doPhast.csh
 
     #	root1 == chrom name, file1 == ss file name without .ss suffix
     # Create gsub file
     # start with parameters from previous human 8way
     cat > template << 'EOF'
 #LOOP
 doPhast.csh $(root1) $(file1) 13 .01 .23
 #ENDLOOP
 'EOF'
     #	happy emacs
 
     # Create parasol batch and run it
     pushd /san/sanvol1/scratch/panTro2/multiz8way/cons
     ls -1 ss/chr*/chr*.ss | sed 's/.ss$//' > \
         /cluster/data/panTro2/bed/multiz8way/cons/run.cons/in.list
     popd
 
     gensub2 in.list single template jobList
     grep chr7 jobList > jobList.chr7
     para make jobList.chr7
         # 346 jobs
     # using para make because jobs are so small they crash 
 
     # Note: 1/3 of jobs crashed, but succeeded on retry
     # These jobs are probably too small/too fast.
     # Consider increasing the size of the msa_split even further
     # in future runs.
 
     para time
 #CPU time in finished jobs:      10586s     176.44m     2.94h    0.12d  0.000 y
 #IO & Wait Time:                  5465s      91.08m     1.52h    0.06d  0.000 y
 #Average job time:                  46s       0.77m     0.01h    0.00d
 #Longest running job:                0s       0.00m     0.00h    0.00d
 #Longest finished job:              66s       1.10m     0.02h    0.00d
 #Submission to last job:           581s       9.68m     0.16h    0.01d
 
     # create Most Conserved track
     ssh kolossus
     cd /san/sanvol1/scratch/panTro2/multiz8way/cons
     #	The sed's and the sort get the file names in chrom,start order
     # (Hiram tricks -- split into columns on [.-/] with 
     #    identifying x,y,z, to allow column sorting and
     #    restoring the filename.  Warning: the sort column
     # will depend on how deep you are in the dir
     find ./bed -name "chr*.bed" | \
         sed -e "s/\// x /g" -e "s/\./ y /g" -e "s/-/ z /g" | \
 	sort -k7,7 -k9,9n | \
 	sed -e "s/ x /\//g" -e "s/ y /\./g" -e "s/ z /-/g" | \
 	xargs cat | \
 	awk '{printf "%s\t%d\t%d\tlod=%d\t%s\n", $1, $2, $3, $5, $5;}' | \
 	/cluster/bin/scripts/lodToBedScore /dev/stdin > mostConserved.bed
     #	~ 1 minute
     cp -p mostConserved.bed /cluster/data/panTro2/bed/multiz8way/cons
 
     # load into database
     ssh hgwdev
     cd /cluster/data/panTro2/bed/multiz8way/cons
     hgLoadBed -strict panTro2 phastConsElements8way mostConserved.bed
     # compare with previous tracks
     hgsql panTro2 -e "select count(*) from phastConsElements8way"
         # 1596760 elements 
     hgsql hg18 -e "select count(*) from phastConsElements17way"
         # 1601903
     # Try for 5% overall cov, and 70% CDS cov 
     featureBits panTro2 -enrichment refGene:cds phastConsElements8way
         # refGene:cds 1.020%, phastConsElements8way 4.999%, both 0.731%, cover 71.68%, enrich 14.34x
     # experiments
     featureBits hg18 -chrom=chr7 -enrichment refGene:cds phastConsElements17way
 # refGene:cds 0.883%, phastConsElements17way 4.838%, both 0.621%, cover 70.31%, enrich 14.53x
     featureBits panTro2 -chrom=chr7 -enrichment refGene:cds phastConsElements8way
     # 13 .01 .23
 # refGene:cds 0.861%, phastConsElements8way 4.914%, both 0.626%, cover 72.71%, enrich 14.80x
     # 13 .007 .23
 # refGene:cds 0.861%, phastConsElements8way 4.780%, both 0.626%, cover 72.69%, enrich 15.21x
     # 12 .007 .23
 # refGene:cds 0.861%, phastConsElements8way 4.599%, both 0.623%, cover 72.34%, enrich 15.73x
     # 12 .006 .23
 # refGene:cds 0.861%, phastConsElements8way 4.553%, both 0.623%, cover 72.35%, enrich 15.89x
     # 12 .006 .24
 # refGene:cds 0.861%, phastConsElements8way 4.657%, both 0.630%, cover 73.11%, enrich 15.70x
     # 12 .006 .25
 # refGene:cds 0.861%, phastConsElements8way 4.765%, both 0.636%, cover 73.90%, enrich 15.51x
     # 12 .004 .27
 # refGene:cds 0.861%, phastConsElements8way 4.845%, both 0.647%, cover 75.14%, enrich 15.51x
     # 12 .006 .27
 # refGene:cds 0.861%, phastConsElements8way 4.962%, both 0.647%, cover 75.17%, enrich 15.15x
     # 13 .007 .27
 # refGene:cds 0.861%, phastConsElements8way 5.241%, both 0.651%, cover 75.57%, enrich 14.42x
     # 12 .008 .28
 # refGene:cds 0.861%, phastConsElements8way 5.152%, both 0.653%, cover 75.82%, enrich 14.72x
     # 12 .007 .27
 #                     phastConsElements8way 5.011%, both 0.648%, cover 75.19%, enrich 15.01x
     # 11 .005 .27
 # refGene:cds 0.861%, phastConsElements8way 4.673%, both 0.644%, cover 74.78%, enrich 16.00x
     # 14 .008 .28
 # refGene:cds 0.861%, phastConsElements8way 5.629%, both 0.659%, cover 76.57%, enrich 13.60x
     # 10,  .11
 # refGene:cds 0.861%, phastConsElements8way 7.936%, both 0.671%, cover 77.92%, enrich 9.82x
     # 12, .20
 # refGene:cds 0.861%, phastConsElements8way 8.316%, both 0.668%, cover 77.56%, enrich 9.33x
     # 12, .11
 # refGene:cds 0.861%, phastConsElements8way 7.304%, both 0.666%, cover 77.34%, enrich 10.59x
     # 14, .11
 # refGene:cds 0.861%, phastConsElements8way 7.936%, both 0.671%, cover 77.92%, enrich 9.82x
 
     # check conservation using consEntropy, at Adam's recommendation
     # first, generate second model file
     set rho = .23
     tree_doctor --scale $rho 
 
     # create wiggle files
     ssh pk
     cd /san/sanvol1/scratch/panTro2/multiz8way/cons/
     # sort by chromName, chromStart so that items are in numerical order 
     #  for wigEncode
     find ./pp -name "chr*.pp" | \
         sed -e "s/\// x /g" -e "s/\./ y /g" -e "s/-/ z /g" | \
 	sort -k7,7 -k9,9n | \
 	sed -e "s/ x /\//g" -e "s/ y /\./g" -e "s/ z /-/g" | \
 	xargs cat | \
         nice wigEncode stdin phastCons8way.wig phastCons8way.wib
     # half an hour or so
     mv phastCons8way.* /cluster/data/panTro2/bed/multiz8way/cons
 
     ssh hgwdev
     cd /cluster/data/panTro2/bed/multiz8way/cons
     ln -s /cluster/data/panTro2/bed/multiz8way/cons/phastCons8way.wib \
         /gbdb/panTro2/multiz8way/phastCons8way.wib
     hgLoadWiggle -pathPrefix=/gbdb/panTro2/multiz8way panTro2 \
         phastCons8way phastCons8way.wig
 
     # add Mark D's codon framing
     ssh kkstore04
     cd /cluster/data/panTro2/bed/multiz8way
     mkdir frames
     cd frames
     cp /cluster/data/hg17/bed/multiz17way/frames/mkMafFrames .
     cp /cluster/data/hg17/bed/multiz17way/frames/Makefile.genes .
     #edit Makefile to correct species names and set and sanDir, prune out all
     # but getGenes target, as the remainder is below.
 
     ssh hgwdev
     cd /cluster/data/panTro2/bed/multiz8way/frames
     make getGenes >&! getGenes.log &
         
     ssh kkstore04
     cd /cluster/data/panTro2/bed/multiz8way/frames
     nice tcsh # easy way to get process niced
     (cat  ../maf/*.maf | genePredToMafFrames panTro2 stdin stdout panTro2 genes/panTro2.gp.gz rheMac2 genes/rheMac2.gp.gz canFam2 genes/canFam2.gp.gz danRer4 genes/danRer4.gp.gz galGal2 genes/galGal2.gp.gz hg18 genes/hg18.gp.gz mm8 genes/mm8.gp.gz | gzip > multiz8way.mafFrames.gz)>&log&
  
     ssh hgwdev
     cd /cluster/data/panTro2/bed/multiz8way/frames
     hgLoadMafFrames panTro2 multiz8wayFrames multiz8way.mafFrames.gz >&log&
 
 ##########################################################################
 # NSCAN track - (2006-07-28 markd)
     cd /cluster/data/panTro2/bed/nscan/
 
     # obtained NSCAN predictions from michael brent's group
     # at WUSTL
     wget -nv -r -np http://ardor.wustl.edu/jeltje/panTro2/chr_gtf
     wget -nv -r -np http://ardor.wustl.edu/jeltje/panTro2/chr_ptx
     # clean up and rename downloaded directorys:
     mv ardor.wustl.edu/jeltje/panTro2/chr_gtf .
     mv ardor.wustl.edu/jeltje/panTro2/chr_ptx .
     rm -rf ardor.wustl.edu
     rm chr_*/index.html*
     gzip chr_*/*
     chmod a-w chr_**.gz
 
     # load tracks.  Note that these have *utr features, rather than
     # exon features.  currently ldHgGene creates separate genePred exons
     # for these.
     ldHgGene -bin -gtf -genePredExt panTro2 nscanGene chr_gtf/chr*.gtf.gz
     # add .a suffix to match transcript id
 
     # add .a suffix to match transcript id
     hgPepPred -suffix=.a panTro2 generic nscanPep chr_ptx/chr*.fa.gz
     rm *.tab
 
     # update trackDb; need a panTro2-specific page to describe informants
     chimp/panTro2/nscanGene.html    (from  CanFam1_hg17_nscan_10_18_2005/description.html )
     chimp/panTro2/trackDb.ra
 
 ##########################################################################
 # AUGUSTUS ab initio track  (DONE, mario, 1/30/2007)
     ssh hgwdev
     mkdir /cluster/data/panTro2/bed/augustus
     cd /cluster/data/panTro2/bed/augustus
 
     # get the program AUGUSTUS
     wget http://augustus.gobics.de/binaries/augustus.2.0.1.src.tar.gz
     tar xzf augustus.2.0.src.tar.gz
     # compile the binary if necessary
     cd augustus/src
     make augustus
 
     # create output directory
     cd /cluster/data/panTro2/bed/augustus
     mkdir out
 
     # create file with sequences and their sizes by modifying chrom.sizes
     cat ../../chrom.sizes | perl -e 'while(<>){s/chr([0-9a-zA-Z]+)/\/cluster\/data\/panTro2\/$1\/chr$1.fa.masked/; print;}' > seq.lst
 
     # create the job list
     augustus/scripts/createAugustusJoblist.pl --sequences seq.lst --chunksize 5300000 --overlap 300000 --command "/cluster/data/panTro2/bed/augustus/augustus/src/augustus --AUGUSTUS_CONFIG_PATH=/cluster/data/panTro2/bed/augustus/augustus/config --species=human --sample=100 --/augustus/verbosity=0" --outputdir /cluster/data/panTro2/bed/augustus/out/ --joblist job.lst
 
     para try
     para check
     para push
 
 # CPU time in finished jobs:    3785702s   63095.03m  1051.58h   43.82d  0.120 y
 # IO & Wait Time:                 38307s     638.45m    10.64h    0.44d  0.001 y
 # Average job time:                5785s      96.42m     1.61h    0.07d
 # Longest running job:                0s       0.00m     0.00h    0.00d
 # Longest finished job:            8767s     146.12m     2.44h    0.10d
 # Submission to last job:         20341s     339.02m     5.65h    0.24d
 
     # check the error files
     cat out/*.err
     # In this case no error was reported.
     rm out/*.err
 
     cat out/*.gff | augustus/scripts/join_aug_pred.pl > augustus.pep.gff
     augustus/scripts/getAnnoFasta.pl augustus.pep.gff
     cat augustus.pep.gff | egrep "CDS|codon"> augustus.gff
 
     # load into database
     ssh hgwdev
     cd /cluster/data/panTro2/bed/augustus/
     ldHgGene -bin panTro2 augustus augustus.gff
     # 33218 gene predictions
 
     hgPepPred panTro2 generic augustusPep augustus.pep.aa
     featureBits panTro2 augustus
     # 32820465 bases of 2909485072 (1.128%) in intersection
 #######################################################################
 # BLASTZ/CHAIN/NET HORSE AND CHIMP (DONE 2/21/07 Fan)
     ssh kkstore05
     mkdir /cluster/data/equCab1/bed/blastz.panTro2.2007-02-19
     cd /cluster/data/equCab1/bed/blastz.panTro2.2007-02-19
     cat << '_EOF_' > DEF
 # Horse vs. Chimp
 
 BLASTZ_M=50
 
 # TARGET: Horse equCab1
 SEQ1_DIR=/san/sanvol1/scratch/equCab1/equCab1.2bit
 SEQ1_LEN=/san/sanvol1/scratch/equCab1/chrom.sizes       
 # Maximum number of scaffolds that can be lumped together
 SEQ1_LIMIT=500     
 SEQ1_CHUNK=30000000
 SEQ1_LAP=10000
 
 # QUERY: Chimp panTro2
 SEQ2_DIR=/scratch/hg/panTro2/panTro2.2bit
 SEQ2_LEN=/cluster/data/panTro2/chrom.sizes 
 SEQ2_CHUNK=10000000
 SEQ2_LAP=0
 
 BASE=/cluster/data/equCab1/bed/blastz.panTro2.2007-02-19
 TMPDIR=/scratch/tmp
 '_EOF_'
     # << this line keeps emacs coloring happy
 doBlastzChainNet.pl DEF -chainMinScore=3000 -chainLinearGap=medium \
       -bigClusterHub pk \
       -blastzOutRoot /cluster/bluearc/equCab1PanTro2 >& do.log &
     tail -f do.log
     ln -s blastz.panTro2.2007-02-19 /cluster/data/equCab1/bed/blastz.panTro2
     nice featureBits equCab1 -chrom=chr1 chainPanTro2Link
 # 121391547 bases of 177498097 (68.390%) in intersection
 
     bash
     time nice -n 19 featureBits equCab1 chainPanTro2Link \
 	> fb.equCab1.chainPanTro2Link.txt 2>&1
     # 1593914101 bases of 2421923695 (65.812%) in intersection
 
     ssh kkstore04
     mkdir /cluster/data/panTro2/bed/blastz.equCab1.swap
     cd /cluster/data/panTro2/bed/blastz.equCab1.swap
     bash
     time doBlastzChainNet.pl \
 	/cluster/data/equCab1/bed/blastz.panTro2.2007-02-19/DEF \
 	-chainMinScore=3000 -chainLinearGap=medium \
 	-verbose=2 -swap -bigClusterHub=pk \
 	-blastzOutRoot /cluster/bluearc/equCab1PanTro2 > swap.log 2>&1 &
 
     # The -blastzOutRoot /cluster/bluearc/equCab1PanTro2 option above
     # should not be there.  It caused an error message during cleaning up.
     # Don't specify it for -swap next time.
 
     ssh hgwdev
     cd /cluster/data/panTro2/bed/blastz.equCab1.swap
     bash
     time nice -n 19 featureBits panTro2 chainEquCab1Link \
 	> fb.panTro2.chainEquCab1Link.txt 2>&1
     # 1636782924 bases of 2909485072 (56.257%) in intersection
 
 #########################################################################
 # Blastz Marmoset calJac1 (DONE - 2007-11-15 - Hiram)
     ssh kkstore04
     screen # use screen to control this job
     mkdir /cluster/data/panTro2/bed/blastzCalJac1.2007-11-13
     cd /cluster/data/panTro2/bed/blastzCalJac1.2007-11-13
 
     cat << '_EOF_' > DEF
 # Chimp vs marmoset
 BLASTZ_M=50
 
 # TARGET: Chimp PanTro2
 SEQ1_DIR=/scratch/hg/panTro2/panTro2.2bit
 SEQ1_LEN=/cluster/data/panTro2/chrom.sizes
 SEQ1_CHUNK=10000000
 SEQ1_LAP=10000
 
 # QUERY: Marmoset calJac1
 SEQ2_DIR=/cluster/bluearc/scratch/data/calJac1/calJac1.2bit
 SEQ2_LEN=/cluster/data/calJac1/chrom.sizes
 SEQ2_CHUNK=10000000
 SEQ2_LAP=0
 
 BASE=/cluster/data/panTro2/bed/blastzCalJac1.2007-11-13
 TMPDIR=/scratch/tmp
 '_EOF_'
     # << happy emacs
 
     time nice -n +19 doBlastzChainNet.pl -verbose=2 `pwd`/DEF \
 	-chainMinScore=3000 -chainLinearGap=medium \
 	-bigClusterHub=pk > blastz.log 2>&1 &
     #	real    1039m19.171s
     cat fb.panTro2.chainCalJac1Link.txt
     #	2220169777 bases of 2909485072 (76.308%) in intersection
 
     mkdir /cluster/data/calJac1/bed/blastz.panTro2.swap
     cd /cluster/data/calJac1/bed/blastz.panTro2.swap
     time nice -n +19 doBlastzChainNet.pl -verbose=2 \
 	/cluster/data/panTro2/bed/blastzCalJac1.2007-11-13/DEF \
 	-chainMinScore=3000 -chainLinearGap=medium \
 	-swap -bigClusterHub=pk > swap.log 2>&1 &
     #	real    320m14.293s
     cat fb.calJac1.chainPanTro2Link.txt
     #	2264115411 bases of 2929139385 (77.296%) in intersection
 
 ###########################################################################
 # Blastz Orangutan ponAbe2 (DONE - 2007-11-16 - Hiram)
     ssh kkstore04
     screen # use screen to control this job
     mkdir /cluster/data/panTro2/bed/blastzPonAbe2.2007-11-13
     cd /cluster/data/panTro2/bed/blastzPonAbe2.2007-11-13
 
     cat << '_EOF_' > DEF
 # Chimp vs marmoset
 BLASTZ_M=50
 
 # TARGET: Chimp PanTro2
 SEQ1_DIR=/scratch/hg/panTro2/panTro2.2bit
 SEQ1_LEN=/cluster/data/panTro2/chrom.sizes
 SEQ1_CHUNK=10000000
 SEQ1_LAP=10000
 
 # QUERY: Orangutan ponAbe2
 SEQ2_DIR=/cluster/bluearc/scratch/data/ponAbe2/ponAbe2.2bit
 SEQ2_LEN=/cluster/data/ponAbe2/chrom.sizes
 SEQ2_CHUNK=10000000
 SEQ2_LAP=0
 
 BASE=/cluster/data/panTro2/bed/blastzPonAbe2.2007-11-13
 TMPDIR=/scratch/tmp
 '_EOF_'
     # << happy emacs
 
     time nice -n +19 doBlastzChainNet.pl -verbose=2 `pwd`/DEF \
 	-chainMinScore=3000 -chainLinearGap=medium \
 	-bigClusterHub=kk > blastz.log 2>&1 &
     #	real    388m20.443s
 # Completed: 129160 of 129168 jobs
 # Crashed: 8 jobs
 # CPU time in finished jobs:   23694314s  394905.23m  6581.75h  274.24d  0.751 y
 # IO & Wait Time:               2536376s   42272.93m   704.55h   29.36d  0.080 y
 # Average job time:                 203s       3.38m     0.06h    0.00d
 # Longest finished job:            5665s      94.42m     1.57h    0.07d
 # Submission to last job:        170709s    2845.15m    47.42h    1.98d
 
     # some jobs failed due to memory limits on kk nodes
     #	finished them manually on kolossus and hgwdev, then continuing
     time nice -n +19 doBlastzChainNet.pl -verbose=2 `pwd`/DEF \
 	-chainMinScore=3000 -chainLinearGap=medium \
 	-continue=cat -bigClusterHub=kk > cat.log 2>&1 &
     #	real    754m45.014s
     cat fb.panTro2.chainPonAbe2Link.txt
     #	2648603020 bases of 2909485072 (91.033%) in intersection
 
     mkdir /cluster/data/ponAbe2/bed/blastz.panTro2.swap
     cd /cluster/data/ponAbe2/bed/blastz.panTro2.swap
     time nice -n +19 doBlastzChainNet.pl -verbose=2 \
 	/cluster/data/panTro2/bed/blastzPonAbe2.2007-11-13/DEF \
 	-chainMinScore=3000 -chainLinearGap=medium \
 	-swap -bigClusterHub=kk > swap.log 2>&1 &
     #	real    160m51.639s
     cat fb.ponAbe2.chainPanTro2Link.txt
     #	2766615040 bases of 3093572278 (89.431%) in intersection
     #	real    2597m24.771s
 
 ###########################################################################
 ############################################################################
 # TRANSMAP vertebrate.2008-05-20 build  (2008-05-24 markd)
 
 vertebrate-wide transMap alignments were built  Tracks are created and loaded
 by a single Makefile. This is available from:
    svn+ssh://hgwdev.cse.ucsc.edu/projects/compbio/usr/markd/svn/projs/transMap/tags/vertebrate.2008-05-20
 
 see doc/builds.txt for specific details.
 ############################################################################
 ############################################################################
 # TRANSMAP vertebrate.2008-06-07 build  (2008-06-30 markd)
 
 vertebrate-wide transMap alignments were built  Tracks are created and loaded
 by a single Makefile. This is available from:
    svn+ssh://hgwdev.cse.ucsc.edu/projects/compbio/usr/markd/svn/projs/transMap/tags/vertebrate.2008-06-30
 
 see doc/builds.txt for specific details.
 ############################################################################
 
 ################################################
 # AUTOMATE UPSTREAM FILE CREATION (2008-10-15 markd)
 update genbank.conf:
 panTro2.upstreamGeneTbl = refGene
 
 ############################################################################
 # SWAP HG19 Chain/Net (DONE - 2009-04-25 - Hiram)
     #	running the swap
     ssh swarm
     mkdir /hive/data/genomes/panTro2/bed/blastz.hg19.swap
     cd /hive/data/genomes/panTro2/bed/blastz.hg19.swap
     time nice -n +19 doBlastzChainNet.pl -verbose=2 \
 	-swap /hive/data/genomes/hg19/bed/lastzPanTro2.2009-03-19/DEF \
 	-noLoadChainSplit -chainMinScore=5000 -chainLinearGap=medium \
 	-workhorse=hgwdev -smallClusterHub=swarm -bigClusterHub=swarm \
 	> swap.log 2>&1 &
     #	real    723m41.377s
     cat fb.panTro2.chainHg19Link.txt 
     #	2761343871 bases of 2909485072 (94.908%) in intersection
 ############################################################################
 ############################################################################
 # TRANSMAP vertebrate.2009-07-01 build  (2009-07-21 markd)
 
 vertebrate-wide transMap alignments were built  Tracks are created and loaded
 by a single Makefile. This is available from:
    svn+ssh://hgwdev.cse.ucsc.edu/projects/compbio/usr/markd/svn/projs/transMap/tags/vertebrate.2009-07-01
 
 see doc/builds.txt for specific details.
 ############################################################################