src/hg/makeDb/doc/ailMel1.txt 1.14

1.14 2010/04/30 21:17:30 chinhli
run doRecipBest.pl for felCatV17e
Index: src/hg/makeDb/doc/ailMel1.txt
===================================================================
RCS file: /projects/compbio/cvsroot/kent/src/hg/makeDb/doc/ailMel1.txt,v
retrieving revision 1.13
retrieving revision 1.14
diff -b -B -U 1000000 -r1.13 -r1.14
--- src/hg/makeDb/doc/ailMel1.txt	20 Apr 2010 00:06:17 -0000	1.13
+++ src/hg/makeDb/doc/ailMel1.txt	30 Apr 2010 21:17:30 -0000	1.14
@@ -1,578 +1,583 @@
 # for emacs: -*- mode: sh; -*-
 
 #	$Id$
 
 # Panda sequence: http://panda.genomics.org.cn/page/panda/download.jsp
 #	Ailuropoda melanoleuca
 
 ##########################################################################
 # Download sequence (DONE - 2009-12-17 - Hiram)
     mkdir /hive/data/genomes/ailMel1
     cd /hive/data/genomes/ailMel1
     mkdir genbank
     cd genbank
 
     wget --timestamping -r --cut-dirs=6 --level=0 -nH -x \
 	--no-remove-listing -np \
 "ftp.ncbi.nlm.nih.gov:genbank/genomes/Eukaryotes/vertebrates_mammals/Ailuropoda_melanoleuca/AilMel_1.0/*"
 
     #	first pass attempt at sequence from public.genomics.org.cn
     #	2009-12-17
     #	NCBI has cleaned the sequence and removed a couple of scaffolds that
     #	are cloning vector sequence
     mkdir download
     cd download
 
 wget --timestamping \
 'ftp://public.genomics.org.cn/BGI/panda/genome_sequence/panda.scafSeq.gapFilled.noMito.GeneScaffold.fa.gz'
 
 wget --timestamping \
 'ftp://public.genomics.org.cn/BGI/panda/genome_sequence/panda_contig100.mapping_scaffolds_nucmer.list.clean.tab.sort.gff.gz'
 
 wget --timestamping \
 'ftp://panda.genomics.org.cn/pub/panda/panda.scafSeq.gapFilled.noMito.GeneScaffold.fa.gz'
 
 
     faCount panda.scafSeq.gapFilled.noMito.GeneScaffold.fa.gz > faCount.txt
 
     grep "^scaff" faCount.txt | awk '{print $1}' > scaffold.list
 
     faSomeRecords panda.scafSeq.gapFilled.noMito.GeneScaffold.fa.gz \
 	scaffold.list stdout | gzip > scaffolds.fa.gz
 
     hgFakeAgp scaffolds.fa.gz stdout | sed -e "s/D/W/" > fake.agp
 
 ##########################################################################
 # Initial genome build (DONE - 2009-12-17 - Hiram)
     cd /hive/data/genomes/ailMel1
 
     cat << '_EOF_' > ailMel1.config.ra
 # Config parameters for makeGenomeDb.pl:
 db ailMel1
 clade mammal
 genomeCladePriority 19
 scientificName Ailuropoda melanoleuca
 commonName Panda
 assemblyDate Dec. 2009
 assemblyLabel BGI-Shenzhen AilMel 1.0 Dec. 2009
 orderKey 220
 mitoAcc NC_009492
 fastaFiles /hive/data/genomes/ailMel1/genbank/Primary_Assembly/unplaced_scaffolds/FASTA/unplaced.scaf.fa.gz
 agpFiles /hive/data/genomes/ailMel1/genbank/Primary_Assembly/unplaced_scaffolds/AGP/unplaced.scaf.agp.gz
 # qualFiles none
 dbDbSpeciesDir panda
 taxId 9646
 
     makeGenomeDb.pl -workhorse=hgwdev -stop=seq ailMel1.config.ra > seq.out 2>&1
 #	real    4m4.783s
     makeGenomeDb.pl -start=agp -stop=agp ailMel1.config.ra > agp.out 2>&1
 #	real    0m20.968s
     makeGenomeDb.pl -start=db -stop=db ailMel1.config.ra > db.out 2>&1
 #	real    6m28.168s
     makeGenomeDb.pl -start=dbDb -stop=dbDb ailMel1.config.ra > dbDb.out 2>&1
     makeGenomeDb.pl -start=trackDb -stop=trackDb ailMel1.config.ra > trackDb.out 2>&1
 
 ##########################################################################
 # running repeat masker (DONE - 2010-02-02 - Hiram)
     mkdir /hive/data/genomes/ailMel1/bed/repeatMasker
     cd /hive/data/genomes/ailMel1/bed/repeatMasker
     doRepeatMasker.pl -buildDir=`pwd` -noSplit -bigClusterHub=swarm \
 	-workhorse=hgwdev ailMel1 > do.log 2>&1
     cat faSize.rmsk.txt
 # 2299509015 bases (54196184 N's 2245312831 real 1371210941 upper 874101890
 #	lower) in 81467 sequences in 1 files
 # %38.01 masked total, %38.93 masked real
 
 ##########################################################################
 # running simple repeat (DONE - 2010-02-02 - Hiram)
     mkdir /hive/data/genomes/ailMel1/bed/simpleRepeat
     cd /hive/data/genomes/ailMel1/bed/simpleRepeat
     time doSimpleRepeat.pl -buildDir=`pwd` -smallClusterHub=swarm \
 	-workhorse=hgwdev ailMel1 > do.log 2>&1 &
 #	real    103m30.880s
     cat fb.simpleRepeat 
 #	25064758 bases of 2245312831 (1.116%) in intersection
 
     cd /hive/data/genomes/ailMel1
     twoBitMask ailMel1.rmsk.2bit \
 	-add bed/simpleRepeat/trfMask.bed ailMel1.2bit
     #	you can safely ignore the warning about fields >= 13
 
     twoBitToFa ailMel1.2bit stdout | faSize stdin > faSize.ailMel1.2bit.txt
     cat faSize.ailMel1.2bit.txt
 # 2299509015 bases (54196184 N's 2245312831 real 1370617891 upper
 #	874694940 lower) in 81467 sequences in 1 files
 #	%38.04 masked total, %38.96 masked real
 
     rm /gbdb/ailMel1/ailMel1.2bit
     ln -s `pwd`/ailMel1.2bit /gbdb/ailMel1/ailMel1.2bit
 
 ##########################################################################
 #  BLATSERVERS ENTRY (DONE - 2009-12-23 - 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 ("ailMel1", "blat12", "17802", "1", "0"); \
 	INSERT INTO blatServers (db, host, port, isTrans, canPcr) \
 	VALUES ("ailMel1", "blat12", "17803", "0", "1");' \
 	    hgcentraltest
     #	test it with some sequence
 
 #########################################################################
 # MAKE 11.OOC FILES FOR BLAT (DONE - 2010-02-03 - Hiram)
     ssh kolossus
     # numerator is ailMel1 gapless bases as reported by faSize 
     # denominator is hg17 gapless bases as reported by featureBits,
     # 1024 is threshold used for human -repMatch:
     calc \( 2245312831 / 2897310462 \) \* 1024
 #	( 2245312831 / 2897310462 ) * 1024 = 793.563675
     # ==> use -repMatch=750 according to size scaled down from 1024 for human.
     #	and rounded down to nearest 50
     cd /hive/data/genomes/ailMel1
     blat ailMel1.2bit /dev/null /dev/null -tileSize=11 \
       -makeOoc=jkStuff/ailMel1.11.ooc -repMatch=750
     #	Wrote 25467 overused 11-mers to jkStuff/ailMel1.11.ooc
 
     mkdir /hive/data/staging/data/ailMel
     cp -p ailMel1.2bit chrom.sizes jkStuff/ailMel1.11.ooc \
 	/hive/data/staging/data/ailMel1
 
 
 #########################################################################
 # genbank run (DONE - 2010-02-04 - Hiram)
     ssh hgwdev
     cd $HOME/kent/src/hg/makeDb/genbank
     # edit etc/genbank.conf to add this section just after canFam2:
 # ailMel1 (panda)
 ailMel1.serverGenome = /hive/data/genomes/ailMel1/ailMel1.2bit
 ailMel1.clusterGenome = /scratch/data/ailMel1/ailMel1.2bit
 ailMel1.lift = no
 ailMel1.ooc = /scratch/data/ailMel1/ailMel1.11.ooc
 ailMel1.refseq.mrna.native.pslCDnaFilter  = ${ordered.refseq.mrna.native.pslCDnaFilter}
 ailMel1.refseq.mrna.xeno.pslCDnaFilter    = ${ordered.refseq.mrna.xeno.pslCDnaFilter}
 ailMel1.genbank.mrna.native.pslCDnaFilter = ${ordered.genbank.mrna.native.pslCDnaFilter}
 ailMel1.genbank.mrna.xeno.pslCDnaFilter   = ${ordered.genbank.mrna.xeno.pslCDnaFilter}
 ailMel1.genbank.est.native.pslCDnaFilter  = ${ordered.genbank.est.native.pslCDnaFilter}
 ailMel1.refseq.mrna.native.load = yes
 ailMel1.refseq.mrna.xeno.load = yes
 ailMel1.genbank.mrna.xeno.load = yes
 ailMel1.downloadDir = ailMel1
 ailMel1.perChromTables = no
 
     # Edit src/lib/gbGenome.c to add new species.  With these four lines:
 # static char *ailMelNames[] = {"Ailuropoda melanoleuca", "Canis familiaris",
 #                                "Canis sp.", "Canis lupus familiaris",
 #                                "Canis lupus", NULL};
 #   ... later ...
 # {"ailMel", ailMelNames},
 
     cvs ci -m "adding ailMel - Panda with Dog specifications" src/lib/gbGenome.c
     make install-server
 
     ssh genbank
     screen  # control this business with a screen since it takes a while
     cd /cluster/data/genbank
     time nice -n +19 bin/gbAlignStep -initial  ailMel1 &
     #	real    652m27.871s
     #	var/build/logs/2010.04.09-15:37:11.ailMel1.initalign.log
 
     ssh hgwdev
     cd /cluster/data/genbank
     time ./bin/gbDbLoadStep -drop -initialLoad ailMel1 &
     #	logFile: var/dbload/hgwdev/logs/2010.04.12-10:13:40.dbload.log
     #	real    30m57.886s
 
     # enable daily alignment and update of hgwdev
     cd ~/kent/src/hg/makeDb/genbank
     cvsup
     # add oryCun2 to:
         etc/align.dbs
         etc/hgwdev.dbs
     cvs ci -m "Adding ailMel1 - Panda - Ailuropoda melanoleuca" \
 	etc/align.dbs etc/hgwdev.dbs
     make etc-update
     #	done - 2010-04-12 - Hiram
 
 ############################################################################
 # reset position to RHO location as found from blat of hg19 RHO gene
     hgsql -e \
 'update dbDb set defaultPos="GL192818.1:558576-566855" where name="ailMel1";' \
 	hgcentraltest
 
 ############################################################################
 # LASTZ Opossum monDom5 (DONE - 2010-02-05 - Hiram)
     #	original alignment
     cd /hive/data/genomes/monDom5/bed/lastzAilMel1.2010-02-04
     cat fb.monDom5.chainAilMel1Link.txt 
     #	223510659 bases of 3501660299 (6.383%) in intersection
 
     #	and for the swap
     mkdir /hive/data/genomes/ailMel1/bed/blastz.monDom5.swap
     cd /hive/data/genomes/ailMel1/bed/blastz.monDom5.swap
     time nice -n +19 doBlastzChainNet.pl -verbose=2 \
 	/hive/data/genomes/monDom5/bed/lastzAilMel1.2010-02-04/DEF \
 	-swap -noLoadChainSplit -syntenicNet \
 	-workhorse=hgwdev -smallClusterHub=memk -bigClusterHub=swarm \
 	-chainMinScore=5000 -chainLinearGap=loose > swap.log 2>&1 &
     #	real    69m35.464s
     cat fb.ailMel1.chainMonDom5Link.txt 
     #	211209682 bases of 2245312831 (9.407%) in intersection
 
 #########################################################################
 # LASTZ Mouse Mm9 (DONE - 2010-02-05 - Hiram)
     #	original alignment
     cd /hive/data/genomes/mm9/bed/lastzAilMel1.2010-02-04
     cat fb.mm9.chainAilMel1Link.txt 
     #	749595031 bases of 2620346127 (28.607%) in intersection
 
     #	and for the swap
     mkdir /hive/data/genomes/ailMel1/bed/blastz.mm9.swap
     cd /hive/data/genomes/ailMel1/bed/blastz.mm9.swap
     time doBlastzChainNet.pl -verbose=2 \
 	/hive/data/genomes/mm9/bed/lastzAilMel1.2010-02-04/DEF \
 	-swap -noLoadChainSplit -bigClusterHub=swarm -smallClusterHub=memk \
 	-workhorse=hgwdev \
 	-chainMinScore=3000 -chainLinearGap=medium > swap.log 2>&1 &
     #	real    54m57.140s
     cat fb.ailMel1.chainMm9Link.txt 
     #	739076250 bases of 2245312831 (32.916%) in intersection
 
 #######################################################################
 # felCatV17e Cat BLASTZ/CHAIN/NET (working  - 2010-03-22 - Chin)
     screen # use a screen to manage this multi-day job
     mkdir /hive/data/genomes/ailMel1/bed/lastzFelCatV17e.2010-03-22
     cd /hive/data/genomes/ailMel1/bed/lastzFelCatV17e.2010-03-22
 
     cat << '_EOF_' > DEF
 # human vs. cat
 # maximum M allowed with lastz is only 254
 BLASTZ_M=254
 
 # QUERY: Panda ailMel1
 SEQ2_DIR=/scratch/data/ailMel1/ailMel1.2bit
 SEQ2_LEN=/scratch/data/ailMel1/chrom.sizes
 SEQ2_CHUNK=10000000
 SEQ2_LAP=10000
 SEQ2_LIMIT=300
 
 # TARGET: Cat (felCatV17e)
 SEQ1_DIR=/scratch/data/felCatV17e/felCatV17e.2bit
 SEQ1_LEN=/scratch/data/felCatV17e/chrom.sizes
 SEQ1_LIMIT=50
 SEQ1_CHUNK=20000000
 SEQ1_LAP=10000
 
 BASE=/hive/data/genomes/ailMel1/bed/lastzFelCatV17e.2010-03-22
 TMPDIR=/scratch/tmp
 '_EOF_'
     # << this line keeps emacs coloring happy
 
     time nice -n +19 doBlastzChainNet.pl -verbose=2 \
         `pwd`/DEF \
         -syntenicNet -noDbNameCheck \
         -chainMinScore=3000 -chainLinearGap=medium \
         -workhorse=hgwdev -smallClusterHub=encodek -bigClusterHub=swarm \
         > do.log 2>&1 &
     #  Command failed:
     #     doClusterRun.csh
     # manually run the doClusterRun.csh  03/28
     time ./doClusterRun.csh > doClusterRun.log 2>&1 &
     # real    93m42.604s
     #  check with counting file number in psl directory to make sure that
     #  all jobs were completed.
     # manually start with "cat" step
     time nice -n +19 doBlastzChainNet.pl -verbose=2 \
         `pwd`/DEF \
         -syntenicNet -continue=cat -noDbNameCheck \
         -chainMinScore=3000 -chainLinearGap=medium \
         -workhorse=hgwdev -smallClusterHub=encodek -bigClusterHub=swarm \
         > cat.log 2>&1 &
     # real    346m8.431s
     cat fb.felCatV17e.chainAilMel1Link.txt
     # 1503647735 bases of 1990635005 (75.536%) in intersection
