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

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


# Meleagris gallopavo (turkey) -- Turkey Genome Consortium Turkey_2.01
(2009-12-10)

ATE:	10-DEC-2009
ORGANISM:	Meleagris gallopavo
TAXID:	9103
ASSEMBLY LONG NAME:	Turkey_2.01
ASSEMBLY SHORT NAME:	Turkey_2.01
ASSEMBLY SUBMITTER:	Turkey Genome Consortium
ASSEMBLY TYPE:	Haploid
NUMBER OF ASSEMBLY-UNITS:	1
Assembly Accession:	GCA_000146605.1


# file template copied from susScr2.txt


# Meleagris gallopavo (NCBI Project ID: 10805, Accession: GCA_000146605.1)
#    by Turkey Genome Consortium
#    sequence:
# ftp://ftp.ncbi.nlm.nih.gov/genbank/genomes/Eukaryotes/vertebrates_other/Meleagris_gallopavo/Turkey_2.01/

#	http://www.ncbi.nlm.nih.gov/Traces/wgs/?val=ADDD00


##########################################################################
# Download sequence (DONE - 2010-10-22 - Chin)
    mkdir /hive/data/genomes/melGal1
    cd /hive/data/genomes/melGal1
    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_other/Meleagris_gallopavo/Turkey_2.01/*"
    # FINISHED --2010-10-22 08:28:53--
    # Downloaded: 243 files, 611M in 9m 19s (1.09 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
    #   2784748484 bases (1600136831 N's 1184611653 real 1184611653 upper
    #   0 lower) in 27 sequences in 27 files
    faSize ucscChr/chr*.fa.gz
    #  1061800384 bases (125894715 N's 935905669 real 935905669 upper 
    #  0 lower) in 5890 sequences in 33 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#^>.*|gb|#>chrUn_#; s#|.*##"  -e "s/\.1//"  \
        | gzip > ucscChr/chrUn.fa.gz
   # about 5858 sequences in the unplaced
    zcat ucscChr/chrUn.fa.gz | grep "^>" | wc
   #    5858    5858   93728

    # Check them with faSize 
    faSize Primary_Assembly/unplaced_scaffolds/FASTA/unplaced.scaf.fa.gz
    #  21513314 bases (2885592 N's 18627722 real 18627722 upper 
    #  0 lower) in 5858 sequences in 1 files
    #  75747883 bases (59104875 N's 16643008 real 16643008 upper 0 lower)
    #     in 1190 sequences in 1 files
    faSize ucscChr/chrUn.fa.gz
    # 21513314 bases (2885592 N's 18627722 real 18627722 upper 
    # 0 lower) in 5858 sequences in 1 files
   
    # N50
    mkdir N50
    export S=Primary_Assembly/assembled_chromosomes
    cut -f2 ${S}/chr2acc  | while read ACC
    do
    	C=`grep "${ACC}" ${S}/chr2acc | cut -f1`
    	echo ">chr${C}" > ucscChr/chr${C}.fa
    	zcat ucscChr/chr${C}.fa.gz  >> N50/melGal1.fa
    done
    zcat ucscChr/chrUn.fa.gz  >> N50/melGal1.fa
    faCount  N50/melGal1.fa | awk ' /^chr/ {print $1, $2}' > N50/chrom.sizes
    n50.pl N50/chrom.sizes
#       reading: N50/chrom.sizes
#       contig count: 5890, total size: 1061800384, one half size: 530900192
# cumulative    N50 count       contig  contig size
502447651       4       chrZ    81011772
530900192 one half size
577312103       5       chr4    74864452

     rm N50/melGal1.fa

#########################################################################
# Initial makeGenomeDb.pl (DONE - 2010-11-03 - Chin)
    cd /hive/data/genomes/melGal1
    cat << '_EOF_' > melGal1.config.ra
# Config parameters for makeGenomeDb.pl:
db melGal1
clade vertebrate
genomeCladePriority 55
scientificName Meleagris gallopavo
commonName Turkey
assemblyDate Dec. 2009
assemblyLabel TGC (NCBI Project ID: 10805, Accession: GCA_000146605.1)
assemblyShortLabel TGC Turkey_2.01 
orderKey 436
mitoAcc NC_010195
fastaFiles /hive/data/genomes/melGal1/genbank/ucscChr/chr*.fa.gz
agpFiles /hive/data/genomes/melGal1/genbank/ucscChr/chr*.agp.gz
# qualFiles none
dbDbSpeciesDir turkey
taxId 9103
'_EOF_'
    # << happy emacs
    time makeGenomeDb.pl -noGoldGapSplit -workhorse=hgwdev melGal1.config.ra \
	> makeGenomeDb.log 2>&1 &
    # real    7m58.244s
    #	add the trackDb entries to the source tree, and the 2bit link:
    ln -s `pwd`/melGal1.unmasked.2bit /gbdb/melGal1/melGal1.2bit

    #  Per instructions in makeGenomeDb.log:
    # cd ~/kent/src/hg/makeDb/trackDb
    # edit makefile to add melGal1 to DBS.
    # git add turkey/melGal1/*.{ra,html}
    # git commit -m "Added melGal1 to DBS." makefile
    # git commit -m "Initial descriptions for melGal1."
    # turkey/melGal1/*.{ra,html}
    # git pull; git push
    # Run make update DBS=melGal1 and make alpha when done.
    # (optional) Clean up /cluster/data/melGal1/TemporaryTrackDbCheckout


#########################################################################
# RepeatMasker (DONE - 2010-11-04 - Chin)
    mkdir /hive/data/genomes/melGal1/bed/repeatMasker
    cd /hive/data/genomes/melGal1/bed/repeatMasker

    time nice -n +19 doRepeatMasker.pl -buildDir=`pwd` \
	-workhorse=hgwdev -bigClusterHub=swarm -noSplit melGal1 > do.log 2>&1 &
    #   real    36m27.415s
    cat faSize.rmsk.txt
    # 1061817101 bases (125894715 N's 935922386 real 874501528 upper 
    # 61420858 lower) in 5891 sequences in 1 files
    # %5.78 masked total, %6.56 masked real

#########################################################################
# simpleRepeats ( DONE - 2010-11-05 - Chin)
    mkdir /hive/data/genomes/melGal1/bed/simpleRepeat
    cd /hive/data/genomes/melGal1/bed/simpleRepeat

    time nice -n +19 doSimpleRepeat.pl -buildDir=`pwd` -workhorse=hgwdev \
	-bigClusterHub=swarm -smallClusterHub=swarm melGal1 > do.log 2>&1 &
    #   real    7m12.044s
    cat fb.simpleRepeat 
    # 4997546 bases of 935922386 (0.534%) in intersection

    #	add to the repeatMasker
    cd /hive/data/genomes/melGal1
    twoBitMask melGal1.rmsk.2bit -add bed/simpleRepeat/trfMask.bed melGal1.2bit
    #	safe to ignore warnings about >=13 fields
    twoBitToFa melGal1.2bit stdout | faSize stdin > melGal1.2bit.faSize.txt
    cat melGal1.2bit.faSize.txt
    # 1061817101 bases (125894715 N's 935922386 real 874132579 upper 
    # 61789807 lower) in 5891 sequences in 1 files
    # %5.82 masked total, %6.60 masked real

#########################################################################
# Marking *all* gaps - they are all in the AGP file
#	(DONE 2010-11-05 - Chin)
    mkdir /hive/data/genomes/melGal1/bed/allGaps
    cd /hive/data/genomes/melGal1/bed/allGaps

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

    hgsql -e "select count(*) from gap;" melGal1
    #	 146748

########################################################################
# Create kluster run files (DONE 2010-11-05 - Chin)
    # numerator is melGal1 gapless bases "real" as reported by: 
    featureBits -noRandom -noHap melGal1 gap
    # 123009123 bases of 917294664 (13.410%) in intersection
    #     1600136831 bases of 1184628269 (135.075%) in interddsection

    # 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 \( 917294664 / 2861349177 \) \* 1024
    #  ( 917294664 / 2861349177 ) * 1024 = 328.275117
    # ==> use -repMatch=400 according to size scaled down from 1024 for human.
    #	and rounded down to nearest 50
    cd /hive/data/genomes/melGal1
    blat melGal1.2bit \
	 /dev/null /dev/null -tileSize=11 -makeOoc=jkStuff/melGal1.11.ooc \
	-repMatch=300 &
    #  Wrote 22129 overused 11-mers to jkStuff/melGal1.11.ooc
    mkdir /hive/data/staging/data/melGal1
    cp -p melGal1.2bit jkStuff/melGal1.11.ooc /hive/data/staging/data/melGal1
    cp -p chrom.sizes /hive/data/staging/data/melGal1
    #	check non-bridged gaps to see what the typical size is:
    hgsql -N \
	-e 'select * from gap where bridge="no" order by size;' melGal1 \
	| sort -k7,7nr
    #   most gaps have size > 100,000
    #	decide on a minimum gap for this break
    gapToLift -verbose=2 -minGap=20000 melGal1 jkStuff/nonBridged.lft \
	-bedFile=jkStuff/nonBridged.bed
    cp -p jkStuff/nonBridged.lft \
	/hive/data/staging/data/melGal1/melGal1.nonBridged.lft
    # ask cluster-admin to copy (evry time if any file chsnged)
    #    /hive/data/staging/data/melGal1 directory to cluster nodes
    #    /scratch/data/melGal1

########################################################################
# GENBANK AUTO UPDATE (DONE - 201w1-02-16 - Chin)
# 2011-10-07 Note: Re-configure per Mark's comments:
# genbank should not have native refseq enabled for this genome, as there are
# no native refseqs. Xeno refseqs should be enabled. Please fix configuration
# file and do a drop and reload.
    ssh hgwdev
    cd $HOME/kent/src/hg/makeDb/genbank
    git pull

    # edit etc/genbank.conf to add melGal1 just after chicken galGal2

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

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

# Edit src/lib/gbGenome.c to add new species.  With these two lines:
# static char *melGalNames[] = {"Meleagris gallopavo", NULL};
#   ... later ...
#    {"melGal", melGal1Names},
#  gbGenome.c is  in
#  /cluster/home/chinhli/kent/src/hg/makeDb/genbank/src/lib
# make and checkin

    make install-server
    git add src/lib/gbGenome.c
    git commit -m "adding melGal1 Turkey" src/lib/gbGenome.c
    git pull
    git push	    

    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 melGal1 &
    #   logFile: var/build/logs/2011.02.16-09:21:59.melGal1.initalign.log
    #   real     109m5.431s

    # load database when finished
    ssh hgwdev
    cd /cluster/data/genbank
    time nice -n +19 ./bin/gbDbLoadStep -drop -initialLoad melGal1 &
    #    logFile: var/dbload/hgwdev/logs/2011.02.16-12:45:15.dbload.log
    #   real    23m10.935s
    # Re-run load step 2011-10-20:
    #  logFile: var/dbload/hgwdev/logs/2011.10.20-10:10:19.dbload.log
    #   real  53m49.262s

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

#########################################################################
# reset default position as galGal3 Chicken (DONE 2011-03-08  - Chin)
#
    hgsql -e \
'update dbDb set defaultPos="chr5:57710001-57780000" where name="melGal1";' \
	hgcentraltest

############################################################################
# ctgPos2 track - showing clone sequence locations on chromosomes
#	(DONE 2011-03-08 - Chin)
# NOTE - create melGal1 entry in all.joiner since this is a new species
    mkdir /hive/data/genomes/melGal1/bed/ctgPos2
    cd /hive/data/genomes/melGal1/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/^N$/);
    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 +x 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 melGal1 ctgPos2 $HOME/kent/src/hg/lib/ctgPos2.sql ctgPos2.tab

############################################################################
# running cpgIsland business (DONE -2011-03-08 - Chin)
    mkdir /hive/data/genomes/melGal1/bed/cpgIsland
    cd /hive/data/genomes/melGal1/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 ../../melGal1.2bit:$C stdout \
	| maskOutFa stdin hard hardMaskedFa/${C}.fa
done

    ssh swarm
    cd /hive/data/genomes/melGal1/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:
    # para time
    # 5891 jobs in batch
    # 130198 jobs (including everybody's) in Parasol queue or running.
    # Checking finished jobs
    # Completed: 5891 of 5891 jobs
    # CPU time in finished jobs:   86s       1.44m     0.02h    0.00d  0.000 y
    # IO & Wait Time:           19261s     321.01m     5.35h    0.22d  0.001 y
    # Average job time:             3s       0.05m     0.00h    0.00d
    # Longest finished job:        26s       0.43m     0.01h    0.00d
    # Submission to last job:     369s       6.15m     0.10h    0.00d


    # 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/melGal1/bed/cpgIsland
    hgLoadBed melGal1 cpgIslandExt -tab \
      -sqlTable=$HOME/kent/src/hg/lib/cpgIslandExt.sql cpgIsland.bed
    # Loaded 13372 elements of size 10
    # Sorted
    # Creating table definition for cpgIslandExt
    # Saving bed.tab
    # Loading melGal1

    #	cleanup
    rm -fr hardMaskedFa

#########################################################################
# all.joiner update, downloads and in pushQ - ( DONE 2011-03-08 - Chin)
    cd $HOME/kent/src/hg/makeDb/schema
    # fixup all.joiner until this is a clean output
    # add melGal after galGal
    joinerCheck -database=melGal1 -all all.joiner

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

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

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

##########################################################################
#  BLATSERVERS ENTRY (DONE 2011-03-23 - Chin)
#       After getting a blat server assigned by the Blat Server Gods,
#   
# blat13
# trans port 17828
# untrans port 17829

    
    hgsql -e 'INSERT INTO blatServers (db, host, port, isTrans, canPcr) \
        VALUES ("melGal1", "blat13", "17828", "1", "0"); \
        INSERT INTO blatServers (db, host, port, isTrans, canPcr) \
        VALUES ("melGal1", "blat13", "17829", "0", "1");' \
            hgcentraltest
    #   test it with some sequence
#########################################################################
#  SWAP melGal1 lastz (DONE - 2011-03-24 - Chin)
    #   original alignment
    cat fb.galGal3.chainMelGal1Link.txt
    # 924755061 bases of 1042591351 (88.698%) in intersection

    #   and the swap
    mkdir /hive/data/genomes/melGal1/bed/blastz.galGal3.swap
    cd /hive/data/genomes/melGal1/bed/blastz.galGal3.swap
    time nice -n +19 doBlastzChainNet.pl -verbose=2 \
        /hive/data/genomes/galGal3/bed/lastzMelGal1.2011-03-24/DEF \
        -swap -noLoadChainSplit -syntenicNet \
        -workhorse=hgwdev -smallClusterHub=memk -bigClusterHub=swarm \
        -chainMinScore=3000 -chainLinearGap=medium > swap.log 2>&1 &
    # real    38m4.406s
    cat fb.melGal1.chainGalGal3Link.txt
    # 874579772 bases of 935922386 (93.446%) in intersection

#########################################################################
# BLASTZ/CHAIN/NET  Turkey vs Zebra Fish taeGut1 (DONE - 2011-03-28 - Chin)
    mkdir /hive/data/genomes/melGal1/bed/lastzTaeGut1.2011-03-28
    cd /hive/data/genomes/melGal1/bed/lastzTaeGut1.2011-03-28

    cat << '_EOF_' > DEF
# Zebra Finch vs Turkey

# TARGET: Turkey melGal1
SEQ1_DIR=/scratch/data/melGal1/melGal1.2bit
SEQ1_LEN=/scratch/data/melGal1/chrom.sizes
SEQ1_CHUNK=10000000
SEQ1_LAP=10000

# QUERY: Zebra Finch  taeGut1 - single chunk big enough to run entire chrom
SEQ2_DIR=/scratch/data/taeGut1/taeGut1.2bit
SEQ2_LEN=/scratch/data/taeGut1/chrom.sizes
SEQ2_CHUNK=10000000
SEQ2_LAP=0

BASE=/hive/data/genomes/melGal1/bed/lastzTaeGut1.2011-03-28
TMPDIR=/scratch/tmp
'_EOF_'
    # << happy emacs

    time nice -n +19 doBlastzChainNet.pl -verbose=2 \
        `pwd`/DEF \
        -syntenicNet \
        -noLoadChainSplit \
        -workhorse=hgwdev -smallClusterHub=memk -bigClusterHub=swarm \
        -chainMinScore=5000 -chainLinearGap=loose > do.log 2>&1
    # Elapsed time: 400m22s
    cat fb.melGal1.chainTaeGut1Link.txt
    # 671105409 bases of 935922386 (71.705%) in intersection
    cd  /hive/data/genomes/melGal1/bed
    ln -s lastzTaeGut1.2011-03-28 lastz.taeGut1

    #   running the swap 
    mkdir /hive/data/genomes/taeGut1/bed/blastz.melGal1.swap
    cd /hive/data/genomes/taeGut1/bed/blastz.melGal1.swap
    time nice -n +19 doBlastzChainNet.pl -verbose=2 \
        /hive/data/genomes/melGal1/bed/lastzTaeGut1.2011-03-28/DEF \
        -swap \
        -noLoadChainSplit \
        -workhorse=hgwdev -smallClusterHub=memk -bigClusterHub=swarm \
        -chainMinScore=5000 -chainLinearGap=loose > swap.log 2>&1
    #   real   36m7.017s 
    cat fb.taeGut1.chainMelGal1Link.txt
    #   795001886 bases of 1222864691 (65.011%) in intersection
    cd /hive/data/genomes/taeGut1/bed
    ln -s blastz.melGal1.swap lastz.melGal1

#########################################################################
#  SWAP hg19 lastz (DONE - 201i1-03-28 - Chin)
    #   original alignment
    cat fb.hg19.chainMelGal1Link.txt
    #   76647912 bases of 2897316137 (2.645%) in intersection

    #   and the swap
    mkdir /hive/data/genomes/melGal1/bed/blastz.hg19.swap
    cd /hive/data/genomes/melGal1/bed/blastz.hg19.swap
    time nice -n +19 doBlastzChainNet.pl -verbose=2 \
        /hive/data/genomes/hg19/bed/lastzMelGal1.2011-03-28/DEF \
        -swap \
        -noLoadChainSplit \
        -workhorse=hgwdev -smallClusterHub=memk -bigClusterHub=swarm \
        -chainMinScore=5000 -chainLinearGap=loose > swap.log 2>&1
    #   real    6m51.280s
    cat fb.melGal1.chainHg19Link.txt
    #   62120143 bases of 935922386 (6.637%) in intersection
    cd /hive/data/genomes/melGal1/bed
    ln -s blastz.hg19.swap lastz.hg19

#########################################################################
#  SWAP mm9 lastz (DONE - 201i1-03-28 - Chin)
    #   original alignment
    cat fb.mm9.chainMelGal1Link.txt
    #   62597891 bases of 2620346127 (2.389%) in intersection

    #   running the swap 
    mkdir /hive/data/genomes/melGal1/bed/blastz.mm9.swap
    cd /hive/data/genomes/melGal1/bed/blastz.mm9.swap
    time nice -n +19 doBlastzChainNet.pl -verbose=2 \
        /hive/data/genomes/mm9/bed/lastzMelGal1.2011-03-30/DEF \
        -swap \
        -noLoadChainSplit \
        -workhorse=hgwdev -smallClusterHub=memk -bigClusterHub=swarm \
        -chainMinScore=5000 -chainLinearGap=loose > swap.log 2>&1
    #   real    6m49.871s
    cat fb.melGal1.chainMm9Link.txt
    #   50649368 bases of 935922386 (5.412%) in intersection
    cd /hive/data/genomes/melGal1/bed
    ln -s blastz.mm9.swap lastz.mm9

#########################################################################
#  SWAP monDom5  lastz (DONE -  2011-03-30 - Chin)
    #   original alignment
    cat fb.monDom5.chainMelGal1Link.txt
    #   79077534 bases of 3501660299 (2.258%) in intersection

    #   running the swap 
    mkdir /hive/data/genomes/melGal1/bed/blastz.monDom5.swap
    cd /hive/data/genomes/melGal1/bed/blastz.monDom5.swap
    time nice -n +19 doBlastzChainNet.pl -verbose=2 \
        /hive/data/genomes/monDom5/bed/lastzMelGal1.2011-03-30/DEF \
        -swap \
        -noLoadChainSplit \
        -workhorse=hgwdev -smallClusterHub=memk -bigClusterHub=swarm \
        -chainMinScore=5000 -chainLinearGap=loose > swap.log 2>&1
    #   Elapsed time: 10m15s
    cat fb.melGal1.chainMonDom5Link.txt
    #   64933013 bases of 935922386 (6.938%) in intersection
    cd /hive/data/genomes/melGal1/bed
    ln -s lastz.monDom5.swap lastz.monDom5

#########################################################################
# BLASTZ/CHAIN/NET  Turkey vs Xenopus xenTro2 (DONE - 2011-04-02 - Chin)
    mkdir /hive/data/genomes/melGal1/bed/lastzXenTro2.2011-04-02
    cd /hive/data/genomes/melGal1/bed/lastzXenTro2.2011-04-02

    cat << '_EOF_' > DEF
# Xenopus vs Turkey

# TARGET: Turkey melGal1
SEQ1_DIR=/scratch/data/melGal1/melGal1.2bit
SEQ1_LEN=/scratch/data/melGal1/chrom.sizes
SEQ1_CHUNK=10000000
SEQ1_LAP=10000

# QUERY: Xenopus  xenTro2 - single chunk big enough to run entire chrom
SEQ2_DIR=/scratch/data/xenTro2/xenTro2.2bit
SEQ2_LEN=/scratch/data/xenTro2/chrom.sizes
SEQ2_CHUNK=20000000
SEQ2_LAP=0
SEQ2_LIMIT=100

BASE=/hive/data/genomes/melGal1/bed/lastzXenTro2.2011-04-02
TMPDIR=/scratch/tmp
'_EOF_'
    # << happy emacs 

    time nice -n +19 doBlastzChainNet.pl -verbose=2 \
        `pwd`/DEF \
        -noLoadChainSplit \
        -workhorse=hgwdev -smallClusterHub=memk -bigClusterHub=swarm \
        -chainMinScore=5000 -chainLinearGap=loose > do.log 2>&1
    # real    288m3.884s
    cat fb.melGal1.chainXenTro2Link.txt
    # 36288270 bases of 935922386 (3.877%) in intersection
    cd  /hive/data/genomes/melGal1/bed
    ln -s lastzXenTro2.2011-04-02 lastz.xenTro2
    
    #   running the swap 
    mkdir /hive/data/genomes/xenTro2/bed/blastz.melGal1.swap
    cd /hive/data/genomes/xenTro2/bed/blastz.melGal1.swap
    time nice -n +19 doBlastzChainNet.pl -verbose=2 \
        /hive/data/genomes/melGal1/bed/lastzXenTro2.2011-04-02/DEF \
        -swap \
        -noLoadChainSplit \
        -workhorse=hgwdev -smallClusterHub=memk -bigClusterHub=swarm \
        -chainMinScore=5000 -chainLinearGap=loose > swap.log 2>&1
    #   real    3m40.138s
    cat fb.xenTro2.chainMelGal1Link.txt
    #   41967849 bases of 1359412157 (3.087%) in intersection
    cd /hive/data/genomes/xenTro2/bed
    ln -s blastz.melGal1.swap lastz.melGal1

#######################################################################
# lastz swap from danRer7 (DONE - 2011-04-02 - Chin)
    # original alignment to danRer7
    cd /hive/data/genomes/danRer7/bed/lastzMelGal1.2011-04-01
    cat fb.danRer7.chainMelGal1Link.txt
    #   35772051 bases of 1409770109 (2.537%) in intersection
    #   running the swap 
    mkdir /hive/data/genomes/melGal1/bed/blastz.danRer7.swap
    cd /hive/data/genomes/melGal1/bed/blastz.danRer7.swap
    time nice -n +19 doBlastzChainNet.pl -verbose=2 \
        /hive/data/genomes/danRer7/bed/lastzMelGal1.2011-04-01/DEF \
        -swap \
        -noLoadChainSplit \
        -workhorse=hgwdev -smallClusterHub=memk -bigClusterHub=swarm \
        -chainMinScore=5000 -chainLinearGap=loose > swap.log 2>&1
    #   real    3m28.672s
    cat fb.melGal1.chainDanRer7Link.txt
    #   25141503 bases of 935922386 (2.686%) in intersection
    cd /hive/data/genomes/melGal1/bed
    ln -s blastz.danRer7.swap lastz.danRer7

#######################################################################
# lastz swap from rn4 (DONE - 2011-04-04 - Chin)
    # original alignment to rn4
    cd /hive/data/genomes/rn4/bed/lastzMelGal1.2011-04-04
    cat fb.rn4.chainMelGal1Link.txt
    #   60383476 bases of 2571531505 (2.348%) in intersection
    #   running the swap 
    mkdir /hive/data/genomes/melGal1/bed/blastz.rn4.swap
    cd /hive/data/genomes/melGal1/bed/blastz.rn4.swap
    time nice -n +19 doBlastzChainNet.pl -verbose=2 \
        /hive/data/genomes/rn4/bed/lastzMelGal1.2011-04-04/DEF \
        -swap \
        -noLoadChainSplit \
        -workhorse=hgwdev -smallClusterHub=memk -bigClusterHub=swarm \
        -chainMinScore=5000 -chainLinearGap=loose > swap.log 2>&1
    #   real    6m49.871s
    cat fb.melGal1.chainRn4Link.txt
    #   48636842 bases of 935922386 (5.197%) in intersection
    cd /hive/data/genomes/melGal1/bed
    ln -s blastz.rn4.swap lastz.rn4

############################################################################
# HUMAN (hg18) PROTEINS TRACK (DONE 2011-04-26 braney)
    # bash  if not using bash shell already
    cd /cluster/data/melGal1
    mkdir /cluster/data/melGal1/blastDb

    awk '{if ($2 > 1000000) print $1}' chrom.sizes > 1meg.lst
    twoBitToFa -seqList=1meg.lst  melGal1.2bit temp.fa
    faSplit gap temp.fa 1000000 blastDb/x -lift=blastDb.lft
#  1140 pieces of 1140 written
    rm temp.fa 1meg.lst

    awk '{if ($2 <= 1000000) print $1}' chrom.sizes > less1meg.lst
    twoBitToFa -seqList=less1meg.lst  melGal1.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
# 1163

    mkdir -p /cluster/data/melGal1/bed/tblastn.hg18KG
    cd /cluster/data/melGal1/bed/tblastn.hg18KG
    echo  ../../blastDb/*.nsq | xargs ls -S | sed "s/\.nsq//"  > query.lst
    wc -l query.lst
# 1163 query.lst

   # we want around 250000 jobs
   calc `wc /cluster/data/hg18/bed/blat.hg18KG/hg18KG.psl | awk '{print $1}'`/\(250000/`wc query.lst | awk '{print $1}'`\)
# 36727/(250000/1163) = 170.854004

   mkdir -p kgfa
   split -l 171 /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
# 215  215 2795 kg.lst

   mkdir -p blastOut
   for i in `cat kg.lst`; do  mkdir blastOut/`basename $i .fa`; done
   tcsh
   cd /cluster/data/melGal1/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/melGal1/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/melGal1/bed/tblastn.hg18KG
    gensub2 query.lst kg.lst blastGsub blastSpec
    para create blastSpec
#    para try, check, push, check etc.

    para time

# Completed: 250045 of 250045 jobs
# CPU time in finished jobs:    5922436s   98707.26m  1645.12h   68.55d  0.188 y
# IO & Wait Time:               1789013s   29816.89m   496.95h   20.71d  0.057 y
# Average job time:                  31s       0.51m     0.01h    0.00d
# Longest finished job:             330s       5.50m     0.09h    0.00d
# Submission to last job:         12017s     200.28m     3.34h    0.14d

    ssh swarm
    cd /cluster/data/melGal1/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: 215 of 215 jobs
# CPU time in finished jobs:       1662s      27.70m     0.46h    0.02d  0.000 y
# IO & Wait Time:                 85567s    1426.12m    23.77h    0.99d  0.003 y
# Average job time:                 406s       6.76m     0.11h    0.00d
# Longest finished job:             510s       8.50m     0.14h    0.01d
# Submission to last job:           537s       8.95m     0.15h    0.01d

    cd /cluster/data/melGal1/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: 35263 failed: 0 errors: 0

    # load table 
    ssh hgwdev
    cd /cluster/data/melGal1/bed/tblastn.hg18KG
    hgLoadPsl melGal1 blastHg18KG.psl

    # check coverage
    featureBits melGal1 blastHg18KG 
# 18362443 bases of 935922386 (1.962%) in intersection

    featureBits melGal1 blastHg18KG ensGene  -enrichment
# blastHg18KG 1.962%, ensGene 2.680%, both 1.787%, cover 91.07%, enrich 33.99x

    rm -rf blastOut
#end tblastn
############################################################################
# N-SCAN gene predictions (nscanGene) - (2011-07-07 markd)

    # obtained NSCAN predictions from Jeltje van Baren <jeltje@mbari.org>
    cd /hive/data/genomes/melGal1/bed/nscan
    cp -r /cluster/home/jeltje/hive/melGal1/for_browser/{melGal1.gtf,	melGal1.prot.fa,readme.txt} .
    gzip melGal1.* 
    chmod a-w *

    # load track
    gtfToGenePred -genePredExt melGal1.gtf.gz stdout | hgLoadGenePred -genePredExt melGal1 nscanGene stdin
    hgPepPred melGal1 generic nscanPep melGal1.prot.fa.gz
    rm *.tab

    # turkey/melGal1/trackDb.ra, add:
    track nscanGene override
    informant Nscan-EST on melGal1 with taeGut1 as informant and galGal3/taeGut1 parameter file. Not pseudogene masked. Used transmapped UCSC gene track as EST track.

    # verify top-level search spec, should produce no results and an error for egrep:
    hgsql -Ne 'select name from nscanGene' melGal1 | egrep -v '^chr[0-9a-zA-Z_]+\.([0-9]+|pasa)((\.[0-9a-z]+)?\.[0-9a-z]+)?$' |head

############################################################################
# lastz  Turkey vs Frog xenTro3 (DONE - 2011-09-22 - Hiram)
    # This was repeated below with different lastz parameters that
    # give a bit higher coverage
    mkdir /hive/data/genomes/melGal1/bed/lastzXenTro3.2011-09-22
    cd /hive/data/genomes/melGal1/bed/lastzXenTro3.2011-09-22

    cat << '_EOF_' > DEF
# Turkey vs. Frog

# TARGET: Turkey melGal1
SEQ1_DIR=/scratch/data/melGal1/melGal1.2bit
SEQ1_LEN=/scratch/data/melGal1/chrom.sizes
SEQ1_CHUNK=10000000
SEQ1_LAP=10000

# QUERY: Frog xenTro3
SEQ2_DIR=/scratch/data/xenTro3/xenTro3.2bit
SEQ2_LEN=/scratch/data/xenTro3/chrom.sizes
SEQ2_CHUNK=10000000
SEQ2_LAP=0
SEQ2_LIMIT=30

BASE=/hive/data/genomes/melGal1/bed/lastzXenTro3.2011-09-22
TMPDIR=/scratch/tmp
'_EOF_'
    # << happy emacs 

    time nice -n +19 doBlastzChainNet.pl -verbose=2 \
        `pwd`/DEF \
        -workhorse=hgwdev -smallClusterHub=memk -bigClusterHub=swarm \
        -chainMinScore=5000 -chainLinearGap=loose > do.log 2>&1 &
    #	real    648m6.498s
    cat fb.melGal1.chainXenTro3Link.txt 
    #	35993896 bases of 935922386 (3.846%) in intersection

    cd  /hive/data/genomes/melGal1/bed
    ln -s lastzXenTro3.2011-09-22 lastz.xenTro3
    
    #   running the swap 
    mkdir /hive/data/genomes/xenTro3/bed/blastz.melGal1.swap
    cd /hive/data/genomes/xenTro3/bed/blastz.melGal1.swap
    time nice -n +19 doBlastzChainNet.pl -verbose=2 \
        /hive/data/genomes/melGal1/bed/lastzXenTro3.2011-09-22/DEF \
        -swap \
        -workhorse=hgwdev -smallClusterHub=memk -bigClusterHub=swarm \
        -chainMinScore=5000 -chainLinearGap=loose > swap.log 2>&1
    #	real    3m29.407s

    cat fb.xenTro3.chainMelGal1Link.txt 
    #	42384698 bases of 1358334882 (3.120%) in intersection

    cd /hive/data/genomes/xenTro3/bed
    ln -s blastz.melGal1.swap lastz.melGal1

#######################################################################
# lastz  Turkey vs Frog xenTro3 (DONE - 2011-09-27 - Hiram)
    mkdir /hive/data/genomes/melGal1/bed/lastzXenTro3.2011-09-27
    cd /hive/data/genomes/melGal1/bed/lastzXenTro3.2011-09-27

    cat << '_EOF_' > DEF
# Turkey vs. Frog
BLASTZ_H=2000
BLASTZ_Y=3400
BLASTZ_L=8000
BLASTZ_K=2200
BLASTZ_Q=/scratch/data/blastz/HoxD55.q

# TARGET: Turkey melGal1
SEQ1_DIR=/scratch/data/melGal1/melGal1.2bit
SEQ1_LEN=/scratch/data/melGal1/chrom.sizes
SEQ1_CHUNK=10000000
SEQ1_LAP=10000

# QUERY: Frog xenTro3
SEQ2_DIR=/scratch/data/xenTro3/xenTro3.2bit
SEQ2_LEN=/scratch/data/xenTro3/chrom.sizes
SEQ2_CHUNK=10000000
SEQ2_LAP=0
SEQ2_LIMIT=30

BASE=/hive/data/genomes/melGal1/bed/lastzXenTro3.2011-09-27
TMPDIR=/scratch/tmp
'_EOF_'
    # << happy emacs

    time nice -n +19 doBlastzChainNet.pl -verbose=2 \
        `pwd`/DEF \
        -workhorse=hgwdev -smallClusterHub=memk -bigClusterHub=swarm \
        -chainMinScore=5000 -chainLinearGap=loose > do.log 2>&1 &
    #	real    537m35.283s
    cat fb.melGal1.chainXenTro3Link.txt
    #	49081729 bases of 935922386 (5.244%) in intersection

    cd  /hive/data/genomes/melGal1/bed
    ln -s lastzXenTro3.2011-09-27 lastz.xenTro3

    #   running the swap
    mkdir /hive/data/genomes/xenTro3/bed/blastz.melGal1.swap
    cd /hive/data/genomes/xenTro3/bed/blastz.melGal1.swap
    time nice -n +19 doBlastzChainNet.pl -verbose=2 \
        /hive/data/genomes/melGal1/bed/lastzXenTro3.2011-09-27/DEF \
        -swap \
        -workhorse=hgwdev -smallClusterHub=memk -bigClusterHub=swarm \
        -chainMinScore=5000 -chainLinearGap=loose > swap.log 2>&1
    #	real    4m9.023s

    cat fb.xenTro3.chainMelGal1Link.txt
    #	57087653 bases of 1358334882 (4.203%) in intersection

    cd /hive/data/genomes/xenTro3/bed
    ln -s blastz.melGal1.swap lastz.melGal1

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