src/hg/makeDb/doc/susScr2.txt 1.3
1.3 2010/03/31 21:13:19 hiram
susScr2 done and in the pushQ
Index: src/hg/makeDb/doc/susScr2.txt
===================================================================
RCS file: /projects/compbio/cvsroot/kent/src/hg/makeDb/doc/susScr2.txt,v
retrieving revision 1.2
retrieving revision 1.3
diff -b -B -U 4 -r1.2 -r1.3
--- src/hg/makeDb/doc/susScr2.txt 29 Mar 2010 16:11:05 -0000 1.2
+++ src/hg/makeDb/doc/susScr2.txt 31 Mar 2010 21:13:19 -0000 1.3
@@ -5,9 +5,9 @@
# Sus scrofa - SGSC Sscrofa9.2 NCBI project 10718, CM000812
# ftp://ftp.ncbi.nlm.nih.gov:genbank/genomes/Eukaryotes/vertebrates_mammals/Sus_scrofa/Sscrofa9.2/
##########################################################################
-# Download sequence (DONE - 2010-03-03 Chin)
+# Download sequence (DONE - 2010-03-25 - Hiram)
mkdir /hive/data/genomes/susScr2
cd /hive/data/genomes/susScr2
mkdir genbank
cd genbank
@@ -47,9 +47,9 @@
# 2262484801 bases (31203023 N's 2231281778 real 2231281778 upper
# 0 lower) in 19 sequences in 19 files
#########################################################################
-# Initial makeGenomeDb.pl (DONE - 2009-11-06 - Hiram)
+# Initial makeGenomeDb.pl (DONE - 2010-03-25 - Hiram)
cd /hive/data/genomes/susScr2
cat << '_EOF_' > susScr2.config.ra
# Config parameters for makeGenomeDb.pl:
db susScr2
@@ -164,9 +164,9 @@
# 196728
# == 100179 + 96549
########################################################################
-# Create kluster run files (DONE - 2009-11-09 - Hiram)
+# Create kluster run files (DONE - 2010-03-26 - Hiram)
# numerator is susScr2 gapless bases "real" as reported by:
# featureBits -noRandom -noHap susScr2 gap
# denominator is hg19 gapless bases as reported by:
# featureBits -noRandom -noHap hg19 gap
@@ -193,9 +193,9 @@
cp -p jkStuff/nonBridged.lft \
/hive/data/staging/data/susScr2/susScr2.nonBridged.lft
########################################################################
-# GENBANK AUTO UPDATE (DONE - 2009-11-09 - Hiram)
+# GENBANK AUTO UPDATE (DONE - 2010-03-26 - Hiram)
ssh hgwdev
cd ~/kent/src/hg/makeDb/genbank
cvsup
@@ -230,26 +230,29 @@
ssh hgwdev
cd /cluster/data/genbank
time nice -n +19 ./bin/gbDbLoadStep -drop -initialLoad susScr2
# logFile: var/dbload/hgwdev/logs/2010.03.26-15:38:17.dbload.log
-XXX - running Fri Mar 26 15:38:45 PDT 2010
- # real 31m29.282s
+ # real 68m
# enable daily alignment and update of hgwdev
cd ~/kent/src/hg/makeDb/genbank
cvsup
# add susScr2 to:
etc/align.dbs
etc/hgwdev.dbs
- cvs ci -m "Added susScr2 - Tetraodon Nigirividis" \
- etc/align.dbs etc/hgwdev.dbs
+ cvs ci -m "Added susScr2 - Pig" etc/align.dbs etc/hgwdev.dbs
make etc-update
+ # DONE 2010-03-31
#########################################################################
# reset position to RHO location as found from blat of hg19 RHO gene
+# (DONE - 2010-03-31 - Hiram)
hgsql -e \
'update dbDb set defaultPos="chr13:57394166-57402412" where name="susScr2";' \
hgcentraltest
+ # and make this the default genome for Pig
+ hgsql -e 'update defaultDb set name="susScr2" where name="susScr1";' \
+ hgcentraltest
############################################################################
# ctgPos2 track - showing clone sequence locations on chromosomes
# (DONE - 2010-03-26 - Hiram)
@@ -295,9 +298,9 @@
hgLoadSqlTab susScr2 ctgPos2 $HOME/kent/src/hg/lib/ctgPos2.sql ctgPos2.tab
############################################################################
-# susScr2 Pig BLASTZ/CHAIN/NET (WORKING - 2010-03-27 - Hiram)
+# susScr2 Pig BLASTZ/CHAIN/NET (DONE - 2010-03-27 - Hiram)
screen # use a screen to manage this multi-day job
mkdir /hive/data/genomes/susScr2/bed/lastzBosTau4.2010-03-27
cd /hive/data/genomes/susScr2/bed/lastzBosTau4.2010-03-27
@@ -337,9 +340,114 @@
netChainSubset -verbose=0 noClass.net susScr2.bosTau4.all.chain.gz stdout \
| chainStitchId stdin stdout | gzip -c > susScr2.bosTau4.over.chain.gz
- # real 2377m21.473s
+ # and, finish the rest of netChains.csh manually, the netToAxt step
+ # and the axtToMaf step, log is in axtChain/finiChains.log
+
+ # after done with netChains - continuing with load:
+ time nice -n +19 doBlastzChainNet.pl -verbose=2 \
+ `pwd`/DEF \
+ -continue=load -noLoadChainSplit -syntenicNet \
+ -workhorse=hgwdev -smallClusterHub=encodek -bigClusterHub=swarm \
+ -chainMinScore=3000 -chainLinearGap=medium > load.log 2>&1 &
+ # this didn't work either due to memory limits with hgLoadChain
+ # add the following to loadUp.csh
+# limit at 160 Gb
+limit datasize 163840m
+limit vmemoryuse 163840m
+ # and finish it manually (7h39m)
+ cat fb.susScr2.chainBosTau4Link.txt
+ # 1478903080 bases of 2231298548 (66.280%) in intersection
+ # then continuing:
+ time nice -n +19 doBlastzChainNet.pl -verbose=2 \
+ `pwd`/DEF \
+ -continue=download -noLoadChainSplit -syntenicNet \
+ -workhorse=hgwdev -smallClusterHub=encodek -bigClusterHub=swarm \
+ -chainMinScore=3000 -chainLinearGap=medium > download.log 2>&1 &
+XXX - running Tue Mar 30 13:18:09 PDT 2010
+ # creating a bigWig graph to see the chain pileups:
+ cd /hive/data/genomes/susScr2/bed/lastzBosTau4.2010-03-27/axtChain
+ zcat susScr2.bosTau4.all.chain.gz | grep "^chain " \
+ | awk '{printf "%s\t%d\t%d\t%s\t%s\t%s\n", $3, $6, $7, $8, $2, $5}' \
+ > all.bed
+ # find the largest score:
+ ave -col=5 all.bed
+Q1 5300.000000
+median 8265.000000
+Q3 14079.000000
+average 12729.401671
+min 3000.000000
+max 1506828099.000000
+count 144157201
+total 1835034915327.000000
+standard deviation 419128.061072
+ # normalize the scores to 0-1000:
+ awk '
+{printf "%s\t%d\t%d\t%s\t%d\t%s\n", $1,$2,$3,$4,(1000*$5)/1506828099, $6}' \
+ all.bed | sort -k1,1 -k2,2n > all.chain.bed
+ bedGraphToBigWig all.chain.overlap.wigVar ../../../chrom.sizes all.chain.bw
+ bigWigInfo all.chain.bw
+version: 3
+isCompressed: yes
+isSwapped: 0
+primaryDataSize: 188,237,695
+primaryIndexSize: 1,267,372
+zoomLevels: 10
+chromCount: 20
+basesCovered: 2,255,615,700
+mean: 28.797958
+min: 1.000000
+max: 14674.000000
+std: 202.526527
+ ln -s `pwd`/all.chain.bw /gbdb/susScr2/bbi/bosTau4ChainPileUp.bw
+
+ hgsql susScr2 -e 'drop table if exists bosTau4ChainPileUp; \
+ create table bosTau4ChainPileUp (fileName varchar(255) not null); \
+ insert into bosTau4ChainPileUp values
+ ("/gbdb/susScr2/bbi/bosTau4ChainPileUp.bw");'
+
+ # and the swap
+ mkdir /hive/data/genomes/bosTau4/bed/blastz.susScr2.swap
+ cd /hive/data/genomes/bosTau4/bed/blastz.susScr2.swap
+ time nice -n +19 doBlastzChainNet.pl -verbose=2 \
+ /hive/data/genomes/susScr2/bed/lastzBosTau4.2010-03-27/DEF \
+ -swap -noLoadChainSplit -syntenicNet \
+ -workhorse=hgwdev -smallClusterHub=encodek -bigClusterHub=pk \
+ -chainMinScore=3000 -chainLinearGap=medium > swap.log 2>&1 &
+ # most interesting, this failed on the first step chainSwap
+ # but the failure didn't make it stop, it continued and produced
+ # zero results to the end. Running manually:
+ cd /hive/data/genomes/bosTau4/bed/blastz.susScr2.swap/axtChain
+export sizeG=188743680
+ulimit -d $sizeG
+ulimit -v $sizeG
+
+chainSwap /hive/data/genomes/susScr2/bed/lastzBosTau4.2010-03-27/axtChain/susScr2.bosTau4.all.chain.gz stdout \
+| nice chainSort stdin stdout | nice gzip -c > bosTau4.susScr2.all.chain.gz
+
+ # it also runs out of memory in the lift over file creation:
+ export sizeG=188743680
+ ulimit -d $sizeG
+ ulimit -v $sizeG
+
+ netChainSubset -verbose=0 noClass.net bosTau4.susScr2.all.chain.gz stdout \
+ | chainStitchId stdin stdout | gzip -c > bosTau4.susScr2.over.chain.gz
+
+ # and netChains.csh is finished manually with this added:
+# memory limit 160 Gb
+limit datasize 163840m
+limit vmemoryuse 163840m
+
+ # manually run the loadUp.csh with this added:
+# memory limit 160 Gb
+limit datasize 163840m
+limit vmemoryuse 163840m
+ # real 498m5.861s
+
+ cat fb.bosTau4.chainSusScr2Link.txt
+ # 1383557633 bases of 2731830700 (50.646%) in intersection
+
#########################################################################
# SWAP mm9 lastz (DONE - 2010-03-27 - Hiram)
# original alignment
cd /hive/data/genomes/mm9/bed/lastzSusScr2.2010-03-26
@@ -395,5 +503,116 @@
cat fb.susScr2.chainMonDom5Link.txt
# 182834643 bases of 2231298548 (8.194%) in intersection
+############################################################################
+# running cpgIsland business (DONE - 2010-03-31 - Hiram)
+ mkdir /hive/data/genomes/susScr2/bed/cpgIsland
+ cd /hive/data/genomes/susScr2/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 ../../susScr2.2bit:$C stdout \
+ | maskOutFa stdin hard hardMaskedFa/${C}.fa
+done
+
+ ssh swarm
+ cd /hive/data/genomes/susScr2/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
+
+ gensub2 chr.list single template jobList
+ para create jobList
+ para try
+ para check ... etc
+ para time
+# Completed: 20 of 20 jobs
+# CPU time in finished jobs: 168s 2.79m 0.05h 0.00d 0.000 y
+# IO & Wait Time: 86s 1.44m 0.02h 0.00d 0.000 y
+# Average job time: 13s 0.21m 0.00h 0.00d
+# Longest finished job: 29s 0.48m 0.01h 0.00d
+# Submission to last job: 31s 0.52m 0.01h 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/susScr2/bed/cpgIsland
+ hgLoadBed susScr2 cpgIslandExt -tab \
+ -sqlTable=$HOME/kent/src/hg/lib/cpgIslandExt.sql cpgIsland.bed
+ # Loaded 38778 elements of size 10
+
+ # cleanup
+ rm -fr hardMaskedFa
+
#########################################################################
+# all.joiner update, downloads and in pushQ - (DONE - 2010-03-31 - Hiram)
+ cd $HOME/kent/src/hg/makeDb/schema
+ # fixup all.joiner until this is a clean output
+ joinerCheck -database=susScr2 -all all.joiner
+
+ mkdir /hive/data/genomes/susScr2/goldenPath
+ cd /hive/data/genomes/susScr2/goldenPath
+ makeDownloads.pl susScr2 > do.log 2>&1
+
+ # now ready for pushQ entry
+ mkdir /hive/data/genomes/susScr2/pushQ
+ cd /hive/data/genomes/susScr2/pushQ
+ makePushQSql.pl susScr2 > susScr2.pushQ.sql 2> stderr.out
+ # check for errors in stderr.out, some are OK, e.g.:
+# WARNING: susScr2 does not have seq
+# WARNING: susScr2 does not have extFile
+
+# WARNING: Could not tell (from trackDb, all.joiner and hardcoded lists of
+# supporting and genbank tables) which tracks to assign these tables to:
+# bosTau4ChainPileUp
+
+ # copy it to hgwbeta
+ scp -p susScr2.pushQ.sql hgwbeta:/tmp
+ ssh hgwbeta
+ cd /tmp
+ hgsql qapushq < susScr2.pushQ.sql
+ # in that pushQ entry walk through each entry and see if the
+ # sizes will set properly
+
+############################################################################