+    # run rbest 2010-04-30
+    time doRecipBest.pl -workhorse=hgwdev -buildDir=`pwd` \
+        felCatV17e ailMel1 > rbest.log 2>&1 &
+    #   real    234m49.556s
+
 
 
     mkdir /hive/data/genomes/felCatV17e/bed/blastz.ailMel1.swap
     cd /hive/data/genomes/felCatV17e/bed/blastz.ailMel1.swap
     time nice -n +19 doBlastzChainNet.pl -verbose=2 \
         /hive/data/genomes/ailMel1/bed/lastzFelCatV17e.2010-03-22/DEF \
         -swap -syntenicNet -noDbNameCheck \
         -workhorse=hgwdev -smallClusterHub=encodek -bigClusterHub=swarm \
         -chainMinScore=3000 -chainLinearGap=medium > swap.log 2>&1 &
     # real    155m57.566s
 
     # need to change name back due to the seq1 and seq2 switch
     # more TODO  XXXX 03/29
     # so should be able to do
     #	cat fb.felCatV17e.chainAilMel1Link.txt
     # without 
     # cd /hive/data/genomes/ailMel1/bed/blastz.felCatV17e.swap
     cat fb.ailMel1.chainFelCatV17eLink.txt
     # 1507273252 bases of 2245312831 (67.130%) in intersection
 
 #########################################################################
 # all.joiner update, downloads and in pushQ - (DONE - 2010-04-01 - Hiram)
     cd $HOME/kent/src/hg/makeDb/schema
     # fixup all.joiner until this is a clean output
     joinerCheck -database=ailMel1 -all all.joiner
 
     mkdir /hive/data/genomes/ailMel1/goldenPath
     cd /hive/data/genomes/ailMel1/goldenPath
     makeDownloads.pl ailMel1 > do.log 2>&1
 
     #	now ready for pushQ entry
     mkdir /hive/data/genomes/ailMel1/pushQ
     cd /hive/data/genomes/ailMel1/pushQ
     makePushQSql.pl ailMel1 > ailMel1.pushQ.sql 2> stderr.out
     #	check for errors in stderr.out, some are OK, e.g.:
 # WARNING: ailMel1 does not have seq
 # WARNING: ailMel1 does not have extFile
 
 # WARNING: Could not tell (from trackDb, all.joiner and hardcoded lists of
 # supporting and genbank tables) which tracks to assign these tables to:
 #  bosTau4ChainPileUp
 
     #	copy it to hgwbeta
     scp -p ailMel1.pushQ.sql hgwbeta:/tmp
     ssh hgwbeta
     cd /tmp
     hgsql qapushq < ailMel1.pushQ.sql
     #	in that pushQ entry walk through each entry and see if the
     #	sizes will set properly
 
 ############################################################################
 # Marking *all* gaps - they are not all in the AGP file
 #	(DONE - 2010-03-25 - Hiram)
     mkdir /hive/data/genomes/ailMel1/bed/allGaps
     cd /hive/data/genomes/ailMel1/bed/allGaps
     time nice -n +19 findMotif -motif=gattaca -verbose=4 \
 	-strand=+ ../../ailMel1.2bit > findMotif.txt 2>&1
     #	real    1m12.153s
     grep "^#GAP " findMotif.txt | sed -e "s/^#GAP //" > allGaps.bed
     featureBits ailMel1 -not gap -bed=notGap.bed
     #	2245312831 bases of 2245312831 (100.000%) in intersection
     featureBits ailMel1 allGaps.bed notGap.bed -bed=new.gaps.bed
     #	0 bases of 2245312831 (0.000%) in intersection
     #	real    84m58.382s
     #	no new gaps
 
 #########################################################################
 # lastz Human hg19 swap (DONE - 2010-02-05 - Hiram)
     # original alignment
     cd /hive/data/genomes/hg19/bed/lastzAilMel1.2010-02-04
     cat fb.hg19.chainAilMel1Link.txt 
     #	1453400264 bases of 2897316137 (50.164%) in intersection
 
     #	and the swap
     mkdir /hive/data/genomes/ailMel1/bed/blastz.hg19.swap
     cd /hive/data/genomes/ailMel1/bed/blastz.hg19.swap
     time nice -n +19 doBlastzChainNet.pl -verbose=2 \
 	/hive/data/genomes/hg19/bed/lastzAilMel1.2010-02-04/DEF \
 	-swap -noLoadChainSplit -syntenicNet \
 	-workhorse=hgwdev -smallClusterHub=memk -bigClusterHub=swarm \
 	-chainMinScore=3000 -chainLinearGap=medium > swap.log 2>&1 &
     real    124m14.393s
     cat fb.ailMel1.chainHg19Link.txt 
     #	1411953704 bases of 2245312831 (62.884%) in intersection
 
 #########################################################################
 # lastz Dog canFam2 swap (DONE - 2010-02-05 - Hiram)
     # original alignment
     cd /hive/data/genomes/canFam2/bed/lastzAilMel1.2010-02-04
     cat fb.canFam2.chainAilMel1Link.txt 
     #	1791212709 bases of 2384996543 (75.103%) in intersection
 
     #	and the swap
     mkdir /hive/data/genomes/ailMel1/bed/blastz.canFam2.swap
     cd /hive/data/genomes/ailMel1/bed/blastz.canFam2.swap
     time doBlastzChainNet.pl -verbose=2 \
 	/hive/data/genomes/canFam2/bed/lastzAilMel1.2010-02-04/DEF \
 	-swap -noLoadChainSplit -workhorse=hgwdev -bigClusterHub=pk \
 	-smallClusterHub=memk -chainMinScore=3000 -chainLinearGap=medium \
 	> swap.log 2>&1 &
     #	real    128m41.005s
      cat fb.ailMel1.chainCanFam2Link.txt 
     #	1788107935 bases of 2245312831 (79.637%) in intersection
 
 #########################################################################
 # lastz Rat rn4 swap (DONE - 2010-02-18 - Hiram)
     # original alignment
     cd /hive/data/genomes/rn4/bed/lastzAilMel1.2010-02-04
     cat fb.rn4.chainAilMel1Link.txt 
     #	708197812 bases of 2571531505 (27.540%) in intersection
 
     #	and the swap
     mkdir /hive/data/genomes/ailMel1/bed/blastz.rn4.swap
     cd /hive/data/genomes/ailMel1/bed/blastz.rn4.swap
     time nice -n +19 doBlastzChainNet.pl -verbose=2 \
 	/hive/data/genomes/rn4/bed/lastzAilMel1.2010-02-04/DEF \
 	-swap -syntenicNet \
 	-workhorse=hgwdev -smallClusterHub=memk -bigClusterHub=swarm \
 	-chainMinScore=3000 -chainLinearGap=medium > swap.log 2>&1 &
     #	real    73m37.751s
     cat fb.ailMel1.chainRn4Link.txt 
     #	695366144 bases of 2245312831 (30.970%) in intersection
 
 #########################################################################
 # HUMAN (hg18) PROTEINS TRACK (DONE braney 2010-04-19)
     # bash  if not using bash shell already
 
     cd /cluster/data/ailMel1
     mkdir /cluster/data/ailMel1/blastDb
 
     awk '{if ($2 > 1000000) print $1}' chrom.sizes > 1meg.lst
     twoBitToFa -seqList=1meg.lst  ailMel1.2bit temp.fa
     faSplit gap temp.fa 1000000 blastDb/x -lift=blastDb.lft
     rm temp.fa 1meg.lst
 
     awk '{if ($2 <= 1000000) print $1}' chrom.sizes > less1meg.lst
     twoBitToFa -seqList=less1meg.lst  ailMel1.2bit temp.fa
     faSplit about temp.fa 1000000 blastDb/y 
     rm temp.fa less1meg.lst
 
     cd blastDb
     for i in *.fa
     do
 	/hive/data/outside/blast229/formatdb -i $i -p F
     done
     rm *.fa
     ls *.nsq | wc -l
 # 2757
 
     mkdir -p /cluster/data/ailMel1/bed/tblastn.hg18KG
     cd /cluster/data/ailMel1/bed/tblastn.hg18KG
     echo  ../../blastDb/*.nsq | xargs ls -S | sed "s/\.nsq//"  > query.lst
     wc -l query.lst
 # 2757 query.lst
 
    # we want around 350000 jobs
    calc `wc /cluster/data/hg18/bed/blat.hg18KG/hg18KG.psl | awk '{print $1}'`/\(350000/`wc query.lst | awk '{print $1}'`\)
 
 # 36727/(350000/2757) = 289.303826
 
    mkdir -p kgfa
    split -l 289 /cluster/data/hg18/bed/blat.hg18KG/hg18KG.psl  kgfa/kg
    cd kgfa
    for i in *; do 
      nice pslxToFa $i $i.fa; 
      rm $i; 
    done
    cd ..
    ls -1S kgfa/*.fa > kg.lst
    wc kg.lst
 #  128  128 1664 kg.lst
 
    mkdir -p blastOut
    for i in `cat kg.lst`; do  mkdir blastOut/`basename $i .fa`; done
    tcsh
    cd /cluster/data/ailMel1/bed/tblastn.hg18KG
    cat << '_EOF_' > blastGsub
 #LOOP
 blastSome $(path1) {check in line $(path2)} {check out exists blastOut/$(root2)/q.$(root1).psl }
 #ENDLOOP
 '_EOF_'
 
    cat << '_EOF_' > blastSome
 #!/bin/sh
 BLASTMAT=/hive/data/outside/blast229/data
 export BLASTMAT
 g=`basename $2`
 f=/tmp/`basename $3`.$g
 for eVal in 0.01 0.001 0.0001 0.00001 0.000001 1E-09 1E-11
 do
 if /hive/data/outside/blast229/blastall -M BLOSUM80 -m 0 -F no -e $eVal -p tblastn -d $1 -i $2 -o $f.8
 then
         mv $f.8 $f.1
         break;
 fi
 done
 if test -f  $f.1
 then
     if /cluster/bin/i386/blastToPsl $f.1 $f.2
     then
 	liftUp -nosort -type=".psl" -nohead $f.3 /cluster/data/ailMel1/blastDb.lft carry $f.2
         liftUp -nosort -type=".psl" -pslQ -nohead $3.tmp /cluster/data/hg18/bed/blat.hg18KG/protein.lft warn $f.3
         if pslCheck -prot $3.tmp
         then                  
             mv $3.tmp $3     
             rm -f $f.1 $f.2 $f.3 $f.4
         fi
         exit 0               
     fi                      
 fi                         
 rm -f $f.1 $f.2 $3.tmp $f.8 $f.3 $f.4
 exit 1
 '_EOF_'
     # << happy emacs
     chmod +x blastSome
     exit 
     
     ssh swarm
     cd /cluster/data/ailMel1/bed/tblastn.hg18KG
     gensub2 query.lst kg.lst blastGsub blastSpec
     para create blastSpec
 #    para try, check, push, check etc.
 
     para time
 # Completed: 352896 of 352896 jobs
 # CPU time in finished jobs:   13166354s  219439.24m  3657.32h  152.39d  0.418 y
 # IO & Wait Time:               1995282s   33254.69m   554.24h   23.09d  0.063 y
 # Average job time:                  43s       0.72m     0.01h    0.00d
 # Longest finished job:             207s       3.45m     0.06h    0.00d
 # Submission to last job:         15726s     262.10m     4.37h    0.18d
 
     ssh swarm
     cd /cluster/data/ailMel1/bed/tblastn.hg18KG
     mkdir chainRun
     cd chainRun
     tcsh
     cat << '_EOF_' > chainGsub
 #LOOP
 chainOne $(path1)
 #ENDLOOP
 '_EOF_'
 
     cat << '_EOF_' > chainOne
 (cd $1; cat q.*.psl | simpleChain -prot -outPsl -maxGap=150000 stdin ../c.`basename $1`.psl)
 '_EOF_'
     chmod +x chainOne
     ls -1dS ../blastOut/kg?? > chain.lst
     gensub2 chain.lst single chainGsub chainSpec
     # do the cluster run for chaining
     para create chainSpec
     para try, check, push, check etc.
 
 # Completed: 128 of 128 jobs
 #CPU time in finished jobs:       1401s      23.35m     0.39h    0.02d  0.000 y
 #IO & Wait Time:                 47407s     790.12m    13.17h    0.55d  0.002 y
 #Average job time:                 381s       6.36m     0.11h    0.00d
 #Longest finished job:             466s       7.77m     0.13h    0.01d
 #Submission to last job:           470s       7.83m     0.13h    0.01d
 
     cd /cluster/data/ailMel1/bed/tblastn.hg18KG/blastOut
     for i in kg??
     do
        cat c.$i.psl | awk "(\$13 - \$12)/\$11 > 0.6 {print}" > c60.$i.psl
        sort -rn c60.$i.psl | pslUniq stdin u.$i.psl
        awk "((\$1 / \$11) ) > 0.60 { print   }" c60.$i.psl > m60.$i.psl
        echo $i
     done
     sort u.*.psl m60* | uniq | sort -T /tmp -k 14,14 -k 16,16n -k 17,17n > ../blastHg18KG.psl
     cd ..
     pslCheck blastHg18KG.psl
 # checked: 64037 failed: 0 errors: 0
 
     # load table 
     ssh hgwdev
     cd /cluster/data/ailMel1/bed/tblastn.hg18KG
     hgLoadPsl ailMel1 blastHg18KG.psl
 
     # check coverage
     featureBits ailMel1 blastHg18KG 
 # 32308509 bases of 2245312831 (1.439%) in intersection
 
     featureBits ailMel1 blastHg18KG xenoRefGene  -enrichment
 # blastHg18KG 1.439%, xenoRefGene 2.186%, both 1.137%, cover 79.05%, enrich  36.17x
     rm -rf blastOut
 #end tblastn