######################################### # # 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