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

#	Caenorhabditis remanei
#	From Wash U GSC
#	http://genome.wustl.edu/genome.cgi?GENOME=Caenorhabditis+remanei
#	http://www.ncbi.nlm.nih.gov/Traces/wgs/?val=AAGD02

#  $Id: caeRem3.txt,v 1.6 2008/08/15 21:53:32 hiram Exp $

#########################################################################
# DOWNLOAD SEQUENCE (DONE - 2008-05-22 - Hiram)
    ssh kkstore01
    mkdir /cluster/store10/caeRem3
    ln -s /cluster/store10/caeRem3 /cluster/data/caeRem3

    mkdir /cluster/data/caeRem3/downloads
    cd /cluster/data/caeRem3/downloads
    wget --timestamping \
ftp://genome.wustl.edu/pub/organism/Invertebrates/Caenorhabditis_remanei/assembly/Caenorhabditis_remanei-15.0.1/ASSEMBLY

    wget --timestamping \
"ftp://genome.wustl.edu/pub/organism/Invertebrates/Caenorhabditis_remanei/assembly/Caenorhabditis_remanei-15.0.1/output/README" \
	-O README.caeRem3

    wget --timestamping -m -np -nd \
"ftp://genome.wustl.edu/pub/organism/Invertebrates/Caenorhabditis_remanei/assembly/Caenorhabditis_remanei-15.0.1/output/"

##########################################################################
## Create chrUn.agp file (DONE - 2008-05-22 - Hiram)
    ssh kkstore01
    cd /cluster/data/caeRem3/downloads
    faSize contigs.fa.gz
# 138406388 bases (185 N's 138406203 real 138406203 upper 0 lower)
#	in 12680 sequences in 1 files
    faSize supercontigs.fa.gz
# 145442736 bases (7036533 N's 138406203 real 138406203 upper 0 lower)
#	in 3670 sequences in 1 files

    #	note the chrUnSize in the following, calculate from the faSize
    #	on the supercontigs:
    # 145442736 + ( (3670 - 1) * 1000 ) = 149111736
    cat << '_EOF_' > mkSuperLift.pl
#!/usr/bin/env perl

use strict;
use warnings;

my $start = 1;
my $end = 1;
my $itemCount = 1;
my $agpItemCount = 1;
my $curContig = "";
my $firstTime = 1;
my $scaffoldGapSize = 1000;
my $scafEnd = 1;
my $superStart = 0;
my $superEnd = 1;
my $superLength = 0;
my $chrUnSize = 149111736;
my $chrUnName = "chrUn";
my $ctgFragCount = 1;

open (CT,">caeRem3.ctgPos2.tab") or die "Can not write to caeRem3.ctgPos2.tab";

open (SC,">caeRem3.CtgScaf.agp") or die "Can not write to caeRem3.CtgScaf.agp";

open (AG,">caeRem3.chrUn.agp") or die "Can not write to caeRem3.chrUn.agp";

open (GL,">caeRem3.gold.tab") or die "Can not write to caeRem3.gold.tab";

open (FH,'zcat supercontigs.agp.gz|') or
	die "can not open zcat supercontigs.agp.gz";
while (my $line=<FH>) {
    chomp $line;
    if ($line =~ m/fragment/) {
	my ($name, $gStart, $gEnd, $gapCounter, $gapWN, $gapLen, $frag, $yesNo) =
	    split('\s+',$line);
	$end = $start + $gapLen - 1;
	$scafEnd += $gapLen;
	printf SC "chrUn\t%d\t%d\t%d\t%s\t%d\t%s\t%s\n",
		    $start, $end, $agpItemCount, $gapWN, $gapLen, $frag, $yesNo;
	printf AG "chrUn\t%d\t%d\t%d\t%s\t%d\t%s\t%s\n",
		    $start, $end, $agpItemCount++, $gapWN, $gapLen, $frag, $yesNo;
	$start = $end + 1;
    } else {
    my ($ctgName, $ctgStart, $ctgEnd, $ctgCounter, $ctgWN, $name, $cStart, $cEnd, $strand) =
	split('\s+',$line);
	my $ctgLen = $ctgEnd - $ctgStart + 1;
	my $cLen = $cEnd - $cStart + 1;
	die "not matching start, end:\n$line" if ($ctgLen != $cLen);
	if (!$firstTime) {
	    if ($ctgName ne $curContig) {
		$superLength = $superEnd - $superStart;
		$end = $start + $scaffoldGapSize - 1;
		printf SC "chrUn\t%d\t%d\t%d\tN\t%d\tscaffold\tno\n",
		    $start, $end, $agpItemCount, $scaffoldGapSize;
		printf AG "chrUn\t%d\t%d\t%d\tN\t%d\tscaffold\tno\n",
		    $start, $end, $agpItemCount++, $scaffoldGapSize;
		printf "%d\t%s\t%d\t%s\t%d\n",
			$superStart, $curContig, $superLength, $chrUnName,
				$chrUnSize;
		printf GL "%s\t%d\t%d\t%d\tW\t%s\t0\t%d\t+\n",
			$chrUnName, $superStart, $superEnd, $itemCount++,
			$curContig, $superLength;
		$start = $end + 1;
		$scafEnd = $cStart - 1;
		$superStart = $start - 1;
		$ctgFragCount = 1;
	    }
	} else {
	    $firstTime = 0;
	    $scafEnd = 0;
	}
	$scafEnd += $ctgLen;
	my $fragStart = $scafEnd - $ctgLen + 1;
	$end = $start + $ctgLen - 1;
	$superEnd = $end;
	my $ctgFragName = sprintf("%s.%d", $ctgName, $ctgFragCount++);
	my $ctgFragStart = 1;
	my $ctgFragEnd = $scafEnd - $fragStart + 1;
	printf SC "chrUn\t%d\t%d\t%d\t%s\t%s\t%d\t%d\t%s\n",
		$start, $end, $agpItemCount, $ctgWN, $ctgFragName,
		$ctgFragStart, $ctgFragEnd, $strand;
	printf AG "chrUn\t%d\t%d\t%d\t%s\t%s\t%d\t%d\t%s\n",
		$start, $end, $agpItemCount++, $ctgWN, $ctgName,
		$fragStart, $scafEnd, $strand;
	printf CT "%s\t%d\tchrUn\t%d\t%d\tW\n",
                $name, $ctgLen, $start-1, $end;

	$start = $end + 1;
	$curContig = $ctgName;
    }
}
close (FH);
close (AG);
close (SC);
$superLength = $superEnd - $superStart;
printf "%d\t%s\t%d\t%s\t%d\n",
    $superStart, $curContig, $superLength, $chrUnName, $chrUnSize;
