src/hg/makeDb/doc/aplCal1.txt 1.4
1.4 2009/05/27 23:14:27 galt
more work on sea hare
Index: src/hg/makeDb/doc/aplCal1.txt
===================================================================
RCS file: /projects/compbio/cvsroot/kent/src/hg/makeDb/doc/aplCal1.txt,v
retrieving revision 1.3
retrieving revision 1.4
diff -b -B -U 1000000 -r1.3 -r1.4
--- src/hg/makeDb/doc/aplCal1.txt 14 May 2009 00:35:33 -0000 1.3
+++ src/hg/makeDb/doc/aplCal1.txt 27 May 2009 23:14:27 -0000 1.4
@@ -1,115 +1,244 @@
#########################################
#
# aplCal1 = Aplysia californica
#
#
# 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
wget -O /usr/local/apache/htdocs/images/Aplysia_californica.jpg \
'http://upload.wikimedia.org/wikipedia/commons/thumb/e/ef/Aplysia_californica.jpg/250px-Aplysia_californica.jpg'
# 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).
-ln -s /cluster/data/aplCal1/aplCal1.2bit /gbdb/aplCal1/
+############################################################################
+# 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
+