#########################################
#
# aplCal1 = Aplysia californica 
#
# (Galt)

# set up main genome directory

ssh hgwdev
cd /hive/data/genomes
mkdir aplCal1
cd aplCal1

mkdir download
cd download

# get sequence from BROAD
# http://www.broad.mit.edu/ftp/pub/assemblies/invertebrates/aplysia/

    cat << '_EOF_' > fetch.sh
#!/bin/sh

wget --timestamping -r -np -l 2 -nd -L
'ftp://ftp.broad.mit.edu/pub/assemblies/invertebrates/aplysia'

'_EOF_'
    # << happy emacs

chmod +x fetch.sh
./fetch.sh


# fixup ids for super fasta and quals
#  change ">scaffold_0.1-1784514 (Draft_v1)" 
#    to   ">scaffold_0"
#  sed 's/\(>scaffold_[0-9]*\)[.].*/\1/'

    cat << '_EOF_' > fix.csh
#!/bin/tcsh

gunzip -c assembly_supers.fasta.gz | sed 's/\(>scaffold_[0-9]*\)[.].*/\1/' \
  | gzip -c > assembly_supers.fasta.fix.gz

gunzip -c assembly_supers.qual.gz | sed 's/\(>scaffold_[0-9]*\)[.].*/\1/' \
  | gzip -c > assembly_supers.qual.fix.gz

'_EOF_'
    # << happy emacs

chmod +x fix.csh
./fix.csh

# totalKinds is something I made to add up the numbers of each kind
# it just reveals that there are too many unplaced fragments
# to consider making them all scaffolds, and as a chrUn it would
# be too large.  We are going to just leave these out of the 
# assembly like Broad did.  JK says ok.
totalKinds assembly.unplaced 2 > unplaced.types
cat unplaced.types
#unplaced 4270994
#low_quality 945330
#other 1781
#deliberate 273497
#vector_or_host 63932
#multiple_instance 27153

# back to the main directory
cd /hive/data/genomes/aplCal1

# Run automation to make the basic genome

    cat << '_EOF_' > aplCal1.config.ra
# Config parameters for makeGenomeDb.pl:
db aplCal1
scientificName Aplysia californica
commonName Sea Hare
assemblyDate Sept. 2008
assemblyLabel Broad Institute v. 1.0
orderKey 825
clade other
genomeCladePriority 19
mitoAcc NC_005827
fastaFiles /hive/data/genomes/aplCal1/download/assembly_supers.fasta.fix.gz
agpFiles /hive/data/genomes/aplCal1/download/assembly.agp
qualFiles /hive/data/genomes/aplCal1/download/assembly_supers.qual.fix.gz
dbDbSpeciesDir seaHare
taxId 6500
'_EOF_'
    # << happy emacs


time makeGenomeDb.pl aplCal1.config.ra > & makeGenomeDb.pl.out &
# took less than 10 minutes

featureBits -countGaps aplCal1 gap
#96577943 bases of 715806041 (13.492%) in intersection

cat chrom.sizes | awk '{sum+=$2;print sum,$0}'
#715806041 scaffold_8765 5001
# same total

# Organism Image
pushd $HOME/browser/images
wget -O Aplysia_californica.jpg \
 'http://upload.wikimedia.org/wikipedia/commons/thumb/e/ef/Aplysia_californica.jpg/250px-Aplysia_californica.jpg'
#  check this .jpg file into the browser doc source tree images directory
cvs add -kb Aplysia_californica.jpg 
cvs commit -m 'organism image for Aplysia' Aplysia_californica.jpg 
cp -p Aplysia_californica.jpg  /usr/local/apache/htdocs/images/
popd

# Edit and check-in templates for description.html, gold.html, gap.html, aplCal1/trackDb.ra

# repeat mask
time doRepeatMasker.pl aplCal1 > & doRepeatMasker.pl.out &

cat bed/RepeatMasker.2009-05-13/faSize.rmsk.txt 
#715806041 bases (96577949 N's 619228092 real 559834485 upper 59393607 lower)
#in 8767 sequences in 1 files