printf GL "%s\t%d\t%d\t%d\tW\t%s\t0\t%d\t+\n",
	$chrUnName, $superStart, $superEnd, $itemCount++,
	$curContig, $superLength;
close (GL);
'_EOF_'
    # << happy emacs
    chmod +x mkSuperLift.pl
    ./mkSuperLift.pl > caeRem3.chrUn.lift

    qaToQac contigs.fa.qual.gz stdout \
        | qacAgpLift caeRem3.CtgScaf.agp stdin chrUn.qac

###########################################################################
## Initial genome build (DONE - 2008-05-22 - Hiram)
    ssh kkstore01
    cd /cluster/data/caeRem3
    cat << '_EOF_' > caeRem3.config.ra
# Config parameters for makeGenomeDb.pl:
db caeRem3
clade worm
genomeCladePriority 10
scientificName Caenorhabditis remanei
commonName C. remanei
assemblyDate Jul. 2007
assemblyLabel Washington University School of Medicine GSC C. remanei 15.0.1
orderKey 878
mitoAcc none
fastaFiles /cluster/data/caeRem3/downloads/supercontigs.fa.gz
agpFiles /cluster/data/caeRem3/downloads/caeRem3.chrUn.agp
qualFiles /cluster/data/caeRem3/downloads/chrUn.qac
dbDbSpeciesDir worm
'_EOF_'
    # << happy emacs

    time nice -n +19 makeGenomeDb.pl caeRem3.config.ra \
	> makeGenomeDb.out 2>&1
    #	real    4m14.263s

    ssh hgwdev
    cd /cluster/data/caeRem3/downloads
    sed -e "s/agpFrag/chrUn_gold/" $HOME/kent/src/hg/lib/agpFrag.sql \
	> chrUn_gold.sql
    # edit that .sql file to add the bin column
    #	 bin smallint(6) NOT NULL default '0',
    #	or, copy the one from caeRem2:
    cp /cluster/data/caeRem2/downloads/chrUn_gold.sql .

    hgLoadBed -sqlTable=chrUn_gold.sql caeRem3 chrUn_gold caeRem3.gold.tab
    hgLoadSqlTab caeRem3 ctgPos2 ~/kent/src/hg/lib/ctgPos2.sql \
	caeRem3.ctgPos2.tab

###########################################################################
## Repeat Masker (DONE - 2008-05-22 - Hiram)
    ssh kkstore01
    mkdir /cluster/data/caeRem3/bed/repeatMasker
    cd /cluster/data/caeRem3/bed/repeatMasker
    time nice -n +19 doRepeatMasker.pl caeRem3 \
	-buildDir=/cluster/data/caeRem3/bed/repeatMasker > do.log 2>&1
    #	real    15m16.812s
    cd /cluster/data/caeRem3
    twoBitToFa caeRem3.rmsk.2bit stdout | faSize stdin \
	> caeRem3.rmsk.2bit.faSize
    cat caeRem3.rmsk.2bit.faSize
    #	149111736 bases (10705533 N's 138406203 real
    #	137196473 upper 1209730 lower) in 1 sequences in 1 files
    #	%0.81 masked total, %0.87 masked real

    #	we are not going to use this sequence, too little is masked,
    #	instead, using the windowMasker business next

    #	ln -s caeRem3.rmsk.2bit caeRem3.2bit

    #	mkdir /san/sanvol1/scratch/worms/caeRem3 
    #	cp -p caeRem3.2bit /san/sanvol1/scratch/worms/caeRem3
    #	cp -p chrom.sizes /san/sanvol1/scratch/worms/caeRem3
    #	cp -p downloads/caeRem3.chrUn.lift /san/sanvol1/scratch/worms/caeRem3

