5002bc264773860c46470525d2f8802f2c45679a
hiram
  Wed Sep 11 09:37:10 2019 -0700
liftOver to regenRn0 and regenRn1 refs #21579

diff --git src/hg/makeDb/doc/rn6.txt src/hg/makeDb/doc/rn6.txt
index c5a770f..c2165bc 100644
--- src/hg/makeDb/doc/rn6.txt
+++ src/hg/makeDb/doc/rn6.txt
@@ -1,4383 +1,4444 @@
 # for emacs: -*- mode: sh; -*-
 #
 #	the above keeps emacs happy while working with this text document
 
 # This file describes how we made the browser database on the
 # Rattus norvegicus genome, March 2012 update (Rnor5.0) from Baylor.
 
 #	http://www.ncbi.nlm.nih.gov/bioproject/16219
 #	http://www.ncbi.nlm.nih.gov/genome/73
 #	http://www.ncbi.nlm.nih.gov/Traces/wgs/?val=AABR07
 
 #	Genome Coverage : 3x BAC; 6x WGS ABI Sanger reads 10.0x PcaBio
 #	chrMt: NC_001665.2
 
 #	DATE:   01-Jul-2014
 #	ORGANISM:       Rattus norvegicus
 #	TAXID:  10116
 #	ASSEMBLY LONG NAME:     Rnor_6.0
 #	ASSEMBLY SHORT NAME:    Rnor_6.0
 #	ASSEMBLY SUBMITTER:     Rat Genome Sequencing Consortium
 #	ASSEMBLY TYPE:  Haploid
 #	NUMBER OF ASSEMBLY-UNITS:       2
 #	ASSEMBLY ACCESSION:     GCA_000001895.4
 #	##Below is a 2 column list with assembly-unit id and name.
 #	##The Primary Assembly unit is listed first.
 #	GCA_000000225.4 Primary Assembly
 #	GCA_000002045.1 non-nuclear
 #	FTP-RELEASE DATE: 03-Jul-2014
 
 #########################################################################
 ## Download sequence (DONE - 2014-07-07 - Hiram)
     mkdir -p /hive/data/genomes/rn6/genbank
     cd /hive/data/genomes/rn6/genbank
 
     time rsync -a -P \
 rsync://ftp.ncbi.nlm.nih.gov/genbank/genomes/Eukaryotes/vertebrates_mammals/Rattus_norvegicus/Rnor_6.0/ ./
 
     # sent 4565 bytes  received 2064907254 bytes  10350435.18 bytes/sec
     # total size is 2064638723  speedup is 1.00
 
     # real    3m19.289s
 
     # measure total sequence:
 
     faSize Primary_Assembly/assembled_chromosomes/FASTA/chr*.fa.gz \
       Primary_Assembly/unlocalized_scaffolds/FASTA/chr*.fa.gz \
       Primary_Assembly/unplaced_scaffolds/FASTA/unplaced.scaf.fa.gz \
       non-nuclear/assembled_chromosomes/FASTA/chrMT.fa.gz
 
 # 2870184193 bases (140322104 N's 2729862089 real 2729862089 upper 0 lower)
 #   in 955 sequences in 46 files
 # Total size: mean 3005428.5 sd 21316687.8
 #   min 230 (gi|661406602|gb|AABR07024268.1|)
 #   max 282763074 (gi|661902784|gb|CM000072.5|) median 7105
 
 #########################################################################
 # fixup names for UCSC standards (DONE - 2014-07-07 - Hiram)
     mkdir /hive/data/genomes/rn6/ucsc
     cd /hive/data/genomes/rn6/ucsc
 
     # can re-use the scripts from hg38 build
     cp -p ../../hg38/ucsc/ucscCompositeAgp.pl .
     time ./ucscCompositeAgp.pl
     # real    0m53.149s
 
     cp -p ../../hg38/ucsc/unlocalized.pl .
     # fixed this script to be more generic so it should work on any
     # assembly in the future as long as the fasta headers maintain the
     # observed structre
     time ./unlocalized.pl
     # real    0m2.097s
 
     cp -p ../../hg38/ucsc/unplaced.pl .
     # fixed this script to be more generic so it should work on any
     # assembly in the future as long as the fasta headers maintain the
     # observed structre
     time ./unplaced.pl
     # real    0m0.494s
 
     zcat ../genbank/non-nuclear/assembled_chromosomes/FASTA/chrMT.fa.gz \
        | sed -e 's/^>.*/>chrM/' > chrM.fa
     zcat ../genbank/non-nuclear/assembled_chromosomes/AGP/chrMT.comp.agp.gz \
       | grep AY172581 | sed -e 's/^AY172581.1/chrM/;' > chrM.agp
 
     # more efficient to have everything zipped:
     time gzip *.fa
     # real    16m8.206s
 
     time gzip *.agp
     #  real    0m1.725s
 
     #	verify all the sequence is still here after all this rigamarole:
     time faSize *.fa.gz
 # 2870184193 bases (140322104 N's 2729862089 real 2729862089 upper 0 lower)
 #    in 955 sequences in 46 files
 # Total size: mean 3005428.5 sd 21316687.8 min 230 (chrUn_AABR07024268v1)
 #    max 282763074 (chr1) median 7105
 
     # real    1m3.343s
 
     # verify same numbers as was in the original files measured above:
 
 # 2870184193 bases (140322104 N's 2729862089 real 2729862089 upper 0 lower)
 #   in 955 sequences in 46 files
 # Total size: mean 3005428.5 sd 21316687.8
 #   min 230 (gi|661406602|gb|AABR07024268.1|)
 #   max 282763074 (gi|661902784|gb|CM000072.5|) median 7105
 
 #########################################################################
 # Create .ra file and run makeGenomeDb.pl (DONE - Hiram - 2014-07-07)
     cd /hive/data/genomes/rn6
     cat << '_EOF_' >rn6.config.ra
 # Config parameters for makeGenomeDb.pl:
 db rn6
 clade mammal
 scientificName Rattus norvegicus
 commonName Rat
 assemblyDate Jul. 2014
 assemblyLabel RGSC Rnor_6.0
 assemblyShortLabel RGSC 6.0
 orderKey 1558
 # chrMt NC_001665.2 included in the genbank release
 mitoAcc none
 fastaFiles /hive/data/genomes/rn6/ucsc/chr*.fa.gz
 agpFiles /hive/data/genomes/rn6/ucsc/chr*.agp.gz
 # qualFiles none
 dbDbSpeciesDir rat
 photoCreditURL http://www.genome.gov/10005141
 photoCreditName NHGRI Press Photos
 ncbiGenomeId 73
 ncbiAssemblyId 382928
 ncbiAssemblyName Rnor_6.0
 ncbiBioProject 10629
 genBankAccessionID GCA_000001895.4
 taxId 10116
 # http://www.ncbi.nlm.nih.gov/Traces/wgs/?val=AABR07
 '_EOF_'
     # << happy emacs
 
     #	run agp step first to verify fasta and agp files agree
     time makeGenomeDb.pl -stop=agp rn6.config.ra > agp.log 2>&1
 # ERROR: duplicate sequence found in rn6.unmasked.2bit
 # echo ERROR: duplicate sequence found in rn6.unmasked.2bit
     twoBitDup rn6.unmasked.2bit
 # chrUn_AABR07022993v1 and chrUn_AABR07022518v1 are identical
 # chrUn_AABR07023006v1 and chrUn_AABR07022518v1 are identical
     # remove chrUn_AABR07022518v1 chrUn_AABR07023006v1
     cd ucsc
     echo "chrUn_AABR07022518v1
 chrUn_AABR07023006v1" > dup.list
     faSomeRecords  -exclude chrUn.fa.gz dup.list chrUn.clean.fa
     mv chrUn.fa.gz chrUn.fa.gz.dups
     mv chrUn.clean.fa chrUn.fa
     gzip chrUn.fa
     zcat chrUn.agp.gz | egrep -v "chrUn_AABR07022518v1|chrUn_AABR07023006v1" \
       > chrUn.agp
     mv chrUn.agp.gz chrUn.agp.gz.dups
     gzip chrUn.agp
 
     cd ..
     rm -fr bed jkStuff rn6.unmasked.2bit
     mv agp.log agp.log.0
     time makeGenomeDb.pl -stop=agp rn6.config.ra > agp.log 2>&1
     # verify end of agp.log indictates:
 # All AGP and FASTA entries agree - both files are valid
 
     # the non-dup sequence is now:
     twoBitToFa rn6.unmasked.2bit stdout | faSize stdin
     #   2870182909 bases (140322104 N's 2729860805 real 2729860805 upper
     #    0 lower) in 953 sequences in 1 files
 
     # continue with the build, don't let this run without a workhorse
     # specified, it can end up on a machine that can't run the wigToBigWig
     # for the gc5Base track:
     time makeGenomeDb.pl -workhorse=hgwdev -fileServer=hgwdev \
       -continue=db rn6.config.ra > db.log 2>&1
     #	real    23m10.146s
 
 #########################################################################
 # fixup search rule for assembly track/gold table (DONE - 2014-07-09 - Hiram)
    export maxLen=`hgsql -N -e 'select frag from gold;' rn6 | awk '{print length($0)}' | sort -run | head -1`
 
    echo $maxLen
    #  14
 
 export C=1
 while [ $C -le $maxLen ];
 do
   echo -n " $C: "
   hgsql -N -e 'select frag from gold;' rn6 | sort -u \
     | awk '{ print substr($0,'$C',1) }' | sort -u | xargs echo | sed -e 's/ //g'
   C=`echo $C | awk '{print $1+1}'`
 done
 #   1: A
 #   2: ACY
 #   3: 012B
 #   4: 012345789R
 #   5: 0123456789
 #   6: 0123456789
 #   7: 0123456789
 #   8: 0123456789
 #   9: .0123456789
 #  10: 0123456789
 #  11: 0123456789
 #  12: 0123456789
 #  13: .
 #  14: 1
 
     # since we have the composite fragment IDs, there are a variety of
     # dot versions:
     hgsql -N -e 'select frag from gold;' rn6 | sed -e 's/.*\.//' \
        | sort | uniq -c  | sort -rn
 #  73744 1
 #    574 5
 #    540 4
 #    390 6
 # ... etc ...
 
     # hence, add to trackDb/rat/rn6/trackDb.ra
 
 searchTable gold
 searchMethod prefix
 searchType bed
 shortCircuit 1
 termRegex A[ACY][B012][R0-9][0-9]+(\.[0-9]+)*
 query select chrom,chromStart,chromEnd,frag from %s where frag like '%s%%'
 searchPriority 8
 
   # test pattern:
    hgsql -N -e 'select frag from gold;' rn6 | wc -l
    # 76195
    hgsql -N -e 'select frag from gold;' rn6 \
      | egrep -e 'A[ACY][B012][R0-9][0-9]+(\.[0-9]+)*' | wc -l
    # 76195
    hgsql -N -e 'select frag from gold;' rn6 \
      | egrep -v -e 'A[ACY][B012][R0-9][0-9]+(\.[0-9]+)*' | wc -l
    # 0
 
 ############################################################################
 # create ucscToINSDC name mapping (DONE - 2014-04-11 - Hiram)
     mkdir /hive/data/genomes/rn6/bed/ucscToINSDC
     cd /hive/data/genomes/rn6/bed/ucscToINSDC
 
     # this script has been maturing over time, it is close to complete.
     # to find a latest copy of it:
     # ls -ogrt /hive/data/genomes/*/bed/ucscToINSDC/translateNames.sh
 
     cp -p /hive/data/genomes/poeRet1/bed/ucscToINSDC/translateNames.sh .
     # to see chrM accession:
     hgsql -e 'select * from gold;' rn6 | grep chrM
     # use that accession here:
     ./translateNames.sh AY172581.1
     # eliminate the duplicates, and needs to be sorted to work with join
     egrep -v "chrUn_AABR07022518v1|chrUn_AABR07023006v1" ucscToINSDC.txt \
         | sort > ucscToINSDC.tab
 
     awk '{printf "%s\t0\t%d\n", $1,$2}' ../../chrom.sizes | sort \
         > name.coordinate.tab
 
     join name.coordinate.tab ucscToINSDC.tab | tr '[ ]' '[\t]' > ucscToINSDC.bed
 
     cut -f1 ucscToINSDC.bed | awk '{print length($0)}' | sort -n | tail -1
 # 27
 
     # use the 27 in this sed:
     sed -e "s/21/27/" $HOME/kent/src/hg/lib/ucscToINSDC.sql \
         | hgLoadSqlTab rn6 ucscToINSDC stdin ucscToINSDC.bed
     checkTableCoords rn6 ucscToINSDC
     # should cover all bases
     featureBits -countGaps rn6 ucscToINSDC
     # 2870182909 bases of 2870182909 (100.000%) in intersection
 
 ##############################################################################
 # running repeat masker (DONE - 2014-07-07,09 - Hiram)
     mkdir /hive/data/genomes/rn6/bed/repeatMasker
     cd /hive/data/genomes/rn6/bed/repeatMasker
     time doRepeatMasker.pl -buildDir=`pwd` \
 	-bigClusterHub=ku -dbHost=hgwdev -workhorse=hgwdev \
 	-smallClusterHub=ku rn6 > do.log 2>&1 &
     #  real    1921m43.379s
 
     cat faSize.rmsk.txt
     # 2870182909 bases (140322104 N's 2729860805 real 1551246793 upper
     #    1178614012 lower) in 953 sequences in 1 files
     # Total size: mean 3011734.4 sd 21338622.2 min 230 (chrUn_AABR07024268v1)
     #    max 282763074 (chr1) median 7294
     # %41.06 masked total, %43.17 masked real
 
     grep -i versi do.log
     # RepeatMasker version open-4.0.5
     #    January 31 2015 (open-4-0-5) version of RepeatMasker
 
     featureBits -countGaps rn6 rmsk
     #	1178980732 bases of 2870182909 (41.077%) in intersection
     # why is it different than the faSize above ?
     # because rmsk masks out some N's as well as bases, the faSize count above
     #	separates out the N's from the bases, it doesn't count lower case N's
 
 ##########################################################################
 # running simple repeat (DONE - 2014-07-07 - Hiram)
     mkdir /hive/data/genomes/rn6/bed/simpleRepeat
     cd /hive/data/genomes/rn6/bed/simpleRepeat
     time doSimpleRepeat.pl -buildDir=`pwd` -bigClusterHub=ku \
 	-dbHost=hgwdev -workhorse=hgwdev -smallClusterHub=ku rn6 > do.log 2>&1 &
     #	real    8m49.125s
 
     cat fb.simpleRepeat
     #	106804757 bases of 2731211732 (3.911%) in intersection
 
     # add the TRF mask to the rmsk sequence:
     #	it masks more sequence
     cd /hive/data/genomes/rn6
     twoBitMask rn6.rmsk.2bit \
 	-add bed/simpleRepeat/trfMask.bed rn6.2bit
     #	you can safely ignore the warning about fields >= 13
 
     twoBitToFa rn6.2bit stdout | faSize stdin > faSize.rn6.2bit.txt
     cat faSize.rn6.2bit.txt
     #   2870182909 bases (140322104 N's 2729860805 real 1547693783 upper
     #   1182167022 lower) in 953 sequences in 1 files
     #   Total size: mean 3011734.4 sd 21338622.2 min 230 (chrUn_AABR07024268v1)
     #       max 282763074 (chr1) median 7294
     #   %41.19 masked total, %43.31 masked real
 
     #	replace the previous symLink which goes to the unmasked 2bit
     rm /gbdb/rn6/rn6.2bit
     ln -s `pwd`/rn6.2bit /gbdb/rn6/rn6.2bit
 
 #########################################################################
 # Verify all gaps are marked (DONE - 2014-07-08 - Hiram)
     mkdir /hive/data/genomes/rn6/bed/gap
     cd /hive/data/genomes/rn6/bed/gap
     time nice -n +19 findMotif -motif=gattaca -verbose=4 \
 	-strand=+ ../../rn6.unmasked.2bit > findMotif.txt 2>&1
     #	real    0m40.856s
     grep "^#GAP " findMotif.txt | sed -e "s/^#GAP //" > allGaps.bed
     time featureBits -countGaps rn6 -not gap -bed=notGap.bed
     #   2731211732 bases of 2870182909 (95.158%) in intersection
     #   real    0m15.822s
 
     time featureBits -countGaps rn6 allGaps.bed notGap.bed -bed=new.gaps.bed
     # 1350927 bases of 2870182909 (0.047%) in intersection
     # real    0m56.160s
     awk '{print $3-$2}' new.gaps.bed | ave stdin | grep total
     #  total 1350927.000000
     awk '{print $3-$2,$0}' new.gaps.bed | sort -nr | head
 # 10377 chr2      172413175       172423552       chr2.4447
 # 7000 chr8       119186607       119193607       chr8.2918
 # 6344 chrX       123693517       123699861       chrX.3045
 
 
     #	what is the highest index in the existing gap table:
     hgsql -N -e "select ix from gap;" rn6 | sort -n | tail -1
     #	14766
     cat << '_EOF_' > mkGap.pl
 #!/bin/env perl
 
 use strict;
 use warnings;
 
 my $ix=`hgsql -N -e "select ix from gap;" rn6 | sort -n | tail -1`;
 chomp $ix;
 
 open (FH,"<new.gaps.bed") or die "can not read new.gaps.bed";
 while (my $line = <FH>) {
     my ($chrom, $chromStart, $chromEnd, $rest) = split('\s+', $line);
     ++$ix;
     printf "%s\t%d\t%d\t%d\tN\t%d\tother\tyes\n", $chrom, $chromStart,
         $chromEnd, $ix, $chromEnd-$chromStart;
 }
 close (FH);
 '_EOF_'
     # << happy emacs
     chmod +x ./mkGap.pl
     ./mkGap.pl > other.bed
     featureBits -countGaps rn6 other.bed
     #	1350927 bases of 2870182909 (0.047%) in intersection
     wc -l other.bed
     #	75068
     # verify no mistake here:
     featureBits -countGaps rn6 gap other.bed
     #	0 bases of 2870182909 (0.000%) in intersection
 
     hgLoadBed -sqlTable=$HOME/kent/src/hg/lib/gap.sql \
 	-noLoad rn6 otherGap other.bed
     #	starting with this many
     hgsql -e "select count(*) from gap;" rn6
     #	70664
     hgsql rn6 -e 'load data local infile "bed.tab" into table gap;'
     #	result count:
     hgsql -e "select count(*) from gap;" rn6
     #	145732
     calc 70664 \+ 75068
     #	70664 + 75068 = 145732.000000
 
     # verify we aren't adding gaps where gaps already exist
     # this would output errors if that were true:
     gapToLift -minGap=1 rn6 nonBridged.lift -bedFile=nonBridged.bed
     # see example in danRer7.txt when problems arise
     # this can happen and is not an error:
 # WARNING: overlapping gap at chrY:692274-692374(scaffold) and chrY:692374-692375(other)
     # the scaffold following the (scaffold) gap begins with a single N
     # how that can happen I have no idea
 
     # there are bridged and non-bridged gaps here:
     hgsql -N -e "select bridge from gap;" rn6 | sort | uniq -c
     #     440 no
     #  145292 yes
 
 #########################################################################
 # cytoBandIdeo - (DONE - 2014-07-14 - Hiram)
     mkdir /hive/data/genomes/rn6/bed/cytoBand
     cd /hive/data/genomes/rn6/bed/cytoBand
     makeCytoBandIdeo.csh rn6
 
 #########################################################################
 ##  CYTOBAND - ideogram track (DONE - 2017-03-15 - Hiram)
     mkdir /hive/data/genomes/rn6/bed/cytoBand
     cd /hive/data/genomes/rn6/bed/cytoBand
     # fetch the ideogram file:
     rsync -a -P rsync://ftp.ncbi.nlm.nih.gov/pub/gdp/ideogram_10116_GCF_000000225.4_NA_V1 ./
 
     # Create bed file
     $HOME/kent/src/utils/ncbi/createNcbiCytoBand.pl \
        ideogram_10116_GCF_000000225.4_NA_V1
 
     # add in the default bands for all other chroms:
     egrep "_|chrM" firstTrack/rn6.cytoBand \
       | awk '{printf "%s\t%s\t%d\t%d\t%s\n", $1,$1,$2,$3,$4}' \
         | sed -e 's/^chr[nUXY0-9]\+_//; s/v1_random/.1/;' \
          | awk '{printf "%s\t%d\t%d\t\t%s\n", $2,$3,$4,$5}' >> cytoBand.bed
 
     ## can now verify before load:
     $HOME/kent/src/utils/ncbi/cytoBandVerify.pl
     #	everything checks out OK on 953 chroms
 
     # these chrom names are too long for the SQL keys:
  cut -f1 cytoBand.bed | sort -u | awk '{print length($0)}' | sort -rn | head -1
 #   27
     sed -e 's/12/27/;' $HOME/kent/src/hg/lib/cytoBand.sql \
         > cytoBand.sql
 
     # Load the bed file
     hgLoadBed -tab -noBin -sqlTable=cytoBand.sql rn6 cytoBand cytoBand.bed
     #   Read 1179 elements of size 5 from cytoBand.bed
 
     # Make cytoBandIdeo track for ideogram gif on hgTracks page.
     # The cytoBandIdeo is just a replicate of the cytoBand track.
     hgsql -e "CREATE TABLE cytoBandIdeo SELECT * FROM cytoBand;" rn6
 
     checkTableCoords rn6 cytoBand
     checkTableCoords rn6 cytoBandIdeo
 
     # should cover %100
     featureBits -countGaps rn6 cytoBand
     # 2870182909 bases of 2870182909 (100.000%) in intersection
     featureBits -countGaps rn6 cytoBandIdeo
     # 2870182909 bases of 2870182909 (100.000%) in intersection
 
 ##########################################################################
 ## WINDOWMASKER (DONE - 2014-07-09 - Hiram)
     mkdir /hive/data/genomes/rn6/bed/windowMasker
     cd /hive/data/genomes/rn6/bed/windowMasker
     time nice -n +19 doWindowMasker.pl -buildDir=`pwd` -workhorse=hgwdev \
 	-dbHost=hgwdev rn6 > do.log 2>&1 &
     #	real    205m59.097s
 
     # Masking statistics
     cat fb.rn6.windowmaskerSdust.clean.txt
     #  949668059 bases of 2870182909 (33.087%) in intersection
     cat fb.rn6.rmsk.windowmaskerSdust.txt
     #  712388228 bases of 2870182909 (24.820%) in intersection
 
 #########################################################################
 # MASK SEQUENCE WITH WM+TRF
     # not running this since RM + TRF is plenty of masking
 #    cd /hive/data/genomes/rn6
 #    twoBitMask -add bed/windowMasker/rn6.cleanWMSdust.2bit \
 #	bed/simpleRepeat/trfMask.bed rn6.2bit
     #	safe to ignore the warnings about BED file with >=13 fields
 #    twoBitToFa rn6.2bit stdout | faSize stdin > faSize.rn6.txt
 #    cat faSize.rn6.txt
 
     #	create symlink to gbdb
 #    ssh hgwdev
 #    rm /gbdb/rn6/rn6.2bit
 #    ln -s `pwd`/rn6.2bit /gbdb/rn6/rn6.2bit
 
 #########################################################################
 # PREPARE LINEAGE SPECIFIC REPEAT FILES FOR BLASTZ (TBD - 2012-03-23 - Hiram)
     ssh ku
     mkdir /hive/data/genomes/rn6/bed/linSpecRep
     cd /hive/data/genomes/rn6/bed/linSpecRep
 
     # split the RM output by chromosome name into separate files
     mkdir rmsk dateRepeats
     head -3 ../repeatMasker/rn6.sorted.fa.out > rmsk.header.txt
     headRest 3 ../repeatMasker/rn6.sorted.fa.out \
 	| splitFileByColumn -ending=.out -col=5 -head=rmsk.header.txt stdin rmsk
 
     ls -1S rmsk/* > rmOut.list
     wc -l rmOut.list
     #	904 rmOut.list
     wc -l ../../chrom.sizes
     # 953 ../../chrom.sizes
     #evidently there are 49 segments without RM annotations, some are very small
 
     cat << '_EOF_' > mkLSR
 #!/bin/csh -fe
 rm -f dateRepeats/$1_homo-sapiens_mus-musculus
 /scratch/data/RepeatMasker140131/DateRepeats \
     $1 -query rat -comp human -comp mouse
 mv $1_homo-sapiens_mus-musculus dateRepeats
 '_EOF_'
     #	<< happy emacs
     chmod +x mkLSR
 
     cat << '_EOF_' > template
 #LOOP
 ./mkLSR $(path1) {check out line+ dateRepeats/$(file1)_homo-sapiens_mus-musculus}
 #ENDLOOP
 '_EOF_'
     #	<< happy emacs
 
     gensub2 rmOut.list single template jobList
     para create jobList
     para try ... check ... push ... etc...
     para time
 # Completed: 904 of 904 jobs
 # CPU time in finished jobs:      21844s     364.06m     6.07h    0.25d  0.001 y
 # IO & Wait Time:                 20183s     336.39m     5.61h    0.23d  0.001 y
 # Average job time:                  46s       0.77m     0.01h    0.00d
 # Longest finished job:              96s       1.60m     0.03h    0.00d
 # Submission to last job:           682s      11.37m     0.19h    0.01d
 
     mkdir notInHuman notInMouse
     for F in dateRepeats/chr*.out_homo-sapiens*
     do
 	B=`basename ${F}`
 	B=${B/.out*/}
 	echo $B
         /cluster/bin/scripts/extractRepeats 1 ${F} > \
 		notInHuman/${B}.out.spec
         /cluster/bin/scripts/extractRepeats 2 ${F} > \
 		notInMouse/${B}.out.spec
     done
 
     #	Verify that these two things are actually different
     #	To check identical
     find ./notInHuman ./notInMouse -name "*.out.spec" | \
 	while read FN; do echo `cat ${FN} | sum -r` ${FN}; done \
 	| sort -k1,1n | sort -t"/" -k3,3 > check.same
     # some of them are the same, but not all:
     sed -e 's#./notInHuman/##; s#./notInMouse/##' check.same \
 	| sort | uniq -c | sort -rn | less
     # you will see a count of two at the beginning, but it becomes one soon
 
     #	Copy to data/staging for cluster replication
     mkdir /hive/data/staging/data/rn6
     rsync -a -P ./notInMouse/ /hive/data/staging/data/rn6/notInMouse/
     rsync -a -P ./notInHuman/ /hive/data/staging/data/rn6/notInHuman/
 
     # We also need the nibs for the lastz runs with lineage specific repeats
     mkdir /hive/data/staging/data/rn6/nib
     mkdir /hive/data/genomes/rn6/nib
     cd /hive/data/genomes/rn6
     cut -f1 chrom.sizes | while read C
 do
     twoBitToFa -seq=${C} rn6.2bit stdout | faToNib -softMask stdin nib/${C}.nib
     ls -og nib/$C.nib
 done
     # verify sequence remains the same
     cut -f1 chrom.sizes | while read C
 do
     nibFrag -masked nib/${C}.nib 0 `grep -w ${C} chrom.sizes | cut -f2` + \
 	stdout
 done | faSize stdin
     #   2870182909 bases (140322104 N's 2729860805 real 1547693783 upper
     #      1182167022 lower) in 953 sequences in 1 files
     #   Total size: mean 3011734.4 sd 21338622.2
     #      min 230 (nib/chrUn_AABR07024268v1.nib:0-230)
     #      max 282763074 (nib/chr1.nib:0-282763074) median 7294
     #   %41.19 masked total, %43.31 masked real
 
     #   Compare to original:
     cat faSize.rn6.2bit.txt
     #   2870182909 bases (140322104 N's 2729860805 real 1547693783 upper
     #      1182167022 lower) in 953 sequences in 1 files
     #   Total size: mean 3011734.4 sd 21338622.2
     #      min 230 (chrUn_AABR07024268v1) max 282763074 (chr1) median 7294
     #   %41.19 masked total, %43.31 masked real
 
     #	Copy to data/genomes staging for cluster replication
     rsync -a -P ./nib/ /hive/data/staging/data/rn6/nib/
 
 #########################################################################
 # cpgIslands - (DONE - 2014-07-09 - Hiram)
     mkdir /hive/data/genomes/rn6/bed/cpgIslands
     cd /hive/data/genomes/rn6/bed/cpgIslands
     time (doCpgIslands.pl -workhorse=hgwdev -buildDir=`pwd` \
        -dbHost=hgwdev -smallClusterHub=ku -bigClusterHub=ku rn6) > do.log 2>&1
     #   real    2m32.227s
 
     cat fb.rn6.cpgIslandExt.txt
     #   11404910 bases of 2729860805 (0.418%) in intersection
 
 #############################################################################
 # CPG Islands Unmasked track (DONE - 2014-07-08 - Hiram)
 
     mkdir /hive/data/genomes/rn6/bed/cpgIslandsUnmasked
     cd /hive/data/genomes/rn6/bed/cpgIslandsUnmasked
     time (doCpgIslands.pl -buildDir=`pwd` -bigClusterHub=ku \
       -tableName=cpgIslandExtUnmasked -dbHost=hgwdev -smallClusterHub=ku \
          -workhorse=hgwdev \
          -maskedSeq=/hive/data/genomes/rn6/rn6.unmasked.2bit rn6) > do.log 2>&1
     # real    14m0.085s
 
     cat fb.rn6.cpgIslandExtUnmasked.txt
     # 13483527 bases of 2731211732 (0.494%) in intersection
 
 #########################################################################
 # genscan - (DONE - 2014-07-09 - Hiram)
     mkdir /hive/data/genomes/rn6/bed/genscan
     cd /hive/data/genomes/rn6/bed/genscan
     time (doGenscan.pl -bigClusterHub=ku -workhorse=hgwdev \
       -dbHost=hgwdev -buildDir=`pwd` rn6) > do.log 2>&1 &
     # real    137m3.209s
 
     cat fb.rn6.genscan.txt
     #   57944086 bases of 2729860805 (2.123%) in intersection
     cat fb.rn6.genscanSubopt.txt
     #   62473326 bases of 2729860805 (2.289%) in intersection
 
 #########################################################################
 # MAKE 11.OOC FILE FOR BLAT/GENBANK (DONE - 2014-07-09 - Hiram)
     # Use -repMatch=1000, based on size -- for human we use 1024
     # use the "real" number from the faSize measurement,
     # hg19 is 2897316137, calculate the ratio factor for 1024:
     calc \( 2729860805 / 2897316137 \) \* 1024
     #	( 2729860805 / 2897316137 ) * 1024 = 964.816172
 
     # round up to 1000  (rn5 was 950)
 
     cd /hive/data/genomes/rn6
     time blat rn6.2bit /dev/null /dev/null -tileSize=11 \
       -makeOoc=jkStuff/rn6.11.ooc -repMatch=1000
     #   Wrote 27021 overused 11-mers to jkStuff/rn6.11.ooc
     #   rn5 had: Wrote 34513 overused 11-mers to jkStuff/rn6.11.ooc
     #	real    1m10.077s
 
     # there are non-bridged gaps, create lift file needed for genbank
     hgsql -N -e "select bridge from gap;" rn6 | sort | uniq -c
     #       440 no
     #    145292 yes
 
     cd /hive/data/genomes/rn6/jkStuff
     gapToLift rn6 rn6.nonBridged.lift -bedFile=rn6.nonBridged.bed
     # largest non-bridged contig:
     awk '{print $3-$2,$0}' rn6.nonBridged.bed | sort -nr | head
     #   52752077 chr18  35449852        88201929        chr18.06
 
 #########################################################################
 # AUTO UPDATE GENBANK (DONE - 2012-05-04 - Hiram)
     # examine the file:
     /cluster/data/genbank/data/organism.lst
     # for your species to see what counts it has for:
 # organism       mrnaCnt estCnt  refSeqCnt
 # Rattus norvegicus       125869  1103595 17503
 
     # to decide which "native" mrna or ests you want to specify in genbank.conf
 
     ssh hgwdev
     cd $HOME/kent/src/hg/makeDb/genbank
     git pull
     # edit etc/genbank.conf to add:
 
 # rn6 (rat)
 rn6.serverGenome = /hive/data/genomes/rn6/rn6.2bit
 rn6.clusterGenome = /hive/data/genomes/rn6/rn6.2bit
 rn6.ooc = /hive/data/genomes/rn6/jkStuff/rn6.11.ooc
 rn6.lift = /hive/data/genomes/rn6/jkStuff/rn6.nonBridged.lift
 rn6.refseq.mrna.native.pslCDnaFilter  = ${finished.refseq.mrna.native.pslCDnaFilter}
 rn6.refseq.mrna.xeno.pslCDnaFilter    = ${finished.refseq.mrna.xeno.pslCDnaFilter}
 rn6.genbank.mrna.native.pslCDnaFilter = ${finished.genbank.mrna.native.pslCDnaFilter}
 rn6.genbank.mrna.xeno.pslCDnaFilter   = ${finished.genbank.mrna.xeno.pslCDnaFilter}
 rn6.genbank.est.native.pslCDnaFilter  = ${finished.genbank.est.native.pslCDnaFilter}
 rn6.downloadDir = rn6
 rn6.refseq.mrna.xeno.load  = yes
 rn6.refseq.mrna.xeno.loadDesc = yes
 rn6.genbank.mrna.xeno.load  = yes
 rn6.perChromTables = no
 rn6.mgc = yes
 # rn6.upstreamGeneTbl = ensGene
 # rn6.upstreamMaf = multiz13way
 # /hive/data/genomes/rn6/bed/multiz13way/species.list.txt
 
     # end of section added to etc/genbank.conf
     git commit -m "adding rn6 rat refs #13578" etc/genbank.conf
     git push
     make etc-update
 
     ssh hgwdev			# used to do this on "genbank" machine
     screen -S rn6           # long running job managed in screen
     cd /cluster/data/genbank
     time nice -n +19 ./bin/gbAlignStep -initial rn6 &
     #	var/build/logs/2014.07.09-12:07:40.rn6.initalign.log
     #   real    902m50.425s
 
     # load database when finished
     ssh hgwdev
     cd /cluster/data/genbank
     time nice -n +19 ./bin/gbDbLoadStep -drop -initialLoad rn6 &
     #	logFile: var/dbload/hgwdev/logs/2014.07.10-18:27:18.rn6.dbload.log
     #   real    68m41.070s
 
     # enable daily alignment and update of hgwdev (DONE - 2012-02-09 - Hiram)
     cd ~/kent/src/hg/makeDb/genbank
     git pull
     # add rn6 to:
         etc/align.dbs etc/hgwdev.dbs
     git commit -m "Added rn6 refs #13578" etc/align.dbs etc/hgwdev.dbs
     git push
     make etc-update
 
 #########################################################################
 # setup search rule for assembly track (DONE - 2014-07-08 - Hiram)
    export maxLen=`hgsql -N -e 'select frag from gold;' rn6 | awk '{print length($0)}' | sort -run | head -1`
 
    echo $maxLen
    #  14
 
 export C=1
 while [ $C -le $maxLen ];
 do
   echo -n " $C: "
   hgsql -N -e 'select frag from gold;' rn6 | sort -u \
     | awk '{ print substr($0,'$C',1) }' | sort -u | xargs echo | sed -e 's/ //g'
   C=`echo $C | awk '{print $1+1}'`
 done
 
 #   1: A
 #   2: ACY
 #   3: 012B
 #   4: 012345789R
 #   5: 0123456789
 #   6: 0123456789
 #   7: 0123456789
 #   8: 0123456789
 #   9: .0123456789
 #  10: 0123456789
 #  11: 0123456789
 #  12: 0123456789
 #  13: .
 #  14: 1
 
 searchTable gold
 searchMethod prefix
 searchType bed
 shortCircuit 1
 termRegex A[ACY][012B][0-9R][0-9]+(\.1)*
 query select chrom,chromStart,chromEnd,frag from %s where frag like '%s%%'
 searchPriority 8
 
   # test pattern:
    hgsql -N -e 'select frag from gold;' rn6 | wc -l
    # 76195
    hgsql -N -e 'select frag from gold;' rn6 \
      | egrep -e 'A[ACY][B012][R0-9][0-9]+(\.[0-9]+)*' | wc -l
    # 76195
    hgsql -N -e 'select frag from gold;' rn6 | sed -e 's/\.[0-9]+$//' \
 	| egrep -e 'A[ACY][B012][R0-9][0-9]+(\.[0-9]+)*' | wc -l
    # 76195
    hgsql -N -e 'select frag from gold;' rn6 \
      | egrep -v -e 'A[ACY][B012][R0-9][0-9]+(\.[0-9]+)*' | wc -l
    # 0
 
 ##############################################################################
 # construct liftOver from rn5 to rn6 (DONE - 2014-07-14 - Hiram)
     # documentation for this step is in rn5.txt
 
 #########################################################################
 ## Default position to same as rn5 found via lift over data n rn5
 ##  (DONE - 2014-07-14 - Hiram)
     ssh hgwdev
     hgsql -e 'update dbDb set defaultPos="chr1:80608553-80639261"
 	where name="rn6";' hgcentraltest
 
 ##############################################################################
 # construct downloads files (DONE - 2014-07-14 - Hiram)
 #   before starting downloads, the joinerCheck should be clean
 #   after rn6 is added to all.joiner:
     joinerCheck -keys -database=rn6 all.joiner
 
     cd /hive/data/genomes/rn6
     makeDownloads.pl -dbHost=hgwdev -workhorse=hgwdev rn6 \
 	> downloads.log 2>&1
     # real    28m50.409s
 
     # examine the goldenPath/*/README.txt files to verify the text
 
 ##############################################################################
 # ready for first pushQ entry (DONE - 2014-07-14 - Hiram)
     mkdir /hive/data/genomes/rn6/pushQ
     cd /hive/data/genomes/rn6/pushQ
     makePushQSql.pl rn6 > rn6.sql 2> stderr.out
     # real    1m52.294s
     # some errors are legitimate and OK:
      head stderr.out
 # WARNING: hgwdev does not have /gbdb/rn6/wib/gc5Base.wib
 # WARNING: hgwdev does not have /gbdb/rn6/wib/quality.wib
 # WARNING: hgwdev does not have /gbdb/rn6/bbi/qualityBw/quality.bw
 # WARNING: rn6 does not have seq
 # WARNING: rn6 does not have extFile
 
     scp -p rn6.sql qateam@hgwbeta:/tmp
     ssh qateam@hgwbeta './bin/x86_64/hgsql qapushq < /tmp/rn6.sql'
 
 ##############################################################################
 ## blat server turned on (DONE - 2014-01-13 - Hiram)
 #	After getting a blat server assigned by the Blat Server Gods,
     ssh hgwdev
 
     hgsql -e 'INSERT INTO blatServers (db, host, port, isTrans, canPcr) \
 	VALUES ("rn6", "blat4b", "17784", "1", "0"); \
 	INSERT INTO blatServers (db, host, port, isTrans, canPcr) \
 	VALUES ("rn6", "blat4b", "17785", "0", "1");' \
 	    hgcentraltest
     #	test it with some sequence
 
 ############################################################################
 # LIFTOVER TO Rn5 (DONE - 2014-07-14 - Hiram )
     mkdir /hive/data/genomes/rn6/bed/blat.rn5.2014-07-25
     cd /hive/data/genomes/rn6/bed/blat.rn5.2014-07-25
     # -debug run to create run dir, preview scripts...
     doSameSpeciesLiftOver.pl -ooc=/hive/data/genomes/rn6/jkStuff/rn6.11.ooc \
         -debug rn6 rn5
     # Real run:
     time doSameSpeciesLiftOver.pl \
 	-bigClusterHub=ku -workhorse=hgwdev -dbHost=hgwdev -buildDir=`pwd` \
       -ooc=/hive/data/genomes/rn6/jkStuff/rn6.11.ooc rn6 rn5 > do.log 2>&1
     # real    73m17.426s
 
     # test with sequence
 #############################################################################
 # UCSC to RefSeq name correspondence (DONE - 2015-04-15 - Hiram)
 
     mkdir /hive/data/genomes/rn6/bed/ucscToRefSeq
     cd /hive/data/genomes/rn6/bed/ucscToRefSeq
 
     rsync -avPL \
   rsync://ftp.ncbi.nlm.nih.gov/genomes/genbank/vertebrate_mammalian/Rattus_norvegicus/all_assembly_versions/GCA_000001895.4_Rnor_6.0/GCA_000001895.4_Rnor_6.0_assembly_report.txt ./
 
     # this assembly_report has "UCSC-style-name" in column 10
     # but it does not name everything
 
     # duplicate contigs removed:
       1 AABR07023006.1
       1 AABR07022518.1
 
     # columns 5 and 7 are the INSDC and RefSeq names
     # eliminate the two duplicate contigs
     grep -v "^#" GCA_000001895.4_Rnor_6.0_assembly_report.txt \
       | egrep -v "AABR07023006.1|AABR07022518.1" \
       | awk -F'\t' '{printf "%s\t%s\n", $5,$7}' | sort > insdc.refSeq.tab
 
     hgsql -N -e 'select name,chrom,chromStart,chromEnd from ucscToINSDC;' \
       rn6 | sort > insdc.ucsc.tab
 
     join insdc.ucsc.tab insdc.refSeq.tab | tr '[ ]' '[\t]' \
        | cut -f2- > ucsc.refSeq.tab
 
     # when working perfectly, all these tab files have the same line count:
     wc -l *.tab
 # 953 insdc.refSeq.tab
 # 953 insdc.ucsc.tab
 # 953 ucsc.refSeq.tab
 
     export chrSize=`cut -f1 ucsc.refSeq.tab | awk '{print length($0)}' | sort -n | tail -1`
     sed -e "s/21/$chrSize/" $HOME/kent/src/hg/lib/ucscToINSDC.sql \
        | sed -e 's/INSDC/RefSeq/g;' > ucscToRefSeq.sql
     hgLoadSqlTab rn6 ucscToRefSeq ./ucscToRefSeq.sql ucsc.refSeq.tab
 
     checkTableCoords  rn6 -table=ucscToRefSeq
     # should cover %100 all bases:
     featureBits -countGaps rn6 ucscToRefSeq
     # 2870182909 bases of 2870182909 (100.000%) in intersection
 
 #########################################################################
 
 # download and load ncbiGene track
 
 db=rn6
 mkdir  /cluster/data/genomes/$db/bed/ncbiGene
 cd  /cluster/data/genomes/$db/bed/ncbiGene
 
 ftpFile=ftp://ftp.ncbi.nlm.nih.gov/genomes/R_norvegicus/GFF/ref_Rnor_6.0_top_level.gff3.gz
 gff3File=`basename $ftpFile`
 
 echo "select * from ucscToRefSeq" | hgsql $db | tail -n +2 | awk '{print 0, $4, $3, $1, $3}' > refSeqToUcsc.lft
 rm -f $gff3File
 wget $ftpFile
 
 /cluster/home/braney/bin/x86_64/gff3ToGenePred -useName -warnAndContinue -attrsOut=attrs -bad=bad.gp $gff3File stdout 2> convertErr.txt | liftUp -type=.gp -extGenePred lift.gp refSeqToUcsc.lft warn  stdin 2> liftErr.txt
 wc -l lift.gp
 # 64457 lift.gp
 wc -l bad.gp
 # 0
 
 tawk '{print $1}'  attrs | sort | uniq > meta
 wc -l meta
 # 63278 meta
 for i in product Dbxref gene gbkey
 do
     echo $i
     tawk -v attr=$i '$2==attr {print $1,$3}' attrs | sort | uniq | join -t $'\t' /dev/stdin meta > out
     mv out meta
 done
 wc -l meta
 # 62775 meta
 
 cat curated.gp predicted.gp | awk '{print $1}' | sort -u > tmp1
 cat meta | awk '{print $1}' | sort -u > tmp2
 join -v 1 tmp1 tmp2 | wc -l
 # 0
 
 egrep "^N(M|R)|^YP" lift.gp > curated.gp
 egrep "^X(M|R)" lift.gp > predicted.gp
 
 wc -l curated.gp predicted.gp
 #    17961 curated.gp
 # 43265 predicted.gp
 # 61226 total
 
 grep dropping convertErr.txt | wc -l
 #    0
 
 awk '/isn/ {print $1}' liftErr.txt | sort -u
 # nothing
 
 hgLoadGenePred -genePredExt $db ncbiRefCurated curated.gp
 hgLoadGenePred -genePredExt $db ncbiRefPredicted predicted.gp
 hgLoadSqlTab $db ncbiRefLink $kent/src/hg/lib/ncbiRefLink.sql meta
 
 hgsql -e 'INSERT INTO trackVersion \
     (db, name, who, version, updateTime, comment, source, dateReference)
     VALUES("rn6", "ncbiRefSeq", "braney", "105", now(),
     "http://www.ncbi.nlm.nih.gov/genome/annotation_euk/Rattus_norvegicus/105/",
     "ftp://ftp.ncbi.nlm.nih.gov/genomes/R_norvegicus",
     "7 August 2014" );' hgFixed
 
 ############################################################################
 # SWAP LASTZ human/hg19 vs. rat/rn6 - (DONE - 2015-06-08 - Hiram)
     # the original alignment
     cd /hive/data/genomes/hg19/bed/lastzRn6.2015-06-08
 
     cat fb.hg19.chainRn6Link.txt
     # 924289029 bases of 2897316137 (31.902%) in intersection
 
     # and for the swap:
     mkdir /hive/data/genomes/rn6/bed/blastz.hg19.swap
     cd /hive/data/genomes/rn6/bed/blastz.hg19.swap
 
     time (doBlastzChainNet.pl -verbose=2 \
       /hive/data/genomes/hg19/bed/lastzRn6.2015-06-08/DEF \
         -swap -chainMinScore=3000 -chainLinearGap=medium \
           -workhorse=hgwdev -smallClusterHub=ku -bigClusterHub=ku \
             -syntenicNet) > swap.log 2>&1
     #  real    77m36.480s
 
     cat fb.rn6.chainHg19Link.txt
     # 949880616 bases of 2729860805 (34.796%) in intersection
 
     time (doRecipBest.pl -buildDir=`pwd` rn6 hg19) > rbest.log 2>&1
     # real    32m19.569s
 
 #########################################################################
 # GENEID GENE PREDICTIONS (DONE - 2015-06-26 - Hiram)
     ssh hgwdev
     mkdir /hive/data/genomes/rn6/bed/geneid
     cd /hive/data/genomes/rn6/bed/geneid
     wget --timestamping \
 http://genome.crg.es/genepredictions/R.norvegicus/rn6/geneid_v1.4/rn6.geneid.prot
     wget --timestamping \
 http://genome.crg.es/genepredictions/R.norvegicus/rn6/geneid_v1.4/rn6.geneid.gtf
     ldHgGene -gtf -genePredExt rn6 geneid rn6.geneid.gtf
 
     #  Read 41652 transcripts in 307250 lines in 1 files
     #     41652 groups 883 seqs 1 sources 3 feature types
     #  41652 gene predictions
 
     featureBits -countGaps rn6 geneid
 # 42028722 bases of 2730871774 (1.539%) in intersection
     featureBits -countGaps rn4 geneid
 # 40255245 bases of 2834127293 (1.420%) in intersection
 
 ##########################################################################
 # SGP GENES (DONE - 2015-07-30 - Hiram)
     mkdir /hive/data/genomes/rn6/bed/sgpGene
     cd /hive/data/genomes/rn6/bed/sgpGene
     wget --timestamping \
 http://genome.crg.es/genepredictions/R.norvegicus/rn6/SGP2/hg38/00README
     wget --timestamping \
 http://genome.crg.es/genepredictions/R.norvegicus/rn6/SGP2/hg38/rn6.sgp2.gtf
     wget --timestamping \
 http://genome.crg.es/genepredictions/R.norvegicus/rn6/SGP2/hg38/rn6.sgp2.gff3
 
     ldHgGene -gtf -genePredExt rn6 sgpGene rn6.sgp2.gtf
     # Read 37981 transcripts in 306550 lines in 1 files
     #   37981 groups 468 seqs 1 sources 3 feature types
     # 37981 gene predictions
 
     featureBits -enrichment rn6 refGene:CDS sgpGene
 # refGene:CDS 0.963%, sgpGene 1.399%, both 0.821%, cover 85.27%, enrich 60.93x
 
 ###########################################################################
 # lastz zebrafish danRer10 (DONE - 2015-09-18,14 - Hiram)
     # establish a screen to control this job with a name to indicate what it is
     screen -S rn6DanRer10
     mkdir /hive/data/genomes/rn6/bed/lastzDanRer10.2015-09-18
     cd /hive/data/genomes/rn6/bed/lastzDanRer10.2015-09-18
 
     printf "%s\n" \
 '# Rat vs. zebrafish
 BLASTZ=/cluster/bin/penn/lastz-distrib-1.03.66/bin/lastz
 BLASTZ_Y=3400
 BLASTZ_L=6000
 BLASTZ_K=2200
 BLASTZ_Q=/scratch/data/blastz/HoxD55.q
 
 # TARGET: Rat Rn6
 SEQ1_DIR=/hive/data/genomes/rn6/rn6.2bit
 SEQ1_LEN=/hive/data/genomes/rn6/chrom.sizes
 SEQ1_CHUNK=20000000
 SEQ1_LAP=10000
 SEQ1_LIMIT=15
 
 # QUERY: zebrafish danRer10
 SEQ2_DIR=/hive/data/genomes/danRer10/danRer10.2bit
 SEQ2_LEN=/hive/data/genomes/danRer10/chrom.sizes
 SEQ2_CHUNK=20000000
 SEQ2_LAP=0
 SEQ2_LIMIT=100
 
 BASE=/hive/data/genomes/rn6/bed/lastzDanRer10.2015-09-18
 TMPDIR=/dev/shm' > DEF
 
     # adjust the SEQ2_LIMIT with -stop=partition to get a reasonable
     #	number of jobs, 50,000 to something under 100,000
     # when not present, SEQ2_LIMIT is a default 100
     time (doBlastzChainNet.pl -verbose=2 `pwd`/DEF \
         -workhorse=hgwdev -smallClusterHub=ku -bigClusterHub=ku \
         -chainMinScore=5000 -chainLinearGap=loose) > do.log 2>&1
     #	real    272m42.068s
 
     cat fb.rn6.chainDanRer10Link.txt
     #	88237547 bases of 2729860805 (3.232%) in intersection
 
     time (doBlastzChainNet.pl -verbose=2 `pwd`/DEF -syntenicNet \
         -workhorse=hgwdev -smallClusterHub=ku -bigClusterHub=ku \
         -chainMinScore=5000 -chainLinearGap=loose \
           -continue=syntenicNet) > synNet.log 2>&1 &
     # real    0m50.415s
 
     time (doRecipBest.pl -buildDir=`pwd` rn6 danRer10) > rbest.log 2>&1 &
     #  real    7m47.791s
 
     #	and for the swap
     mkdir /hive/data/genomes/danRer10/bed/blastz.rn6.swap
     cd /hive/data/genomes/danRer10/bed/blastz.rn6.swap
     time (doBlastzChainNet.pl -verbose=2 \
 	/hive/data/genomes/rn6/bed/lastzDanRer10.2015-09-18/DEF \
         -workhorse=hgwdev -smallClusterHub=ku -bigClusterHub=ku \
         -swap -chainMinScore=5000 -chainLinearGap=loose) > swap.log 2>&1
     #	real    20m40.003s
 
     cat  fb.danRer10.chainRn6Link.txt
     #	79649796 bases of 1369683683 (5.815%) in intersection
 
     time (doRecipBest.pl -buildDir=`pwd` danRer10 rn6) > rbest.log 2>&1
     # real    7m59.071s
 
     time (doBlastzChainNet.pl -verbose=2 -syntenicNet \
 	/hive/data/genomes/rn6/bed/lastzDanRer10.2015-09-18/DEF \
         -workhorse=hgwdev -smallClusterHub=ku -bigClusterHub=ku \
         -swap -chainMinScore=5000 -chainLinearGap=loose \
             -continue=syntenicNet) > synNet.log 2>&1 &
     # real    0m41.726s
 
 #########################################################################
 # cloneEnds (DONE - 2016-05-26 - Hiram)
 
     mkdir /hive/data/genomes/rn6/bed/cloneEnds
     cd /hive/data/genomes/rn6/bed/cloneEnds
 
     # fetch the NCBI INSDC name correspondence file:
     rsync -L -a -P rsync://ftp.ncbi.nlm.nih.gov/genomes/refseq/vertebrate_mammalian/Rattus_norvegicus/all_assembly_versions/GCF_000001895.5_Rnor_6.0/GCF_000001895.5_Rnor_6.0_assembly_report.txt ./
     rsync -L -a -P rsync://ftp.ncbi.nlm.nih.gov/genomes/refseq/vertebrate_mammalian/Rattus_norvegicus/all_assembly_versions/GCF_000001895.5_Rnor_6.0/GCF_000001895.5_Rnor_6.0_genomic.fna.gz ./
     faToTwoBit GCF_000001895.5_Rnor_6.0_genomic.fna.gz \
 	GCF_000001895.5_Rnor_6.0.2bit
     twoBitDup -keyList=stdout GCF_000001895.5_Rnor_6.0.2bit \
 	| sort > GCF_000001895.5_Rnor_6.0.idKeys.txt
 
     # fetch the clone reports
     mkdir reports
     rsync -a -P \
 rsync://ftp.ncbi.nih.gov/repository/clone/reports/Rattus_norvegicus/*.GCF_000001895.5.102.*.gff \
        ./reports/
 
     # this procedure to correspond the names is a work in progress
     # turns out the RefSeq assembly has three identical contigs, two dups
     # to be removed, the same ones removed for rn6 build::
     cat remove.list
 NW_007906462.1
 NW_007906649.1
 
     faSomeRecords -exclude GCF_000001895.5_Rnor_6.0_genomic.fna.gz \
       remove.list stdout | gzip -c > noDups.GCF_000001895.5_Rnor_6.0.fa.gz
 
     faToTwoBit noDups.GCF_000001895.5_Rnor_6.0.fa.gz \
       noDups.GCF_000001895.5_Rnor_6.0.2bit
     twoBitDup noDups.GCF_000001895.5_Rnor_6.0.2bit
 
     $HOME/kent/src/utils/twoBitDup/twoBitDup -keyList=stdout \
           ../../rn6.2bit | sort -k1,1 > rn6.idKeys.txt
     $HOME/kent/src/utils/twoBitDup/twoBitDup -keyList=stdout \
           noDups.GCF_000001895.5_Rnor_6.0.2bit | sort -k1,1 \
               > GCF_000001895.5_Rnor_6.0.idKeys.txt
 
     join GCF_000001895.5_Rnor_6.0.idKeys.txt rn6.idKeys.txt \
          | tr '[ ]' '[\t]' | cut -f2- > refSeq.ucscName.tab
 
     # establish full library list:
     ls reports/*.GCF_000001895.5.102.*.gff | sed -e 's#reports/##' \
        | cut -d"." -f1 | sort -u > library.list.txt
     # only four libraries:
     cat library.list.txt
 CH230
 RNB1
 RNB2
 RP32
 
     # most unfortunate, the GFF files use a secondary accession identifier,
     # they do not use the identifier used in the assembly itself.
     cd /hive/data/genomes/rn6/bed/cloneEnds/reports
     grep -h -v "^#" *.gff | cut -f1 | sed -e 's/gpp.//; s/|$//;' \
         | sed '/^$/d' | sort -u | sort > ../gpc.gps.identifiers.txt
     cd /hive/data/genomes/rn6/bed/cloneEnds
     # there are only 75 identifiers:
     wc -l gpc.gps.identifiers.txt
 # 75 gpc.gps.identifiers.txt
     # fetch corresponding identifiers from Entrez:
 export PATH=/hive/data/outside/ncbi/EDirect/edirect:$PATH
 export EMAIL="hiram@soe.ucsc.edu"
 
 for ID in `cat  gpc.gps.identifiers.txt`
 do
   esearch -db nucleotide -query "${ID}" | efetch -format xml \
      | xtract -pattern Seq-id_other -element Textseq-id_accession,Textseq-id_version | awk -v id=$ID '{printf "%s\t%s.%s\n", id, $1,$2}'
 done > gpcGps.to.INSDC.tab
 
     # real    3m30.775s
     # verify have all:
     wc -l gpcGps.to.INSDC.tab gpc.gps.identifiers.txt
 #    75 gpcGps.to.INSDC.tab
 #    75 gpc.gps.identifiers.txt
 
     # rework that refSeq.ucscName.tab file:
     join GCF_000001895.5_Rnor_6.0.idKeys.txt rn6.idKeys.txt \
          | tr '[ ]' '[\t]' | cut -f2- > INSDC.ucscName.tab
     join -1 2 <(sort -k2 gpcGps.to.INSDC.tab) <(sort INSDC.ucscName.tab) \
        | tr '[ ]' '[\t]' | cut -f2- > refSeq.ucscName.tab
 
     # a script to scan the GFF files, with the refSeq.ucscName.tab
     # name correspondence to construct bed files
 
     printf '#!/usr/bin/env perl
 
 use strict;
 use warnings;
 
 my $argc = scalar(@ARGV);
 
 if ($argc < 1) {
   printf STDERR "usage: ./rn6.pl <report.gff> [moreReports.gff]\n";
   exit 255;
 }
 
 my %%refSeqToUcsc;   # key is refSeq name, value is UCSC chrom name
 open (FH, "<refSeq.ucscName.tab") or die "can not read refSeq.ucscName.tab";
 while (my $line = <FH>) {
   chomp $line;
   my ($refSeq, $ucsc) = split("\t", $line);
   $refSeqToUcsc{$refSeq} = $ucsc;
 }
 close (FH);
 
 my %%chromSizes;    # key is UCSC chrom name, key is chrom size
 open (FH, "</hive/data/genomes/rn6/chrom.sizes") or die "can not read rn6/chrom.sizes";
 while (my $line = <FH>) {
   chomp $line;
   my ($chr, $size) = split("\t", $line);
   $chromSizes{$chr} = $size;
 }
 close (FH);
 
 while (my $file = shift) {
 my %%starts;   # key is parent ID, value is start end coordinates start,end
 my %%ends;	# key is parent ID, value is end end coordinates start,end
 my %%parents;	# key is parent ID, value is 1 to signify exists
 my %%endNames;   # key is parent ID, value is the Name of the parent clone_insert
 
 printf STDERR "# processing $file\n";
 
 open (FH, "<$file") or die "can not read $file";
 while (my $line = <FH>) {
   chomp $line;
   next if ($line=~ m/^#/);
   my @a = split("\t", $line);
   next if (scalar(@a) < 1);
   my $contig = $a[0];
   $contig =~ s/gpp//;
   $contig =~ s/\|/g/;
   my $ucscChr = $refSeqToUcsc{$contig};
   if (!defined($ucscChr)) {
     printf STDERR "# ERR: contig not in refSeqToUcsc: %%s\\n", $contig;
     next;
   }
   next if (! exists($chromSizes{$ucscChr}));
   my $chromSize = $chromSizes{$ucscChr};
   my $chromStart = $a[3] - 1;
   my $chromEnd = $a[4];
   if ($chromStart > $chromSize) {
     printf STDERR "# warning chromStart over size $ucscChr $chromStart $chromEnd\\n";
     $chromStart = $chromSize-1;
   }
   if ($chromEnd > $chromSize) {
     my $overRun = $chromEnd - $chromSize;
     printf STDERR "# warning chromEnd over size by $overRun -> $ucscChr $chromStart $chromEnd\\n";
     $chromEnd = $chromSize;
   }
   my $id="notFound";
   my $name="notFound";
   my $parent="notFound";
   my @b = split(";", $a[8]);
   for (my $i = 0; $i < scalar(@b); ++$i) {
      my ($tag, $value) = split("=", $b[$i]);
      if ($tag eq "ID") {
         $id = $value;
         if ($id !~ m/-/) {
           if (exists($parents{$id})) {
             printf STDERR "# WARN: duplicate parent: $id";
           } else {
             $parents{$id} = $ucscChr;
           }
         }
      } elsif ($tag eq "Parent") {
         $parent = $value;
      } elsif ($tag eq "Name") {
         $name = $value;
      }
   }
   my $type="notFound";
   my $insertType = $a[2];
   if ($insertType =~ m/clone_insert_start/) {
      $type = "start";
      if ($parent eq "notFound") {
        printf STDERR "# ERR: can not find parent for start $name Ttype $id\\n";
      } else {
        if (!exists($parents{$parent})) {
          printf STDERR "# ERR: start found $name  with no parent $parent declared\\n";
        } elsif (exists($starts{$parent})) {
          printf STDERR "# ERR: duplicate start for $parent\\n";
        } elsif ($ucscChr eq $parents{$parent}) {
          $starts{$parent} = sprintf("%%s\t%%s", $chromStart, $chromEnd);
        } else {
          printf STDERR "# ERR: start on different chrom $ucscChr than parent $parent $parents{$parent}\\n";
        }
      }
   } elsif ($insertType =~ m/clone_insert_end/) {
      $type = "end";
      if ($parent eq "notFound") {
        printf STDERR "# ERR: can not find parent for end $name Ttype $id\\n";
      } else {
        if (!exists($parents{$parent})) {
          printf STDERR "# ERR: end found $name  with no parent $parent declared\\n";
        } elsif (exists($ends{$parent})) {
          printf STDERR "# ERR: duplicate end for $parent\\n";
        } elsif ($ucscChr eq $parents{$parent}) {
          $ends{$parent} = sprintf("%%s\t%%s", $chromStart, $chromEnd);
        } else {
          printf STDERR "# ERR: end on different chrom $ucscChr than parent $parent $parents{$parent}\\n";
        }
      }
   } elsif ($insertType =~ m/clone_insert/) {
      $type = "insert";
      $endNames{$id} = $name;
   }
   $name =~ s/gi\|//g;
   $id =~ s/gi\|//g;
   printf STDERR "%%s\t%%d\t%%d\t%%s_%%s_%%s\t0\t%%s\\n", $ucscChr, $chromStart, $chromEnd, $name, $type, $id, $a[6];
 }	# while (my $line = <FH>)
 
 close (FH);
 
 foreach my $parent (keys %%parents) {
   if (! exists($starts{$parent}) ) {
     printf STDERR "# ERR: no start for $parent\\n";
   } elsif (! exists($ends{$parent}) ) {
     printf STDERR "# ERR: no end for $parent\\n";
   } else {
     my $strand = "+";
     my $chrStart = 0;
     my $chrEnd = 0;
     my $blockStart = 0;
     my ($sStart, $sEnd) = split("\t", $starts{$parent});
     my ($eStart, $eEnd) = split("\t", $ends{$parent});
     my $startSize = $sEnd - $sStart;
     my $endSize = $eEnd - $eStart;
     if ($eStart < $sStart) {
       $chrStart = $eStart;
       $chrEnd = $sEnd;
       $blockStart = $sStart - $chrStart;
       $strand = "-";
       $startSize = $eEnd - $eStart;
       $endSize = $sEnd - $sStart;
     } else {
       $chrStart = $sStart;
       $chrEnd = $eEnd;
       $blockStart = $eStart - $chrStart;
     }
     if ($startSize > $blockStart) {
       printf STDERR "# startSize > blockStart $endNames{$parent}\\n";
     } else {
       printf "%%s\t%%d\t%%d\t%%s\t0\t%%s\t%%d\t%%d\t0\t2\t%%d,%%d\t0,%%d\\n", $parents{$parent}, $chrStart, $chrEnd, $endNames{$parent}, $strand, $chrStart, $chrEnd, $startSize, $endSize, $blockStart;
     }
   }
 }
 }
 ' > rn6.pl
 
     chmod +x rn6.pl
 
     # process GFF files into bed files into separateLibs/ directory
 for L in `cat library.list.txt`
 do
    export destDir="separateLibs/${L}"
    echo "working: ${L}" 1>&1
    mkdir -p "${destDir}"
    ./rn6.pl reports/${L}.GCF_000001895.5.102.*.gff \
        2> ${destDir}/tmp.bed6 | sort -k1,1 -k2,2n > ${destDir}/rn6.${L}.bed
    sort -k1,1 -k2,2n ${destDir}/tmp.bed6 > ${destDir}/rn6.${L}.items.bed6
 done
 
     # for human clone ends we used only the libraries with more than
     #  20,000 clone ends.  Only lib RP32 here has a small number: 106
     # so use them all
     # use only those libraries with more than 20,000 clone ends
     wc -l separateLibs/*/*.bed | sort -n | grep -v total | awk '$1 > 100' \
         | sed -e 's#.*separateLibs/##; s#/.*##' > libs.over100.list
 
     # filter out bad ends, length must be <= median size times three
     cat libs.over100.list | while read D
 do
    if [ ! -s separateLibs/${D}/lengths.txt ]; then
       awk '{print $3-$2}' separateLibs/${D}/rn6.${D}.bed \
         > separateLibs/${D}/lengths.txt
    fi
    median3X=`ave separateLibs/${D}/lengths.txt | grep median | awk '{printf "%d", $2*3}'`
    awk '($3-$2) < '$median3X'' separateLibs/${D}/rn6.${D}.bed > separateLibs/${D}/rn6.median3X.bed
    awk '($3-$2) >= '$median3X'' separateLibs/${D}/rn6.${D}.bed > separateLibs/${D}/rn6.badMap.bed
    before=`cat separateLibs/${D}/rn6.${D}.bed | wc -l`
    after=`cat separateLibs/${D}/rn6.median3X.bed | wc -l`
    dropped=`echo $before $after | awk '{print $1-$2}'`
    perCent=`echo $dropped $before | awk '{printf "%.2f", 100*'$dropped/$before'}'`
    printf "#\t$D $before - $after = $dropped -> %% $perCent dropped\n"
 done
 
 #       RP32 106 - 105 = 1 -> % 0.94 dropped
 #       CH230 83579 - 78516 = 5063 -> % 6.06 dropped
 #       RNB2 96590 - 89108 = 7482 -> % 7.75 dropped
 #       RNB1 134066 - 124987 = 9079 -> % 6.77 dropped
 
    wc -l separateLibs/*/*.median3X.bed
 #     78516 separateLibs/CH230/rn6.median3X.bed
 #    124987 separateLibs/RNB1/rn6.median3X.bed
 #     89108 separateLibs/RNB2/rn6.median3X.bed
 #       105 separateLibs/RP32/rn6.median3X.bed
 #    292716 total
 
    # loading the median3X files
 for L in `cat libs.over100.list`
 do
     echo $L 1>&2
     hgLoadBed -type=bed12 rn6 cloneEnd${L} \
        separateLibs/${L}/rn6.median3X.bed \
         > separateLibs/loadBed.${L}.log 2>&1
 done
 
    # loading the dropped ends:
    mkdir /hive/data/genomes/rn6/bed/cloneEnds/droppedTooBig
    cd /hive/data/genomes/rn6/bed/cloneEnds/droppedTooBig
    # link them to here
    cat ../libs.over100.list | while read L
 do
   ln -s ../separateLibs/${L}/rn6.badMap.bed ${L}.badMap.bed
 done
   wc -l *.bed
 #    5063 CH230.badMap.bed
 #    9079 RNB1.badMap.bed
 #    7482 RNB2.badMap.bed
 #       1 RP32.badMap.bed
 #   21625 total
 
   # then load
   hgLoadBed -type=bed12 rn6 cloneEndbadEnds *.badMap.bed
   # Read 21625 elements of size 12 from CH230.badMap.bed
 
     # construct multiple mapped ends:
    cd /hive/data/genomes/rn6/bed/cloneEnds
 for L in `cat libs.over100.list`
 do
     cat separateLibs/${L}/rn6.median3X.bed
 done | sort -k4 > allEnds.bed
     wc -l allEnds.bed
 #    292716 allEnds.bed
 
     cut -f4 allEnds.bed | sort | uniq -c | sort -rn > allEnds.names.count.txt
 
     awk '$1 > 1' allEnds.names.count.txt | awk '{print $2}' \
        | sort > multiples.names.txt
 
     join -o "2.1,2.2,2.3,2.4,2.5,2.6,2.7,2.8,2.9,2.10,2.11,2.12" \
        -2 4 multiples.names.txt allEnds.bed | tr '[ ]' '[\t]' \
           | sort -k1,1 -k2,2n  > allEnds.multiple.locations.bed
 
     hgLoadBed -type=bed12 rn6 cloneEndmultipleMaps \
         allEnds.multiple.locations.bed > load.multipleMaps.log 2>&1
 
     awk '$6 == "+"' allEnds.bed | sort -k1,1 -k2,2n \
       | bedItemOverlapCount rn6 stdin > allEnds.forward.bedGraph
 
     awk '$6 == "-"' allEnds.bed | sort -k1,1 -k2,2n \
       | bedItemOverlapCount rn6 stdin > allEnds.reverse.bedGraph
     wc -l *.bedGraph
 #    272228 allEnds.forward.bedGraph
 #    270927 allEnds.reverse.bedGraph
 #    543155 total
 
     bedGraphToBigWig allEnds.forward.bedGraph \
        /hive/data/genomes/rn6/chrom.sizes \
          cloneEndcoverageForward.bw
 
     bedGraphToBigWig allEnds.reverse.bedGraph \
        /hive/data/genomes/rn6/chrom.sizes \
           cloneEndcoverageReverse.bw
 
     bigWigInfo cloneEndcoverageForward.bw
 version: 4
 isCompressed: yes
 isSwapped: 0
 primaryDataSize: 1,359,215
 primaryIndexSize: 21,052
 zoomLevels: 7
 chromCount: 69
 basesCovered: 2,781,456,290
 mean: 8.347160
 min: 1.000000
 max: 141.000000
 std: 4.341902
     bigWigInfo cloneEndcoverageReverse.bw
 version: 4
 isCompressed: yes
 isSwapped: 0
 primaryDataSize: 1,353,539
 primaryIndexSize: 20,988
 zoomLevels: 7
 chromCount: 62
 basesCovered: 2,781,413,397
 mean: 8.322005
 min: 1.000000
 max: 168.000000
 std: 4.457940
 
     mkdir /gbdb/rn6/bbi/cloneEnd
     ln -s `pwd`/cloneEndcoverageForward.bw /gbdb/rn6/bbi/cloneEnd
     ln -s `pwd`/cloneEndcoverageReverse.bw /gbdb/rn6/bbi/cloneEnd
 
     hgBbiDbLink rn6 cloneEndcoverageForward \
         /gbdb/rn6/bbi/cloneEnd/cloneEndcoverageForward.bw
     hgBbiDbLink rn6 cloneEndcoverageReverse \
         /gbdb/rn6/bbi/cloneEnd/cloneEndcoverageReverse.bw
 
     ### Fixup the scores to indicate how many multiple mappings as mentioned
     ### in the hg19 bacEnds description page: one mapping: score = 1000
     ### multiple mappings: score = 1500/count
     ### the sort | uniq -c | awk does this score calculation with the name
     ###   in column 1
     ### The join puts the existing table together with those scores
     ### DONE - 2015-06-18 - Hiram
 
     mkdir /hive/data/genomes/rn6/bed/cloneEnds/addCounts
     cd /hive/data/genomes/rn6/bed/cloneEnds/addCounts
     mkdir score withScore noScore
     for table in cloneEndCH230 cloneEndRNB1 cloneEndRNB2 cloneEndRP32 \
 cloneEndbadEnds cloneEndmultipleMaps
 do
   hgsql -N -e "select name from $table;" rn6 | sort | uniq -c |
       awk '{ if (1 == $1) {printf "%s\t1000\n", $2} else {printf "%s\t%d\n", $2, 1500/$1} }' \
          | sort > score/rn6.$table.score.tab
   hgsql -N -e "select * from $table order by name;" rn6 \
       | sort -k5 > noScore/rn6.$table.tab
   join -1 5 noScore/rn6.$table.tab score/rn6.$table.score.tab \
     | tr '[ ]' '[\t]' | awk '{printf "%d\t%s\t%d\t%d\t%s\t%d\t%s\t%d\t%d\t%d\t%d\t%s\t%s\n", $2,$3,$4,$5,$1,$14,$7,$8,$9,$10,$11,$12,$13}' \
     | sort -k2,2 -k3,3n > withScore/rn6.$table.withScore.tab
   hgsql -e "delete from $table;" rn6
   hgsql -e "load data local infile \"withScore/rn6.$table.withScore.tab\" into table $table;" rn6
 done
 
 ##############################################################################
 # ensemblToRgd (TBD - 2016-08-17 - jcasper)
 # We would like to have otherOrgs links from other species be able to link to
 # RGD pages, but this requires RGD identifiers.  Long term solution involves
 # incorporating that data from GenBank, medium term involves building an RGD
 # Genes track, but short term (this) is to cook up a table that gives RGD
 # IDs for the Ensembl transcript IDs that we have on rn6.
 
 screen -S ensToRgd
 BASEDIR=/hive/data/genomes/rn6/bed/ensToRgd
 mkdir -p ${BASEDIR}
 cd ${BASEDIR}
 
 mkdir rgdDownload
 cd rgdDownload
 wget ftp://ftp.rgd.mcw.edu/pub/data_release/GENES_RAT.txt
 cd ${BASEDIR}
 grep -v '^#' rgdDownload/GENES_RAT.txt | tail -n +2 | tawk '$38 ~ /[A-Z]/ {print $38, $1}' > ensGAndRgd.tab
 perl -ne '($genes, $rgd) = split; foreach $gene (split /;/, $genes) {print "$gene\t$rgd\n";}' < ensGAndRgd.tab | sort > ensGToRgd.tab
 # Halfway there - have Ensembl gene identifier translation, but need
 # transcript identifier translation.  NB: RGD doesn't use version
 # numbers for the Ensembl gene IDs.
 hgsql rn6 -Ne 'select substring_index(gene, ".", 1), transcript from ensGtp' | sort > ensGToEnsT.tab
 # Workaround for join being awkward on tabs:
 TAB=$(printf "\t")
 join -t "$TAB" ensGToEnsT.tab ensGToRgd.tab | tawk '{print $2, $3}' > ensToRgd.tab
 # Use a generic name/value table structure to load the data
 hgLoadSqlTab rn6 ensToRgd ${HOME}/kent/src/hg/lib/knownTo.sql ensToRgd.tab
 
 # Edit hgGene/hgGeneData otherOrgs.ra files for hg38 and mm10 to grab IDs from
 # the ensToRgd table in rn6 and use them to link to RGD.
 
 ##############################################################################
 # lastz X. tropicalis xenTro7 (DONE - 2017-01-11 - Hiram)
     # establish a screen to control this job with a name to indicate what it is
     screen -S rn6XenTro7
     mkdir /hive/data/genomes/rn6/bed/lastzXenTro7.2017-01-11
     cd /hive/data/genomes/rn6/bed/lastzXenTro7.2017-01-11
 
     printf "%s\n" \
 '# Rat vs. X. tropicalis
 BLASTZ=/cluster/bin/penn/lastz-distrib-1.03.66/bin/lastz
 BLASTZ_Y=3400
 BLASTZ_L=6000
 BLASTZ_K=2200
 BLASTZ_Q=/scratch/data/blastz/HoxD55.q
 
 # TARGET: Rat Rn6
 SEQ1_DIR=/hive/data/genomes/rn6/rn6.2bit
 SEQ1_LEN=/hive/data/genomes/rn6/chrom.sizes
 SEQ1_CHUNK=20000000
 SEQ1_LAP=10000
 SEQ1_LIMIT=15
 
 # QUERY: X. tropicalis xenTro7
 SEQ2_DIR=/hive/data/genomes/xenTro7/xenTro7.2bit
 SEQ2_LEN=/hive/data/genomes/xenTro7/chrom.sizes
 SEQ2_CHUNK=20000000
 SEQ2_LAP=0
 SEQ2_LIMIT=20
 
 BASE=/hive/data/genomes/rn6/bed/lastzXenTro7.2017-01-11
 TMPDIR=/dev/shm' > DEF
 
     # adjust the SEQ2_LIMIT with -stop=partition to get a reasonable
     #	number of jobs, 50,000 to something under 100,000
     # when not present, SEQ2_LIMIT is a default 100
     time (doBlastzChainNet.pl -verbose=2 `pwd`/DEF \
         -workhorse=hgwdev -smallClusterHub=ku -bigClusterHub=ku \
         -chainMinScore=5000 -chainLinearGap=loose) > do.log 2>&1
     #	real    1754m46.659s  (much contention for ku resources)
 
     cat fb.rn6.chainXenTro7Link.txt
     #	112036484 bases of 2729860805 (4.104%) in intersection
 
     time (doBlastzChainNet.pl -verbose=2 `pwd`/DEF \
         -workhorse=hgwdev -smallClusterHub=ku -bigClusterHub=ku \
         -continue=syntenicNet -chainMinScore=5000 -chainLinearGap=loose \
 	-syntenicNet) > synNet.log 2>&1 &
     # real    2m23.289s
 
     time (doRecipBest.pl -workhorse=hgwdev -buildDir=`pwd` rn6 xenTro7) \
        > rbest.log 2>&1 &
     #  real    1079m9.975s
 
     #	and for the swap
     mkdir /hive/data/genomes/xenTro7/bed/blastz.rn6.swap
     cd /hive/data/genomes/xenTro7/bed/blastz.rn6.swap
     time (doBlastzChainNet.pl -verbose=2 \
 	/hive/data/genomes/rn6/bed/lastzXenTro7.2017-01-11/DEF \
         -workhorse=hgwdev -smallClusterHub=ku -bigClusterHub=ku \
         -swap -chainMinScore=5000 -chainLinearGap=loose) > swap.log 2>&1
     #	real    66m47.360s
 
     cat  fb.xenTro7.chainRn6Link.txt
     #	108092395 bases of 1365936747 (7.913%) in intersection
 
     time (doRecipBest.pl -workhorse=hgwdev -buildDir=`pwd` xenTro7 rn6) \
          > rbest.log 2>&1
     # real    1027m44.180s
 
     time (doBlastzChainNet.pl -verbose=2 -syntenicNet \
 	/hive/data/genomes/rn6/bed/lastzXenTro7.2017-01-11/DEF \
         -workhorse=hgwdev -smallClusterHub=ku -bigClusterHub=ku \
         -swap -chainMinScore=5000 -chainLinearGap=loose \
             -continue=syntenicNet) > synNet.log 2>&1 &
     #  real    2m23.304s
 
 ##############################################################################
 # lastz chimp panTro5 (DONE - 2017-01-11 - Hiram)
     # establish a screen to control this job with a name to indicate what it is
     screen -S panTro5Rn6
     mkdir /hive/data/genomes/rn6/bed/lastzPanTro5.2017-01-11
     cd /hive/data/genomes/rn6/bed/lastzPanTro5.2017-01-11
 
     printf "%s\n" \
 '# Rat vs. chimp
 BLASTZ=/cluster/bin/penn/lastz-distrib-1.03.66/bin/lastz
 BLASTZ_Y=3400
 BLASTZ_L=6000
 BLASTZ_K=2200
 BLASTZ_Q=/scratch/data/blastz/HoxD55.q
 
 # TARGET: Rat Rn6
 SEQ1_DIR=/hive/data/genomes/rn6/rn6.2bit
 SEQ1_LEN=/hive/data/genomes/rn6/chrom.sizes
 SEQ1_CHUNK=20000000
 SEQ1_LAP=10000
 SEQ1_LIMIT=15
 
 # QUERY: chimp panTro5
 SEQ2_DIR=/hive/data/genomes/panTro5/panTro5.2bit
 SEQ2_LEN=/hive/data/genomes/panTro5/chrom.sizes
 SEQ2_CHUNK=20000000
 SEQ2_LAP=0
 SEQ2_LIMIT=200
 
 BASE=/hive/data/genomes/rn6/bed/lastzPanTro5.2017-01-11
 TMPDIR=/dev/shm' > DEF
 
     # adjust the SEQ2_LIMIT with -stop=partition to get a reasonable
     #	number of jobs, 50,000 to something under 100,000
     # when not present, SEQ2_LIMIT is a default 100
     time (doBlastzChainNet.pl -verbose=2 `pwd`/DEF \
         -workhorse=hgwdev -smallClusterHub=ku -bigClusterHub=ku \
         -chainMinScore=3000 -chainLinearGap=medium) > do.log 2>&1
     #	real    2703m15.450s
 
     cat fb.rn6.chainPanTro5Link.txt
     #	1114217604 bases of 2729860805 (40.816%) in intersection
 
     time (doBlastzChainNet.pl -verbose=2 `pwd`/DEF -syntenicNet \
         -workhorse=hgwdev -smallClusterHub=ku -bigClusterHub=ku \
         -chainMinScore=3000 -chainLinearGap=medium \
           -continue=syntenicNet) > synNet.log 2>&1 &
     #  real    20m31.634s
 
     time (doRecipBest.pl -workhorse=hgwdev -buildDir=`pwd` rn6 panTro5) \
            > rbest.log 2>&1 &
     #  real    1753m33.467s
 
     #	and for the swap
     mkdir /hive/data/genomes/panTro5/bed/blastz.rn6.swap
     cd /hive/data/genomes/panTro5/bed/blastz.rn6.swap
     time (doBlastzChainNet.pl -verbose=2 \
 	/hive/data/genomes/rn6/bed/lastzPanTro5.2017-01-11/DEF \
         -workhorse=hgwdev -smallClusterHub=ku -bigClusterHub=ku \
         -swap -chainMinScore=3000 -chainLinearGap=medium) > swap.log 2>&1
     #	real    185m16.607s
 
     cat  fb.panTro5.chainRn6Link.txt
     #	1107582383 bases of 3132620660 (35.356%) in intersection
 
     time (doRecipBest.pl -workhorse=hgwdev -buildDir=`pwd` panTro5 rn6) \
          > rbest.log 2>&1
     # real    1411m22.158s
 
     time (doBlastzChainNet.pl -verbose=2 -syntenicNet \
 	/hive/data/genomes/rn6/bed/lastzPanTro5.2017-01-11/DEF \
         -workhorse=hgwdev -smallClusterHub=ku -bigClusterHub=ku \
         -swap -chainMinScore=3000 -chainLinearGap=medium \
            -continue=syntenicNet) > synNet.log 2>&1 &
     # real    21m41.393s
 
 ##############################################################################
 # lastz opossum monDom5 (DONE - 2017-01-11 - Hiram)
     # establish a screen to control this job with a name to indicate what it is
     screen -S monDom5Rn6
     mkdir /hive/data/genomes/rn6/bed/lastzMonDom5.2017-01-11
     cd /hive/data/genomes/rn6/bed/lastzMonDom5.2017-01-11
 
     printf "%s\n" \
 '# Rat vs. opossum
 BLASTZ=/cluster/bin/penn/lastz-distrib-1.03.66/bin/lastz
 BLASTZ_M=254
 BLASTZ_Y=3400
 BLASTZ_L=6000
 BLASTZ_K=2200
 BLASTZ_Q=/scratch/data/blastz/HoxD55.q
 
 # TARGET: Rat Rn6
 SEQ1_DIR=/hive/data/genomes/rn6/rn6.2bit
 SEQ1_LEN=/hive/data/genomes/rn6/chrom.sizes
 SEQ1_CHUNK=20000000
 SEQ1_LAP=10000
 SEQ1_LIMIT=5
 
 # QUERY: opossum monDom5
 SEQ2_DIR=/hive/data/genomes/monDom5/monDom5.2bit
 SEQ2_LEN=/hive/data/genomes/monDom5/chrom.sizes
 SEQ2_CHUNK=20000000
 SEQ2_LAP=0
 SEQ2_LIMIT=30
 
 BASE=/hive/data/genomes/rn6/bed/lastzMonDom5.2017-01-11
 TMPDIR=/dev/shm' > DEF
 
     # adjust the SEQ2_LIMIT with -stop=partition to get a reasonable
     #	number of jobs, 50,000 to something under 100,000
     # when not present, SEQ2_LIMIT is a default 100
     time (doBlastzChainNet.pl -verbose=2 `pwd`/DEF \
         -workhorse=hgwdev -smallClusterHub=ku -bigClusterHub=ku \
         -chainMinScore=5000 -chainLinearGap=loose) > do.log 2>&1
     #	real    2258m54.815s  (much contention for ku resources)
 
     cat fb.rn6.chainMonDom5Link.txt
     #	263494827 bases of 2729860805 (9.652%) in intersection
 
     time (doBlastzChainNet.pl -verbose=2 `pwd`/DEF -syntenicNet \
         -workhorse=hgwdev -smallClusterHub=ku -bigClusterHub=ku \
         -continue=syntenicNet -chainMinScore=5000 \
           -chainLinearGap=loose) > synNet.log 2>&1 &
     # real    4m16.551s
 
     time (doRecipBest.pl -workhorse=hgwdev -buildDir=`pwd` rn6 monDom5) \
            > rbest.log 2>&1 &
     #  real    2596m10.639s
 
     #	and for the swap
     mkdir /hive/data/genomes/monDom5/bed/blastz.rn6.swap
     cd /hive/data/genomes/monDom5/bed/blastz.rn6.swap
     time (doBlastzChainNet.pl -verbose=2 \
 	/hive/data/genomes/rn6/bed/lastzMonDom5.2017-01-11/DEF \
         -workhorse=hgwdev -smallClusterHub=ku -bigClusterHub=ku \
         -swap -chainMinScore=5000 -chainLinearGap=loose) > swap.log 2>&1
     #	real    41m5.953s
 
     cat  fb.monDom5.chainRn6Link.txt
     #	251790387 bases of 3501660299 (7.191%) in intersection
 
     time (doRecipBest.pl -workhorse=hgwdev -buildDir=`pwd` monDom5 rn6) \
          > rbest.log 2>&1
     # real    2521m59.900s
 
     time (doBlastzChainNet.pl -verbose=2 -syntenicNet \
 	/hive/data/genomes/rn6/bed/lastzMonDom5.2017-01-11/DEF \
         -workhorse=hgwdev -smallClusterHub=ku -bigClusterHub=ku \
         -swap -chainMinScore=5000 -chainLinearGap=loose \
           -continue=syntenicNet) > synNet.log 2>&1 &
     # real    5m25.250s
 
 ##############################################################################
 # lastz rhesus rheMac8 (DONE - 2017-01-13 - Hiram)
     # establish a screen to control this job with a name to indicate what it is
     screen -S rheMac8Rn6
     mkdir /hive/data/genomes/rn6/bed/lastzRheMac8.2017-01-13
     cd /hive/data/genomes/rn6/bed/lastzRheMac8.2017-01-13
 
     printf "%s\n" \
 '# Rat vs. rhesus
 BLASTZ=/cluster/bin/penn/lastz-distrib-1.03.66/bin/lastz
 BLASTZ_Y=3400
 BLASTZ_L=6000
 BLASTZ_K=2200
 BLASTZ_Q=/scratch/data/blastz/HoxD55.q
 
 # TARGET: Rat Rn6
 SEQ1_DIR=/hive/data/genomes/rn6/rn6.2bit
 SEQ1_LEN=/hive/data/genomes/rn6/chrom.sizes
 SEQ1_CHUNK=20000000
 SEQ1_LAP=10000
 SEQ1_LIMIT=15
 
 # QUERY: rhesus rheMac8
 SEQ2_DIR=/hive/data/genomes/rheMac8/rheMac8.2bit
 SEQ2_LEN=/hive/data/genomes/rheMac8/chrom.sizes
 SEQ2_CHUNK=20000000
 SEQ2_LAP=0
 SEQ2_LIMIT=1000
 
 BASE=/hive/data/genomes/rn6/bed/lastzRheMac8.2017-01-13
 TMPDIR=/dev/shm' > DEF
 
     # adjust the SEQ2_LIMIT with -stop=partition to get a reasonable
     #	number of jobs, 50,000 to something under 100,000
     # when not present, SEQ2_LIMIT is a default 100
     time (doBlastzChainNet.pl -verbose=2 `pwd`/DEF \
         -workhorse=hgwdev -smallClusterHub=ku -bigClusterHub=ku \
         -chainMinScore=3000 -chainLinearGap=medium) > do.log 2>&1
     #	real    1196m51.891s
 
     cat fb.rn6.chainRheMac8Link.txt
     #	1096456425 bases of 2729860805 (40.165%) in intersection
 
     time (doBlastzChainNet.pl -verbose=2 `pwd`/DEF -syntenicNet \
         -workhorse=hgwdev -smallClusterHub=ku -bigClusterHub=ku \
         -chainMinScore=3000 -chainLinearGap=medium \
            -continue=syntenicNet) > synNet.log 2>&1 &
     # real    19m21.810s
 
     time (doRecipBest.pl -workhorse=hgwdev -buildDir=`pwd` rn6 rheMac8) \
            > rbest.log 2>&1 &
     #  real    1569m42.488s
 
     #	and for the swap
     mkdir /hive/data/genomes/rheMac8/bed/blastz.rn6.swap
     cd /hive/data/genomes/rheMac8/bed/blastz.rn6.swap
     time (doBlastzChainNet.pl -verbose=2 \
 	/hive/data/genomes/rn6/bed/lastzRheMac8.2017-01-13/DEF \
         -workhorse=hgwdev -smallClusterHub=ku -bigClusterHub=ku \
         -swap -chainMinScore=3000 -chainLinearGap=medium) > swap.log 2>&1
     #	real    205m19.557s
     # lost swap.log
 
     cat  fb.rheMac8.chainRn6Link.txt
     #	1051607110 bases of 3142093174 (33.468%) in intersection
 
     time (doRecipBest.pl -workhorse=hgwdev -buildDir=`pwd` rheMac8 rn6) \
          > rbest.log 2>&1
     # real    1601m35.848s
 
     time (doBlastzChainNet.pl -verbose=2 -syntenicNet \
 	/hive/data/genomes/rn6/bed/lastzRheMac8.2017-01-13/DEF \
         -workhorse=hgwdev -smallClusterHub=ku -bigClusterHub=ku \
         -swap -chainMinScore=3000 -chainLinearGap=medium \
             -continue=syntenicNet) > swap.log 2>&1 &
     # real    19m14.673s
 
 ##############################################################################
 # lastz guinea pig cavPor3 (DONE - 2017-01-13 - Hiram)
     # establish a screen to control this job with a name to indicate what it is
     screen -S cavPor3Rn6
     mkdir /hive/data/genomes/rn6/bed/lastzCavPor3.2017-01-13
     cd /hive/data/genomes/rn6/bed/lastzCavPor3.2017-01-13
 
     printf "%s\n" \
 '# Rat vs. guinea pig
 BLASTZ=/cluster/bin/penn/lastz-distrib-1.03.66/bin/lastz
 
 # TARGET: Rat Rn6
 SEQ1_DIR=/hive/data/genomes/rn6/rn6.2bit
 SEQ1_LEN=/hive/data/genomes/rn6/chrom.sizes
 SEQ1_CHUNK=20000000
 SEQ1_LAP=10000
 SEQ1_LIMIT=15
 
 # QUERY: guinea pig cavPor3
 SEQ2_DIR=/hive/data/genomes/cavPor3/cavPor3.2bit
 SEQ2_LEN=/hive/data/genomes/cavPor3/chrom.sizes
 SEQ2_CHUNK=20000000
 SEQ2_LAP=0
 SEQ2_LIMIT=15
 
 BASE=/hive/data/genomes/rn6/bed/lastzCavPor3.2017-01-13
 TMPDIR=/dev/shm' > DEF
 
     # adjust the SEQ2_LIMIT with -stop=partition to get a reasonable
     #	number of jobs, 50,000 to something under 100,000
     # when not present, SEQ2_LIMIT is a default 100
     time (doBlastzChainNet.pl -verbose=2 `pwd`/DEF \
         -workhorse=hgwdev -smallClusterHub=ku -bigClusterHub=ku \
         -chainMinScore=3000 -chainLinearGap=medium) > do.log 2>&1
     #	real    742m21.264s
 
     cat fb.rn6.chainCavPor3Link.txt
     #	759825883 bases of 2729860805 (27.834%) in intersection
 
     time (doBlastzChainNet.pl -verbose=2 `pwd`/DEF -syntenicNet \
         -workhorse=hgwdev -smallClusterHub=ku -bigClusterHub=ku \
         -chainMinScore=3000 -chainLinearGap=medium \
           -continue=syntenicNet) > synNet.log 2>&1 &
     # real    12m11.706s
 
     time (doRecipBest.pl -workhorse=hgwdev -buildDir=`pwd` rn6 cavPor3) \
            > rbest.log 2>&1 &
     #  real    535m40.395s
 
     #	and for the swap
     mkdir /hive/data/genomes/cavPor3/bed/blastz.rn6.swap
     cd /hive/data/genomes/cavPor3/bed/blastz.rn6.swap
     time (doBlastzChainNet.pl -verbose=2 \
 	/hive/data/genomes/rn6/bed/lastzCavPor3.2017-01-13/DEF \
         -workhorse=hgwdev -smallClusterHub=ku -bigClusterHub=ku \
         -swap -chainMinScore=3000 -chainLinearGap=medium) > swap.log 2>&1
     #	real    54m21.781s
 
     cat  fb.cavPor3.chainRn6Link.txt
     #	755393931 bases of 2663369733 (28.362%) in intersection
 
     time (doRecipBest.pl -workhorse=hgwdev -buildDir=`pwd` cavPor3 rn6) \
          > rbest.log 2>&1
     # real    564m23.551s
 
     time (doBlastzChainNet.pl -verbose=2 -syntenicNet \
 	/hive/data/genomes/rn6/bed/lastzCavPor3.2017-01-13/DEF \
         -workhorse=hgwdev -smallClusterHub=ku -bigClusterHub=ku \
         -swap -chainMinScore=3000 -chainLinearGap=medium \
           -continue=syntenicNet) > synNet.log 2>&1 &
     # real    12m6.911s
 
 ##############################################################################
 # lastz dog canFam3 (DONE - 2017-01-13 - Hiram)
     # establish a screen to control this job with a name to indicate what it is
     screen -S canFam3Rn6
     mkdir /hive/data/genomes/rn6/bed/lastzCanFam3.2017-01-13
     cd /hive/data/genomes/rn6/bed/lastzCanFam3.2017-01-13
 
     printf "%s\n" \
 '# Rat vs. dog
 BLASTZ=/cluster/bin/penn/lastz-distrib-1.03.66/bin/lastz
 
 # TARGET: Rat Rn6
 SEQ1_DIR=/hive/data/genomes/rn6/rn6.2bit
 SEQ1_LEN=/hive/data/genomes/rn6/chrom.sizes
 SEQ1_CHUNK=20000000
 SEQ1_LAP=10000
 SEQ1_LIMIT=15
 
 # QUERY: dog canFam3
 SEQ2_DIR=/hive/data/genomes/canFam3/canFam3.2bit
 SEQ2_LEN=/hive/data/genomes/canFam3/chrom.sizes
 SEQ2_CHUNK=20000000
 SEQ2_LAP=0
 SEQ2_LIMIT=15
 
 BASE=/hive/data/genomes/rn6/bed/lastzCanFam3.2017-01-13
 TMPDIR=/dev/shm' > DEF
 
     # adjust the SEQ2_LIMIT with -stop=partition to get a reasonable
     #	number of jobs, 50,000 to something under 100,000
     # when not present, SEQ2_LIMIT is a default 100
     time (doBlastzChainNet.pl -verbose=2 `pwd`/DEF \
         -workhorse=hgwdev -smallClusterHub=ku -bigClusterHub=ku \
         -chainMinScore=3000 -chainLinearGap=medium) > do.log 2>&1
     #	real    694m15.665s
 
     cat fb.rn6.chainCanFam3Link.txt
     #	779099966 bases of 2729860805 (28.540%) in intersection
 
     time (doBlastzChainNet.pl -verbose=2 `pwd`/DEF -syntenicNet \
         -workhorse=hgwdev -smallClusterHub=ku -bigClusterHub=ku \
         -chainMinScore=3000 -chainLinearGap=medium \
           -continue=syntenicNet) > synNet.log 2>&1 &
     # real    12m27.782s
 
     time (doRecipBest.pl -workhorse=hgwdev -buildDir=`pwd` rn6 canFam3) \
            > rbest.log 2>&1 &
     #  real    471m35.128s
 
     #	and for the swap
     mkdir /hive/data/genomes/canFam3/bed/blastz.rn6.swap
     cd /hive/data/genomes/canFam3/bed/blastz.rn6.swap
     time (doBlastzChainNet.pl -verbose=2 \
 	/hive/data/genomes/rn6/bed/lastzCanFam3.2017-01-13/DEF \
         -workhorse=hgwdev -smallClusterHub=ku -bigClusterHub=ku \
         -swap -chainMinScore=3000 -chainLinearGap=medium) > swap.log 2>&1
     #	real    47m17.732s
 
     cat  fb.canFam3.chainRn6Link.txt
     #	738483963 bases of 2392715236 (30.864%) in intersection
 
     time (doRecipBest.pl -workhorse=hgwdev -buildDir=`pwd` canFam3 rn6) \
          > rbest.log 2>&1
     # real    357m54.846s
 
     time (doBlastzChainNet.pl -verbose=2 -syntenicNet \
 	/hive/data/genomes/rn6/bed/lastzCanFam3.2017-01-13/DEF \
         -workhorse=hgwdev -smallClusterHub=ku -bigClusterHub=ku \
         -swap -chainMinScore=3000 -chainLinearGap=medium \
           -continue=syntenicNet) > synNet.log 2>&1 &
     # real    12m37.764s
 
 ##############################################################################
 # lastz cow bosTau8 (DONE - 2017-01-16 - Hiram)
     # establish a screen to control this job with a name to indicate what it is
     screen -S bosTau8Rn6
     mkdir /hive/data/genomes/rn6/bed/lastzBosTau8.2017-01-16
     cd /hive/data/genomes/rn6/bed/lastzBosTau8.2017-01-16
 
     printf "%s\n" \
 '# Rat vs. cow
 BLASTZ=/cluster/bin/penn/lastz-distrib-1.03.66/bin/lastz
 
 # TARGET: Rat Rn6
 SEQ1_DIR=/hive/data/genomes/rn6/rn6.2bit
 SEQ1_LEN=/hive/data/genomes/rn6/chrom.sizes
 SEQ1_CHUNK=20000000
 SEQ1_LAP=10000
 SEQ1_LIMIT=15
 
 # QUERY: cow bosTau8
 SEQ2_DIR=/hive/data/genomes/bosTau8/bosTau8.2bit
 SEQ2_LEN=/hive/data/genomes/bosTau8/chrom.sizes
 SEQ2_CHUNK=20000000
 SEQ2_LAP=0
 SEQ2_LIMIT=15
 
 BASE=/hive/data/genomes/rn6/bed/lastzBosTau8.2017-01-16
 TMPDIR=/dev/shm' > DEF
 
     # adjust the SEQ2_LIMIT with -stop=partition to get a reasonable
     #	number of jobs, 50,000 to something under 100,000
     # when not present, SEQ2_LIMIT is a default 100
     time (doBlastzChainNet.pl -verbose=2 `pwd`/DEF \
         -workhorse=hgwdev -smallClusterHub=ku -bigClusterHub=ku \
         -chainMinScore=3000 -chainLinearGap=medium) > do.log 2>&1
     #	real    379m56.864s
 
     cat fb.rn6.chainBosTau8Link.txt
     #	707100078 bases of 2729860805 (25.902%) in intersection
 
     time (doBlastzChainNet.pl -verbose=2 `pwd`/DEF -syntenicNet \
         -workhorse=hgwdev -smallClusterHub=ku -bigClusterHub=ku \
         -chainMinScore=3000 -chainLinearGap=medium \
             -continue=syntenicNet) > synNet.log 2>&1 &
     # real    12m37.764s
 
     time (doRecipBest.pl -workhorse=hgwdev -buildDir=`pwd` rn6 bosTau8) \
            > rbest.log 2>&1 &
     #  real    688m13.841s
 
     #	and for the swap
     mkdir /hive/data/genomes/bosTau8/bed/blastz.rn6.swap
     cd /hive/data/genomes/bosTau8/bed/blastz.rn6.swap
     time (doBlastzChainNet.pl -verbose=2 \
 	/hive/data/genomes/rn6/bed/lastzBosTau8.2017-01-16/DEF \
         -workhorse=hgwdev -smallClusterHub=ku -bigClusterHub=ku \
         -swap -chainMinScore=3000 -chainLinearGap=medium) > swap.log 2>&1
     #	real    53m51.901s
 
     cat  fb.bosTau8.chainRn6Link.txt
     #	671724400 bases of 2649307237 (25.355%) in intersection
 
     time (doRecipBest.pl -workhorse=hgwdev -buildDir=`pwd` bosTau8 rn6) \
          > rbest.log 2>&1
     # real    580m38.055s
 
     time (doBlastzChainNet.pl -verbose=2 -syntenicNet \
 	/hive/data/genomes/rn6/bed/lastzBosTau8.2017-01-16/DEF \
         -workhorse=hgwdev -smallClusterHub=ku -bigClusterHub=ku \
         -swap -chainMinScore=3000 -chainLinearGap=medium \
           -continue=syntenicNet) > synNet.log 2>&1 &
     #  real    11m19.902s
 
 ##############################################################################
 # lastz panda ailMel1 (DONE - 2017-01-16 - Hiram)
     # establish a screen to control this job with a name to indicate what it is
     screen -S ailMel1Rn6
     mkdir /hive/data/genomes/rn6/bed/lastzAilMel1.2017-01-16
     cd /hive/data/genomes/rn6/bed/lastzAilMel1.2017-01-16
 
     printf "%s\n" \
 '# Rat vs. panda
 BLASTZ=/cluster/bin/penn/lastz-distrib-1.03.66/bin/lastz
 
 # TARGET: Rat Rn6
 SEQ1_DIR=/hive/data/genomes/rn6/rn6.2bit
 SEQ1_LEN=/hive/data/genomes/rn6/chrom.sizes
 SEQ1_CHUNK=20000000
 SEQ1_LAP=10000
 SEQ1_LIMIT=15
 
 # QUERY: panda ailMel1
 SEQ2_DIR=/hive/data/genomes/ailMel1/ailMel1.2bit
 SEQ2_LEN=/hive/data/genomes/ailMel1/chrom.sizes
 SEQ2_CHUNK=20000000
 SEQ2_LAP=0
 SEQ2_LIMIT=400
 
 BASE=/hive/data/genomes/rn6/bed/lastzAilMel1.2017-01-16
 TMPDIR=/dev/shm' > DEF
 
     # adjust the SEQ2_LIMIT with -stop=partition to get a reasonable
     #	number of jobs, 50,000 to something under 100,000
     # when not present, SEQ2_LIMIT is a default 100
     time (doBlastzChainNet.pl -verbose=2 `pwd`/DEF \
         -workhorse=hgwdev -smallClusterHub=ku -bigClusterHub=ku \
         -chainMinScore=3000 -chainLinearGap=medium) > do.log 2>&1
     #	real    309m30.719s
 
     cat fb.rn6.chainAilMel1Link.txt
     #	829980217 bases of 2729860805 (30.404%) in intersection
 
     time (doBlastzChainNet.pl -verbose=2 `pwd`/DEF -syntenicNet \
         -workhorse=hgwdev -smallClusterHub=ku -bigClusterHub=ku \
         -chainMinScore=3000 -chainLinearGap=medium \
            -continue=syntenicNet) > synNet.log 2>&1 &
     #  real    12m19.496s
 
     time (doRecipBest.pl -workhorse=hgwdev -buildDir=`pwd` rn6 ailMel1) \
            > rbest.log 2>&1 &
     #  real    405m31.270s
 
     #	and for the swap
     mkdir /hive/data/genomes/ailMel1/bed/blastz.rn6.swap
     cd /hive/data/genomes/ailMel1/bed/blastz.rn6.swap
     time (doBlastzChainNet.pl -verbose=2 \
 	/hive/data/genomes/rn6/bed/lastzAilMel1.2017-01-16/DEF \
         -workhorse=hgwdev -smallClusterHub=ku -bigClusterHub=ku \
         -swap -chainMinScore=3000 -chainLinearGap=medium) > swap.log 2>&1
     #	real    53m19.878s
 
     cat  fb.ailMel1.chainRn6Link.txt
     #	780846350 bases of 2245312831 (34.777%) in intersection
 
     time (doRecipBest.pl -workhorse=hgwdev -buildDir=`pwd` ailMel1 rn6) \
          > rbest.log 2>&1
     # real    324m22.903s
 
     time (doBlastzChainNet.pl -verbose=2 -syntenicNet \
 	/hive/data/genomes/rn6/bed/lastzAilMel1.2017-01-16/DEF \
         -workhorse=hgwdev -smallClusterHub=ku -bigClusterHub=ku \
         -swap -chainMinScore=3000 -chainLinearGap=medium \
           -continue=syntenicNet) > synNet.log 2>&1 &
     # real    14m19.950s
 
 ##############################################################################
 # lastz prairie vole micOch1 (DONE - 2017-01-17 - Hiram)
     # establish a screen to control this job with a name to indicate what it is
     screen -S micOch1Rn6
     mkdir /hive/data/genomes/rn6/bed/lastzMicOch1.2017-01-17
     cd /hive/data/genomes/rn6/bed/lastzMicOch1.2017-01-17
 
     printf "%s\n" \
 '# Rat vs. prairie vole
 BLASTZ=/cluster/bin/penn/lastz-distrib-1.03.66/bin/lastz
 
 # TARGET: Rat Rn6
 SEQ1_DIR=/hive/data/genomes/rn6/rn6.2bit
 SEQ1_LEN=/hive/data/genomes/rn6/chrom.sizes
 SEQ1_CHUNK=20000000
 SEQ1_LAP=10000
 SEQ1_LIMIT=15
 
 # QUERY: prairie vole micOch1
 SEQ2_DIR=/hive/data/genomes/micOch1/micOch1.2bit
 SEQ2_LEN=/hive/data/genomes/micOch1/chrom.sizes
 SEQ2_CHUNK=20000000
 SEQ2_LAP=0
 SEQ2_LIMIT=20
 
 BASE=/hive/data/genomes/rn6/bed/lastzMicOch1.2017-01-17
 TMPDIR=/dev/shm' > DEF
 
     # adjust the SEQ2_LIMIT with -stop=partition to get a reasonable
     #	number of jobs, 50,000 to something under 100,000
     # when not present, SEQ2_LIMIT is a default 100
     time (doBlastzChainNet.pl -verbose=2 `pwd`/DEF \
         -workhorse=hgwdev -smallClusterHub=ku -bigClusterHub=ku \
         -chainMinScore=3000 -chainLinearGap=medium) > do.log 2>&1
     #	real    1009m6.065s
 
     cat fb.rn6.chainMicOch1Link.txt
     #	1561627848 bases of 2729860805 (57.205%) in intersection
 
     time (doBlastzChainNet.pl -verbose=2 `pwd`/DEF -syntenicNet \
         -workhorse=hgwdev -smallClusterHub=ku -bigClusterHub=ku \
         -chainMinScore=3000 -chainLinearGap=medium \
           -continue=syntenicNet) > synNet.log 2>&1 &
     # real    26m19.328s
 
     time (doRecipBest.pl -workhorse=hgwdev -buildDir=`pwd` rn6 micOch1) \
            > rbest.log 2>&1 &
     #  real    1213m36.203s
 
     #	and for the swap
     mkdir /hive/data/genomes/micOch1/bed/blastz.rn6.swap
     cd /hive/data/genomes/micOch1/bed/blastz.rn6.swap
     time (doBlastzChainNet.pl -verbose=2 \
 	/hive/data/genomes/rn6/bed/lastzMicOch1.2017-01-17/DEF \
         -workhorse=hgwdev -smallClusterHub=ku -bigClusterHub=ku \
         -swap -chainMinScore=3000 -chainLinearGap=medium) > swap.log 2>&1
     #	real    138m3.307s
 
     cat  fb.micOch1.chainRn6Link.txt
     #	1475851323 bases of 2104321675 (70.134%) in intersection
 
     time (doRecipBest.pl -workhorse=hgwdev -buildDir=`pwd` micOch1 rn6) \
          > rbest.log 2>&1
     # real    1096m12.914s
 
     time (doBlastzChainNet.pl -verbose=2 -syntenicNet \
 	/hive/data/genomes/rn6/bed/lastzMicOch1.2017-01-17/DEF \
         -workhorse=hgwdev -smallClusterHub=ku -bigClusterHub=ku \
         -swap -chainMinScore=3000 -chainLinearGap=medium \
             -continue=syntenicNet) > synNet.log 2>&1 &
     # real    24m30.723s
 
 ##############################################################################
 # lastz rabbit oryCun2 (DONE - 2017-01-17 - Hiram)
     # establish a screen to control this job with a name to indicate what it is
     screen -S oryCun2Rn6
     mkdir /hive/data/genomes/rn6/bed/lastzOryCun2.2017-01-17
     cd /hive/data/genomes/rn6/bed/lastzOryCun2.2017-01-17
 
     printf "%s\n" \
 '# Rat vs. rabbit
 BLASTZ=/cluster/bin/penn/lastz-distrib-1.03.66/bin/lastz
 
 # TARGET: Rat Rn6
 SEQ1_DIR=/hive/data/genomes/rn6/rn6.2bit
 SEQ1_LEN=/hive/data/genomes/rn6/chrom.sizes
 SEQ1_CHUNK=20000000
 SEQ1_LAP=10000
 SEQ1_LIMIT=15
 
 # QUERY: rabbit oryCun2
 SEQ2_DIR=/hive/data/genomes/oryCun2/oryCun2.2bit
 SEQ2_LEN=/hive/data/genomes/oryCun2/chrom.sizes
 SEQ2_CHUNK=20000000
 SEQ2_LAP=0
 SEQ2_LIMIT=20
 
 BASE=/hive/data/genomes/rn6/bed/lastzOryCun2.2017-01-17
 TMPDIR=/dev/shm' > DEF
 
     # adjust the SEQ2_LIMIT with -stop=partition to get a reasonable
     #	number of jobs, 50,000 to something under 100,000
     # when not present, SEQ2_LIMIT is a default 100
     time (doBlastzChainNet.pl -verbose=2 `pwd`/DEF \
         -workhorse=hgwdev -smallClusterHub=ku -bigClusterHub=ku \
         -chainMinScore=3000 -chainLinearGap=medium) > do.log 2>&1
     #	real    722m23.231s
 
     cat fb.rn6.chainOryCun2Link.txt
     #	673020370 bases of 2729860805 (24.654%) in intersection
 
     time (doBlastzChainNet.pl -verbose=2 `pwd`/DEF -syntenicNet \
         -workhorse=hgwdev -smallClusterHub=ku -bigClusterHub=ku \
         -chainMinScore=3000 -chainLinearGap=medium \
            -continue=syntenicNet) > synNet.log 2>&1 &
     # real    10m40.444s
 
     time (doRecipBest.pl -workhorse=hgwdev -buildDir=`pwd` rn6 oryCun2) \
            > rbest.log 2>&1 &
     #  real    701m28.973s
 
     #	and for the swap
     mkdir /hive/data/genomes/oryCun2/bed/blastz.rn6.swap
     cd /hive/data/genomes/oryCun2/bed/blastz.rn6.swap
     time (doBlastzChainNet.pl -verbose=2 \
 	/hive/data/genomes/rn6/bed/lastzOryCun2.2017-01-17/DEF \
         -workhorse=hgwdev -smallClusterHub=ku -bigClusterHub=ku \
         -swap -chainMinScore=3000 -chainLinearGap=medium) > swap.log 2>&1
     #	real    45m34.526s
 
     cat  fb.oryCun2.chainRn6Link.txt
     #	650509122 bases of 2604023284 (24.981%) in intersection
 
     time (doRecipBest.pl -workhorse=hgwdev -buildDir=`pwd` oryCun2 rn6) \
          > rbest.log 2>&1
     # real    641m54.166s
 
     time (doBlastzChainNet.pl -verbose=2 -syntenicNet \
 	/hive/data/genomes/rn6/bed/lastzOryCun2.2017-01-17/DEF \
         -workhorse=hgwdev -smallClusterHub=ku -bigClusterHub=ku \
         -swap -chainMinScore=3000 -chainLinearGap=medium \
           -continue=syntenicNet) > synNet.log 2>&1 &
     # real    10m40.444s
 
 ##############################################################################
 # lastz cat felCat8 (DONE - 2017-01-17 - Hiram)
     # establish a screen to control this job with a name to indicate what it is
     screen -S felCat8Rn6
     mkdir /hive/data/genomes/rn6/bed/lastzFelCat8.2017-01-17
     cd /hive/data/genomes/rn6/bed/lastzFelCat8.2017-01-17
 
     printf "%s\n" \
 '# Rat vs. cat
 BLASTZ=/cluster/bin/penn/lastz-distrib-1.03.66/bin/lastz
 
 # TARGET: Rat Rn6
 SEQ1_DIR=/hive/data/genomes/rn6/rn6.2bit
 SEQ1_LEN=/hive/data/genomes/rn6/chrom.sizes
 SEQ1_CHUNK=20000000
 SEQ1_LAP=10000
 SEQ1_LIMIT=15
 
 # QUERY: cat felCat8
 SEQ2_DIR=/hive/data/genomes/felCat8/felCat8.2bit
 SEQ2_LEN=/hive/data/genomes/felCat8/chrom.sizes
 SEQ2_CHUNK=20000000
 SEQ2_LAP=0
 SEQ2_LIMIT=1000
 
 BASE=/hive/data/genomes/rn6/bed/lastzFelCat8.2017-01-17
 TMPDIR=/dev/shm' > DEF
 
     # adjust the SEQ2_LIMIT with -stop=partition to get a reasonable
     #	number of jobs, 50,000 to something under 100,000
     # when not present, SEQ2_LIMIT is a default 100
     time (doBlastzChainNet.pl -verbose=2 `pwd`/DEF \
         -workhorse=hgwdev -smallClusterHub=ku -bigClusterHub=ku \
         -chainMinScore=3000 -chainLinearGap=medium) > do.log 2>&1
     #	real    830m52.261s
 
     cat fb.rn6.chainFelCat8Link.txt
     #	805848215 bases of 2729860805 (29.520%) in intersection
 
     time (doBlastzChainNet.pl -verbose=2 `pwd`/DEF -syntenicNet \
         -workhorse=hgwdev -smallClusterHub=ku -bigClusterHub=ku \
         -chainMinScore=3000 -chainLinearGap=medium \
         -continue=syntenicNet) > synNet.log 2>&1 &
     # real    13m49.129s
 
     time (doRecipBest.pl -workhorse=hgwdev -buildDir=`pwd` rn6 felCat8) \
            > rbest.log 2>&1 &
     #  real    648m55.901s
 
     #	and for the swap
     mkdir /hive/data/genomes/felCat8/bed/blastz.rn6.swap
     cd /hive/data/genomes/felCat8/bed/blastz.rn6.swap
     time (doBlastzChainNet.pl -verbose=2 \
 	/hive/data/genomes/rn6/bed/lastzFelCat8.2017-01-17/DEF \
         -workhorse=hgwdev -smallClusterHub=ku -bigClusterHub=ku \
         -swap -chainMinScore=3000 -chainLinearGap=medium) > swap.log 2>&1
     #	real    112m39.614s
 
     cat  fb.felCat8.chainRn6Link.txt
     #	771419833 bases of 2599716822 (29.673%) in intersection
 
     time (doRecipBest.pl -workhorse=hgwdev -buildDir=`pwd` felCat8 rn6) \
          > rbest.log 2>&1
     # real    577m31.945s
 
     time (doBlastzChainNet.pl -verbose=2 -syntenicNet \
 	/hive/data/genomes/rn6/bed/lastzFelCat8.2017-01-17/DEF \
         -workhorse=hgwdev -smallClusterHub=ku -bigClusterHub=ku \
         -swap -chainMinScore=3000 -chainLinearGap=medium \
           -continue=syntenicNet) > synNet.log 2>&1 &
     # real    13m51.296s
 
 ##############################################################################
 # lastz turkey melGal5 (DONE - 2017-01-17 - Hiram)
     # establish a screen to control this job with a name to indicate what it is
     screen -S melGal5Rn6
     mkdir /hive/data/genomes/rn6/bed/lastzMelGal5.2017-01-17
     cd /hive/data/genomes/rn6/bed/lastzMelGal5.2017-01-17
 
     printf "%s\n" \
 '# Rat vs. turkey
 BLASTZ=/cluster/bin/penn/lastz-distrib-1.03.66/bin/lastz
 
 # TARGET: Rat Rn6
 SEQ1_DIR=/hive/data/genomes/rn6/rn6.2bit
 SEQ1_LEN=/hive/data/genomes/rn6/chrom.sizes
 SEQ1_CHUNK=20000000
 SEQ1_LAP=10000
 SEQ1_LIMIT=15
 
 # QUERY: turkey melGal5
 SEQ2_DIR=/hive/data/genomes/melGal5/melGal5.2bit
 SEQ2_LEN=/hive/data/genomes/melGal5/chrom.sizes
 SEQ2_CHUNK=20000000
 SEQ2_LAP=0
 SEQ2_LIMIT=600
 
 BASE=/hive/data/genomes/rn6/bed/lastzMelGal5.2017-01-17
 TMPDIR=/dev/shm' > DEF
 
     # adjust the SEQ2_LIMIT with -stop=partition to get a reasonable
     #	number of jobs, 50,000 to something under 100,000
     # when not present, SEQ2_LIMIT is a default 100
     time (doBlastzChainNet.pl -verbose=2 `pwd`/DEF \
         -workhorse=hgwdev -smallClusterHub=ku -bigClusterHub=ku \
         -chainMinScore=5000 -chainLinearGap=loose) > do.log 2>&1
     #	real    504m54.176s
 
     cat fb.rn6.chainMelGal5Link.txt
     #	66617232 bases of 2729860805 (2.440%) in intersection
 
     time (doBlastzChainNet.pl -verbose=2 `pwd`/DEF -syntenicNet \
         -workhorse=hgwdev -smallClusterHub=ku -bigClusterHub=ku \
         -chainMinScore=5000 -chainLinearGap=loose \
           -continue=syntenicNet) > synNet.log 2>&1 &
     # real    0m52.732s
 
     time (doRecipBest.pl -workhorse=hgwdev -buildDir=`pwd` rn6 melGal5) \
            > rbest.log 2>&1 &
     #  real    312m57.146s
     # lost the log
 
     #	and for the swap
     mkdir /hive/data/genomes/melGal5/bed/blastz.rn6.swap
     cd /hive/data/genomes/melGal5/bed/blastz.rn6.swap
     time (doBlastzChainNet.pl -verbose=2 \
 	/hive/data/genomes/rn6/bed/lastzMelGal5.2017-01-17/DEF \
         -workhorse=hgwdev -smallClusterHub=ku -bigClusterHub=ku \
         -swap -chainMinScore=5000 -chainLinearGap=loose) > swap.log 2>&1
     #	real    22m31.971s
 
     cat  fb.melGal5.chainRn6Link.txt
     #	53668356 bases of 1093044709 (4.910%) in intersection
 
     time (doRecipBest.pl -workhorse=hgwdev -buildDir=`pwd` melGal5 rn6) \
          > rbest.log 2>&1
     # real    243m51.163s
 
     time (doBlastzChainNet.pl -verbose=2 -syntenicNet \
 	/hive/data/genomes/rn6/bed/lastzMelGal5.2017-01-17/DEF \
         -workhorse=hgwdev -smallClusterHub=ku -bigClusterHub=ku \
         -swap -chainMinScore=5000 -chainLinearGap=loose \
            -continue=syntenicNet) > synNet.log 2>&1 &
     # real    0m45.406s
 
 ##############################################################################
 # lastz Elephant shark calMil1 (DONE - 2018-01-18 - Hiram)
     # establish a screen to control this job with a name to indicate what it is
     screen -S calMil1Rn6
     mkdir /hive/data/genomes/rn6/bed/lastzCalMil1.2018-01-18
     cd /hive/data/genomes/rn6/bed/lastzCalMil1.2018-01-18
 
     printf "%s\n" \
 '# Rat vs. Elephant shark
 BLASTZ=/cluster/bin/penn/lastz-distrib-1.03.66/bin/lastz
 
 # TARGET: Rat Rn6
 SEQ1_DIR=/hive/data/genomes/rn6/rn6.2bit
 SEQ1_LEN=/hive/data/genomes/rn6/chrom.sizes
 SEQ1_CHUNK=20000000
 SEQ1_LAP=10000
 SEQ1_LIMIT=15
 
 # QUERY: Elephant shark calMil1
 SEQ2_DIR=/hive/data/genomes/calMil1/calMil1.2bit
 SEQ2_LEN=/hive/data/genomes/calMil1/chrom.sizes
 SEQ2_CHUNK=20000000
 SEQ2_LAP=0
 SEQ2_LIMIT=50
 
 BASE=/hive/data/genomes/rn6/bed/lastzCalMil1.2018-01-18
 TMPDIR=/dev/shm' > DEF
 
     # adjust the SEQ2_LIMIT with -stop=partition to get a reasonable
     #	number of jobs, 50,000 to something under 100,000
     # when not present, SEQ2_LIMIT is a default 100
     time (doBlastzChainNet.pl -verbose=2 `pwd`/DEF \
         -workhorse=hgwdev -smallClusterHub=ku -bigClusterHub=ku \
         -chainMinScore=5000 -chainLinearGap=loose) > do.log 2>&1
     #	real    579m52.297s
 
     cat fb.rn6.chainCalMil1Link.txt
     #	45434042 bases of 2729860805 (1.664%) in intersection
 
     time (doBlastzChainNet.pl -verbose=2 `pwd`/DEF -syntenicNet \
         -workhorse=hgwdev -smallClusterHub=ku -bigClusterHub=ku \
         -chainMinScore=5000 -chainLinearGap=loose \
            -continue=syntenicNet) > synNet.log 2>&1 &
     # real    0m57.126s
 
     time (doRecipBest.pl -workhorse=hgwdev -buildDir=`pwd` rn6 calMil1) \
            > rbest.log 2>&1 &
     #  real    193m26.630s
 
     #	and for the swap
     mkdir /hive/data/genomes/calMil1/bed/blastz.rn6.swap
     cd /hive/data/genomes/calMil1/bed/blastz.rn6.swap
     time (doBlastzChainNet.pl -verbose=2 \
 	/hive/data/genomes/rn6/bed/lastzCalMil1.2018-01-18/DEF \
         -workhorse=hgwdev -smallClusterHub=ku -bigClusterHub=ku \
         -swap -chainMinScore=5000 -chainLinearGap=loose) > swap.log 2>&1
     #	real    17m57.104s
 
     cat  fb.calMil1.chainRn6Link.txt
     #	40370300 bases of 936953458 (4.309%) in intersection
 
     time (doRecipBest.pl -workhorse=hgwdev -buildDir=`pwd` calMil1 rn6) \
          > rbest.log 2>&1
     # real    163m30.392s
 
     time (doBlastzChainNet.pl -verbose=2 -syntenicNet \
 	/hive/data/genomes/rn6/bed/lastzCalMil1.2018-01-18/DEF \
         -workhorse=hgwdev -smallClusterHub=ku -bigClusterHub=ku \
         -swap -chainMinScore=5000 -chainLinearGap=loose \
             -continue=syntenicNet) > synNet.log 2>&1 &
     # real    0m51.155s
 
 ##############################################################################
 # lastz tarsier tarSyr2 (DONE - 2017-01-18 - Hiram)
     # establish a screen to control this job with a name to indicate what it is
     screen -S tarSyr2Rn6
     mkdir /hive/data/genomes/rn6/bed/lastzTarSyr2.2017-01-18
     cd /hive/data/genomes/rn6/bed/lastzTarSyr2.2017-01-18
 
     printf "%s\n" \
 '# Rat vs. tarsier
 BLASTZ=/cluster/bin/penn/lastz-distrib-1.03.66/bin/lastz
 BLASTZ_Y=3400
 BLASTZ_L=6000
 BLASTZ_K=2200
 BLASTZ_Q=/scratch/data/blastz/HoxD55.q
 
 # TARGET: Rat Rn6
 SEQ1_DIR=/hive/data/genomes/rn6/rn6.2bit
 SEQ1_LEN=/hive/data/genomes/rn6/chrom.sizes
 SEQ1_CHUNK=20000000
 SEQ1_LAP=10000
 SEQ1_LIMIT=15
 
 # QUERY: tarsier tarSyr2
 SEQ2_DIR=/hive/data/genomes/tarSyr2/tarSyr2.2bit
 SEQ2_LEN=/hive/data/genomes/tarSyr2/chrom.sizes
 SEQ2_CHUNK=20000000
 SEQ2_LAP=0
 SEQ2_LIMIT=1000
 
 BASE=/hive/data/genomes/rn6/bed/lastzTarSyr2.2017-01-18
 TMPDIR=/dev/shm' > DEF
 
     # adjust the SEQ2_LIMIT with -stop=partition to get a reasonable
     #	number of jobs, 50,000 to something under 100,000
     # when not present, SEQ2_LIMIT is a default 100
     time (doBlastzChainNet.pl -verbose=2 `pwd`/DEF \
         -workhorse=hgwdev -smallClusterHub=ku -bigClusterHub=ku \
         -chainMinScore=3000 -chainLinearGap=medium) > do.log 2>&1
     #	real    803m10.196s
 
     cat fb.rn6.chainTarSyr2Link.txt
     #	1051238729 bases of 2729860805 (38.509%) in intersection
 
     time (doBlastzChainNet.pl -verbose=2 `pwd`/DEF -syntenicNet \
         -workhorse=hgwdev -smallClusterHub=ku -bigClusterHub=ku \
         -chainMinScore=3000 -chainLinearGap=medium \
            -continue=syntenicNet) > synNet.log 2>&1 &
     # real    15m30.483s
 
     time (doRecipBest.pl -workhorse=hgwdev -buildDir=`pwd` rn6 tarSyr2) \
            > rbest.log 2>&1 &
     #  real    1075m22.381s
 
     #	and for the swap
     mkdir /hive/data/genomes/tarSyr2/bed/blastz.rn6.swap
     cd /hive/data/genomes/tarSyr2/bed/blastz.rn6.swap
     time (doBlastzChainNet.pl -verbose=2 \
 	/hive/data/genomes/rn6/bed/lastzTarSyr2.2017-01-18/DEF \
         -workhorse=hgwdev -smallClusterHub=ku -bigClusterHub=ku \
         -swap -chainMinScore=3000 -chainLinearGap=medium) > swap.log 2>&1
     #	real    306m57.784s
 
     cat  fb.tarSyr2.chainRn6Link.txt
     #	1110521625 bases of 3405755564 (32.607%) in intersection
 
     time (doRecipBest.pl -workhorse=hgwdev -buildDir=`pwd` tarSyr2 rn6) \
          > rbest.log 2>&1
     # real    1247m45.887s
 
 ##############################################################################
 # lastz platypus ornAna2 (DONE - 2017-01-18 - Hiram)
     # establish a screen to control this job with a name to indicate what it is
     screen -S ornAna2Rn6
     mkdir /hive/data/genomes/rn6/bed/lastzOrnAna2.2017-01-18
     cd /hive/data/genomes/rn6/bed/lastzOrnAna2.2017-01-18
 
     printf "%s\n" \
 '# Rat vs. platypus
 BLASTZ=/cluster/bin/penn/lastz-distrib-1.03.66/bin/lastz
 BLASTZ_M=254
 BLASTZ_Y=3400
 BLASTZ_L=6000
 BLASTZ_K=2200
 BLASTZ_Q=/scratch/data/blastz/HoxD55.q
 
 # TARGET: Rat Rn6
 SEQ1_DIR=/hive/data/genomes/rn6/rn6.2bit
 SEQ1_LEN=/hive/data/genomes/rn6/chrom.sizes
 SEQ1_CHUNK=20000000
 SEQ1_LAP=10000
 SEQ1_LIMIT=5
 
 # QUERY: platypus ornAna2
 SEQ2_DIR=/hive/data/genomes/ornAna2/ornAna2.2bit
 SEQ2_LEN=/hive/data/genomes/ornAna2/chrom.sizes
 SEQ2_CHUNK=20000000
 SEQ2_LAP=0
 SEQ2_LIMIT=1000
 
 BASE=/hive/data/genomes/rn6/bed/lastzOrnAna2.2017-01-18
 TMPDIR=/dev/shm' > DEF
 
     # adjust the SEQ2_LIMIT with -stop=partition to get a reasonable
     #	number of jobs, 50,000 to something under 100,000
     # when not present, SEQ2_LIMIT is a default 100
     time (doBlastzChainNet.pl -verbose=2 `pwd`/DEF \
         -workhorse=hgwdev -smallClusterHub=ku -bigClusterHub=ku \
         -chainMinScore=5000 -chainLinearGap=loose) > do.log 2>&1
     #	real    636m10.747s
 
     cat fb.rn6.chainOrnAna2Link.txt
     #	174410156 bases of 2729860805 (6.389%) in intersection
 
     time (doBlastzChainNet.pl -verbose=2 `pwd`/DEF -syntenicNet \
         -workhorse=hgwdev -smallClusterHub=ku -bigClusterHub=ku \
         -chainMinScore=5000 -chainLinearGap=loose \
           -continue=syntenicNet) > synNet.log 2>&1 &
     # real    2m21.696s
 
     time (doRecipBest.pl -workhorse=hgwdev -buildDir=`pwd` rn6 ornAna2) \
            > rbest.log 2>&1 &
     #  real    801m25.537s
 
     #	and for the swap
     mkdir /hive/data/genomes/ornAna2/bed/blastz.rn6.swap
     cd /hive/data/genomes/ornAna2/bed/blastz.rn6.swap
     time (doBlastzChainNet.pl -verbose=2 \
 	/hive/data/genomes/rn6/bed/lastzOrnAna2.2017-01-18/DEF \
         -workhorse=hgwdev -smallClusterHub=ku -bigClusterHub=ku \
         -swap -chainMinScore=5000 -chainLinearGap=loose) > swap.log 2>&1
     #	real    94m30.003s
 
     cat  fb.ornAna2.chainRn6Link.txt
     #	159391765 bases of 1842252119 (8.652%) in intersection
 
     time (doRecipBest.pl -workhorse=hgwdev -buildDir=`pwd` ornAna2 rn6) \
          > rbest.log 2>&1 &
     # real    794m22.755s
 
     time (doBlastzChainNet.pl -verbose=2 -syntenicNet \
 	/hive/data/genomes/rn6/bed/lastzOrnAna2.2017-01-18/DEF \
         -workhorse=hgwdev -smallClusterHub=ku -bigClusterHub=ku \
         -swap -chainMinScore=5000 -chainLinearGap=loose \
            -continue=syntenicNet) > synNet.log 2>&1 &
     # real    2m59.537s
 
 ##############################################################################
 ## 20-Way Multiz (DONE - 2017-01-20 - Hiram)
     ssh hgwdev
     mkdir /hive/data/genomes/rn6/bed/multiz20way
     cd /hive/data/genomes/rn6/bed/multiz20way
 
     # from the 191-way in the source tree, select out the 5 used here:
     /cluster/bin/phast/tree_doctor \
         --prune-all-but ailMel1,bosTau8,calMil1,canFam3,cavPor3,danRer10,\
 felCat8,galGal5,hg38,melGal5,micOch1,mm10,monDom5,ornAna2,oryCun2,panTro5,\
 rheMac8,rn6,tarSyr2,xenTro7 \
         /cluster/home/hiram/kent/src/hg/utils/phyloTrees/191way.nh \
           > rn6.20way.nh.0
     cat rn6.20way.nh.0
 # (((((((((((hg38:0.006550,panTro5:0.006840):0.029424,
 # rheMac8:0.043601):0.073920,tarSyr2:0.111400):0.036014,(((micOch1:0.140000,
 # (mm10:0.084509,rn6:0.091589):0.047773):0.183142,cavPor3:0.175779):0.025746,
 # oryCun2:0.215690):0.015313):0.020593,(bosTau8:0.189080,(felCat8:0.098612,
 # (canFam3:0.052458,ailMel1:0.080000):0.020000):0.060572):0.032898):0.258392,
 # monDom5:0.340786):0.071664,ornAna2:0.456592):0.109504,(galGal5:0.041254,
 # melGal5:0.085718):0.518188):0.050000,xenTro7:0.347944):0.211354,
 # danRer10:1.201072):0.100000,calMil1:0.750000);
 
     # using TreeGraph2 tree editor on the Mac, rearrange to get rn6
     # at the top:
 
     #	what that looks like:
  ~/kent/src/hg/utils/phyloTrees/asciiTree.pl rn6.20way.nh | sed -e 's/^/# /;'
 # ((((((((((((rn6:0.091589,
 #            mm10:0.084509):0.047773,
 #           micOch1:0.14):0.183142,
 #          cavPor3:0.175779):0.025746,
 #         oryCun2:0.21569):0.015313,
 #        (((hg38:0.00655,
 #          panTro5:0.00684):0.029424,
 #         rheMac8:0.043601):0.07392,
 #        tarSyr2:0.1114):0.036014):0.020593,
 #       (((canFam3:0.052458,
 #         ailMel1:0.08):0.02,
 #        felCat8:0.098612):0.060572,
 #       bosTau8:0.18908):0.032898):0.258392,
 #      monDom5:0.340786):0.071664,
 #     ornAna2:0.456592):0.109504,
 #    (galGal5:0.041254,
 #    melGal5:0.085718):0.518188):0.05,
 #   xenTro7:0.347944):0.211354,
 #  danRer10:1.201072):0.1,
 # calMil1:0.75);
 
     # extract species list from that .nh file
     sed 's/[a-z][a-z]*_//g; s/:[0-9\.][0-9\.]*//g; s/;//; /^ *$/d' \
         rn6.20way.nh | xargs echo | sed 's/ //g; s/,/ /g' \
         | sed 's/[()]//g; s/,/ /g' | tr '[ ]' '[\n]' > species.list.txt
 
     # construct db to name translation list:
     cat species.list.txt | while read DB
 do
 hgsql -N -e "select name,organism from dbDb where name=\"${DB}\";" hgcentraltest
 done | sed -e "s/\t/->/; s/ /_/g;" | sed -e 's/$/;/' | sed -e 's/\./_/g' \
         | sed -e 's/-nosed/_nosed/; s/-eating/_eating/;' > db.to.name.txt
 
     # construct a common name .nh file:
     /cluster/bin/phast/tree_doctor --rename \
     "`cat db.to.name.txt`" rn6.20way.nh | sed -e 's/00*)/)/g; s/00*,/,/g' \
        | $HOME/kent/src/hg/utils/phyloTrees/asciiTree.pl /dev/stdin \
          > rn6.20way.commonNames.nh
     cat rn6.20way.commonNames.nh | sed -e 's/^/# /;'
 # ((((((((((((Rat:0.091589,
 #            Mouse:0.084509):0.047773,
 #           Prairie_vole:0.14):0.183142,
 #          Guinea_pig:0.175779):0.025746,
 #         Rabbit:0.21569):0.015313,
 #        (((Human:0.00655,
 #          Chimp:0.00684):0.029424,
 #         Rhesus:0.043601):0.07392,
 #        Tarsier:0.1114):0.036014):0.020593,
 #       (((Dog:0.052458,
 #         Panda:0.08):0.02,
 #        Cat:0.098612):0.060572,
 #       Cow:0.18908):0.032898):0.258392,
 #      Opossum:0.340786):0.071664,
 #     Platypus:0.456592):0.109504,
 #    (Chicken:0.041254,
 #    Turkey:0.085718):0.518188):0.05,
 #   X__tropicalis:0.347944):0.211354,
 #  Zebrafish:1.201072):0.1,
 # Elephant_shark:0.75);
 
 #	Use this specification in the phyloGif tool:
 #	http://genome.ucsc.edu/cgi-bin/phyloGif
 #	to obtain a png image for src/hg/htdocs/images/phylo/rn6_20way.png
 
     ~/kent/src/hg/utils/phyloTrees/asciiTree.pl rn6.20way.nh > t.nh
     ~/kent/src/hg/utils/phyloTrees/scientificNames.sh t.nh \
        | $HOME/kent/src/hg/utils/phyloTrees/asciiTree.pl /dev/stdin \
           > rn6.20way.scientificNames.nh
     rm -f t.nh
     cat rn6.20way.scientificNames.nh | sed -e 's/^/# /;'
 # ((((((((((((Rattus_norvegicus:0.091589,
 #            Mus_musculus:0.084509):0.047773,
 #           Microtus_ochrogaster:0.14):0.183142,
 #          Cavia_porcellus:0.175779):0.025746,
 #         Oryctolagus_cuniculus:0.21569):0.015313,
 #        (((Homo_sapiens:0.00655,
 #          Pan_troglodytes:0.00684):0.029424,
 #         Macaca_mulatta:0.043601):0.07392,
 #        Tarsius_syrichta:0.1114):0.036014):0.020593,
 #       (((Canis_lupus_familiaris:0.052458,
 #         Ailuropoda_melanoleuca:0.08):0.02,
 #        Felis_catus:0.098612):0.060572,
 #       Bos_taurus:0.18908):0.032898):0.258392,
 #      Monodelphis_domestica:0.340786):0.071664,
 #     Ornithorhynchus_anatinus:0.456592):0.109504,
 #    (Gallus_gallus:0.041254,
 #    Meleagris_gallopavo:0.085718):0.518188):0.05,
 #   Xenopus_tropicalis:0.347944):0.211354,
 #  Danio_rerio:1.201072):0.1,
 # Callorhinchus_milii:0.75);
 
     /cluster/bin/phast/all_dists rn6.20way.nh | grep rn6 \
         | sed -e "s/rn6.//" | sort -k2n > 20way.distances.txt
     #	Use this output to create the table below
     cat 20way.distances.txt | sed -e 's/^/# /;'
 # mm10  0.176098
 # micOch1       0.279362
 # cavPor3       0.498283
 # hg38  0.509471
 # panTro5       0.509761
 # tarSyr2       0.510977
 # rheMac8       0.517098
 # canFam3       0.550084
 # oryCun2       0.563940
 # felCat8       0.576238
 # ailMel1       0.577626
 # bosTau8       0.606134
 # monDom5       0.983334
 # ornAna2       1.170804
 # xenTro7       1.221660
 # galGal5       1.383158
 # melGal5       1.427622
 # calMil1       1.935070
 # danRer10      2.286142
 
     printf '#!/usr/bin/env perl
 
 use strict;
 use warnings;
 
 open (FH, "<20way.distances.txt") or
         die "can not read 20way.distances.txt";
 
 my $count = 0;
 while (my $line = <FH>) {
     chomp $line;
     my ($D, $dist) = split('"'"'\\s+'"'"', $line);
     my $chain = "chain" . ucfirst($D);
     my $B="/hive/data/genomes/rn6/bed/lastz.$D/fb.rn6." .
         $chain . "Link.txt";
     my $chainLinkMeasure =
         `awk '"'"'{print \\$5}'"'"' ${B} 2> /dev/null | sed -e "s/(//; s/)//"`;
     chomp $chainLinkMeasure;
     $chainLinkMeasure = 0.0 if (length($chainLinkMeasure) < 1);
     $chainLinkMeasure =~ s/\\%%//;
     my $swapFile="/hive/data/genomes/${D}/bed/lastz.rn6/fb.${D}.chainRn6Link.txt";
     my $swapMeasure = "N/A";
     if ( -s $swapFile ) {
 	$swapMeasure =
 	    `awk '"'"'{print \\$5}'"'"' ${swapFile} 2> /dev/null | sed -e "s/(//; s/)//"`;
 	chomp $swapMeasure;
 	$swapMeasure = 0.0 if (length($swapMeasure) < 1);
 	$swapMeasure =~ s/\\%%//;
     }
     my $orgName=
     `hgsql -N -e "select organism from dbDb where name='"'"'$D'"'"';" hgcentraltest`;
     chomp $orgName;
     if (length($orgName) < 1) {
         $orgName="N/A";
     }
     ++$count;
     printf "# %%02d  %%.4f (%%%% %%06.3f) (%%%% %%06.3f) - %%s %%s\\n", $count, $dist,
         $chainLinkMeasure, $swapMeasure, $orgName, $D;
 }
 close (FH);
 ' > sizeStats.pl
     chmod +x ./sizeStats.pl
     ./sizeStats.pl
 
 #	If you can fill in all the numbers in this table, you are ready for
 #	the multiple alignment procedure
 
 #       featureBits chainLink measures
 #               chainLink
 #  N distance  on rn6  on other     other species
 # 01  0.1761 (% 71.015) (% 70.886) - Mouse mm10
 # 02  0.2794 (% 57.205) (% 70.134) - Prairie vole micOch1
 # 03  0.4983 (% 27.834) (% 28.362) - Guinea pig cavPor3
 # 04  0.5095 (% 34.899) (% 31.077) - Human hg38
 # 05  0.5098 (% 40.816) (% 35.356) - Chimp panTro5
 # 06  0.5110 (% 38.509) (% 32.607) - Tarsier tarSyr2
 # 07  0.5171 (% 40.165) (% 33.468) - Rhesus rheMac8
 # 08  0.5501 (% 28.540) (% 30.864) - Dog canFam3
 # 09  0.5639 (% 24.654) (% 24.981) - Rabbit oryCun2
 # 10  0.5762 (% 29.520) (% 29.673) - Cat felCat8
 # 11  0.5776 (% 30.404) (% 34.777) - Panda ailMel1
 # 12  0.6061 (% 25.902) (% 25.355) - Cow bosTau8
 # 13  0.9833 (% 09.652) (% 07.191) - Opossum monDom5
 # 14  1.1708 (% 06.389) (% 08.652) - Platypus ornAna2
 # 15  1.2217 (% 04.104) (% 07.913) - X. tropicalis xenTro7
 # 16  1.3832 (% 03.983) (% 07.853) - Chicken galGal5
 # 17  1.4276 (% 02.440) (% 04.910) - Turkey melGal5
 # 18  1.9351 (% 01.664) (% 04.309) - Elephant shark calMil1
 # 19  2.2861 (% 03.232) (% 05.815) - Zebrafish danRer10
 
 # None of this concern for distances matters in building the first step, the
 # maf files.  The distances will be better calibrated later.
 
     # create species list and stripped down tree for autoMZ
     sed 's/[a-z][a-z]*_//g; s/:[0-9\.][0-9\.]*//g; s/;//; /^ *$/d' \
 	rn6.20way.nh | xargs echo | sed 's/ //g; s/,/ /g' > tree.nh
 
     sed 's/[()]//g; s/,/ /g' tree.nh > species.list
 # rn6 mm10 micOch1 cavPor3 oryCun2 hg38 panTro5 rheMac8 tarSyr2 canFam3
 #   ailMel1 felCat8 bosTau8 monDom5 ornAna2 galGal5 melGal5 xenTro7
 #   danRer10 calMil1
 
     # survey N50 for each
     for db in `cat species.list`
 do
 printf "# %7s:\t" "${db}"
 n50.pl /hive/data/genomes/$db/chrom.sizes 2>&1 | tail -1
 done | sort -k4n
 
 # monDom5:      1817563546      3       chr3    527952102
 # galGal5:      630658223       5       chrZ    82310166
 # xenTro7:      774070486       5       KB021657        124127367
 # melGal5:      591651097       6       chr5    59006440
 #     rn6:      1538522439      8       chr7    145729302
 # felCat8:      1440507093      8       chrB4   142431058
 #    hg38:      1692529807      9       chr8    145138636
 #    mm10:      1442871972      9       chr10   130694993
 # oryCun2:      1389219202      9       chr8    111795807
 # rheMac8:      1622316098      9       chr8    144306982
 # panTro5:      1725767609      10      chr10   135926727
 # bosTau8:      1366174162      11      chr9    105708250
 # danRer10:     738680380       12      chr17   53345113
 # micOch1:      1154595365      14      chrX    61813736
 # canFam3:      1241874881      15      chr13   63241923
 # cavPor3:      1384780737      28      scaffold_27     27942054
 # calMil1:      489260987       60      KI635914        4521921
 # ornAna2:      998675786       182     chrUn_DS181001v1        991605
 # ailMel1:      1150550801      521     GL192858.1      1281781
 # tarSyr2:      1727089816      2384    KE946086v1      401181
 
 
     #	bash shell syntax here ...
     cd /hive/data/genomes/rn6/bed/multiz20way
     export H=/hive/data/genomes/rn6/bed
     mkdir mafLinks
     # good assemblies can use syntenic net:
     # monDom5 rn6 felCat8 hg38 mm10 oryCun2 rheMac8
     # panTro5 bosTau8 micOch1 canFam3
 
     for G in monDom5 felCat8 hg38 mm10 oryCun2 cavPor3 rheMac8 \
              panTro5 bosTau8 micOch1 canFam3
     do
       mkdir mafLinks/$G
       echo ln -s ${H}/lastz.$G/axtChain/rn6.${G}.synNet.maf.gz ./mafLinks/$G
       ln -s ${H}/lastz.$G/axtChain/rn6.${G}.synNet.maf.gz ./mafLinks/$G
     done
 
     # distant species, use unfiltered mafNet
     for G in ornAna2 calMil1 danRer10 galGal5 melGal5 xenTro7
     do
       mkdir mafLinks/$G
       echo ln -s ${H}/lastz.$G/mafNet/rn6.${G}.net.maf.gz ./mafLinks/$G
       ln -s ${H}/lastz.$G/mafNet/rn6.${G}.net.maf.gz ./mafLinks/$G
     done
 
     # lower quality assemblies using recip best net:
     #  ailMel1 tarSyr2
     for G in ailMel1 tarSyr2
     do
       mkdir mafLinks/$G
       echo ln -s ${H}/lastz.$G/mafRBestNet/rn6.${G}.rbest.maf.gz ./mafLinks/$G
       ln -s ${H}/lastz.$G/mafRBestNet/rn6.${G}.rbest.maf.gz ./mafLinks/$G
     done
 
     # verify the symLinks are good:
     ls -ogrtL mafLinks/*/* | sed -e 's/-rw-rw-r-- 1//;' | while read L
 do
    let "n++"
    printf "# %02d\t%s\n" $n "$L"
 done
 # 01    687455612 Feb 17  2015 mafLinks/hg38/rn6.hg38.synNet.maf.gz
 # 02    61987372 Sep 18  2015 mafLinks/danRer10/rn6.danRer10.net.maf.gz
 # 03    1192966029 Apr  8  2016 mafLinks/mm10/rn6.mm10.synNet.maf.gz
 # 04    80158096 May  3  2016 mafLinks/galGal5/rn6.galGal5.net.maf.gz
 # 05    74068999 Jan 12 17:25 mafLinks/xenTro7/rn6.xenTro7.net.maf.gz
 # 06    604285485 Jan 17 03:57 mafLinks/ailMel1/rn6.ailMel1.rbest.maf.gz
 # 07    49852699 Jan 17 18:41 mafLinks/melGal5/rn6.melGal5.net.maf.gz
 # 08    32487411 Jan 18 19:14 mafLinks/calMil1/rn6.calMil1.net.maf.gz
 # 09    126292792 Jan 18 20:27 mafLinks/ornAna2/rn6.ornAna2.net.maf.gz
 # 10    132321627 Jan 19 10:13 mafLinks/monDom5/rn6.monDom5.synNet.maf.gz
 # 11    562609533 Jan 19 10:24 mafLinks/canFam3/rn6.canFam3.synNet.maf.gz
 # 12    535101244 Jan 19 10:24 mafLinks/cavPor3/rn6.cavPor3.synNet.maf.gz
 # 13    749344283 Jan 19 10:31 mafLinks/panTro5/rn6.panTro5.synNet.maf.gz
 # 14    739119670 Jan 19 10:32 mafLinks/rheMac8/rn6.rheMac8.synNet.maf.gz
 # 15    497464725 Jan 19 10:33 mafLinks/bosTau8/rn6.bosTau8.synNet.maf.gz
 # 16    578966556 Jan 19 10:39 mafLinks/felCat8/rn6.felCat8.synNet.maf.gz
 # 17    1039380892 Jan 19 10:52 mafLinks/micOch1/rn6.micOch1.synNet.maf.gz
 # 18    463610104 Jan 19 11:06 mafLinks/oryCun2/rn6.oryCun2.synNet.maf.gz
 # 19    724617931 Jan 20 02:02 mafLinks/tarSyr2/rn6.tarSyr2.rbest.maf.gz
 
     # and note that it is the correct count 20-way minus 1 == 19
 
     # split the maf files into a set of hashed named files
     # this hash named split keeps the same chr/contig names in the same
     # named hash file.
     mkdir /hive/data/genomes/rn6/bed/multiz20way/mafSplit
     cd /hive/data/genomes/rn6/bed/multiz20way/mafSplit
     time for D in `sed -e "s/rn6 //" ../species.list`
 do
     echo "${D}"
     mkdir $D
     cd $D
     echo "mafSplit -byTarget -useHashedName=8 /dev/null . ../../mafLinks/${D}/*.maf.gz"
     mafSplit -byTarget -useHashedName=8 /dev/null . \
 	../../mafLinks/${D}/*.maf.gz
     cd ..
 done
     # real    6m51.805s
 
     # construct a list of all possible maf file names.
     # they do not all exist in each of the species directories
     find . -type f | wc -l
     # 2173
     find . -type f | grep ".maf$" | xargs -L 1 basename | sort -u > maf.list
     wc -l maf.list
     # 227 maf.list
 
     mkdir /hive/data/genomes/rn6/bed/multiz20way/splitRun
     cd /hive/data/genomes/rn6/bed/multiz20way/splitRun
     mkdir maf run
     cd run
     mkdir penn
     cp -p /cluster/bin/penn/multiz.2009-01-21_patched/multiz penn
     cp -p /cluster/bin/penn/multiz.2009-01-21_patched/maf_project penn
     cp -p /cluster/bin/penn/multiz.2009-01-21_patched/autoMZ penn
 
     #	verify the db and pairs settings are correct
     printf '#!/bin/csh -ef
 set db = rn6
 set c = $1
 set result = $2
 set run = `/bin/pwd`
 set tmp = /dev/shm/$db/multiz.$c
 set pairs = /hive/data/genomes/rn6/bed/multiz20way/mafSplit
 /bin/rm -fr $tmp
 /bin/mkdir -p $tmp
 /bin/cp -p ../../tree.nh ../../species.list $tmp
 pushd $tmp > /dev/null
 foreach s (`/bin/sed -e "s/$db //" species.list`)
     set in = $pairs/$s/$c
     set out = $db.$s.sing.maf
     if (-e $in.gz) then
         /bin/zcat $in.gz > $out
         if (! -s $out) then
             echo "##maf version=1 scoring=autoMZ" > $out
         endif
     else if (-e $in) then
         /bin/ln -s $in $out
     else
         echo "##maf version=1 scoring=autoMZ" > $out
     endif
 end
 set path = ($run/penn $path); rehash
 $run/penn/autoMZ + T=$tmp E=$db "`cat tree.nh`" $db.*.sing.maf $c \
         > /dev/null
 popd > /dev/null
 /bin/rm -f $result
 /bin/cp -p $tmp/$c $result
 /bin/rm -fr $tmp
 ' > autoMultiz.csh
 
     chmod +x autoMultiz.csh
 
     printf '#LOOP
 ./autoMultiz.csh $(file1) {check out line+ /hive/data/genomes/rn6/bed/multiz20way/splitRun/maf/$(root1).maf}
 #ENDLOOP
 ' > template
 
     ln -s ../../mafSplit/maf.list maf.list
     ssh ku
     cd /hive/data/genomes/rn6/bed/multiz20way/splitRun/run
     gensub2 maf.list single template jobList
     para create jobList
     para try ... check ... push ... etc...
 # Completed: 227 of 227 jobs
 # CPU time in finished jobs:     458979s    7649.65m   127.49h    5.31d  0.015 y
 # IO & Wait Time:                  2311s      38.52m     0.64h    0.03d  0.000 y
 # Average job time:                2032s      33.87m     0.56h    0.02d
 # Longest finished job:           41181s     686.35m    11.44h    0.48d
 # Submission to last job:         41324s     688.73m    11.48h    0.48d
 
     # combine into one file  (the 1>&2 redirect sends the echo to stderr)
     cd /hive/data/genomes/rn6/bed/multiz20way
     head -1 splitRun/maf/013.maf > multiz20way.maf
     time for F in splitRun/maf/*.maf
 do
     echo "${F}" 1>&2
     egrep -v "^#" ${F}
 done >> multiz20way.maf
     # real    1m49.082s
 
 
     tail -1 splitRun/maf/013.maf >> multiz20way.maf
 # -rw-rw-r--  1 29010594028 Jan 22 09:08 multiz20way.maf
 
 # Load into database
     ssh hgwdev
     cd /hive/data/genomes/rn6/bed/multiz20way
     mkdir /gbdb/rn6/multiz20way
     ln -s `pwd`/multiz20way.maf /gbdb/rn6/multiz20way
     cd /dev/shm
     time hgLoadMaf rn6 multiz20way
 # Loaded 24754572 mafs in 1 files from /gbdb/rn6/multiz20way
 # real    8m51.299s
 
     time hgLoadMafSummary -verbose=2 -minSize=30000 \
 	-mergeGap=1500 -maxSize=200000 rn6 multiz20waySummary \
 	/gbdb/rn6/multiz20way/multiz20way.maf
 # Created 3877691 summary blocks from 169823109 components and 24754572 mafs from /gbdb/rn6/multiz20way/multiz20way.maf
 # real    12m4.373s
 
 # -rw-rw-r-- 1 1265516402 Jan 22 09:15 multiz20way.tab
 # -rw-rw-r-- 1  182504514 Jan 22 09:37 multiz20waySummary.tab
 
     wc -l multiz20way*.tab
 #   24754572 multiz20way.tab
 #    3877691 multiz20waySummary.tab
 
     rm multiz20way*.tab
 
 ##############################################################################
 # GAP ANNOTATE MULTIZ20WAY MAF AND LOAD TABLES (DONE - 2017-01-22 - Hiram)
     # mafAddIRows has to be run on single chromosome maf files, it does not
     #	function correctly when more than one reference sequence
     #	are in a single file.  Need to split of the maf file into individual
     #   maf files
     mkdir -p /hive/data/genomes/rn6/bed/multiz20way/anno/mafSplit
     cd /hive/data/genomes/rn6/bed/multiz20way/anno/mafSplit
 
     time mafSplit -outDirDepth=2 -byTarget -useFullSequenceName \
         /dev/null . ../../multiz20way.maf
     #   real    10m33.882s
 
     find . -type f | wc -l
     #   624
 
     # check for N.bed files everywhere:
     cd /hive/data/genomes/rn6/bed/multiz20way/anno
     for DB in `cat ../species.list`
 do
     if [ ! -s /hive/data/genomes/${DB}/${DB}.N.bed ]; then
         echo "MISS: ${DB}"
 #         cd /hive/data/genomes/${DB}
 #         twoBitInfo -nBed ${DB}.2bit ${DB}.N.bed
     else
         echo "  OK: ${DB}"
     fi
 done
 
     cd /hive/data/genomes/rn6/bed/multiz20way/anno
     for DB in `cat ../species.list`
 do
     echo "${DB} "
     ln -s  /hive/data/genomes/${DB}/${DB}.N.bed ${DB}.bed
     echo ${DB}.bed  >> nBeds
     ln -s  /hive/data/genomes/${DB}/chrom.sizes ${DB}.len
     echo ${DB}.len  >> sizes
 done
     # make sure they all are successful symLinks:
     ls -ogrtL
 
     screen -S gapAnno      # use a screen to control this longish job
     ssh ku
     cd /hive/data/genomes/rn6/bed/multiz20way/anno
     mkdir result
     find ./mafSplit -type d | sed -e 's#./mafSplit/##' | while read D
 do
     echo mkdir -p result/${D}
     mkdir -p result/${D}
 done
     printf '#LOOP
 mafAddIRows -nBeds=nBeds mafSplit/$(path1) /hive/data/genomes/rn6/rn6.2bit {check out exists+ result/$(path1)}
 #ENDLOOP
 ' > template
     # << happy emacs
 
     find ./mafSplit -type f | sed -e 's#^./mafSplit/##' > maf.list
     gensub2 maf.list single template jobList
     # limit jobs on a node with the ram=32g requirement because they go fast
     para -ram=32g create jobList
     para try ... check ... push ...
 # Completed: 624 of 624 jobs
 # CPU time in finished jobs:       1802s      30.03m     0.50h    0.02d  0.000 y
 # IO & Wait Time:                  1624s      27.07m     0.45h    0.02d  0.000 y
 # Average job time:                   5s       0.09m     0.00h    0.00d
 # Longest finished job:             148s       2.47m     0.04h    0.00d
 # Submission to last job:           177s       2.95m     0.05h    0.00d
 
     # verify all result files have some content, look for 0 size files:
     find ./result -type f -size 0
     # should see none
     # or in this manner:
     find ./result -type f | xargs ls -og | sort -k3nr | tail
 
     # combine into one file  (the 1>&2 redirect sends the echo to stderr)
     head -q -n 1 result/4/0/chrUn_KL568283v1.maf > rn6.20way.maf
     time find ./result -type f | while read F
 do
     echo "${F}" 1>&2
     grep -h -v "^#" ${F}
 done >> rn6.20way.maf
     # real    2m13.643s
 
     #	these maf files do not have the end marker, this does nothing:
     #	tail -q -n 1 result/4/0/chrUn_KL568283v1.maf >> rn6.20way.maf
     # How about an official end marker:
     echo "##eof maf" >> rn6.20way.maf
     # it almost doubles in size
     ls -og
 # -rw-rw-r--  1 48341214170 Jan 22 10:39 rn6.20way.maf
 # -rw-rw-r--  1 29010594028 Jan 22 09:08 multiz20way.maf
 
     du -hsc rn6.20way.maf
     # 46G     rn6.20way.maf
 
     # construct symlinks to get the individual maf files into gbdb:
     rm /gbdb/rn6/multiz20way/multiz20way.maf   # remove previous results
     ln -s `pwd`/rn6.20way.maf /gbdb/rn6/multiz20way/multiz20way.maf
 
     # Load into database
     cd /dev/shm
     time hgLoadMaf -pathPrefix=/gbdb/rn6/multiz20way rn6 multiz20way
     # Loaded 25393063 mafs in 1 files from /gbdb/rn6/multiz20way
     # real    12m43.103s
 
     time hgLoadMafSummary -verbose=2 -minSize=30000 \
 	-mergeGap=1500 -maxSize=200000 rn6 multiz20waySummary \
         /gbdb/rn6/multiz20way/multiz20way.maf
     # Created 3877691 summary blocks from 169823109 components
     #    and 25393063 mafs from /gbdb/rn6/multiz20way/multiz20way.maf
     # real    15m39.900s
 
 
 # -rw-rw-r-- 1 1301989038 Jan 22 10:53 multiz20way.tab
 # -rw-rw-r-- 1  190259896 Jan 22 11:12 multiz20waySummary.tab
 
     rm multiz20way*.tab
 
 ######################################################################
 # MULTIZ20WAY MAF FRAMES (DONE - 2017-01-22 - Hiram)
     ssh hgwdev
     mkdir /hive/data/genomes/rn6/bed/multiz20way/frames
     cd /hive/data/genomes/rn6/bed/multiz20way/frames
 #   survey all the genomes to find out what kinds of gene tracks they have
     printf '#!/bin/csh -fe
 foreach db (`cat ../species.list`)
     printf "# ${db}: "
     set tables = `hgsql $db -N -e "show tables" | egrep "Gene|ncbiRefSeq"`
     foreach table ($tables)
         if ($table == "ensGene" || $table == "refGene" || \
            $table == "ncbiRefSeq" || $table == "mgcGenes" || \
            $table == "knownGene" || $table == "xenoRefGene" ) then
            set count = `hgsql $db -N -e "select count(*) from $table"`
             echo -n "${table}: ${count}, "
         endif
     end
     set orgName = `hgsql hgcentraltest -N -e \
             "select scientificName from dbDb where name='"'"'$db'"'"'"`
     set orgId = `hgsql $db -N -e \
             "select id from organism where name='"'"'$orgName'"'"'"`
     if ($orgId == "") then
         echo "Mrnas: 0"
     else
         set count = `hgsql $db -N -e "select count(*) from gbCdnaInfo where organism=$orgId"`
         echo "Mrnas: ${count}"
     endif
 end
 ' > showGenes.csh
     # << happy emacs
     chmod +x ./showGenes.csh
     time ./showGenes.csh
 # rn6: ensGene: 40459, mgcGenes: 7017, refGene: 19039, xenoRefGene: 197625, Mrnas: 1250600
 # mm10: ensGene: 103734, knownGene: 63759, mgcGenes: 26777, ncbiRefSeq: 106507, refGene: 36803, xenoRefGene: 177505, Mrnas: 5364985
 # micOch1: Mrnas: 19
 # cavPor3: ensGene: 26129, refGene: 488, xenoRefGene: 305485, Mrnas: 21235
 # oryCun2: ensGene: 24964, refGene: 1445, xenoRefGene: 313991, Mrnas: 40491
 # hg38: ensGene: 208239, knownGene: 197782, mgcGenes: 35305, ncbiRefSeq: 159322, refGene: 67251, xenoRefGene: 184734, Mrnas: 11474400
 # panTro5: refGene: 2897, xenoRefGene: 223650, Mrnas: 11351
 # rheMac8: ensGene: 56743, refGene: 6465, xenoRefGene: 214674, Mrnas: 447197
 # tarSyr2: xenoRefGene: 337666, Mrnas: 0
 # canFam3: ensGene: 39074, refGene: 2291, xenoRefGene: 262153, Mrnas: 388033
 # ailMel1: ensGene: 25055, xenoRefGene: 338914, Mrnas: 276
 # felCat8: refGene: 433, xenoRefGene: 225635, Mrnas: 3785
 # bosTau8: refGene: 14630, xenoRefGene: 297152, Mrnas: 1617801
 # monDom5: ensGene: 24882, refGene: 995, xenoRefGene: 251960, Mrnas: 2814
 # ornAna2: refGene: 695, xenoRefGene: 605912, Mrnas: 10740
 # galGal5: ensGene: 48760, refGene: 7541, xenoRefGene: 235493, Mrnas: 638577
 # melGal5: refGene: 109, xenoRefGene: 231886, Mrnas: 17894
 # xenTro7: mgcGenes: 9029, refGene: 8910, xenoRefGene: 265136, Mrnas: 1298763
 # danRer10: ensGene: 58274, mgcGenes: 16947, refGene: 16551, xenoRefGene: 127693, Mrnas: 1535288
 # calMil1: xenoRefGene: 361187, Mrnas: 146271
 
     # from that summary, use these gene sets:
     # ensGene: rn6 cavPor3 oryCun2 rheMac8 canFam3 ailMel1 monDom5 galGal5
     #          danRer10
     # refGene: panTro5 bosTau8 xenTro7
     # knownGene - hg38 mm10
     # none: micOch1 tarSyr2 felCat8 ornAna2 melGal5 calMil1
 
     mkdir genes
     #   1. knownGene: hg38 mm10
     for DB in hg38 mm10
 do
     hgsql -N -e "select name,chrom,strand,txStart,txEnd,cdsStart,cdsEnd,exonCount,exonStarts,exonEnds from knownGene" ${DB} \
       | genePredSingleCover stdin stdout | gzip -2c \
         > genes/${DB}.gp.gz
     printf "# ${DB}: "
     genePredCheck -db=${DB} genes/${DB}.gp.gz
 done
 # hg38: checked: 21375 failed: 0
 # mm10: checked: 21100 failed: 0
 
     #   2. ensGene: rn6 cavPor3 oryCun2 rheMac8 canFam3 ailMel1
     #               monDom5 galGal5 danRer10
 
     for DB in rn6 cavPor3 oryCun2 rheMac8 canFam3 ailMel1 monDom5 galGal5 danRer10
 do
 hgsql -N -e "select name,chrom,strand,txStart,txEnd,cdsStart,cdsEnd,exonCount,exonStarts,exonEnds
 from ensGene" ${DB} \
       | genePredSingleCover stdin stdout | gzip -2c \
         > /dev/shm/${DB}.tmp.gz
     mv /dev/shm/${DB}.tmp.gz genes/$DB.gp.gz
     printf "# ${DB}: "
     genePredCheck -db=${DB} genes/${DB}.gp.gz
 done
 # rn6: checked: 22124 failed: 0
 # cavPor3: checked: 18631 failed: 0
 # oryCun2: checked: 19165 failed: 0
 # rheMac8: checked: 20859 failed: 0
 # canFam3: checked: 19507 failed: 0
 # ailMel1: checked: 19204 failed: 0
 # monDom5: checked: 21033 failed: 0
 # galGal5: checked: 16229 failed: 0
 # danRer10: checked: 25250 failed: 0
 
     #   3. refGene for panTro5 bosTau8 refGene
     for DB in panTro5 bosTau8 xenTro7
 do
 hgsql -N -e "select * from refGene" ${DB} | cut -f2- \
       | genePredSingleCover stdin stdout | gzip -2c \
         > /dev/shm/${DB}.tmp.gz
     mv /dev/shm/${DB}.tmp.gz genes/$DB.gp.gz
     echo -n "# ${DB}: "
     genePredCheck -db=${DB} genes/${DB}.gp.gz
 done
 # panTro5: checked: 2068 failed: 0
 # bosTau8: checked: 13344 failed: 0
 # xenTro7: checked: 8517 failed: 0
 
     # verify counts for genes are reasonable:
     for T in genes/*.gz
 do
     echo -n "# $T: "
     zcat $T | cut -f1 | sort | uniq -c | wc -l
 done
 # genes/ailMel1.gp.gz: 19204
 # genes/bosTau8.gp.gz: 13071
 # genes/canFam3.gp.gz: 19507
 # genes/cavPor3.gp.gz: 18631
 # genes/danRer10.gp.gz: 25250
 # genes/galGal5.gp.gz: 16229
 # genes/hg38.gp.gz: 21375
 # genes/mm10.gp.gz: 21100
 # genes/monDom5.gp.gz: 21033
 # genes/oryCun2.gp.gz: 19165
 # genes/panTro5.gp.gz: 1973
 # genes/rheMac8.gp.gz: 20859
 # genes/rn6.gp.gz: 22124
 # genes/xenTro7.gp.gz: 8338
 
     time (cat ../anno/rn6.20way.maf | genePredToMafFrames rn6 stdin stdout \
           `ls genes | sed -e 's/.gp.gz//;' | xargs echo \
             | sed -e "s#\([a-zA-Z0-9]*\)#\1 genes/\1.gp.gz#g;"` \
 		| gzip > multiz20wayFrames.bed.gz)
     # real    11m43.455s
 
     # verify there are frames on everything, should be 14 species:
     zcat multiz20wayFrames.bed.gz | awk '{print $4}' | sort | uniq -c \
        | sed -e 's/^/# /;'
 #  234094 ailMel1
 #  148818 bosTau8
 #  247686 canFam3
 #  216773 cavPor3
 #  452415 danRer10
 #  403321 galGal5
 #  241462 hg38
 #  232386 mm10
 #  233904 monDom5
 #  214841 oryCun2
 #   20024 panTro5
 #  251019 rheMac8
 #  196117 rn6
 #  198989 xenTro7
 
     #   load the resulting file
     ssh hgwdev
     cd /hive/data/genomes/rn6/bed/multiz20way/frames
     time hgLoadMafFrames rn6 multiz20wayFrames multiz20wayFrames.bed.gz
     #   real    0m37.422s
 
     time featureBits -countGaps rn6 multiz20wayFrames
     # 51905008 bases of 2870182909 (1.808%) in intersection
     # real    0m21.147s
 
     #   enable the trackDb entries:
 # frames multiz20wayFrames
 # irows on
     #   appears to work OK
 
 #########################################################################
 # Phylogenetic tree for 20-way (DONE - 2017-01-23 - Hiram)
     mkdir /hive/data/genomes/rn6/bed/multiz20way/4d
     cd /hive/data/genomes/rn6/bed/multiz20way/4d
 
     # using the ensGene
     hgsql -N -e "select * from ensGene;" rn6 \
       | cut -f2- | genePredSingleCover stdin stdout > rn6.ensGene.gp
     genePredCheck -db=rn6 rn6.ensGene.gp
     # checked: 22124 failed: 0
 
     # the annotated maf is:
     og ../anno/rn6.20way.maf
 # -rw-rw-r-- 1 48341214170 Jan 22 10:39 ../anno/rn6.20way.maf
 
     mkdir annoSplit
     cd annoSplit
     time mafSplit -verbose=2 -outDirDepth=2 -byTarget -useFullSequenceName \
 	/dev/null . ../../anno/rn6.20way.maf
     # real    20m9.735s
 
 
     find . -type f | wc -l
     #   624
     ssh ku
     mkdir /hive/data/genomes/rn6/bed/multiz20way/4d/run
     cd /hive/data/genomes/rn6/bed/multiz20way/4d/run
     mkdir ../mfa
 
     # newer versions of msa_view have a slightly different operation
     # the sed of the gp file inserts the reference species in the chr name
     printf '#!/bin/csh -fe
 set PHASTBIN = /cluster/bin/phast.build/cornellCVS/phast.2010-12-30/bin
 set GP = rn6.ensGene.gp
 set r = "/hive/data/genomes/rn6/bed/multiz20way"
 set c = $1:r
 set infile = $r/4d/annoSplit/$2
 set outDir = $r/4d/mfa/$3:h
 set outfile = $r/4d/mfa/$3
 /bin/mkdir -p $outDir
 cd /dev/shm
 /bin/awk -v C=$c '"'"'$2 == C {print}'"'"' $r/4d/$GP | sed -e "s/\\t$c\\t/\\trn6.$c\\t/" > $c.gp
 set NL=`wc -l $c.gp| gawk '"'"'{print $1}'"'"'`
 echo $NL
 if ("$NL" != "0") then
     $PHASTBIN/msa_view --4d --features $c.gp -i MAF $infile -o SS > $c.ss
     $PHASTBIN/msa_view -i SS --tuple-size 1 $c.ss > $outfile
 else
     echo "" > $outfile
 endif
 /bin/rm -f /dev/shm/$c.gp /dev/shm/$c.ss
 ' > 4d.csh
 
     chmod +x 4d.csh
 
     find ../annoSplit -type f | sed -e "s#../annoSplit/##" > maf.list
     wc -l maf.list
 # 624 maf.list
 
     printf '#LOOP
 4d.csh $(file1) $(path1) {check out line+ ../mfa/$(dir1)/$(dir2)$(root1).mfa}
 #ENDLOOP
 ' > template
 
     gensub2 maf.list single template jobList
     # better to use 32g here to avoid out of memory, and they run quickly
     # anyway, no need to run so many at a time
     para -ram=32g create jobList
     para try ... check
     para time
     # one job ran out of memory, finish on hgwdev:
     time ./4d.csh chr1.maf 0/3/chr1.maf ../mfa/0/3/chr1.mfa
     # real    3m5.094s
     # 17 others have no results
 # Completed: 606 of 624 jobs
 # Crashed: 18 jobs
 # CPU time in finished jobs:       1241s      20.68m     0.34h    0.01d  0.000 y
 # IO & Wait Time:                  1805s      30.09m     0.50h    0.02d  0.000 y
 # Average job time:                   5s       0.08m     0.00h    0.00d
 # Longest finished job:             130s       2.17m     0.04h    0.00d
 # Submission to last job:           170s       2.83m     0.05h    0.00d
 
     # Not all results have contents, or finish successfully, that is OK
     # it is because not all contigs have genes, only gene sequences are measured
 
     # combine mfa files
     ssh hgwdev
     cd /hive/data/genomes/rn6/bed/multiz20way/4d
     # remove the broken empty files, size 0 and size 1:
     find ./mfa -type f -size 0 | xargs rm -f
     # sometimes this doesn't work, don't know why
     find ./mfa -type f -size 1 | xargs rm -f
     # when it doesn't, use this empty list procedure
     find ./mfa -type f | xargs ls -og | awk '$3 < 2' | awk '{print $NF}' \
         > empty.list
     cat empty.list | xargs rm -f
     # see what is left:
     ls -ogrt mfa/*/*/*.mfa | sort -k3nr | wc -l
     #  116
 
     # want comma-less species.list
     time /cluster/bin/phast.build/cornellCVS/phast.2010-12-30/bin/msa_view \
 	--aggregate "`cat ../species.list`" mfa/*/*/*.mfa | sed s/"> "/">"/ \
 	    > 4d.all.mfa
     # real    0m2.356s
 
     # check they are all in there, should be 20:
     grep "^>" 4d.all.mfa | wc -l
     # 20
     grep "^>" 4d.all.mfa | sed -e 's/^/# /;'
 # >rn6
 # >mm10
 # >micOch1
 # >cavPor3
 # >oryCun2
 # >hg38
 # >panTro5
 # >rheMac8
 # >tarSyr2
 # >canFam3
 # >ailMel1
 # >felCat8
 # >bosTau8
 # >monDom5
 # >ornAna2
 # >galGal5
 # >melGal5
 # >xenTro7
 # >danRer10
 # >calMil1
 
     sed 's/[a-z][a-z]*_//g; s/:[0-9\.][0-9\.]*//g; s/;//; /^ *$/d' \
 	../rn6.20way.nh | xargs echo | sed -e 's/ //g' > tree_commas.nh
     # tree_commas.nh looks like:
 #  ((((((((((((rn6,mm10),micOch1),cavPor3),oryCun2),(((hg38,panTro5),rheMac8),
 #    tarSyr2)),(((canFam3,ailMel1),felCat8),bosTau8)),monDom5),ornAna2),
 #    (galGal5,melGal5)),xenTro7),danRer10),calMil1)
 
     # use phyloFit to create tree model (output is phyloFit.mod)
     time /cluster/bin/phast.build/cornellCVS/phast.2010-12-30/bin/phyloFit \
 	    --EM --precision MED --msa-format FASTA --subst-mod REV \
 		--tree tree_commas.nh 4d.all.mfa
     #  real    31m34.749s
 
 
     mv phyloFit.mod all.mod
 
     grep TREE all.mod
 # ((((((((((((rn6:0.0895665,mm10:0.0844801):0.0651926,
 #    micOch1:0.125658):0.158112,cavPor3:0.215578):0.0359592,
 #    oryCun2:0.199508):0.015313,(((hg38:0.00685733,panTro5:0.0068924):0.0246446,
 #    rheMac8:0.0361324):0.0780051,tarSyr2:0.139334):0.0250771):0.0171773,
 #    (((canFam3:0.088804,ailMel1:0.077809):0.0216705,
 #    felCat8:0.0896594):0.0508452,bosTau8:0.167604):0.0287444):0.244728,
 #    monDom5:0.336766):0.0739297,ornAna2:0.491447):0.142595,
 #    (galGal5:0.0589443,melGal5:0.0696907):0.398276):0.146609,
 #    xenTro7:0.856129):0.151588,danRer10:1.03966):0.340195,calMil1:0.340195);
 
     # compare these calculated lengths to the tree extracted from 191way:
     grep TREE all.mod | sed -e 's/TREE: //' \
        | /cluster/bin/phast/all_dists /dev/stdin | grep rn6 \
           | sed -e "s/rn6.//;"  | sort > new.dists
     /cluster/bin/phast/all_dists ../rn6.20way.nh | grep rn6 \
         | sed -e "s/rn6.//;" | sort > old.dists
      # printing out the 'new', the 'old' the 'difference' and percent difference
     join new.dists old.dists | awk '{
   printf "#\t%s\t%8.6f\t%8.6f\t%8.6f\t%8.6f\n", $1, $2, $3, $2-$3, 100*($2-$3)/$3 }' \
       | sort -k3n
 #       mm10    0.174047        0.176098        -0.002051       -1.164692
 #       micOch1 0.280417        0.279362        0.001055        0.377646
 #       hg38    0.498727        0.509471        -0.010744       -2.108854
 #       panTro5 0.498762        0.509761        -0.010999       -2.157678
 #       rheMac8 0.503358        0.517098        -0.013740       -2.657137
 #       cavPor3 0.528449        0.498283        0.030166        6.053989
 #       tarSyr2 0.528554        0.510977        0.017577        3.439881
 #       oryCun2 0.548338        0.563940        -0.015602       -2.766606
 #       felCat8 0.550570        0.576238        -0.025668       -4.454409
 #       ailMel1 0.560390        0.577626        -0.017236       -2.983938
 #       canFam3 0.571385        0.550084        0.021301        3.872318
 #       bosTau8 0.577669        0.606134        -0.028465       -4.696156
 #       monDom5 0.962815        0.983334        -0.020519       -2.086677
 #       ornAna2 1.191425        1.170804        0.020621        1.761268
 #       galGal5 1.299794        1.383158        -0.083364       -6.027077
 #       melGal5 1.310540        1.427622        -0.117082       -8.201191
 #       calMil1 1.821160        1.935070        -0.113910       -5.886609
 #       xenTro7 1.845311        1.221660        0.623651        51.049474
 #       danRer10  2.180430      2.286142        -0.105712       -4.624035
 
 #########################################################################
 # phastCons 20-way (DONE - 2017-01-22 - Hiram)
     # split 20way mafs into 10M chunks and generate sufficient statistics
     # files for # phastCons
     ssh ku
     mkdir -p /hive/data/genomes/rn6/bed/multiz20way/cons/SS
     cd /hive/data/genomes/rn6/bed/multiz20way/cons/SS
     mkdir result done
 
     printf '#!/bin/csh -ef
 set d = $1
 set c = $2
 set doneDir = done/$d
 set MAF = /hive/data/genomes/rn6/bed/multiz20way/anno/result/$d/$c.maf
 set WINDOWS = /hive/data/genomes/rn6/bed/multiz20way/cons/SS/result/$d/$c
 set WC = `cat $MAF | wc -l`
 set NL = `grep "^#" $MAF | wc -l`
 if ( -s $3 ) then
     exit 0
 endif
 if ( -s $3.running ) then
     exit 0
 endif
 
 /bin/mkdir -p $doneDir
 /bin/date >> $3.running
 
 /bin/rm -fr $WINDOWS
 /bin/mkdir -p $WINDOWS
 pushd $WINDOWS > /dev/null
 if ( $WC != $NL ) then
 /cluster/bin/phast.build/cornellCVS/phast.2010-12-30/bin/msa_split \\
     $MAF -i MAF -o SS -r $WINDOWS/$c -w 10000000,0 -I 1000 -B 5000
 endif
 popd > /dev/null
 /bin/date >> $3
 /bin/rm -f $3.running
 ' > mkSS.csh
 
     chmod +x mkSS.csh
 
     printf '#LOOP
 mkSS.csh $(dir1) $(root1) {check out line+ done/$(dir1)/$(root1)}
 #ENDLOOP
 ' > template
 
     find ../../anno/result -type f | sed -e "s#../../anno/result/##" > maf.list
     wc -l maf.list
 # 624 maf.list
 
     ssh ku
     cd /hive/data/genomes/rn6/bed/multiz20way/cons/SS
 
     gensub2 maf.list single template jobList
     # beware overwhelming the cluster with these quick high I/O jobs
     para create jobList
     para try ... check ... etc
     para push
 # Completed: 624 of 624 jobs
 # CPU time in finished jobs:       5547s      92.45m     1.54h    0.06d  0.000 y
 # IO & Wait Time:                  1900s      31.67m     0.53h    0.02d  0.000 y
 # Average job time:                  12s       0.20m     0.00h    0.00d
 # Longest finished job:             681s      11.35m     0.19h    0.01d
 # Submission to last job:           733s      12.22m     0.20h    0.01d
 
     find ./result -type f | wc -l
     # 638
 
     # Run phastCons
     #	This job is I/O intensive in its output files, beware where this
     #	takes place or do not run too many at once.
     ssh ku
     mkdir -p /hive/data/genomes/rn6/bed/multiz20way/cons/run.cons
     cd /hive/data/genomes/rn6/bed/multiz20way/cons/run.cons
 
     #	This is setup for multiple runs based on subsets, but only running
     #   the 'all' subset here.
     #   It triggers off of the current working directory
     #	$cwd:t which is the "grp" in this script.  Running:
     #	all and vertebrates
 
     printf '#!/bin/csh -fe
 set PHASTBIN = /cluster/bin/phast.build/cornellCVS/phast.2010-12-30/bin
 set c = $1
 set d = $2
 set f = $3
 set len = $4
 set cov = $5
 set rho = $6
 set grp = $cwd:t
 set cons = /hive/data/genomes/rn6/bed/multiz20way/cons
 set tmp = $cons/tmp/${d}_${c}
 mkdir -p $tmp
 set ssSrc = $cons/SS/result
 set useGrp = "$grp.mod"
 if (-s $cons/$grp/$grp.non-inf) then
   ln -s $cons/$grp/$grp.mod $tmp
   ln -s $cons/$grp/$grp.non-inf $tmp
   ln -s $ssSrc/$d/$f $tmp
 else
   ln -s $ssSrc/$d/$f $tmp
   ln -s $cons/$grp/$grp.mod $tmp
 endif
 pushd $tmp > /dev/null
 if (-s $grp.non-inf) then
   $PHASTBIN/phastCons $f $useGrp \
     --rho $rho --expected-length $len --target-coverage $cov --quiet \\
     --not-informative `cat $grp.non-inf` \\
     --seqname $c --idpref $c --most-conserved $c.bed --score > $c.pp
 else
   $PHASTBIN/phastCons $f $useGrp \\
     --rho $rho --expected-length $len --target-coverage $cov --quiet \\
     --seqname $c --idpref $c --most-conserved $c.bed --score > $c.pp
 endif
 popd > /dev/null
 mkdir -p pp/$d bed/$d
 sleep 4
 touch pp/$d bed/$d
 rm -f pp/$d/$c.pp
 rm -f bed/$d/$c.bed
 mv $tmp/$c.pp pp/$d
 mv $tmp/$c.bed bed/$d
 rm -fr $tmp
 rmdir --ignore-fail-on-non-empty $cons/tmp/$d:h
 ' > doPhast.csh
 
     chmod +x doPhast.csh
 
     #	this template will serve for all runs
     #	root1 == chrom name, file1 == ss file name without .ss suffix
     printf '#LOOP
 ../run.cons/doPhast.csh $(root1) $(dir1) $(file1) 45 0.3 0.3 {check out line+ pp/$(dir1)/$(root1).pp}
 #ENDLOOP
 ' > template
 
     find ../SS/result -type f | sed -e "s#../SS/result/##" > ss.list
     wc -l ss.list
     #	638 ss.list
 
     # Create parasol batch and run it
     # run for all species
     cd /hive/data/genomes/rn6/bed/multiz20way/cons
     mkdir -p all
     cd all
     #	Using the .mod tree
     cp -p ../../4d/all.mod ./all.mod
 
     gensub2 ../run.cons/ss.list single ../run.cons/template jobList
     para -ram=32g create jobList
     para try ... check ...
     para push
 # Completed: 638 of 638 jobs
 # CPU time in finished jobs:      10814s     180.24m     3.00h    0.13d  0.000 y
 # IO & Wait Time:                  4534s      75.56m     1.26h    0.05d  0.000 y
 # Average job time:                  24s       0.40m     0.01h    0.00d
 # Longest finished job:              57s       0.95m     0.02h    0.00d
 # Submission to last job:           101s       1.68m     0.03h    0.00d
 
     # create Most Conserved track
     cd /hive/data/genomes/rn6/bed/multiz20way/cons/all
     time cut -f1 ../../../../chrom.sizes | while read C
 do
     ls -d bed/?/?/${C} 2> /dev/null | while read D
     do
         echo ${D}/${C}*.bed 1>&2
         cat ${D}/${C}*.bed
     done | sort -k1,1 -k2,2n \
     | awk '{printf "%s\t%d\t%d\tlod=%d\t%s\n", "'${C}'", $2, $3, $5, $5;}'
 done > tmpMostConserved.bed
     # real    0m40.602s
 
     time /cluster/bin/scripts/lodToBedScore tmpMostConserved.bed \
          > mostConserved.bed
     # real    0m23.400s
 
     # -rw-rw-r--  1 92173490 Jan 23 09:13 mostConserved.bed
     # load into database
     ssh hgwdev
     cd /hive/data/genomes/rn6/bed/multiz20way/cons/all
     time hgLoadBed rn6 phastConsElements20way mostConserved.bed
     # Read 2626510 elements of size 5 from mostConserved.bed
     # real    0m24.079s
 
     # on human we often try for 5% overall cov, and 70% CDS cov
     # most bets are off here for that goal, these alignments are too few
     #	and too far between
     #	--rho 0.3 --expected-length 45 --target-coverage 0.3
     time featureBits rn6 -enrichment ncbiRefSeq:cds phastConsElements20way
 # ncbiRefSeq:cds 1.217%, phastConsElements20way 3.976%, both 0.819%,
 #   cover 67.27%, enrich 16.92x
     # real    2m33.330s
 
     # Create merged posterier probability file and wiggle track data files
     cd /hive/data/genomes/rn6/bed/multiz20way/cons/all
     mkdir downloads
 
     # the third sed fixes the chrom names, removing the partition extensions
     time (find ./pp -type f | sed -e "s#^./##; s#\.# d #g; s#-# m #;" \
 	| sort -k1,1 -k3,3n | sed -e "s# d #.#g; s# m #-#g;" | xargs cat \
 	| sed -e 's/\.[0-9][0-9]*-[0-9][0-9]* start/ start/' \
         | gzip -c > downloads/phastCons20way.wigFix.gz)
     #	real    21m37.207s
 
 # -rw-rw-r-- 1 2339247044 Jan 23 09:36 phastCons20way.wigFix.gz
 
     # check integrity of data with wigToBigWig
     time (zcat downloads/phastCons20way.wigFix.gz \
 	| wigToBigWig -verbose=2 stdin /hive/data/genomes/rn6/chrom.sizes \
 	    phastCons20way.bw) > bigWig.log 2>&1
     egrep "real|VmPeak" bigWig.log
     # pid=61873: VmPeak:    23109920 kB
     # real    30m44.092s
 
     bigWigInfo phastCons20way.bw | sed -e 's/^/# /;'
 # version: 4
 # isCompressed: yes
 # isSwapped: 0
 # primaryDataSize: 3,494,424,791
 # primaryIndexSize: 78,132,788
 # zoomLevels: 10
 # chromCount: 368
 # basesCovered: 2,078,464,015
 # mean: 0.135364
 # min: 0.000000
 # max: 1.000000
 # std: 0.268417
 
     #	encode those files into wiggle data
     time (zcat downloads/phastCons20way.wigFix.gz \
 	| wigEncode stdin phastCons20way.wig phastCons20way.wib)
     # Converted stdin, upper limit 1.00, lower limit 0.00
     #  real    10m54.962s
 
     du -hsc *.wi?
     # 2.0G    phastCons20way.wib
     # 233M    phastCons20way.wig
 
     # Load gbdb and database with wiggle.
     ln -s `pwd`/phastCons20way.wib /gbdb/rn6/multiz20way/phastCons20way.wib
     time hgLoadWiggle -pathPrefix=/gbdb/rn6/multiz20way \
 	rn6 phastCons20way phastCons20way.wig
     #   real    0m29.148s
 
     # use to set trackDb.ra entries for wiggle min and max
     # and verify table is loaded correctly
 
     wigTableStats.sh rn6 phastCons20way
 # db.table               min max mean       count sumData
 # rn6.phastCons20way   0 1 0.135364 2078464015 2.8135e+08
 #    0.268417 viewLimits=0:1
 
     #  Create histogram to get an overview of all the data
     time hgWiggle -doHistogram -db=rn6 \
 	-hBinSize=0.001 -hBinCount=1000 -hMinVal=0.0 -verbose=2 \
 	    phastCons20way > histogram.data 2>&1
     #	real    1m55.313s
 
     #	create plot of histogram:
 
     printf 'set terminal png small x000000 xffffff xc000ff x66ff66 xffff00 x00ffff font \
 "/usr/share/fonts/default/Type1/n022004l.pfb"
 set size 1.4, 0.8
 set key left box
 set grid noxtics
 set grid ytics
 set title " Rat rn6 Histogram phastCons20way track"
 set xlabel " phastCons20way score"
 set ylabel " Relative Frequency"
 set y2label " Cumulative Relative Frequency (CRF)"
 set y2range [0:1]
 set y2tics
 set yrange [0:0.02]
 
 plot "histogram.data" using 2:5 title " RelFreq" with impulses, \\
         "histogram.data" using 2:7 axes x1y2 title " CRF" with lines
 ' | gnuplot > histo.png
 
     display histo.png &
 
 #########################################################################
 # phyloP for 20-way (DONE - 2016-05-09,11 - Hiram)
     # run phyloP with score=LRT
     ssh ku
     mkdir /cluster/data/rn6/bed/multiz20way/consPhyloP
     cd /cluster/data/rn6/bed/multiz20way/consPhyloP
 
     mkdir run.phyloP
     cd run.phyloP
     # Adjust model file base composition background and rate matrix to be
     # representative of the chromosomes in play
     grep BACKGROUND ../../4d/all.mod | awk '{printf "%0.3f\n", $3 + $4}'
     #	0.536
     /cluster/bin/phast.build/cornellCVS/phast.2010-12-30/bin/modFreqs \
 	../../4d/all.mod 0.536 > all.mod
     # verify, the BACKGROUND should now be paired up:
     grep BACK all.mod
     #   BACKGROUND: 0.232000 0.268000 0.268000 0.232000 
 
     printf '#!/bin/csh -fe
 set PHASTBIN = /cluster/bin/phast.build/cornellCVS/phast.2010-12-30/bin
 set f = $1
 set d = $f:h
 set file1 = $f:t
 set out = $2
 set cName = $f:t:r
 set grp = $cwd:t
 set cons = /hive/data/genomes/rn6/bed/multiz20way/consPhyloP
 set tmp = $cons/tmp/$grp/$f
 /bin/rm -fr $tmp
 /bin/mkdir -p $tmp
 set ssSrc = "/hive/data/genomes/rn6/bed/multiz20way/cons/SS/result/$f"
 set useGrp = "$grp.mod"
 /bin/ln -s $cons/run.phyloP/$grp.mod $tmp
 pushd $tmp > /dev/null
 $PHASTBIN/phyloP --method LRT --mode CONACC --wig-scores --chrom $cName \\
     -i SS $useGrp $ssSrc.ss > $file1.wigFix
 popd > /dev/null
 /bin/mkdir -p $out:h
 sleep 4
 /bin/touch $out:h
 /bin/mv $tmp/$file1.wigFix $out
 /bin/rm -fr $tmp
 /bin/rmdir --ignore-fail-on-non-empty $cons/tmp/$grp/$d
 /bin/rmdir --ignore-fail-on-non-empty $cons/tmp/$grp/$d:h
 /bin/rmdir --ignore-fail-on-non-empty $cons/tmp/$grp
 /bin/rmdir --ignore-fail-on-non-empty $cons/tmp
 ' > doPhyloP.csh
 
     chmod +x doPhyloP.csh
 
     # Create list of chunks
     find ../../cons/SS/result -type f | grep ".ss$" \
 	| sed -e "s/.ss$//; s#^../../cons/SS/result/##" > ss.list
     # make sure the list looks good
     wc -l ss.list
     #	638 ss.list
 
     # Create template file
     #	file1 == $chr/$chunk/file name without .ss suffix
     printf '#LOOP
 ../run.phyloP/doPhyloP.csh $(path1) {check out line+ wigFix/$(dir1)/$(file1).wigFix}
 #ENDLOOP
 ' > template
 
     ######################   Running all species  #######################
     # setup run for all species
     mkdir /hive/data/genomes/rn6/bed/multiz20way/consPhyloP/all
     cd /hive/data/genomes/rn6/bed/multiz20way/consPhyloP/all
     rm -fr wigFix
     mkdir wigFix
 
     gensub2 ../run.phyloP/ss.list single ../run.phyloP/template jobList
     #	beware overwhelming the cluster with these fast running high I/O jobs
     para create jobList
     para try ... check ... push ... etc ...
     para push
     para time > run.time
 # Completed: 638 of 638 jobs
 # CPU time in finished jobs:     277914s    4631.90m    77.20h    3.22d  0.009 y
 # IO & Wait Time:                  8871s     147.85m     2.46h    0.10d  0.000 y
 # Average job time:                 450s       7.49m     0.12h    0.01d
 # Longest finished job:            1737s      28.95m     0.48h    0.02d
 # Submission to last job:          3476s      57.93m     0.97h    0.04d
 
     mkdir downloads
 
     time (find ./wigFix -type f | sed -e "s#^./##; s#\.# d #g; s#-# m #;" \
 	| sort -k1,1 -k3,3n | sed -e "s# d #.#g; s# m #-#g;" | xargs cat \
 	| gzip -c > downloads/phyloP20way.wigFix.gz)
     #   real    32m19.962s
 
     # check integrity of data with wigToBigWig
     time (zcat downloads/phyloP20way.wigFix.gz \
 	| wigToBigWig -verbose=2 stdin /hive/data/genomes/rn6/chrom.sizes \
 	phyloP20way.bw) > bigWig.log 2>&1
 XXX - running - Mon Jan 23 19:45:33 PST 2017
     egrep "real|VmPeak" bigWig.log
     # pid=25924: VmPeak:    23109924 kB
     # real    30m51.433s
 
     bigWigInfo phyloP20way.bw | sed -e 's/^/# /;'
 # version: 4
 # isCompressed: yes
 # isSwapped: 0
 # primaryDataSize: 4,248,984,727
 # primaryIndexSize: 78,132,788
 # zoomLevels: 10
 # chromCount: 368
 # basesCovered: 2,078,464,015
 # mean: 0.108077
 # min: -8.546000
 # max: 3.488000
 # std: 0.723885
 
     #	encode those files into wiggle data
     time (zcat downloads/phyloP20way.wigFix.gz \
 	| wigEncode stdin phyloP20way.wig phyloP20way.wib)
     # Converted stdin, upper limit 3.49, lower limit -8.55
     # real    12m5.611s
 
 
     du -hsc *.wi?
     # 2.0G    phyloP20way.wib
     # 239M    phyloP20way.wig
 
     # Load gbdb and database with wiggle.
     ln -s `pwd`/phyloP20way.wib /gbdb/rn6/multiz20way/phyloP20way.wib
     time hgLoadWiggle -pathPrefix=/gbdb/rn6/multiz20way rn6 \
 	phyloP20way phyloP20way.wig
     # real    0m27.237s
 
     # use to set trackDb.ra entries for wiggle min and max
     # and verify table is loaded correctly
 
     wigTableStats.sh rn6 phyloP20way
 # db.table           min   max     mean      count sumData
 # rn6.phyloP20way -8.546 3.488 0.108077 2078464015 2.24635e+08
 #       stdDev viewLimits
 #     0.723885 viewLimits=-3.51135:3.488
 
     #	that range is: 8.546+3.488 = 12.034 for hBinSize=0.012034
 
     #  Create histogram to get an overview of all the data
     time hgWiggle -doHistogram \
 	-hBinSize=0.012034 -hBinCount=1000 -hMinVal=-8.546 -verbose=2 \
 	    -db=rn6 phyloP20way > histogram.data 2>&1
     # real    1m52.632s
 
     # find the Y range for the 2:5 graph
     grep -v chrom histogram.data | grep "^[0-9]" | ave -col=5 stdin \
       | sed -e 's/^/# /;'
 # Q1 0.000000
 # median 0.000045
 # Q3 0.000649
 # average 0.001104
 # min 0.000000
 # max 0.052600
 # count 906
 # total 0.999990
 # standard deviation 0.003101
 
     # find the X range for the 2:5 graph
     grep "^[0-9]" histogram.data | ave -col=2 stdin \
       | sed -e 's/^/# /;'
 # Q1 -4.695120
 # median -1.969420
 # Q3 0.756282
 # average -1.979713
 # min -8.546000
 # max 3.475970
 # count 906
 # total -1793.619961
 # standard deviation 3.167347
 
     #	create plot of histogram:
     printf 'set terminal png small x000000 xffffff xc000ff x66ff66 xffff00 x00ffff font \
 "/usr/share/fonts/default/Type1/n022004l.pfb"
 set size 1.4, 0.8
 set key left box
 set grid noxtics
 set grid ytics
 set title " Rat rn6 Histogram phyloP20way track"
 set xlabel " phyloP20way score"
 set ylabel " Relative Frequency"
 set y2label " Cumulative Relative Frequency (CRF)"
 set y2range [0:1]
 set y2tics
 set xtics
 set xrange [-2:2]
 set yrange [0:0.053]
 
 plot "histogram.data" using 2:5 title " RelFreq" with impulses, \
         "histogram.data" using 2:7 axes x1y2 title " CRF" with lines
 ' | gnuplot > histo.png
 
     display histo.png &
 
 #############################################################################
 # hgPal downloads (DONE - 2017-01-24 - Hiram)
 #   FASTA from 20-way for ensGene
 
     ssh hgwdev
     screen -S rn6HgPal
     mkdir /hive/data/genomes/rn6/bed/multiz20way/pal
     cd /hive/data/genomes/rn6/bed/multiz20way/pal
     cat ../species.list | tr '[ ]' '[\n]' > order.list
 
     # beware large contig count assemblies, this can take a long time (hours)
     export mz=multiz20way
     export gp=ensGene
     export db=rn6
     export I=0
     export D=0
     mkdir exonAA exonNuc
     printf '#!/bin/sh\n' > $gp.jobs
 
     time for C in `sort -nk2 ../../../chrom.sizes | cut -f1`
     do
         I=`echo $I | awk '{print $1+1}'`
         D=`echo $D | awk '{print $1+1}'`
         dNum=`echo $D | awk '{printf "%03d", int($1/1000)}'`
         mkdir -p exonNuc/${dNum} > /dev/null
         mkdir -p exonAA/${dNum} > /dev/null
 	echo "mafGene -chrom=$C -exons -noTrans $db $mz $gp order.list stdout | gzip -c > exonNuc/${dNum}/$C.exonNuc.fa.gz &"
 	echo "mafGene -chrom=$C -exons $db $mz $gp order.list stdout | gzip -c > exonAA/${dNum}/$C.exonAA.fa.gz &"
         if [ $I -gt 16 ]; then
             echo "date"
             echo "wait"
             I=0
         fi
     done >> $gp.jobs
     # real    0m21.771s
 
     echo "date" >> $gp.jobs
     echo "wait" >> $gp.jobs
 
     chmod +x  $gp.jobs
 
     time (./$gp.jobs) > $gp.jobs.log 2>&1 &
     # real    28m28.221s
 
     export mz=multiz20way
     export gp=ensGene
     time find ./exonAA -type f | grep exonAA.fa.gz | xargs zcat \
      | gzip -c > $gp.$mz.exonAA.fa.gz
     #  real    0m29.606s
 
     time find ./exonNuc -type f | grep exonNuc.fa.gz | xargs zcat \
      | gzip -c > $gp.$mz.exonNuc.fa.gz
     #   real    2m27.672s
 
 # -rw-rw-r-- 1 133052531 Jan 24 08:13 ensGene.multiz20way.exonAA.fa.gz
 # -rw-rw-r-- 1 231318306 Jan 24 08:16 ensGene.multiz20way.exonNuc.fa.gz
 
     export mz=multiz20way
     export gp=ensGene
     export db=rn6
     export pd=/usr/local/apache/htdocs-hgdownload/goldenPath/$db/$mz/alignments
     mkdir -p $pd
     md5sum *.fa.gz > md5sum.txt
     ln -s `pwd`/$gp.$mz.exonAA.fa.gz $pd/$gp.exonAA.fa.gz
     ln -s `pwd`/$gp.$mz.exonNuc.fa.gz $pd/$gp.exonNuc.fa.gz
     ln -s `pwd`/md5sum.txt $pd/
 
     rm -rf exonAA exonNuc
 
 #############################################################################
 # construct download files for 20-way (DONE - 2017-01-24 - Hiram)
     mkdir /usr/local/apache/htdocs-hgdownload/goldenPath/rn6/multiz20way
     mkdir /usr/local/apache/htdocs-hgdownload/goldenPath/rn6/phastCons20way
     mkdir /usr/local/apache/htdocs-hgdownload/goldenPath/rn6/phyloP20way
     mkdir /hive/data/genomes/rn6/bed/multiz20way/downloads
     cd /hive/data/genomes/rn6/bed/multiz20way/downloads
     printf '#!/bin/sh
 
 for db in `cat ../species.list`
 do
   orgName=`hgsql -N -e "select organism from dbDb where name=\\"${db}\\";" hgcentraltest`
   sciName=`hgsql -N -e "select scientificName from dbDb where name=\\"${db}\\";" hgcentraltest`
   genome=`hgsql -N -e "select genome from dbDb where name=\\"${db}\\";" hgcentraltest`
   descr=`hgsql -N -e "select description from dbDb where name=\\"${db}\\";" hgcentraltest`
   printf \"%%-15s %%-24s %%-20s\\n\" \"$genome\" \"$sciName\" \"$descr\"
 done
 ' > docInfo.sh
     chmod +x docInfo.sh
     ./docInfo.sh
 # Rat             Rattus norvegicus        Jul. 2014 (RGSC 6.0/rn6)
 # Mouse           Mus musculus             Dec. 2011 (GRCm38/mm10)
 # Prairie vole    Microtus ochrogaster     Oct 2012 (MicOch1.0/micOch1)
 # Guinea pig      Cavia porcellus          Feb. 2008 (Broad/cavPor3)
 # Rabbit          Oryctolagus cuniculus    Apr. 2009 (Broad/oryCun2)
 # Human           Homo sapiens             Dec. 2013 (GRCh38/hg38)
 # Chimp           Pan troglodytes          May 2016 (Pan_tro 3.0/panTro5)
 # Rhesus          Macaca mulatta           Nov. 2015 (BCM Mmul_8.0.1/rheMac8)
 # Tarsier         Tarsius syrichta         Sep. 2013 (Tarsius_syrichta-2.0.1/tarSyr2)
 # Dog             Canis lupus familiaris   Sep. 2011 (Broad CanFam3.1/canFam3)
 # Panda           Ailuropoda melanoleuca   Dec. 2009 (BGI-Shenzhen 1.0/ailMel1)
 # Cat             Felis catus              Nov. 2014 (ICGSC Felis_catus_8.0/felCat8)
 # Cow             Bos taurus               Jun. 2014 (Bos_taurus_UMD_3.1.1/bosTau8)
 # Opossum         Monodelphis domestica    Oct. 2006 (Broad/monDom5)
 # Platypus        Ornithorhynchus anatinus Feb. 2007 (ASM227v2/ornAna2)
 # Chicken         Gallus gallus            Dec 2015 (Gallus_gallus-5.0/galGal5)
 # Turkey          Meleagris gallopavo      Nov. 2014 (Turkey_5.0/melGal5)
 # X. tropicalis   Xenopus tropicalis       Sep. 2012 (JGI 7.0/xenTro7)
 # Zebrafish       Danio rerio              Sep. 2014 (GRCz10/danRer10)
 # Elephant shark  Callorhinchus milii      Dec. 2013 (Callorhinchus_milii-6.1.3/calMil1)
 
 
     printf '#!/bin/sh
 
 for db in `cat ../species.list`
 do
   orgName=`hgsql -N -e "select organism from dbDb where name=\\"${db}\\";" hgcentraltest`
   sciName=`hgsql -N -e "select scientificName from dbDb where name=\\"${db}\\";" hgcentraltest`
   genome=`hgsql -N -e "select genome from dbDb where name=\\"${db}\\";" hgcentraltest`
   descr=`hgsql -N -e "select description from dbDb where name=\\"${db}\\";" hgcentraltest`
   printf "<tr align=left><td>%%s</td><td>%%s</td><td>\\n    %%s</td><td><a href=\\"../cgi-bin/hgTracks?db=%%s\\"\\n    target=_blank>%%s</a></td><td>%%s</td></tr>\\n" "$genome" "$sciName" "$descr" "${db}" "$descr" "MAF Net"
 done
 ' > descrPageInfo.sh
 
     chmod +x descrPageInfo.sh
     # use this to add the table in the description page
     # it doesn't check the maf type, that needs to be fixed
 
     mkdir multiz20way phastCons20way phyloP20way
     cd multiz20way
     time cp -p ../../anno/rn6.20way.maf .
     #   real    1m13.217s
 
     # -rw-rw-r-- 1 48341214170 Jan 22 10:39 rn6.20way.maf
 
     du -hsc *
     #  46G     rn6.20way.maf
 
     time gzip *.maf
     #   real    69m38.345s
 
     # -rw-rw-r-- 1 7725464919 Jan 22 10:39 rn6.20way.maf.gz
 
     du -hsc *.maf.gz
     #  7.2G    rn6.20way.maf.gz
 
     grep TREE ../../4d/all.mod | awk '{print $NF}' \
       | ~/kent/src/hg/utils/phyloTrees/asciiTree.pl /dev/stdin \
          > rn6.20way.nh
     ~/kent/src/hg/utils/phyloTrees/commonNames.sh rn6.20way.nh \
       | ~/kent/src/hg/utils/phyloTrees/asciiTree.pl /dev/stdin \
          > rn6.20way.commonNames.nh
     ~/kent/src/hg/utils/phyloTrees/scientificNames.sh rn6.20way.nh \
 	| $HOME/kent/src/hg/utils/phyloTrees/asciiTree.pl /dev/stdin \
 	    > rn6.20way.scientificNames.nh
 
     du -hsc *.maf.gz ../../anno/rn6.20way.maf
     # 7.2G    rn6.20way.maf.gz
     # 46G     ../../anno/rn6.20way.maf
 
     ### construct the upstream MAF files for ensGene
     # bash script
 #!/bin/sh
 export geneTbl="ensGene"
 for S in 1000 2000 5000
 do
     echo "making upstream${S}.maf"
     featureBits rn6 ${geneTbl}:upstream:${S} -fa=/dev/null -bed=stdout \
         | perl -wpe 's/_up[^\t]+/\t0/' | sort -k1,1 -k2,2n \
         | /cluster/bin/$MACHTYPE/mafFrags rn6 multiz20way \
                 stdin stdout \
                 -orgs=/hive/data/genomes/rn6/bed/multiz20way/species.list \
         | gzip -c > upstream${S}.${geneTbl}.maf.gz
     echo "done upstream${S}.${geneTbl}.maf.gz"
 done
     # real    18m28.217s
 
     time md5sum *.nh *.maf.gz > md5sum.txt
     #   real    0m26.660s
 
     # obtain the README.txt from galVar1/multiz5way and update for this
     #   situation
 
     ln -s `pwd`/*.nh `pwd`/*.maf.gz `pwd`/*.txt \
         /usr/local/apache/htdocs-hgdownload/goldenPath/rn6/multiz20way
 
     #####################################################################
     cd /hive/data/genomes/rn6/bed/multiz20way/downloads/phastCons20way
 
     ln -s ../../cons/all/downloads/phastCons20way.wigFix.gz \
         ./rn6.phastCons20way.wigFix.gz
     ln -s ../../cons/all/phastCons20way.bw ./rn6.phastCons20way.bw
     ln -s ../../cons/all/all.mod ./rn6.phastCons20way.mod
     time md5sum *.gz *.mod *.bw > md5sum.txt
     #   real    0m20.354s
 
     # obtain the README.txt from galVar1/phastCons5way and update for this
     #   situation
 
     ln -s `pwd`/* \
       /usr/local/apache/htdocs-hgdownload/goldenPath/rn6/phastCons20way
 
     #####################################################################
     cd /hive/data/genomes/rn6/bed/multiz20way/downloads/phyloP20way
 
     ln -s ../../consPhyloP/all/downloads/phyloP20way.wigFix.gz \
         ./rn6.phyloP20way.wigFix.gz
     ln -s ../../consPhyloP/run.phyloP/all.mod rn6.phyloP20way.mod
     ln -s ../../consPhyloP/all/phyloP20way.bw rn6.phyloP20way.bw
 
     time md5sum *.mod *.bw *.gz > md5sum.txt
     #   real    0m29.648s
 
     # obtain the README.txt from galVar1/phyloP5way and update for this
     #   situation
     ln -s `pwd`/* \
       /usr/local/apache/htdocs-hgdownload/goldenPath/rn6/phyloP20way
 
 #############################################################################
 # wiki page for 20-way (DONE - 2017-01-24 - Hiram)
     mkdir /hive/users/hiram/bigWays/rn6.20way
     cd /hive/users/hiram/bigWays
     echo "rn6" > rn6.20way/ordered.list
     awk '{print $1}' /hive/data/genomes/rn6/bed/multiz20way/20way.distances.txt \
        >> rn6.20way/ordered.list
 
     # sizeStats.sh catches up the cached measurements required for data
     # in the tables.  They are usually already mostly done, only new
     # assemblies will have updates.
     ./sizeStats.sh rn6.20way/ordered.list
     # dbDb.sh constructs rn6.20way/Rn6_20-way_conservation_alignment.html
     # may need to add new assembly references to srcReference.list and
     # urlReference.list
     ./dbDb.sh rn6 20way
     # sizeStats.pl constructs rn6.20way/Rn6_20-way_Genome_size_statistics.html
     # this requires entries in coverage.list for new sequences
     ./sizeStats.pl rn6 20way
 
     # defCheck.pl constructs Rn6_20-way_conservation_lastz_parameters.html
     ./defCheck.pl rn6 20way
 
     # this constructs the html pages in rn6.20way/:
 # -rw-rw-r-- 1 3317 May 11 11:08 Rn6_20-way_conservation_alignment.html
 # -rw-rw-r-- 1 4827 May 11 11:08 Rn6_20-way_Genome_size_statistics.html
 # -rw-rw-r-- 1 3308 May 11 11:08 Rn6_20-way_conservation_lastz_parameters.html
     # add those pages to the genomewiki.  Their page names are the
     # names of the .html files without the .html:
 #  Rn6_20-way_conservation_alignment
 #  Rn6_20-way_Genome_size_statistics
 #  Rn6_20-way_conservation_lastz_parameters
 
     # when you view the first one you enter, it will have links to the
     # missing two.
 
 ############################################################################
 # composite chainNet trackDb (DONE - 2017-01-27 - Hiram)
 
     cd ~/kent/src/hg/makeDb/trackDb/rat/rn6
 
     printf 'mm10
 micOch1
 cavPor3
 oryCun2
 ' > glires.list
 
     printf 'hg38
 panTro5
 tarSyr2
 rheMac8
 ' > primates.list
 
     printf 'canFam3
 felCat8
 ailMel1
 bosTau8
 ' > laurasiatheria.list
 
     printf 'monDom5
 ornAna2
 xenTro7
 ' > sarcopterygii.list
 
     printf 'galGal5
 melGal5
 ' > birds.list
 
     printf 'calMil1
 danRer10' > fish.list
 
     $HOME/kent/src/hg/utils/phyloTrees/chainNetCompositeTrackDb.pl \
        vertebrate glires.list primates.list \
           laurasiatheria.list sarcopterygii.list birds.list fish.list \
               > vertebrateChainNet.ra
 
     # edit that result to fixup the shortLabel, longLabel and default
     # visibilities desired on some of the tracks
 
     # insert this line in the trackDb.ra to get it included:
 
     include vertebrateChainNet.ra
 
 ############################################################################
 # crispr 10K track (DONE - 2018-07-02 - Hiram)
     mkdir -p /hive/data/genomes/rn6/bed/crispr.10K
     cd /hive/data/genomes/rn6/bed/crispr.10K
 
     # since this time, the indexing has been built into the primary script
 export PATH=/hive/data/outside/crisprTrack/crispor/tools/usrLocalBin:$PATH
 export TMPDIR=/dev/shm
 
 time (/hive/data/outside/crisprTrack/crispor/tools/crisprAddGenome ucscLocal \
    rn6 --geneTable ensGene \
       --baseDir /hive/data/outside/crisprTrack/crispor/genomes) > index.log 2>&1
 
     # real    109m47.348s
 
     time (/cluster/home/hiram/kent/src/hg/utils/automation/doCrispr.pl \
    -buildDir=`pwd` -bigClusterHub=ku -smallClusterHub=ku \
   -continue=guides -stop=guides -workhorse=hgwdev rn6 ensGene) > guides.log 2>&1
     # real    18m23.512s
 
     # step specScores were run manually to test fixes in the script to avoid
     # the bottleneck
     /cluster/home/hiram/kent/src/hg/utils/automation/doCrispr.pl \
    -buildDir=`pwd` -bigClusterHub=ku -smallClusterHub=ku -debug \
         -continue=specScores -stop=specScores -workhorse=hgwdev rn6 ensGene
 
     /cluster/home/hiram/kent/src/hg/utils/automation/doCrispr.pl -debug \
    -buildDir=`pwd` -bigClusterHub=ku -smallClusterHub=ku \
         -continue=effScores -stop=offTargets -workhorse=hgwdev rn6 ensGene)
 
     /cluster/home/hiram/kent/src/hg/utils/automation/doCrispr.pl -debug \
    -buildDir=`pwd` -bigClusterHub=ku -smallClusterHub=ku \
         -continue=load -stop=load -workhorse=hgwdev rn6 ensGene
 
 ############################################################################
+# LIFTOVER TO regenRn0 (DONE - 2019-06-29 - Hiram)
+    ssh hgwdev
+    mkdir /hive/data/genomes/rn6/bed/blat.regenRn0.2019-06-29
+    cd /hive/data/genomes/rn6/bed/blat.regenRn0.2019-06-29
+    doSameSpeciesLiftOver.pl -verbose=2 \
+        -debug -bigClusterHub=ku -dbHost=hgwdev -workhorse=hgwdev \
+        -ooc=/hive/data/genomes/rn6/jkStuff/rn6.11.ooc \
+         rn6 regenRn0
+    time (doSameSpeciesLiftOver.pl -verbose=2 \
+        -bigClusterHub=ku -dbHost=hgwdev -workhorse=hgwdev \
+        -ooc=/hive/data/genomes/rn6/jkStuff/rn6.11.ooc \
+         rn6 regenRn0) > doLiftOverToRegenRn0.log 2>&1
+    # real    249m30.284s
+
+    # see if the liftOver menus function in the browser from regenRn0 to galGal5
+
+    # would like to see this as a track:
+    chainToPsl rn6ToRegenRn0.over.chain.gz ../../chrom.sizes \
+   /hive/data/genomes/regenRn0/chrom.sizes ../../rn6.2bit \
+   /hive/data/genomes/regenRn0/regenRn0.2bit rn6ToRegenRn0.psl
+
+#########################################################################
+# LIFTOVER TO regenRn1 (DONE - 2019-09-09 - Hiram)
+    ssh hgwdev
+    mkdir /hive/data/genomes/rn6/bed/blat.regenRn1.2019-09-09
+    cd /hive/data/genomes/rn6/bed/blat.regenRn1.2019-09-09
+    doSameSpeciesLiftOver.pl -verbose=2 \
+        -debug -bigClusterHub=ku -dbHost=hgwdev -workhorse=hgwdev \
+        -ooc=/hive/data/genomes/rn6/jkStuff/rn6.11.ooc \
+         rn6 regenRn1
+    time (doSameSpeciesLiftOver.pl -verbose=2 \
+        -bigClusterHub=ku -dbHost=hgwdev -workhorse=hgwdev \
+        -ooc=/hive/data/genomes/rn6/jkStuff/rn6.11.ooc \
+         rn6 regenRn1) > doLiftOverToRegenRn1.log 2>&1
+    # real    252m38.421s
+
+    # see if the liftOver menus function in the browser from rn6 to regenRn1
+
+    # would like to see this as a track:
+    hgLoadChain rn6 chainRegenRn1 rn6ToRegenRn1.over.chain.gz
+    # Loading 237485 chains into rn6.chainRegenRn1
+
+    chainSort rn6ToRegenRn1.over.chain.gz stdout \
+      | chainPreNet stdin \
+        /hive/data/genomes/rn6/chrom.sizes \
+        /hive/data/genomes/regenRn1/chrom.sizes stdout \
+      | chainNet stdin -minSpace=1 /hive/data/genomes/rn6/chrom.sizes \
+        /hive/data/genomes/regenRn1/chrom.sizes stdout /dev/null \
+      | netSyntenic stdin noClass.net
+
+   netClass -verbose=0 -noAr noClass.net rn6 regenRn1 rn6.regenRn1.net
+
+    netFilter -minGap=10 rn6.regenRn1.net \
+      | hgLoadNet -verbose=0 rn6 netRegenRn1 stdin
+
+featureBits rn6 chainRegenRn1Link > fb.rn6.chainRegenRn1Link.txt 2>&1
+
+    cat fb.rn6.chainRegenRn1Link.txt
+    # 2265178557 bases of 2729860805 (82.978%) in intersection
+
+#########################################################################