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

#       $Id: bosTau6.txt,v 1.6 2010/05/06 16:27:44 chinhli Exp $


# Bos taurus (bovine) --  Bos_taurus_UMD_3.1 (2009-11-19)


# file template copied from oviAri1.txt


# Bos taurus (NCBI Project ID: 10709, Accession: GCA_000003055.3) 
#   by International Cow Genomics Consortium (ISGC)
#    assembly] sequence: 
# ftp://ftp.ncbi.nlm.nih.gov/genbank/genomes/Eukaryotes/vertebrates_mammals/Bos_taurus/Bos_taurus_UMD_3.1/

##########################################################################
# Download sequence (DONE - 2011-05-02 Chin)
    mkdir /hive/data/genomes/bosTau6
    cd /hive/data/genomes/bosTau6
    mkdir genbank
    cd genbank
    wget --timestamping -r --cut-dirs=6 --level=0 -nH -x \
        --no-remove-listing -np \
"ftp://ftp.ncbi.nlm.nih.gov/genbank/genomes/Eukaryotes/vertebrates_mammals/Bos_taurus/Bos_taurus_UMD_3.1/*"
    # FINISHED --2011-05-02 15:20:18--
    # Downloaded: 229 files, 1.9G in 11m 53s (2.73 MB/s)
    # Read ASSEMBLY_INFO 

    mkdir ucscChr
    # stay at genbank directory
    #   fixup the accession names to become UCSC chrom names

export S=Primary_Assembly/assembled_chromosomes
cut -f2 ${S}/chr2acc  | while read ACC
do
    C=`grep "${ACC}" ${S}/chr2acc | cut -f1`
    echo "${ACC} -> chr${C}"
    zcat ${S}/AGP/chr${C}.comp.agp.gz \
        | sed -e "s/^${ACC}/chr${C}/" | gzip > ucscChr/chr${C}.agp.gz
done


export S=Primary_Assembly/assembled_chromosomes
cut -f2 ${S}/chr2acc  | while read ACC
do
    C=`grep "${ACC}" ${S}/chr2acc | cut -f1`
    echo "${ACC} -> chr${C}"
    echo ">chr${C}" > ucscChr/chr${C}.fa
    zcat ${S}/FASTA/chr${C}.fa.gz | grep -v "^>" >> ucscChr/chr${C}.fa
    gzip ucscChr/chr${C}.fa &
done
   # Check them with faSize 
    faSize Primary_Assembly/assembled_chromosomes/FASTA/chr*.fa.gz
    #  2660906405 bases (20740230 N's 2640166175 real 2640166175 upper 
    #  0 lower) in 30 sequences in 30 files

    faSize ucscChr/chr*.fa.gz
    #  2660906405 bases (20740230 N's 2640166175 real 2640166175 upper 
    #  0 lower) in 30 sequences in 30 files

    # For unplaced scalfolds, named them as chrUn_xxxxxxxx
    #   where xxxxxx is the original access id as: chrUn_GL340781.1  
    #   The ".1" at the end need to be filter out since
    #   MySQL does not allow "." as part of the table name and 
    #   will casue problems in genbank task step later

export S=Primary_Assembly/unplaced_scaffolds
zcat ${S}/AGP/unplaced.scaf.agp.gz | grep "^#" > ucscChr/chrUn.agp
   # append the gap records
zcat ${S}/AGP/unplaced.scaf.agp.gz | grep -v "^#" \
        | sed -e "s/^/chrUn_/" -e "s/\.1//"  >> ucscChr/chrUn.agp
gzip ucscChr/chrUn.agp &

zcat ${S}/FASTA/unplaced.scaf.fa.gz \
        | sed -e "s#^>.*|tpg|#>chrUn_#; s#|.*##"  -e "s/\.1//"  \
        | gzip > ucscChr/chrUn.fa.gz
   # about 3286 sequences in the unplaced
zcat ucscChr/chrUn.fa.gz | grep "^>" | wc
   #     3286    3286   13144

   # Check them with faSize 
   faSize Primary_Assembly/unplaced_scaffolds/FASTA/unplaced.scaf.fa.gz
   #   9499556 bases (40 N's 9499516 real 9499516 upper 0 
   #   lower) in 3286 sequences in 1 files

   faSize ucscChr/chrUn.fa.gz
   #   9499556 bases (40 N's 9499516 real 9499516 upper 0 
   #   lower) in 3286 sequences in 1 files

#########################################################################
# Initial  database build (DONE - 2011-05-09 - Chin)
    cd /hive/data/genomes/bosTau6
    cat << '_EOF_' > bosTau6.config.ra