###########################################################################
## Window Masker (DONE - 2008-05-22 - Hiram)
    ssh kolossus
    mkdir /cluster/data/caeRem3/bed/windowMasker
    cd /cluster/data/caeRem3/bed/windowMasker
    time nice -n +19 ~/kent/src/hg/utils/automation/doWindowMasker.pl \
	-workhorse kolossus \
	-buildDir=/cluster/data/caeRem3/bed/windowMasker caeRem3 > do.log 2>&1
    #	real    10m40.787s

    twoBitToFa caeRem3.wmsk.sdust.2bit stdout \
	| faSize stdin > caeRem3.wmsk.sdust.2bit.faSize 2>&1
    cat caeRem3.wmsk.sdust.2bit.faSize
    #	149111736 bases (10705533 N's 138406203 real
    #	96266616 upper 42139587 lower) in 1 sequences in 1 files
    #	%28.26 masked total, %30.45 masked real

    ssh hgwdev
    cd /cluster/data/caeRem3/bed/windowMasker
    hgLoadBed caeRem3 windowmaskerSdust windowmasker.sdust.bed.gz
    #	Loaded 1122372 elements of size 3

    #	eliminate the gaps from the masking (WM bug)
    featureBits caeRem3 -not gap -bed=notGap.bed
    #	138406388 bases of 138406388 (100.000%) in intersection
 
    time nice -n +19 featureBits caeRem3 windowmaskerSdust notGap.bed \
        -bed=stdout | gzip -c > cleanWMask.bed.gz
    #	42139616 bases of 138406388 (30.446%) in intersection

    #	reload track to get it clean
    hgLoadBed caeRem3 windowmaskerSdust cleanWMask.bed.gz
    #	Loaded 1075009 elements of size 4
    featureBits caeRem3 windowmaskerSdust
    #	42139616 bases of 138406388 (30.446%) in intersection

###########################################################################
## Simple Repeats (DONE - 2008-05-22 - Hiram)
    ssh titan
    screen # use a screen to manage this job
    mkdir /cluster/data/caeRem3/bed/simpleRepeat
    cd /cluster/data/caeRem3/bed/simpleRepeat
    time twoBitToFa ../../caeRem3.unmasked.2bit stdout \
	| nice -n +19 trfBig -trf=/cluster/bin/i386/trf stdin /dev/null \
	-bedAt=simpleRepeat.bed -tempDir=/scratch/tmp > do.out 2>&1
    #	real    43m39.381s
    awk '$5 <= 12' simpleRepeat.bed > trfMask.bed

    ssh hgwdev
    cd /cluster/data/caeRem3/bed/simpleRepeat
    nice -n +19 hgLoadBed caeRem3 simpleRepeat \
      simpleRepeat.bed -sqlTable=$HOME/kent/src/hg/lib/simpleRepeat.sql
    #	Loaded 39971 elements of size 16
    featureBits caeRem3 simpleRepeat > fb.caeRem3.simpleRepeat.txt 2>&1
    cat fb.caeRem3.simpleRepeat.txt 
    #	5219673 bases of 138406388 (3.771%) in intersection

##############################################################################
## combine trf mask with windowmasker (DONE - 2008-05-23 - Hiram)
    ssh kkstore01
    cd /cluster/data/caeRem3
    zcat bed/windowMasker/cleanWMask.bed.gz \
	| twoBitMask caeRem3.unmasked.2bit stdin \
	    -type=.bed caeRem3.cleanWMSdust.2bit
    cat bed/simpleRepeat/trfMask.bed \
      | twoBitMask -add -type=.bed caeRem3.cleanWMSdust.2bit stdin caeRem3.2bit
    #	can safely ignore the warning about BED file >= 13 fields
    #	check it
    twoBitToFa caeRem3.2bit stdout | faSize stdin
    #	149111736 bases (10705533 N's 138406203 real 96215522 upper
    #	42190681 lower) in 1 sequences in 1 files
    #	%28.29 masked total, %30.48 masked real

    #	create gbdb symlink
    ssh hgwdev
    ln -s /cluster/data/caeRem3/caeRem3.2bit /gbdb/caeRem3

