src/hg/makeDb/doc/loxAfr3.txt 1.2
1.2 2009/07/23 00:53:28 braney
start human proteins track
Index: src/hg/makeDb/doc/loxAfr3.txt
===================================================================
RCS file: /projects/compbio/cvsroot/kent/src/hg/makeDb/doc/loxAfr3.txt,v
retrieving revision 1.1
retrieving revision 1.2
diff -b -B -U 1000000 -r1.1 -r1.2
--- src/hg/makeDb/doc/loxAfr3.txt 21 Jul 2009 20:42:58 -0000 1.1
+++ src/hg/makeDb/doc/loxAfr3.txt 23 Jul 2009 00:53:28 -0000 1.2
@@ -1,143 +1,306 @@
# for emacs: -*- mode: sh; -*-
# This file describes how we made the elephant browser database on
# Broad Institute loxAfr3 (NCBI project 12569, AAGU03000000)
# "$Id$";
#############################################################################
# Download sequence (DONE - 2009-07-15 - Hiram)
mkdir -p /hive/data/genomes/loxAfr3/broad
/hive/data/genomes/loxAfr3/broad
wget --timestamping \
"ftp://ftp.broadinstitute.org/pub/assemblies/mammals/elephant/loxAfr3/*"
# lift quality scores to scaffold coordinates
qaToQac assembly.quals.gz assembly.quals.qac
qacAgpLift assembly.agp assembly.quals.qac loxAfr3.quals.qac
#############################################################################
# Elephant loxAfr3 browser initialization (DONE - 2009-07-15 - Hiram)
cd /hive/data/genomes/loxAfr3
cat << '_EOF_' > loxAft3.config.ra
# Config parameters for makeGenomeDb.pl:
db loxAfr3
clade mammal
genomeCladePriority 35
scientificName Loxodonta africana
commonName Elephant
assemblyDate Jul. 2009
assemblyLabel Broad Institute loxAfr3 (NCBI project 12569, AAGU03000000)
orderKey 340
mitoAcc NC_000934
fastaFiles /hive/data/genomes/loxAfr3/broad/assembly.bases.gz
agpFiles /hive/data/genomes/loxAfr3/broad/assembly.agp
qualFiles /hive/data/genomes/loxAfr3/broad/loxAfr3.quals.qac
dbDbSpeciesDir elephant
taxId 9785
'_EOF_'
# << happy emacs
# run stepwise to verify each step
makeGenomeDb.pl -stop=seq loxAfr3.config.ra > seq.log 2>&1
makeGenomeDb.pl -continue=agp -stop=agp loxAfr3.config.ra > agp.log 2>&1
makeGenomeDb.pl -continue=db -stop=db loxAfr3.config.ra > db.log 2>&1
makeGenomeDb.pl -continue=dbDb -stop=dbDb loxAfr3.config.ra > dbDb.log 2>&1
makeGenomeDb.pl -continue=trackDb loxAfr3.config.ra > trackDb.log 2>&1
#############################################################################
# loxAfr3 repeatMasker (DONE - 2009-07-15 - Hiram)
mkdir /hive/data/genomes/loxAfr3/bed/repeatMasker
cd /hive/data/genomes/loxAfr3/bed/repeatMasker
doRepeatMasker.pl -buildDir=`pwd` loxAfr3 > do.log 2>&1
# about 6 hours
cat faSize.rmsk.txt
# 3196760833 bases (78195493 N's 3118565340 real 1633809371 upper
# 1484755969 lower) in 2353 sequences in 1 files
#############################################################################
# loxAfr3 simpleRepeat (DONE - 2009-07-15 - Hiram)
mkdir /hive/data/genomes/loxAfr3/bed/simpleRepeat
cd /hive/data/genomes/loxAfr3/bed/simpleRepeat
time doSimpleRepeat.pl -buildDir=`pwd` loxAfr3 > do.log 2>&1
# real 23m42.536s
cat fb.simpleRepeat
# 27746420 bases of 3118565340 (0.890%) in intersection
# add to RM after done above:
twoBitMask bed/repeatMasker/loxAfr3.clean.2bit \
-add bed/simpleRepeat/trfMask.bed loxAfr3.2bit
twoBitToFa loxAfr3.2bit stdout | faSize stdin > faSize.loxAfr3.2bit.txt
#############################################################################
# create ooc file and populate /scratch/data (DONE - 2009-07-16 - Hiram)
# repMatch = 1024 * sizeof(loxAfr3)/sizeof(hg19)
# 1102 = 1024 * (3118565340/2897310462)
time blat loxAfr3.2bit \
/dev/null /dev/null -tileSize=11 -makeOoc=jkStuff/loxAfr3.11.ooc \
-repMatch=1100
# Wrote 41026 overused 11-mers to jkStuff/loxAfr3.11.ooc
mkdir /hive/data/staging/data/loxAfr3
cp -p loxAfr3.2bit /hive/data/staging/data/loxAfr3
cp -p jkStuff/loxAfr3.11.ooc /hive/data/staging/data/loxAfr3
cp -p chrom.sizes /hive/data/staging/data/loxAfr3
# request push to kluster nodes
##########################################################################
## GENBANK alignments (DONE - 2007-08-03 - Hiram)
cd $HOME/kent/src/hg/makeDb/genbank/etc
cvs up
# edit genbank.conf and add the following entry just above loxAfr1:
# loxAfr3 (elephant)
loxAfr3.serverGenome = /hive/data/genomes/loxAfr3/loxAfr3.2bit
loxAfr3.clusterGenome = /scratch/data/loxAfr3/loxAfr3.2bit
loxAfr3.ooc = /scratch/data/loxAfr3/loxAfr3.11.ooc
loxAfr3.lift = no
loxAfr3.refseq.mrna.native.pslCDnaFilter = ${lowCover.refseq.mrna.native.pslCDnaFilter}
loxAfr3.refseq.mrna.xeno.pslCDnaFilter = ${lowCover.refseq.mrna.xeno.pslCDnaFilter}
loxAfr3.genbank.mrna.native.pslCDnaFilter = ${lowCover.genbank.mrna.native.pslCDnaFilter}
loxAfr3.genbank.mrna.xeno.pslCDnaFilter = ${lowCover.genbank.mrna.xeno.pslCDnaFilter}
loxAfr3.genbank.est.native.pslCDnaFilter = ${lowCover.genbank.est.native.pslCDnaFilter}
loxAfr3.refseq.mrna.native.load = yes
loxAfr3.refseq.mrna.xeno.load = yes
loxAfr3.genbank.mrna.xeno.load = yes
loxAfr3.genbank.est.native.load = no
loxAfr3.downloadDir = loxAfr3
loxAfr3.perChromTables = no
# after commiting that edit, install thusly:
cd $HOME/kent/src/hg/makeDb/genbank
make etc-update
ssh genbank
screen # use a screen to manage this long lived job
cd /cluster/data/genbank
time nice -n +19 bin/gbAlignStep -initial loxAfr3 &
## logFile: var/build/logs/2009.07.21-10:14:29.loxAfr3.initalign.log
# real 174m1.258s
ssh hgwdev
cd /cluster/data/genbank
time nice -n +19 ./bin/gbDbLoadStep -drop -initialLoad loxAfr3
# var/dbload/hgwdev/logs/2009.07.21-13:22:25.dbload.log
# real 14m0.774s
featureBits loxAfr3 xenoMrna
# 65367813 bases of 3118565340 (2.096%) in intersection
featureBits loxAfr3 xenoRefGene
# 50364043 bases of 3118565340 (1.615%) in intersection
featureBits loxAfr3 all_mrna
# 15550 bases of 3118565340 (0.000%) in intersection
# enable daily alignment and update of hgwdev (DONE - 2009-07-21 - Hiram)
cd ~/kent/src/hg/makeDb/genbank
cvsup
# add mm9 to:
etc/align.dbs
etc/hgwdev.dbs
cvs ci -m "Added loxAfr3 - Loxodonta africana" etc/align.dbs etc/hgwdev.dbs
make etc-update
##########################################################################
+# HUMAN (hg18) PROTEINS TRACK (started 2009-07-22 braney )
+ # bash if not using bash shell already
+
+ cd /cluster/data/loxAfr3
+ mkdir /cluster/data/loxAfr3/blastDb
+
+ awk '{if ($2 > 1000000) print $1}' chrom.sizes > 1meg.lst
+ twoBitToFa -seqList=1meg.lst loxAfr3.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 loxAfr3.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
+# 3934
+
+ mkdir -p /cluster/data/loxAfr3/bed/tblastn.hg18KG
+ cd /cluster/data/loxAfr3/bed/tblastn.hg18KG
+ echo ../../blastDb/*.nsq | xargs ls -S | sed "s/\.nsq//" > query.lst
+ wc -l query.lst
+# 3934 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/3934) = 577.936072
+
+ mkdir -p kgfa
+ split -l 578 /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
+# 64 64 832 kg.lst
+
+ mkdir -p blastOut
+ for i in `cat kg.lst`; do mkdir blastOut/`basename $i .fa`; done
+ tcsh
+ cd /cluster/data/loxAfr3/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/loxAfr3/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/loxAfr3/bed/tblastn.hg18KG
+ gensub2 query.lst kg.lst blastGsub blastSpec
+ para create blastSpec
+# para try, check, push, check etc.
+
+ para time
+# Completed: 251158 of 251158 jobs
+# CPU time in finished jobs: 8587545s 143125.76m 2385.43h 99.39d 0.272 y
+# IO & Wait Time: 885281s 14754.68m 245.91h 10.25d 0.028 y
+# Average job time: 38s 0.63m 0.01h 0.00d
+# Longest finished job: 1167s 19.45m 0.32h 0.01d
+# Submission to last job: 9852s 164.20m 2.74h 0.11d
+
+ ssh swarm
+ cd /cluster/data/loxAfr3/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: 166 of 166 jobs
+# CPU time in finished jobs: 107801s 1796.69m 29.94h 1.25d 0.003 y
+# IO & Wait Time: 23945s 399.08m 6.65h 0.28d 0.001 y
+# Average job time: 794s 13.23m 0.22h 0.01d
+# Longest finished job: 25056s 417.60m 6.96h 0.29d
+# Submission to last job: 25085s 418.08m 6.97h 0.29d
+
+
+ cd /cluster/data/loxAfr3/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: 43305 failed: 0 errors: 0
+
+ # load table
+ ssh hgwdev
+ cd /cluster/data/loxAfr3/bed/tblastn.hg18KG
+ hgLoadPsl loxAfr3 blastHg18KG.psl
+
+ # check coverage
+ featureBits loxAfr3 blastHg18KG
+# 19115900 bases of 1506896106 (1.269%) in intersection
+
+ featureBits loxAfr3 blastHg18KG refGene -enrichment
+# blastHg18KG 1.269%, refGene 1.827%, both 0.659%, cover 51.97%, enrich 28.45x
+
+ rm -rf blastOut
+#end tblastn