# Config parameters for makeGenomeDb.pl:
db bosTau6
clade mammal
genomeCladePriority 31
scientificName Bos taurus
commonName Cow
assemblyDate Nov. 2009
assemblyLabel UMD_3.1 (NCBI project 10708, accession GCA_000003055.3)
assemblyShortLabel Bos_taurus_UMD_3.1
orderKey 236
mitoAcc NC_006853
fastaFiles /hive/data/genomes/bosTau6/genbank/ucscChr/chr*.fa.gz
agpFiles /hive/data/genomes/bosTau6/genbank/ucscChr/chr*.agp.gz
# qualFiles none
dbDbSpeciesDir cow
taxId 9913
'_EOF_'
    # << happy emacs
    time makeGenomeDb.pl -noGoldGapSplit -workhorse=hgwdev bosTau6.config.ra \
	> makeGenomeDb.log 2>&1 &
    # real    21m31.993s

    # real    real    12m42.419s
    #	add the trackDb entries to the source tree, and the 2bit link:
    ln -s `pwd`/bosTau6.unmasked.2bit /gbdb/bosTau6/bosTau6.2bit

    #  Per instructions in makeGenomeDb.log:
#Search for '***' notes in each file in and make corrections (sometimes the
#files used for a previous assembly might make a better template):
#  description.html /cluster/data/bosTau6/html/{trackDb.ra,gap.html,gold.html}
#
#Then copy these files to your ~/kent/src/hg/makeDb/trackDb/cow/bosTau6
# - cd ~/kent/src/hg/makeDb/trackDb
# - edit makefile to add bosTau6 to DBS.
# - git add cow/bosTau6/*.{ra,html}
# - git commit -m "Added bosTau6 to DBS." makefile
# - git commit -m "Initial descriptions for bosTau6." cow/bosTau6/*.{ra,html}
# - git pull; git push
# - Run make update DBS=bosTau6 and make alpha when done.
# - (optional) Clean up /cluster/data/bosTau6/TemporaryTrackDbCheckout



#########################################################################
# RepeatMasker (DONE 2011-05-09 - Chin)
    mkdir /hive/data/genomes/bosTau6/bed/repeatMasker
    cd /hive/data/genomes/bosTau6/bed/repeatMasker
    time nice -n +19 doRepeatMasker.pl -buildDir=`pwd` \
	-workhorse=hgwdev -bigClusterHub=swarm -noSplit bosTau6 > do.log 2>&1 &
    #   real    355m43.294s
    cat faSize.rmsk.txt
    # 2670422299 bases (20740270 N's 2649682029 real 1370807360 upper 
    # 1278874669 lower) in 3317 sequences in 1 files
    # %47.89 masked total, %48.27 masked real

#########################################################################
# simpleRepeats (DONE - 2011-05-09 - Chin)
    mkdir /hive/data/genomes/bosTau6/bed/simpleRepeat
    cd /hive/data/genomes/bosTau6/bed/simpleRepeat

    time nice -n +19 doSimpleRepeat.pl -buildDir=`pwd` -workhorse=hgwdev \
	-dbHost=hgwdev -bigClusterHub=swarm -smallClusterHub=memk \
        bosTau6 > do.log 2>&1 &
    #   real    13m47.346s

    # Batch failed after 4 tries on 
    #   ./TrfRun.csh /hive/data/genomes/bosTau6/TrfPart/071/071.lst.bed
    # Command failed:
    # ssh -x -o 'StrictHostKeyChecking = no' -o 'BatchMode = yes' memk 
    #  nice /hive/data/genomes/bosTau6/bed/simpleRepeat/run.cluster/doTrf.csh
    # 071 failed due to no input TrfPart/071/071.lst.bed
    #  which is chrM and can be ignored.
    # to continue, create an empty TrfPart/071/071.lst.bed by:
    touch /hive/data/genomes/bosTau6/TrfPart/071/071.lst.bed 
    # run a para time > run.time on memk to get that file to exist:
    # Checking finished jobs
    # Completed: 71 of 72 jobs
    # Crashed: 1 jobs
    # CPU time in finished jobs: 12340s     205.67m     3.43h    0.14d  0.000 y
    # IO & Wait Time:              236s       3.93m     0.07h    0.00d  0.000 y
    # Average job time:            177s       2.95m     0.05h    0.00d
    # Longest finished job:        606s      10.10m     0.17h    0.01d
    # Submission to last job:      984s      16.40m     0.27h    0.01d


    #   continue
    time doSimpleRepeat.pl -buildDir=`pwd` -bigClusterHub=swarm \
        -dbHost=hgwdev -workhorse=hgwdev -smallClusterHub=memk \
        -continue=filter bosTau6 > filter.log 2>&1 &
    # real	0m20.430s

    cat fb.simpleRepeat 
    # 32390355 bases of 2649685036 (1.222%) in intersection

    #	add to the repeatMasker
    cd /hive/data/genomes/bosTau6
    twoBitMask bosTau6.rmsk.2bit -add bed/simpleRepeat/trfMask.bed bosTau6.2bit
    #	safe to ignore warnings about >=13 fields
    twoBitToFa bosTau6.2bit stdout | faSize stdin > bosTau6.2bit.faSize.txt
    cat bosTau6.2bit.faSize.txt
    # 2670422299 bases (20740270 N's 2649682029 real 1370313596 upper 
    #    1279368433 lower) in 3317 sequences in 1 files
    # %47.91 masked total, %48.28 masked real

    # double check with featureBits
    featureBits -countGaps bosTau6 gap
    # 20740270 bases of 2670422299 (0.777%) in intersection

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