##############################################################################
## create masked contigs for blastz (DONE - 2008-05-23 - Hiram)
    ssh kkstore01
    mkdir /cluster/data/caeRem3/maskedContigs
    cd /cluster/data/caeRem3/maskedContigs
    #	copy over lft2BitToFa.pl from /cluster/data/caePb2/jkStuff
    cp -p /cluster/data/caePb2/jkStuff/lft2BitToFa.pl ../jkStuff
    ln -s ../downloads/caeRem3.chrUn.lift .

    time nice -n +19 ../jkStuff/lft2BitToFa.pl \
	../caeRem3.2bit caeRem3.chrUn.lift > super.fa
    #	real    1m40.247s
    #	verify it is OK, should have same numbers of real, upper and lower as
    #	measured in the caeRem3.2bit, only N's are different:
    faSize super.fa
    #	145442736 bases (7036533 N's 138406203 real 96215522 upper
    #	42190681 lower) in 3670 sequences in 1 files
    #	%29.01 masked total, %30.48 masked real
    #	in fact, the missing N's should be:
    # (3670 - 1) * 1000  = 149111736 - 145442736 = 10705533 - 7036533 =
    #	3669000


    faToTwoBit super.fa caeRem3.supercontigs.2bit
    twoBitInfo caeRem3.supercontigs.2bit stdout | sort -k2nr \
        > caeRem3.supercontigs.sizes
    #	Copy this 2bit file and sizes over to /san/sanvol1/scratch/worms/caeRem3

    #  the lft2BitToFa.pl script is:
    cat << '_EOF_'
#!/usr/bin/env perl

use strict;
use warnings;
use File::Basename;
    
sub usage()
{
printf "usage: %s <file.2bit> <file.lft> [more_files.lft]\n",  basename($0);
printf "\tfasta output is to stdout, therefore route stdout to result file\n";
exit 255;
}   
    
my $argc = scalar(@ARGV);

usage if ($argc < 2);

my $twoBitFile = shift;
    
while (my $liftFile = shift)
{   
open (FH,"<$liftFile") or die "Can not open $liftFile";
while (my $line=<FH>)
    {
    chomp $line;
    my ($start, $contig, $length, $chrom, $chrom_length) = split('\s',$line);
    my $cmd=sprintf("twoBitToFa $twoBitFile:%s:%d-%d stdout",
        $chrom, $start, $start+$length);
    print `$cmd | sed -e "s#^>.*#>$contig#"`;
    }
close (FH);
}
'_EOF_'
# << happy emacs

#########################################################################
# MAKE 11.OOC FILE FOR BLAT/GENBANK (DONE - 2008-05-23 - Hiram)
    # Use -repMatch=50 (based on size -- for human we use 1024, and 
    # C. remanei size is ~5.1% of human judging by gapless caeRem2 vs. hg18 
    # genome sizes from featureBits.  So we would use 50, but that yields a very
    # high number of tiles to ignore, especially for a small more compact 
    # genome.  Bump that up a bit to be more conservative.
    #	This more conservative setting was done June 4th.

    ssh kkstore01
    cd /cluster/data/caeRem3
    blat caeRem3.2bit /dev/null /dev/null -tileSize=11 \
      -makeOoc=jkStuff/11.ooc -repMatch=100
    #	Wrote 10173 overused 11-mers to 11.ooc
    cp -p jkStuff/11.ooc /san/sanvol1/scratch/worms/caeRem3
    # the san directory should have:
# -rw-rw-r--  1       16 May 13 09:52 chrom.sizes
# -rw-rw-r--  1 59427484 May 13 11:21 caePb2.2bit
# -rw-rw-r--  1 58631157 May 13 11:31 caePb2.supercontigs.2bit
# -rw-rw-r--  1    56860 May 13 11:31 caePb2.supercontigs.sizes
# -rw-rw-r--  1    40700 Jun  4 11:42 11.ooc
    #	copy those files to /cluster/bluearc/scratch/data/caeRem3/
    #	and to the Iservers /iscratch/i/worms/caeRem3/
    #	request an rsync of /scratch/data/ from the cluster admins

#########################################################################
# GENBANK AUTO UPDATE (DONE - 2008-05-23 - Hiram)
    # align with latest genbank process.
    ssh hgwdev
    cd ~/kent/src/hg/makeDb/genbank
    cvsup
    # edit etc/genbank.conf to add caePb2 just before caeRem2

# caeRem3 (C. remanei)
caeRem3.serverGenome = /cluster/data/caeRem3/caeRem3.2bit
caeRem3.clusterGenome = /iscratch/i/worms/caeRem3/caeRem3.2bit
caeRem3.ooc = /iscratch/i/worms/caeRem3/11.ooc
caeRem3.lift = /iscratch/i/worms/caeRem3/caeRem3.chrUn.lift
caeRem3.align.unplacedChroms = chrUn
caeRem3.refseq.mrna.native.pslCDnaFilter  = ${lowCover.refseq.mrna.native.pslCDnaFilter}
caeRem3.refseq.mrna.xeno.pslCDnaFilter    = ${lowCover.refseq.mrna.xeno.pslCDnaFilter}
caeRem3.genbank.mrna.native.pslCDnaFilter = ${lowCover.genbank.mrna.native.pslCDnaFilter}
caeRem3.genbank.mrna.xeno.pslCDnaFilter   = ${lowCover.genbank.mrna.xeno.pslCDnaFilter}
caeRem3.genbank.est.native.pslCDnaFilter  = ${lowCover.genbank.est.native.pslCDnaFilter}
caeRem3.refseq.mrna.native.load = no
caeRem3.refseq.mrna.xeno.load  = yes
caeRem3.refseq.mrna.xeno.loadDesc = yes
caeRem3.genbank.mrna.xeno.load = yes
caeRem3.genbank.est.native.load = yes
caeRem3.genbank.est.native.loadDesc = no
caeRem3.downloadDir = caeRem3
caeRem3.perChromTables = no

    cvs ci -m "Added caeRem3." etc/genbank.conf
    # update /cluster/data/genbank/:
    make etc-update

    ssh genbank
    cd /cluster/data/genbank
    time nice -n +19 bin/gbAlignStep -initial caeRem3 &
    #	var/build/logs/2008.05.23-10:34:13.caeRem3.initalign.log
    #	real    133m52.226s

    # load database when finished
    ssh hgwdev
    cd /cluster/data/genbank
    time nice -n +19 ./bin/gbDbLoadStep -drop -initialLoad caeRem3
    #	logFile: var/dbload/hgwdev/logs/2008.05.23-14:22:08.dbload.log
    #	real    17m45.023s

    # enable daily alignment and update of hgwdev
    cd ~/kent/src/hg/makeDb/genbank
    cvsup
    # add caeRem3 to:
        etc/align.dbs
        etc/hgwdev.dbs
    cvs ci -m "Added caeRem3 - C. remanei" etc/align.dbs etc/hgwdev.dbs
    make etc-update