[hgwdev:aplCal1> featureBits -countGaps aplCal1 rmsk
59403223 bases of 715806041 (8.299%) in intersection


# simple repeat masker trf
time doSimpleRepeat.pl aplCal1 > & doSimpleRepeat.pl.out &

[hgwdev:aplCal1> featureBits -countGaps aplCal1 simpleRepeat
#55409917 bases of 715806041 (7.741%) in intersection


# make final masked .2bit
[hgwdev:aplCal1> twoBitMask aplCal1.rmsk.2bit -add bed/simpleRepeat.2009-05-13/trfMask.bed aplCal1.2bit
#Warning: BED file bed/simpleRepeat.2009-05-13/trfMask.bed has >=13 fields
#which means it might contain block coordinates, but this program uses only the
#first three fields (the entire span -- no support for blocks).



############################################################################
#	prepare cluster data (DONE - 2009-05-19 - Galt)

ssh hgwdev
cd /hive/data/genomes/aplCal1

# create gbdb symlink
rm /gbdb/aplCal1/aplCal1.2bit
ln -s `pwd`/aplCal1.2bit /gbdb/aplCal1/

time blat aplCal1.2bit \
    /dev/null /dev/null -tileSize=11 -makeOoc=11.ooc -repMatch=1024
#Wrote 2216 overused 11-mers to 11.ooc
#29.646u 0.750s 0:33.12 91.7%    0+0k 0+0io 3pf+0w

mkdir /hive/data/staging/data/aplCal1
cp -p aplCal1.2bit /hive/data/staging/data/aplCal1
cp -p 11.ooc /hive/data/staging/data/aplCal1
cp -p chrom.sizes /hive/data/staging/data/aplCal1

# ask admin to sync this directory: /hive/data/staging/data/aplCal1/
#	to the kluster nodes /scratch/data/aplCal1/

############################################################################
# running cpgIsland (DONE - 2009-05-19 - Galt)

ssh hgwdev
cd /hive/data/genomes/aplCal1

mkdir bed/cpgIsland
cd bed/cpgIsland
cvs -d /projects/compbio/cvsroot checkout -P hg3rdParty/cpgIslands
cd hg3rdParty/cpgIslands
# comment out the following two lines if it compiles cleanly
# some day  (there were some other fixups too, adding include lines)
sed -i -e "s#\(extern char\* malloc\)#// \1#" cpg_lh.c
make
#warning: incompatible implicit declaration of built-in function
# ignore the warnings
cd ../../ 
ln -s hg3rdParty/cpgIslands/cpglh.exe

# make hardmasked fasta files
mkdir -p hardMaskedFa
bash
cut -f1 ../../chrom.sizes | while read C
do
echo ${C}
twoBitToFa ../../aplCal1.2bit:$C stdout \
    | maskOutFa stdin hard hardMaskedFa/${C}.fa
done

#exit bash 
exit

cut -f1 ../../chrom.sizes > chr.list
cat << '_EOF_' > template
#LOOP
./runOne $(root1) {check out line results/$(root1).cpg}
#ENDLOOP
'_EOF_'
# << happy emacs

cat << '_EOF_' > runOne
#!/bin/csh -fe
./cpglh.exe hardMaskedFa/$1.fa > /scratch/tmp/$1.$$
mv /scratch/tmp/$1.$$ $2
'_EOF_'
# << happy emacs

mkdir results
chmod +x runOne

pk
cd /hive/data/genomes/aplCal1/bed/cpgIsland

gensub2 chr.list single template jobList
para create jobList
para try
para check
para push
para time
#Completed: 8767 of 8767 jobs
#CPU time in finished jobs:         47s       0.79m     0.01h    0.00d  0.000 y
#IO & Wait Time:                 43996s     733.26m    12.22h    0.51d  0.001 y
#Average job time:                   5s       0.08m     0.00h    0.00d
#Longest finished job:              26s       0.43m     0.01h    0.00d
#Submission to last job:          1161s      19.35m     0.32h    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

# took around 15 minutes

cd /hive/data/genomes/aplCal1/bed/cpgIsland
hgLoadBed aplCal1 cpgIslandExt -tab \
  -sqlTable=$HOME/kent/src/hg/lib/cpgIslandExt.sql cpgIsland.bed

#Reading cpgIsland.bed
#Loaded 23247 elements of size 10
#Sorted
#Creating table definition for cpgIslandExt
#Saving bed.tab
#Loading aplCal1


#############################################################
#
# Since there is no chrUn and no huge gaps,
# MarkD and Hiram say we can skip the unbridged-lift step
#

############################################################################
# AUTO UPDATE GENBANK RUN  (DONE - 2009-05-27,13 - Galt)
    # align with latest genbank process.

    # first time species step
    cd ~/kent/src/hg/makeDb/genbank/src/lib
    vi gbGenome.c
    # add aplCal
    cd ~/kent/src/hg/makeDb/genbank
    make server-install

    cd ~/kent/src/hg/makeDb/genbank
    cvsup
    # edit etc/genbank.conf to add aplCal1 

# aplCal1 - (Sea Hare Broad WGS Aplcal 2.0)
#       Assembly Accession: AASC00000000.2
aplCal1.serverGenome = /hive/data/genomes/aplCal1/aplCal1.2bit
aplCal1.clusterGenome = /scratch/data/aplCal1/aplCal1.2bit
aplCal1.ooc = /scratch/data/aplCal1/11.ooc
aplCal1.refseq.mrna.native.pslCDnaFilter  = ${lowCover.refseq.mrna.native.pslCDnaFilter}
aplCal1.refseq.mrna.xeno.pslCDnaFilter    = ${lowCover.refseq.mrna.xeno.pslCDnaFilter}
aplCal1.genbank.mrna.native.pslCDnaFilter = ${lowCover.genbank.mrna.native.pslCDnaFilter}
aplCal1.genbank.mrna.xeno.pslCDnaFilter   = ${lowCover.genbank.mrna.xeno.pslCDnaFilter}
aplCal1.genbank.est.native.pslCDnaFilter = ${lowCover.genbank.est.native.pslCDnaFilter}
aplCal1.genbank.est.xeno.pslCDnaFilter   = ${lowCover.genbank.est.xeno.pslCDnaFilter}
aplCal1.genbank.est.native.load = yes
aplCal1.genbank.est.xeno.load = no
aplCal1.refseq.mrna.xeno.load  = yes
aplCal1.refseq.mrna.xeno.loadDesc = yes
aplCal1.perChromTables = no
aplCal1.downloadDir = aplCal1
aplCal1.genbank.mrna.blatTargetDb = yes
aplCal1.lift = no

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

    ssh genbank
    screen		#	use a screen to manage this job
    cd /cluster/data/genbank
    time nice +19 bin/gbAlignStep -initial aplCal1 &
    #   real 3:21:30.93
    #   var/build/logs/2009.05.27-20:11:50.aplCal1.initalign.log

    # load database when finished
    ssh hgwdev
    screen	# use screen to manage this long running command
    cd /cluster/data/genbank
    time nice +19 ./bin/gbDbLoadStep -drop -initialLoad aplCal1 &
    # logFile: var/dbload/hgwdev/logs/2009.05.27-23:42:23.dbload.log
    # real   23:26.57

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

# markd 2010-01-19: turns out there are no refseqs, so disable it.
 aplCal1.refseq.mrna.native.load  = no

#########################################################################
#  BLATSERVERS ENTRY (DONE - 2009-05-28 - Galt)
#	After getting a blat server assigned by the Blat Server Gods,
    ssh hgwdev

    hgsql -e 'INSERT INTO blatServers (db, host, port, isTrans, canPcr) \
	VALUES ("aplCal1", "blat13", "17822", "1", "0"); \
	INSERT INTO blatServers (db, host, port, isTrans, canPcr) \
	VALUES ("aplCal1", "blat13", "17823", "0", "1");' \
	    hgcentraltest
    #	test it with some sequence

############################################################################
# Making download files (DONE - 2009-05-28 - Galt)


    cd /hive/data/genomes/aplCal1
    ln -s bed/RepeatMasker.2009-05-13/aplCal1.fa.out .
    cd bed
    ln -s simpleRepeat.2009-05-13 simpleRepeat
    cd ..

    makeDownloads.pl aplCal1 >& downloads.log


###########################################################################
# HUMAN (hg18) PROTEINS TRACK (DONE 2009-06-08 braney )
    # bash  if not using bash shell already

    cd /cluster/data/aplCal1
    mkdir /cluster/data/aplCal1/blastDb

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

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


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

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

# 36727/(50000/656) = 481.858240

   mkdir -p kgfa
   split -l 482 /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
#  77   77 1001 kg.lst

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

    para time
# Completed: 50512 of 50512 jobs
# CPU time in finished jobs:    4656584s   77609.74m  1293.50h   53.90d  0.148 y
# IO & Wait Time:                281459s    4690.98m    78.18h    3.26d  0.009 y
# Average job time:                  98s       1.63m     0.03h    0.00d
# Longest finished job:             409s       6.82m     0.11h    0.00d
# Submission to last job:         10255s     170.92m     2.85h    0.12d

    ssh swarm
    cd /cluster/data/aplCal1/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=12000 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: 77 of 77 jobs
# CPU time in finished jobs:       2940s      49.01m     0.82h    0.03d  0.000 y
# IO & Wait Time:                  5988s      99.79m     1.66h    0.07d  0.000 y
# Average job time:                 116s       1.93m     0.03h    0.00d
# Longest finished job:             297s       4.95m     0.08h    0.00d
# Submission to last job:           455s       7.58m     0.13h    0.01d

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

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

    # check coverage
    featureBits aplCal1 blastHg18KG 
# 6922807 bases of 619228098 (1.118%) in intersection

    featureBits aplCal1 blastHg18KG xenoRefGene  -enrichment
# blastHg18KG 1.118%, xenoRefGene 0.845%, both 0.490%, cover 43.82%, enrich 51.85x

    rm -rf blastOut
#end tblastn


############################################################################
# Creating the pushQ entry (DONE - 2009-06-08 - Galt)
    mkdir /hive/data/genomes/aplCal1/pushQ
    cd /hive/data/genomes/aplCal1/pushQ
    makePushQSql.pl aplCal1 > aplCal1.pushQ.sql
    #	aplCal1 does not have seq
    #	aplCal1 does not have extFile
    ssh hgwbeta hgsql -h hgsqlbeta qapushq < /hive/data/genomes/aplCal1/pushQ/aplCal1.pushQ.sql

    # Then, look at file(s) named in each of the following
    # wiggle tables and added these to the files section of qapushq.aplCal1
    # queue table:
          hgsql aplCal1 -e 'select distinct(file) from gc5Base'
          hgsql aplCal1 -e 'select distinct(file) from quality'

    # added to all.joiner both $gbd and $chainDest

############################################################################
# Creating the chain/net trackDb.ra entry (DONE - 2010-02-04 - Galt)
# To Do #1813 remove aplCal1 <-> hg19 chain/nets (2011-10-07 Chin)
# However, only tables and download files are physically drop/removed
# All data created stay on hive.

    cd /hive/data/genomes/hg19/bed
    ln -s lastzAplCal1.2009-06-08/ lastz.aplCal1
    findScores.pl hg19 aplCal1
#looking in file:
#  /hive/data/genomes/hg19/bed/lastz.aplCal1/axtChain/run/chain.csh
#-scoreScheme=/scratch/data/blastz/HoxD55.q
#matrix 16 91,-90,-25,-100,-90,100,-100,-25,-25,-100,100,-90,-100,-25,-90,91
#-minScore=5000
#-linearGap=loose

    cd kent/src/hg/makeDb/trackDb
    vi trackDb.chainNet.ra 
#cloned fr2 setting,  substituted aplCal1 for fr2 and AplCal1 for Fr2
# checked that these settings match the above output
#matrix 16 91,-90,-25,-100,-90,100,-100,-25,-25,-100,100,-90,-100,-25,-90,91
#chainMinScore 5000 
#chainLinearGap loose

    make DBS=hg19
    cvs commit -m 'adding aplCal1 chain/net' trackDb.chainNet.ra