#########################################################################
# Marking *all* gaps - they are not all in the AGP file
#	(DONE - 2011-05-10 - Chin)
    mkdir /hive/data/genomes/bosTau6/bed/allGaps
    cd /hive/data/genomes/bosTau6/bed/allGaps

    time nice -n +19 findMotif -motif=gattaca -verbose=4 \
	-strand=+ ../../bosTau6.unmasked.2bit > findMotif.txt 2>&1
    #   real    1m10.942s
    grep "^#GAP " findMotif.txt | sed -e "s/^#GAP //" > allGaps.bed
    featureBits bosTau6 -not gap -bed=notGap.bed
    # 2649685036 bases of 2649685036 (100.000%) in intersection
    featureBits bosTau6 allGaps.bed notGap.bed -bed=new.gaps.bed
    # 3007 bases of 2649685036 (0.000%) in intersection
    #	what is the highest index in the existing gap table:
    hgsql -N -e "select ix from gap;" bosTau6 | sort -n | tail -1
    #	17354

    # use tcsh and ctrl-c to create the here doc
    cat << '_EOF_' > mkGap.pl
#!/usr/bin/env perl

use strict;
use warnings;

my $ix=`hgsql -N -e "select ix from gap;" bosTau6 | 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 bosTau6 other.bed
    # 3007 bases of 2649685036 (0.000%) in intersection
    hgLoadBed -sqlTable=$HOME/kent/src/hg/lib/gap.sql \
	-noLoad bosTau6 otherGap other.bed
    #  Loaded 2011 elements of size 8
    #	adding this many:
    wc -l bed.tab
    #   2011 bed.tab 
    #	starting with this many
    hgsql -e "select count(*) from gap;" bosTau6
    #	72454
    hgsql bosTau6 -e 'load data local infile "bed.tab" into table gap;'
    #	result count:
    hgsql -e "select count(*) from gap;" bosTau6
    #	74465
    # == 72454 + 2011


########################################################################
# Create kluster run files (DONE - 2011-05-10 - Chin)
    # numerator is bosTau6 gapless bases "real" as reported by: 
    featureBits -noRandom -noHap bosTau6 gap
    # 20740230 bases of 2640182513 (0.786%) in intersection

    # denominator is hg19 gapless bases as reported by:
    #	featureBits -noRandom -noHap hg19 gap
    #     234344806 bases of 2861349177 (8.190%) in intersection
    # 1024 is threshold used for human -repMatch:
    calc \( 2640182513 / 2861349177 \) \* 1024
    ( 2640182513 / 2861349177 ) * 1024 = 944.850393
    # ==> use -repMatch=900 according to size scaled down from 1024 for human.
    #	and rounded down to nearest 50
    cd /hive/data/genomes/bosTau6
    blat bosTau6.2bit \
	 /dev/null /dev/null -tileSize=11 -makeOoc=jkStuff/bosTau6.11.ooc \
	-repMatch=900 &
    #	Wrote 33615 overused 11-mers to jkStuff/bosTau6.11.ooc
    mkdir /hive/data/staging/data/bosTau6
    cp -p bosTau6.2bit jkStuff/bosTau6.11.ooc /hive/data/staging/data/bosTau6
    cp -p chrom.sizes /hive/data/staging/data/bosTau6
    #	check non-bridged gaps to see what the typical size is:
    hgsql -N \
	-e 'select * from gap where bridge="no" order by size;' bosTau6 \
	| sort -k7,7nr
    #   most non-bridges gaps have size = 100, follow pig's example (5000) 
    #	decide on a minimum gap for this break
    gapToLift -verbose=2 -minGap=100 bosTau6 jkStuff/nonBridged.lft \
	-bedFile=jkStuff/nonBridged.bed
    cp -p jkStuff/nonBridged.lft \
	/hive/data/staging/data/bosTau6/bosTau6.nonBridged.lft
    # ask cluster-admin to copy (evry time if any file chsnged)
    #    /hive/data/staging/data/bosTau6 directory to
    #    /scratch/data/bosTau6 on cluster nodes
    # before porceed to genbank step

########################################################################
# GENBANK AUTO UPDATE (DONE - 2011-05-10 - Chin)
    ssh hgwdev
    cd $HOME/kent/src/hg/makeDb/genbank
    git pull

    # edit etc/genbank.conf to add bosTau6 just before susScr2

# bosTau6 (Cow)
bosTau6.serverGenome = /hive/data/genomes/bosTau6/bosTau6.2bit
bosTau6.clusterGenome = /scratch/data/bosTau6/bosTau6.2bit
bosTau6.ooc = /scratch/data/bosTau6/bosTau6.11.ooc
bosTau6.lift = no
bosTau6.perChromTables = no
bosTau6.refseq.mrna.native.pslCDnaFilter  = ${ordered.refseq.mrna.native.pslCDnaFilter}
bosTau6.refseq.mrna.xeno.pslCDnaFilter    = ${ordered.refseq.mrna.xeno.pslCDnaFilter}
bosTau6.genbank.mrna.native.pslCDnaFilter = ${ordered.genbank.mrna.native.pslCDnaFilter}
bosTau6.genbank.mrna.xeno.pslCDnaFilter   = ${ordered.genbank.mrna.xeno.pslCDnaFilter}
bosTau6.genbank.est.native.pslCDnaFilter  = ${ordered.genbank.est.native.pslCDnaFilter}
bosTau6.genbank.est.xeno.pslCDnaFilter    = ${ordered.genbank.est.xeno.pslCDnaFilter}
bosTau6.downloadDir = bosTau6
bosTau6.refseq.mrna.native.load  = yes
bosTau6.refseq.mrna.xeno.load = yes
bosTau6.refseq.mrna.xeno.loadDesc  = yes


    git add etc/genbank.conf
    git commit -m "Added bosTau6" etc/genbank.conf
    git push
    # update /cluster/data/genbank/:
    make etc-update


# Edit src/lib/gbGenome.c to add new species.  Skipped 

    screen	#  control this business with a screen since it takes a while
    cd /cluster/data/genbank
     time nice -n +19 ./bin/gbAlignStep -initial bosTau6 &
    #   logFile: var/build/logs/2011.05.10-15:10:42.bosTau6.initalign.log
    #   real    736m26.569s
    #   To re-do, rm the dir first:
    #     /cluster/data/genbank/data/aligned/genbank.183.0/bosTau6


    # load database when finished
    ssh hgwdev
    cd /cluster/data/genbank
    time nice -n +19 ./bin/gbDbLoadStep -drop -initialLoad bosTau6 &
    # logFile: var/dbload/hgwdev/logs/2011.05.11-08:24:53.dbload.log
    #   real    32m5.964s



    # enable daily alignment and update of hgwdev
    cd ~/kent/src/hg/makeDb/genbank
    gitpull
    # add bosTau6 to:
    #   etc/align.dbs
    #   etc/hgwdev.dbs
    git add etc/align.dbs
    git add etc/hgwdev.dbs
    git commit -m "Added bosTau6 - Cow" etc/align.dbs etc/hgwdev.dbs
    git push
    make etc-update

#########################################################################
# reset default position (DONE - 2011-05-11 - Chin)
#       Reset default position to an area of chr3 with a number of genes
    hgsql -e \
'update dbDb set defaultPos="chr3:21,592,157-21,596,916" where name="bosTau6";' \
        hgcentraltest


############################################################################
# ctgPos2 track - showing clone sequence locations on chromosomes
#	(DONE 2011-05-11 - Chin)
# NOTE - create bosTau6 entry in all.joiner since this is a new species
    mkdir /hive/data/genomes/bosTau6/bed/ctgPos2
    cd /hive/data/genomes/bosTau6/bed/ctgPos2
    cat << '_EOF_' > agpToCtgPos2.pl
#!/usr/bin/env perl

use warnings;
use strict;

my $argc = scalar(@ARGV);

if ($argc != 1) {
    printf STDERR "usage: zcat your.files.agp.gz | agpToCtgPos2.pl /dev/stdin > ctgPos2.tab\n";
    exit 255;
}

my $agpFile = shift;

open (FH, "<$agpFile") or die "can not read $agpFile";
while (my $line = <FH>) {
    next if ($line =~ m/^#/);
    chomp $line;
    my @a = split('\s+', $line);
    next if ($a[4] =~ m/^U$/);
    my $chrSize = $a[2]-$a[1]+1;
    my $ctgSize = $a[7]-$a[6]+1;
    die "sizes differ $chrSize != $ctgSize\n$line\n" if ($chrSize != $ctgSize);
    printf "%s\t%d\t%s\t%d\t%d\t%s\n", $a[5], $chrSize, $a[0], $a[1]-1, $a[2], $a[4];
}
close (FH);
'_EOF_'
    # << happy emacs

chmod 775 agpToCtgPos2.pl

export S=../../genbank/Primary_Assembly/assembled_chromosomes
cut -f2 ${S}/chr2acc | while read ACC
do
    C=`grep "${ACC}" ${S}/chr2acc | cut -f1`
    zcat ${S}/AGP/chr${C}.agp.gz \
        | sed -e "s/^${ACC}/chr${C}/"
done | ./agpToCtgPos2.pl /dev/stdin > ctgPos2.tab

    hgLoadSqlTab bosTau6 ctgPos2 $HOME/kent/src/hg/lib/ctgPos2.sql ctgPos2.tab
    # add the track ctgPos2 to src/hg/makeDb/trackDb/cow/bosTau6/trackDb.ra
    # at src/makeDb/trackdb do "make update DBS=bosTau6" or/and "make alpha"
    # based on result of
    # hgsql -N -e "select type from ctgPos2;" bosTau6 | sort | uniq
    # W
    # prepare the src/hg/makeDb/trackDb/cow/bosTau6/ctgPos2.html 


############################################################################
# running cpgIsland business (DONE -2011-05-16 - Chin)
    mkdir /hive/data/genomes/bosTau6/bed/cpgIsland
    cd /hive/data/genomes/bosTau6/bed/cpgIsland
    cvs -d /projects/compbio/cvsroot checkout -P hg3rdParty/cpgIslands
    cd hg3rdParty/cpgIslands
    #   needed to fixup this source, adding include to readseq.c:
#include "string.h"
    #   and to cpg_lh.c:
#include "unistd.h"
#include "stdlib.h"
    # and fixing a declaration in cpg_lh.c
    sed -e "s#\(extern char\* malloc\)#// \1#" cpg_lh.c > tmp.c
    mv tmp.c cpg_lh.c
    make
    cd ../../ 
    ln -s hg3rdParty/cpgIslands/cpglh.exe
    mkdir -p hardMaskedFa
    cut -f1 ../../chrom.sizes | while read C
do
    echo ${C}
    twoBitToFa ../../bosTau6.2bit:$C stdout \
        | maskOutFa stdin hard hardMaskedFa/${C}.fa
done

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

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


    chmod 775 runOne
    gensub2 chr.list single template jobList
    para create jobList
    para try
    para check ... etc
    para time
    para problems
    para status
    # then, kick it with para push
    # check it with plb
    # when all are done, para time shows:
    #  Checking finished jobs
    # Completed: 3317 of 3317 jobs
    # CPU time in finished jobs:   197s       3.29m     0.05h    0.00d  0.000 y
    # IO & Wait Time:             9373s     156.21m     2.60h    0.11d  0.000 y
    # Average job time:              3s       0.05m     0.00h    0.00d
    # Longest finished job:         18s       0.30m     0.01h    0.00d
    # Submission to last job:      677s      11.28m     0.19h    0.01d

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

    ssh hgwdev
    cd /hive/data/genomes/bosTau6/bed/cpgIsland
    hgLoadBed bosTau6 cpgIslandExt -tab \
      -sqlTable=$HOME/kent/src/hg/lib/cpgIslandExt.sql cpgIsland.bed
    # Loaded 37638 elements of size 10
    # Sorted
    # Creating table definition for cpgIslandExt
    # Saving bed.tab
    # Loading bosTau6

    #   cleanup
    rm -fr hardMaskedFa


###########################################################################
# ornAna1 Platypus LASTZ/CHAIN/NET (DONE - 2011-05-20 - Chin)
    screen # use screen to control this job
    mkdir /cluster/data/bosTau6/bed/blastzOrnAna1.2011-05-20
    cd /cluster/data/bosTau6/bed/blastzOrnAna1.2011-05-20

    cat << '_EOF_' > DEF
# Cow vs. Platypus
BLASTZ_Y=3400
BLASTZ_L=6000
BLASTZ_K=2200
BLASTZ_Q=/cluster/data/blastz/HoxD55.q
BLASTZ_M=50

# TARGET: Cow bosTau6
SEQ1_DIR=/scratch/data/bosTau6/bosTau6.2bit
SEQ1_LEN=/cluster/data/bosTau6/chrom.sizes
# Maximum number of scaffolds that can be lumped together
SEQ1_LIMIT=300
SEQ1_CHUNK=20000000
SEQ1_LAP=0

# QUERY: Platypus ornAna1
SEQ2_DIR=/scratch/data/ornAna1/ornAna1.2bit
SEQ2_LEN=/cluster/data/ornAna1/chrom.sizes
SEQ2_CHUNK=20000000
SEQ2_LIMIT=300
SEQ2_LAP=0

BASE=/cluster/data/bosTau6/bed/blastzOrnAna1.2011-05-20
TMPDIR=/scratch/tmp
'_EOF_'
    # << happy emacs
    
    time nice -n +19 doBlastzChainNet.pl `pwd`/DEF \
        -syntenicNet -noDbNameCheck \
        -workhorse=hgwdev -smallClusterHub=memk -bigClusterHub=swarm \
        -chainMinScore=5000 -verbose=2 -chainLinearGap=loose > do.log 2>&1 &
    #   real     329m59.325s

    # filter with doRecipBest.pl 
    time doRecipBest.pl -workhorse=hgwdev -buildDir=`pwd` \
        bosTau6 ornAna1 > rbest.log 2>&1
    # real    15m15.120s
    cat fb.bosTau6.chainOrnAna1Link.txt
    #   195910907 bases of 2649682029 (7.394%) in intersection
    cd /hive/data/genomes/bosTau6/bed
    ln -s blastzOrnAna1.2011-05-20 lastz.ornAna1


    mkdir /cluster/data/ornAna1/bed/blastz.bosTau6.swap
    cd /cluster/data/ornAna1/bed/blastz.bosTau6.swap
    time nice -n +19 doBlastzChainNet.pl \
        /cluster/data/bosTau6/bed/blastzOrnAna1.2011-05-20/DEF \
        -chainMinScore=5000 -verbose=2 -smallClusterHub=memk \
        -swap -syntenicNet \
          -chainLinearGap=loose -bigClusterHub=swarm > swap.log 2>&1 &
    #   real       59m45.336
    cat fb.ornAna1.chainBosTau6Link.txt
    #   190401000 bases of 1842236818 (10.335%) in intersection
# TBD for all chain/net:
 *** All done !  Elapsed time: 59m18s
 *** Make sure that goldenPath/ornAna1/vsBosTau6/README.txt is accurate.
 *** Add {chain,net}BosTau6 tracks to trackDb.ra if necessary.

    cd /cluster/data/ornAna1/bed
    ln -s blastz.bosTau6.swap lastz.bosTau6

#####################################################################
# susScr2 Pig BLASTZ/CHAIN/NET (DONE - 2011-05-24 - Chin)
    screen # use a screen to manage this multi-day job
    mkdir /hive/data/genomes/bosTau6/bed/lastzSusScr2.2011-05-24
    cd /hive/data/genomes/bosTau6/bed/lastzSusScr2.2011-05-24
    
    cat << '_EOF_' > DEF
# Pig vs. Cow
BLASTZ_M=50

# TARGET: Cow BosTau6
SEQ1_DIR=/scratch/data/bosTau6/bosTau6.2bit
SEQ1_LEN=/scratch/data/bosTau6/chrom.sizes
SEQ1_CHUNK=10000000
SEQ1_LAP=10000

# QUERY: Pig SusScr2
SEQ2_DIR=/scratch/data/susScr2/susScr2.2bit
SEQ2_LEN=/scratch/data/susScr2/chrom.sizes
SEQ2_CHUNK=10000000
SEQ2_LAP=0

BASE=/hive/data/genomes/bosTau6/bed/lastzSusScr2.2011-05-24
TMPDIR=/scratch/tmp
'_EOF_'
    # << this line keeps emacs coloring happy

    time nice -n +19 doBlastzChainNet.pl -verbose=2 \
        `pwd`/DEF \
        -noLoadChainSplit -syntenicNet \
        -workhorse=hgwdev -smallClusterHub=memk -bigClusterHub=swarm \
        -chainMinScore=3000 -chainLinearGap=medium > do.log 2>&1 &
    #   real    2198m44.401s
    cat fb.bosTau6.chainSusScr2Link.txt
    #   1306738899 bases of 2649682029 (49.317%) in intersection
    cd /hive/data/genomes/bosTau6/bed
    ln -s lastzSusScr2.2011-05-24 lastz.susScr2

    # then swap
    mkdir /hive/data/genomes/susScr2/bed/blastz.bosTau6.swap
    cd /hive/data/genomes/susScr2/bed/blastz.bosTau6.swap
    time nice -n +19 doBlastzChainNet.pl -verbose=2 \
        /hive/data/genomes/bosTau6/bed/lastzSusScr2.2011-05-24/DEF \
        -swap -noLoadChainSplit -syntenicNet \
        -workhorse=hgwdev -smallClusterHub=memk -bigClusterHub=pk \
        -chainMinScore=3000 -chainLinearGap=medium > swap.log 2>&1 &
    #  real    380m53.612s

    cat fb.susScr2.chainBosTau6Link.txt
    # 1475323050 bases of 2231298548 (66.119%) in intersection



#########################################################################
#  SWAP mm9 LASTZ (DONE - 2010-05-17 - Chin)
    #	original alignment done at mm9.txt
    cd /hive/data/genomes/mm9/bed/lastzBosTau6.2011-05-17
    cat fb.mm9.chainBosTau6Link.txt
    # 699351036 bases of 2620346127 (26.689%) in intersection
    # Create link
    cd /hive/data/genomes/mm9/bed
    ln -s  lastzBosTau6.2011-05-17 lastz.bosTau6

    #   and the swap 
    mkdir /hive/data/genomes/bosTau6/bed/blastz.mm9.swap
    cd /hive/data/genomes/bosTau6/bed/blastz.mm9.swap
    time nice -n +19 doBlastzChainNet.pl -verbose=2 \
        /hive/data/genomes/mm9/bed/lastzBosTau6.2011-05-17/DEF \
        -swap -syntenicNet \
        -noLoadChainSplit \
        -workhorse=hgwdev -smallClusterHub=memk -bigClusterHub=swarm \
        -chainMinScore=3000 -chainLinearGap=medium > swap.log 2>&1
    #   real     53m5.237s
    cat fb.bosTau6.chainMm9Link.txt
    # 688894115 bases of 2649682029 (25.999%) in intersection
    cd /hive/data/genomes/bosTau6/bed
    ln -s blastz.mm9.swap lastz.mm9

############################################################################
#  SWAP hg19 LASTZ (DONE 2011-05-16 - Chin)
    #	original alignment
    cd /hive/data/genomes/hg19/bed/lastzBosTau6.2011-05-16
    cat fb.hg19.chainBosTau6Link.txt
    # 1370696434 bases of 2897316137 (47.309%) in intersection

    #   and the swap
    mkdir /hive/data/genomes/bosTau6/bed/blastz.hg19.swap
    cd /hive/data/genomes/bosTau6/bed/blastz.hg19.swap
    time nice -n +19 doBlastzChainNet.pl -verbose=2 \
        /hive/data/genomes/hg19/bed/lastzBosTau6.2011-05-16/DEF \
        -swap -syntenicNet \
        -noLoadChainSplit \
        -workhorse=hgwdev -smallClusterHub=memk -bigClusterHub=swarm \
        -chainMinScore=3000 -chainLinearGap=medium > swap.log 2>&1
    #   real     72m57.889s
    cat fb.bosTau6.chainHg19Link.txt
    #   1336966333 bases of 2649682029 (50.458%) in intersection
    cd /hive/data/genomes/bosTau6/bed
    ln -s blastz.hg19.swap lastz.hg19


############################################################################
#  SWAP canFam2 LASTZ (DONE 2011-05-25 - Chin)
    #   original alignment
    cd /hive/data/genomes/canFam2/bed/lastzBosTau6.2011-05-25
    cat fb.canFam2.chainBosTau6Link.txt
    # 1389712912 bases of 2384996543 (58.269%) in intersection

    #   and the swap 
    mkdir /hive/data/genomes/bosTau6/bed/blastz.canFam2.swap
    cd /hive/data/genomes/bosTau6/bed/blastz.canFam2.swap
    time nice -n +19 doBlastzChainNet.pl -verbose=2 \
        /hive/data/genomes/canFam2/bed/lastzBosTau6.2011-05-25/DEF \
        -swap -syntenicNet  \
        -noLoadChainSplit \
        -workhorse=hgwdev -smallClusterHub=memk -bigClusterHub=swarm \
        -chainMinScore=3000 -chainLinearGap=medium > swap.log 2>&1
    #   real     88m38.555s
    cat fb.bosTau6.chainCanFam2Link.txt
    # 1404315725 bases of 2649682029 (52.999%) in intersection
    cd /hive/data/genomes/bosTau6/bed
    ln -s blastz.canFam2.swap lastz.canFam2

############################################################################
#  SWAP rn4 LASTZ (DONE 2011-05-31 - Chin)
    cd /hive/data/genomes/rn4/bed/blastzBosTau6.2011-05-31
    cat fb.rn4.chainBosTau6Link.txt
    # 660647959 bases of 2571531505 (25.691%) in intersection

    #   and the swap
    mkdir /hive/data/genomes/bosTau6/bed/blastz.rn4.swap
    cd /hive/data/genomes/bosTau6/bed/blastz.rn4.swap
    time nice -n +19 doBlastzChainNet.pl -verbose=2 \
        /hive/data/genomes/rn4/bed/blastzBosTau6.2011-05-31/DEF \
        -swap -noLoadChainSplit -syntenicNet \
        -workhorse=hgwdev -smallClusterHub=memk -bigClusterHub=swarm \
        -chainMinScore=3000 -chainLinearGap=medium > swap.log 2>&1 &
    #   real    49m43.881

    cat fb.bosTau6.chainRn4Link.txt
    # 648232024 bases of 2649682029 (24.465%) in intersection
    cd /hive/data/genomes/bosTau6/bed
    ln -s blastz.rn4.swap lastz.rn4

############################################################################
#  SWAP monDom5 LASTZ (DONE 2011-05-31 - Chin)
    cd /hive/data/genomes/monDom5/bed/lastzBosTau6.2011-05-31
    cat fb.monDom5.chainBosTau6Link.txt
    #   207703303 bases of 3501660299 (5.932%) in intersection

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

    cat fb.bosTau6.chainMonDom5Link.txt
    #   199980646 bases of 2649682029 (7.547%) in intersection
    cd /hive/data/genomes/bosTau6/bed
    ln -s blastz.monDom5.swap lastz.monDom5


#########################################################################
# all.joiner update, downloads and in pushQ - (working 2011-06-14 - Chin)
    cd $HOME/kent/src/hg/makeDb/schema
    # fixup all.joiner until this is a clean output
    joinerCheck -database=bosTau6 -all all.joiner

    mkdir /hive/data/genomes/bosTau6/goldenPath
    cd /hive/data/genomes/bosTau6/goldenPath
    makeDownloads.pl bosTau6 > do.log 2>&1

    #	now ready for pushQ entry
    mkdir /hive/data/genomes/bosTau6/pushQ
    cd /hive/data/genomes/bosTau6/pushQ
    makePushQSql.pl bosTau6 > bosTau6.pushQ.sql 2> stderr.out
    #	check for errors in stderr.out, some are OK, e.g.:
    # WARNING: hgwdev does not have /gbdb/bosTau6/wib/gc5Base.wib
    # WARNING: hgwdev does not have /gbdb/bosTau6/wib/quality.wib
    # WARNING: hgwdev does not have /gbdb/bosTau6/bbi/quality.bw
    # WARNING: bosTau6 does not have seq
    # WARNING: bosTau6 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:
    #   tableList

    #	copy it to hgwbeta
    scp -p bosTau6.pushQ.sql hgwbeta:/tmp
    ssh hgwbeta
    cd /tmp
    hgsql qapushq < bosTau6.pushQ.sql
    #	in that pushQ entry walk through each entry and see if the
    #	sizes will set properly

############################################################################
#  bosTau6 - Cow - Ensembl Genes version 64  (DONE - 2011-10-12 - hiram)
    ssh hgwdev
    cd /hive/data/genomes/bosTau6
    cat << '_EOF_' > bosTau6.ensGene.ra
# required db variable
db bosTau6
# optional nameTranslation, the sed command that will transform
#       Ensemble names to UCSC names.  With quotes just to make sure.
nameTranslation "s/^\([0-9X][0-9]*\)/chr\1/; s/^MT/chrM/; s/^GJ\([0-9]*\).1/chrUn_GJ\1/"
'_EOF_'
#  << happy emacs

    doEnsGeneUpdate.pl  -ensVersion=64 bosTau6.ensGene.ra
    ssh hgwdev
    cd /hive/data/genomes/bosTau6/bed/ensGene.64
    featureBits bosTau6 ensGene
    # 42264149 bases of 2649682029 (1.595%) in intersection

############################################################################
# ucscToEnsembl name translation table (DONE - 2011-10-12 - Hiram)
    mkdir /hive/data/genomes/bosTau6/bed/ucscToEnsembl
    cd /hive/data/genomes/bosTau6/bed/ucscToEnsembl

    cut -f1 ../../chrom.sizes | while read C
do
E=`echo ${C} | sed -e "s/^chrM/MT/; s/chrUn_GJ\([0-9]*\)/GJ\1.1/; s/^chr//;"`
    echo -e "${C}\t${E}"
done > ucscToEnsembl.tab

    # determine longest UCSC name:
    cut -f1 ucscToEnsembl.tab | awk '{print length($0)}' | sort -nr | head -1
    #	14

    # use that size in this SQL definition index:
    cat << '_EOF_' > ucscToEnsembl.sql
# UCSC to Ensembl chr name translation
CREATE TABLE ucscToEnsembl (
    ucsc varchar(255) not null,        # UCSC chromosome name
    ensembl varchar(255) not null,     # Ensembl chromosome name
              #Indices
    PRIMARY KEY(ucsc(14))
);
'_EOF_'

    hgLoadSqlTab bosTau6 ucscToEnsembl ucscToEnsembl.sql ucscToEnsembl.tab

#############################################################################
# LIFTOVER TO bosTau4 (DONE - 2011-10-20 - Chin)

    sh hgwdev
    cd /hive/data/genomes/bosTau6
    ln -s jkStuff/bosTau6.11.ooc 11.ooc
    time nice -n +19 doSameSpeciesLiftOver.pl -verbose=2 \
        -bigClusterHub=swarm -dbHost=hgwdev -workhorse=hgwdev \
         bosTau6 bosTau4 > doLiftOverToBosTau4.log 2>&1 &
    # real    110m26.710s
    pushd /usr/local/apache/htdocs-hgdownload/goldenPath/bosTau6/liftOver/
    md5sum bosTau6ToBosTau4.over.chain.gz >> md5sum.txt
    popd

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