############################################################################
#  BLATSERVERS ENTRY (DONE - 2008-05-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 ("caeRem3", "blat11", "17788", "1", "0"); \
	INSERT INTO blatServers (db, host, port, isTrans, canPcr) \
	VALUES ("caeRem3", "blat11", "17789", "0", "1");' \
	    hgcentraltest
    #	test it with some sequence

############################################################################
## Default position (DONE - 2008-05-23 - Hiram)
    ssh hgwdev
    hgsql -e 'update dbDb set defaultPos="chrUn:119043304-119068959"
	where name="caeRem3";' hgcentraltest

##########################################################################
## Creating downloads (DONE - 2008-05-23 - Hiram)
    #	There is only one chrom, make its trfMaskChrom file exist
    ssh hgwdev
    mkdir /cluster/data/caeRem3/bed/simpleRepeat/trfMaskChrom
    cd /cluster/data/caeRem3/bed/simpleRepeat/trfMaskChrom
    ## symlink didn't work here, the symlink ended up in the tar file
    cp -p ../trfMask.bed ./chrUn.bed
    cd  /cluster/data/caeRem3
    #	there is no RM file
    time nice -n +19 /cluster/bin/scripts/makeDownloads.pl caeRem3 
    #	real    2m10.175s
    ## *!* EDIT THE README.txt FILES *!* in goldenPath/*/README.txt ##

    cp -p bed/windowMasker/cleanWMask.bed.gz \
	goldenPath/bigZips/chromWMSdust.bed.gz

##########################################################################
## Creating pushQ (DONE - 2008-05-23 - Hiram)
    ssh hgwdev
    mkdir /cluster/data/caeRem3/pushQ
    cd /cluster/data/caeRem3/pushQ
    /cluster/bin/scripts/makePushQSql.pl caeRem3 > caeRem3.sql 2> stderr.out

##########################################################################
#  Fixup assembly date after confirmation from John Spieth
#	2008-05-27 - DONE - Hiram
    ssh hgwdev
    hgsql -e 'update dbDb set description="May 2007" where name="caeRem3"' \
	hgcentraltest

#########################################################################
# ELEGANS (ce6) PROTEINS TRACK (DONE -  Hiram - 2008-06-11)
    ssh kkstore01

    #	breaking up this target genome into manageable pieces
    #	using the contigs so we don't have genes spanning non-bridged gaps
    mkdir /cluster/data/caeRem3/blastDb
    cd /cluster/data/caeRem3
    twoBitToFa -noMask maskedContigs/caeRem3.supercontigs.2bit temp.fa
    faSplit gap temp.fa 1000000 blastDb/x -lift=blastDb.lft
    #	3709 pieces of 3709 written
    rm temp.fa
    cd blastDb
    for i in *.fa
    do
	/cluster/bluearc/blast229/formatdb -i $i -p F
    done
    rm *.fa

    ##	copy to san for kluster access
    mkdir -p /san/sanvol1/scratch/worms/caeRem3/blastDb
    cd /san/sanvol1/scratch/worms/caeRem3/blastDb
    rsync -a --progress --stats /cluster/data/caeRem3/blastDb/ .

    ## create the query protein set
    mkdir -p /cluster/data/caeRem3/bed/tblastn.ce6SG
    cd /cluster/data/caeRem3/bed/tblastn.ce6SG
    echo  /san/sanvol1/scratch/worms/caeRem3/blastDb/*.nsq | xargs ls -S \
	| sed "s/\.nsq//"  > query.lst
    wc -l query.lst
    # 3709 query.lst

   # we want around 100000 jobs
   calc `wc /cluster/data/ce6/bed/blat.ce6SG/ce6SG.psl | awk "{print \\\$1}"`/\(100000/`wc query.lst | awk "{print \\\$1}"`\)

    #	23741/(100000/3709) = 880.553690

   mkdir -p /cluster/bluearc/worms/caeRem3/bed/tblastn.ce6SG/sgfa
   split -l 880 /cluster/data/ce6/bed/blat.ce6SG/ce6SG.psl \
	/cluster/bluearc/worms/caeRem3/bed/tblastn.ce6SG/sgfa/sg
   ln -s /cluster/bluearc/worms/caeRem3/bed/tblastn.ce6SG/sgfa sgfa
   cd sgfa
   for i in *; do 
     nice pslxToFa $i $i.fa; 
     rm $i; 
   done
   cd ..
   ls -1S sgfa/*.fa > sg.lst
   mkdir -p /cluster/bluearc/worms/caeRem3/bed/tblastn.ce6SG/blastOut
   ln -s /cluster/bluearc/worms/caeRem3/bed/tblastn.ce6SG/blastOut
   for i in `cat sg.lst`; do  mkdir blastOut/`basename $i .fa`; done
   
   cd /cluster/data/caeRem3/bed/tblastn.ce6SG
   cat << '_EOF_' > template
#LOOP
blastSome $(path1) {check in line $(path2)} {check out exists blastOut/$(root2)/q.$(root1).psl }
#ENDLOOP
'_EOF_'
    # << happy emacs

   cat << '_EOF_' > blastSome
#!/bin/sh
CBIN=/cluster/bin/i386
DB=caeRem3
BLASTMAT=/scratch/data/blast229/data
SCR="/scratch/tmp/${DB}"
g=`/bin/basename $2`
D=`/bin/basename $1`
R=`/bin/basename $3`
inputFile=${SCR}/${D}_${g}
export CBIN DB BLASTMAT SCR g D R
/bin/mkdir -p ${SCR}
/bin/cp -p $1.* ${SCR}
/bin/cp -p $2 ${inputFile}
/bin/cp -p /cluster/store10/${DB}/blastDb.lft ${SCR}/${D}_blastDb.lft
/bin/cp -p /cluster/store3/worm/ce6/bed/blat.ce6SG/protein.lft \
	${SCR}/${D}_protein.lft
f=${SCR}/${R}.$g
for eVal in 0.01 0.001 0.0001 0.00001 0.000001 1E-09 1E-11
do
if /scratch/data/blast229/blastall -M BLOSUM80 -m 0 -F no \
        -e $eVal -p tblastn -d ${SCR}/${D} -i ${inputFile} -o $f.8
then
        /bin/mv $f.8 $f.1
        break;
fi
done
if test -f  $f.1
then
    if ${CBIN}/blastToPsl $f.1 $f.2
    then
        ${CBIN}/liftUp -nosort -type=".psl" -nohead $f.3 \
            ${SCR}/${D}_blastDb.lft carry $f.2 > /dev/null
        ${CBIN}/liftUp -nosort -type=".psl" -pslQ -nohead ${SCR}/${R}.tmp \
            ${SCR}/${D}_protein.lft warn $f.3 > /dev/null
        if ${CBIN}/pslCheck -prot ${SCR}/${R}.tmp
        then
            /bin/rm -f $3
            /bin/mv ${SCR}/${R}.tmp $3
            /bin/rm -f $f.1 $f.2 $f.3 $f.4 ${SCR}/${D}.* ${inputFile} \
		${SCR}/${D}_protein.lft ${SCR}/${D}_blastDb.lft
            /bin/rmdir --ignore-fail-on-non-empty ${SCR}
        fi
        exit 0
    fi
fi
/bin/rm -f $f.1 $f.2 $3.tmp $f.8 $f.3 $f.4 ${SCR}/${D}.* ${inputFile} \
	${SCR}/${D}_protein.lft ${SCR}/${D}_blastDb.lft
/bin/rmdir --ignore-fail-on-non-empty ${SCR}
exit 1
'_EOF_'
    # << happy emacs
    chmod +x blastSome

    ssh kk
    cd /cluster/data/caeRem3/bed/tblastn.ce6SG
    gensub2 query.lst sg.lst template jobList
    para create jobList
#    para try, check, push, check etc.
# Completed: 100142 of 100143 jobs
# Other count: 1 jobs
# CPU time in finished jobs:    9711211s  161853.51m  2697.56h  112.40d  0.308 y
# IO & Wait Time:               1135276s   18921.27m   315.35h   13.14d  0.036 y
# Average job time:                 108s       1.81m     0.03h    0.00d
# Longest finished job:             706s      11.77m     0.20h    0.01d
# Submission to last job:         32977s     549.62m     9.16h    0.38d
    #	the 1 stray job was a tracking error.  It did get finished,
    #	there was a result file and para recover said nothing was left to do

    # do the cluster run for chaining
    ssh kk
    mkdir /cluster/data/caeRem3/bed/tblastn.ce6SG/chainRun
    cd /cluster/data/caeRem3/bed/tblastn.ce6SG/chainRun
    cat << '_EOF_' > template
#LOOP
chainOne $(path1)
#ENDLOOP
'_EOF_'
    # << happy emacs

    cat << '_EOF_' > chainOne
(cd $1; cat q.*.psl | simpleChain -prot -outPsl -maxGap=50000 stdin /cluster/bluearc/worms/caeRem3/bed/tblastn.ce6SG/blastOut/c.`basename $1`.psl)
'_EOF_'
    # << happy emacs
    chmod +x chainOne

    ls -1dS /cluster/bluearc/worms/caeRem3/bed/tblastn.ce6SG/blastOut/sg?? \
	> chain.lst
    gensub2 chain.lst single template jobList
    cd /cluster/data/caeRem3/bed/tblastn.ce6SG/chainRun
    para create jobList
    # para maxNode 30 # there are only 27 jobs
    para try, check, push, check etc.
# Completed: 27 of 27 jobs
# CPU time in finished jobs:        394s       6.57m     0.11h    0.00d  0.000 y
# IO & Wait Time:                  2095s      34.91m     0.58h    0.02d  0.000 y
# Average job time:                  92s       1.54m     0.03h    0.00d
# Longest finished job:             110s       1.83m     0.03h    0.00d
# Submission to last job:           183s       3.05m     0.05h    0.00d

    ssh kkstore01
    cd /cluster/data/caeRem3/bed/tblastn.ce6SG/blastOut
    for i in sg??
    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 -T /scratch/tmp -k 14,14 -k 16,16n -k 17,17n u.*.psl m60* | uniq \
	> /cluster/data/caeRem3/bed/tblastn.ce6SG/blastCe6SG.psl
    cd ..
    pslCheck blastCe6SG.psl
    #	checked: 29004 failed: 0 errors: 0

    # load table 
    ssh hgwdev
    cd /cluster/data/caeRem3/bed/tblastn.ce6SG
    liftUp -nohead -type=.psl stdout ../../downloads/caeRem3.chrUn.lift \
	error blastCe6SG.psl | sort -k14,14 -k16,16n \
	| hgLoadPsl -table=blastCe6SG caeRem3 stdin

    # check coverage
    featureBits caeRem3 blastCe6SG 
    #	20355944 bases of 138406388 (14.707%) in intersection
    featureBits caeRem2 blastCe4SG 
    #	19763359 bases of 146898439 (13.454%) in intersection

    featureBits ce6 sangerGene
    #	28134889 bases of 100281426 (28.056%) in intersection

    ssh kkstore01
    rm -rf /cluster/data/caeRem3/bed/tblastn.ce6SG/blastOut
    rm -rf /cluster/bluearc/worms/caeRem3/bed/tblastn.ce6SG

###########################################################################
# ELEGANS (ce9) PROTEINS TRACK (DONE - 2010-10-07 - Hiram)
    cd /hive/data/genomes/caeRem3
    mv blastDb blastDb.2008-06-17
    mv blastDb.lft blastDb.lft.2008-06-17
    mkdir blastDb
    twoBitToFa caeRem3.unmasked.2bit temp.fa
    faSplit gap temp.fa 1000000 blastDb/x -lift=blastDb.lft
    #	177 pieces of 177 written
    rm temp.fa
    cd blastDb
    for i in *.fa
    do
	/scratch/data/blast-2.2.11/bin/formatdb -i $i -p F
    done
    rm *.fa

    ## create the query protein set
    mkdir -p /hive/data/genomes/caeRem3/bed/tblastn.ce9SG
    cd /hive/data/genomes/caeRem3/bed/tblastn.ce9SG
    echo  /hive/data/genomes/caeRem3/blastDb/*.nsq | xargs ls -S \
	| sed "s/\.nsq//"  > query.lst
    wc -l query.lst
    # 9538 query.lst

   # we want around 50000 jobs
    calc `wc -l /hive/data/genomes/ce9/bed/blat.ce9SG/ce9SG.psl  | cut -d" " -f1`/\(50000/`wc -l query.lst | cut -d" " -f1`\)
    #	28103/(50000/177) = 99.484620

   mkdir sgfa
   split -l 90 /hive/data/genomes/ce9/bed/blat.ce9SG/ce9SG.psl sgfa/sg
   cd sgfa
   for i in *; do 
     nice pslxToFa $i $i.fa; 
     rm $i; 
   done
   cd ..
   ls -1S sgfa/*.fa > sg.lst
   mkdir blastOut
   for i in `cat sg.lst`; do  mkdir blastOut/`basename $i .fa`; done
   
   cat << '_EOF_' > template
#LOOP
blastSome $(path1) {check in line $(path2)} {check out exists blastOut/$(root2)/q.$(root1).psl }
#ENDLOOP
'_EOF_'
    # << happy emacs

   cat << '_EOF_' > blastSome
#!/bin/sh
DB=caeRem3
BLASTMAT=/scratch/data/blast-2.2.11/data
SCR="/scratch/tmp/${DB}"
g=`basename $2`
D=`basename $1`
export BLASTMAT DB SCR g D
mkdir -p ${SCR}
cp -p $1.* ${SCR}
f=${SCR}/`basename $3`.$g
for eVal in 0.01 0.001 0.0001 0.00001 0.000001 1E-09 1E-11
do
if /scratch/data/blast-2.2.11/bin/blastall -M BLOSUM80 -m 0 -F no \
	-e $eVal -p tblastn -d ${SCR}/$D -i $2 -o $f.8
then
        mv $f.8 $f.1
        break;
fi
done
if test -f  $f.1
then
    if /cluster/bin/x86_64/blastToPsl $f.1 $f.2
    then
        liftUp -nosort -type=".psl" -nohead $f.3 \
	    /hive/data/genomes/${DB}/blastDb.lft carry $f.2 > /dev/null
        liftUp -nosort -type=".psl" -pslQ -nohead $3.tmp \
	    /hive/data/genomes/ce9/bed/blat.ce9SG/protein.lft warn $f.3 > /dev/null
        if pslCheck -prot $3.tmp
        then
            mv $3.tmp $3
            rm -f $f.1 $f.2 $f.3 $f.4 ${SCR}/$D.*
            rmdir --ignore-fail-on-non-empty ${SCR}
        fi
        exit 0
    fi
fi
rm -f $f.1 $f.2 $3.tmp $f.8 $f.3 $f.4 ${SCR}/$D.*
rmdir --ignore-fail-on-non-empty ${SCR}
exit 1
'_EOF_'
    # << happy emacs
    chmod +x blastSome

    ssh swarm
    cd /hive/data/genomes/caeRem3/bed/tblastn.ce9SG
    gensub2 query.lst sg.lst template jobList
    para create jobList
    para try, check, push, check etc.
# Completed: 55401 of 55401 jobs
# CPU time in finished jobs:     449440s    7490.67m   124.84h    5.20d  0.014 y
# IO & Wait Time:                171423s    2857.05m    47.62h    1.98d  0.005 y
# Average job time:                  11s       0.19m     0.00h    0.00d
# Longest finished job:              28s       0.47m     0.01h    0.00d
# Submission to last job:          1997s      33.28m     0.55h    0.02d

    # do the cluster run for chaining
    ssh swarm
    mkdir /hive/data/genomes/caeRem3/bed/tblastn.ce9SG/chainRun
    cd /hive/data/genomes/caeRem3/bed/tblastn.ce9SG/chainRun
    cat << '_EOF_' > template
#LOOP
chainOne $(path1) {check out exists+ ../blastOut/c.$(file1).psl}
#ENDLOOP
'_EOF_'
    # << happy emacs

    cat << '_EOF_' > chainOne
#!/bin/csh -fe
cd $1
set b = $1:t
cat q.*.psl | simpleChain -prot -outPsl -maxGap=50000 stdin \
/hive/data/genomes/caeRem3/bed/tblastn.ce9SG/blastOut/c.$b.psl
'_EOF_'
    # << happy emacs
    chmod +x chainOne

    ls -1dS /hive/data/genomes/caeRem3/bed/tblastn.ce9SG/blastOut/sg?? \
	> chain.lst
    gensub2 chain.lst single template jobList
    cd /hive/data/genomes/caeRem3/bed/tblastn.ce9SG/chainRun
    para create jobList
    para -maxJob=30 push
    para try, check, push, check etc.
# Completed: 313 of 313 jobs
# CPU time in finished jobs:       3546s      59.10m     0.99h    0.04d  0.000 y
# IO & Wait Time:                 10163s     169.38m     2.82h    0.12d  0.000 y
# Average job time:                  44s       0.73m     0.01h    0.00d
# Longest finished job:             394s       6.57m     0.11h    0.00d
# Submission to last job:           461s       7.68m     0.13h    0.01d

XXX - ready
    cd /hive/data/genomes/caeRem3/bed/tblastn.ce9SG/blastOut
    for i in sg??
    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 -T /scratch/tmp -k 14,14 -k 16,16n -k 17,17n u.*.psl m60* | uniq \
	> ../blastCe9SG.psl
    cd ..
    pslCheck blastCe9SG.psl
    #	checked: 33684 failed: 0 errors: 0

    # load table 
    ssh hgwdev
    cd /hive/data/genomes/caeRem3/bed/tblastn.ce9SG
    hgLoadPsl caeRem3 blastCe9SG.psl

    # check coverage
    featureBits caeRem3 blastCe9SG
    #	20302540 bases of 138406388 (14.669%) in intersection
    featureBits cb3 blastCe9SG
    #	18490367 bases of 108433446 (17.052%) in intersection
    featureBits caePb2 blastCe9SG
    #	23730009 bases of 170473138 (13.920%) in intersection
    featureBits caeJap3 blastCe9SG
    #	12894398 bases of 154057934 (8.370%) in intersection
    featureBits melHap1 blastCe9SG
    #	4376245 bases of 53017507 (8.254%) in intersection
    featureBits melInc1 blastCe9SG
    #	3882043 bases of 82095019 (4.729%) in intersection
    featureBits bruMal1 blastCe9SG 
    #	4424694 bases of 89235536 (4.958%) in intersection
    featureBits priPac2 blastCe9SG
    #	5436779 bases of 133634773 (4.068%) in intersection
    featureBits haeCon1 blastCe9SG
    #	4990746 bases of 278844984 (1.790%) in intersection
    featureBits ce9 sangerGene
    #	28689552 bases of 100286004 (28.608%) in intersection

    rm -rf blastOut

#########################################################################