06d7be056190c14b85e71bc12523f18ea6815b5e markd Mon Dec 7 00:50:29 2020 -0800 BLAT mmap index support merge with master diff --git src/hg/makeDb/doc/calJac3.txt src/hg/makeDb/doc/calJac3.txt index e676988..a891718 100644 --- src/hg/makeDb/doc/calJac3.txt +++ src/hg/makeDb/doc/calJac3.txt @@ -1,3305 +1,3326 @@ # for emacs: -*- mode: sh; -*- # $Id: calJac3.txt,v 1.16 2010/06/11 16:00:57 chinhli Exp $ # Marmoset sequence: # ftp.ncbi.nlm.nih.gov/genbank/genomes/Eukaryotes/vertebrates_mammals/ # Callithrix_jacchus/Callithrix_jacchus-3.2 # Callithrix jacchus # http://www.hgsc.bcm.tmc.edu/project-species-p-Marmoset.hgsc?pageLocation=Marmoset # http://www.ncbi.nlm.nih.gov/Traces/wgs/?val=ACFV00 ########################################################################## # Download sequence (DONE - 2010-02-04 - Hiram) mkdir /hive/data/genomes/calJac3 cd /hive/data/genomes/calJac3 mkdir genbank cd genbank wget --timestamping -r --cut-dirs=6 --level=0 -nH -x \ --no-remove-listing -np \ "ftp://ftp.ncbi.nlm.nih.gov/genbank/genomes/Eukaryotes/vertebrates_mammals/Callithrix_jacchus/Callithrix_jacchus-3.2/*" mkdir ucscChr cd ucscChr # fixup the accession names to become UCSC chrom names zcat ../Primary_Assembly/unplaced_scaffolds/FASTA/unplaced.scaf.fa.gz \ | sed -e "s/^>.*gb|\([A-Z]*[0-9]*\).1.*/>chrUn_\1/" > chrUn.fa zcat ../Primary_Assembly/unplaced_scaffolds/AGP/unplaced.scaf.agp.gz \ | sed -e "s/^\([A-Z]*[0-9]*\).1/chrUn_\1/" > chrUn.agp find ../Primary_Assembly/unlocalized_scaffolds/FASTA -type f \ | grep "unlocalized.scaf.fa.gz$" | head -1 | while read F do C=`basename ${F}` C=${C/.unlocalized.scaf.fa.gz} zcat "${F}" | sed -e "s/^>.*gb|\([A-Z]*[0-9]*\).1.*/>${C}_\1_random/" done > chr_randoms.fa find ../Primary_Assembly/unlocalized_scaffolds/AGP -type f | grep ".agp.gz$" \ | while read F do C=`basename ${F}` echo -n "${C} " C=${C/.unlocalized.scaf.agp.gz} echo "${C}" zcat "${F}" | sed -e "s/^\([A-Z]*[0-9]*\).1/${C}_\1_random/" done > chr_randoms.agp AC="../Primary_Assembly/assembled_chromosomes" for F in ${AC}/FASTA/chr*.fa.gz do C=`basename $F` C=${C/.fa.gz} echo -n "${C} " H=`zcat "${F}" | head -1` chrN=`echo $H | sed -e "s/.*Callithrix jacchus chromosome /chr/; s/, .*//"` A=`echo $H | sed -e "s/. Callithrix.*//; s/.*gb.//"` echo $chrN $A zcat ${AC}/AGP/${chrN}.comp.agp.gz \ | sed -e "s/^${A}/${chrN}/" > ${chrN}.agp echo ">${chrN}" > ${chrN}.fa zcat ${AC}/FASTA/${chrN}.fa.gz | grep -v "^>" >> ${chrN}.fa done ########################################################################## # Initial genome build (DONE - 2009-12-17 - Hiram) cd /hive/data/genomes/calJac3 cat << '_EOF_' > calJac3.config.ra # Config parameters for makeGenomeDb.pl: db calJac3 clade mammal genomeCladePriority 16 scientificName Callithrix jacchus commonName Marmoset assemblyDate Mar. 2009 assemblyLabel WUSTL 3.2 (GCA_000004665.1) orderKey 40 mitoAcc none fastaFiles /hive/data/genomes/calJac3/genbank/ucscChr/*.fa agpFiles /hive/data/genomes/calJac3/genbank/ucscChr/*.agp # qualFiles none dbDbSpeciesDir marmoset taxId 9483 '_EOF_' makeGenomeDb.pl -workhorse=hgwdev -stop=seq calJac3.config.ra > seq.out 2>&1 # real 4m5.924s makeGenomeDb.pl -continue=agp -stop=agp calJac3.config.ra > agp.out 2>&1 # real 0m20.968s makeGenomeDb.pl -continue=db -stop=db calJac3.config.ra > db.out 2>&1 # real 5m39.181s # XXX - chromInfo doesn't have large enough fields for the name keys # been fixed in later versions of makeGenomeDb.pl makeGenomeDb.pl -continue=dbDb -stop=dbDb calJac3.config.ra > dbDb.out 2>&1 makeGenomeDb.pl -continue=trackDb -stop=trackDb calJac3.config.ra > trackDb.out 2>&1 ########################################################################## # running repeat masker (DONE - 2010-02-02 - Hiram) mkdir /hive/data/genomes/calJac3/bed/repeatMasker cd /hive/data/genomes/calJac3/bed/repeatMasker time doRepeatMasker.pl -buildDir=`pwd` -noSplit \ -bigClusterHub=swarm -dbHost=hgwdev -workhorse=hgwdev \ -smallClusterHub=memk calJac3 > do.log 2>&1 & # real 443m55.891s cat faSize.rmsk.txt # 2914958544 bases (162452744 N's 2752505800 real 1439949519 upper 1312556281 # lower) in 14205 sequences in 1 files # %45.03 masked total, %47.69 masked real ########################################################################## # running simple repeat (DONE - 2010-02-02 - Hiram) mkdir /hive/data/genomes/calJac3/bed/simpleRepeat cd /hive/data/genomes/calJac3/bed/simpleRepeat time doSimpleRepeat.pl -buildDir=`pwd` -bigClusterHub=swarm \ -dbHost=hgwdev -workhorse=hgwdev -smallClusterHub=memk \ calJac3 > do.log 2>&1 & # real 198m33.953s cat fb.simpleRepeat # 66211508 bases of 2752505800 (2.405%) in intersection cd /hive/data/genomes/calJac3 twoBitMask calJac3.rmsk.2bit \ -add bed/simpleRepeat/trfMask.bed calJac3.2bit # you can safely ignore the warning about fields >= 13 twoBitToFa calJac3.2bit stdout | faSize stdin > faSize.calJac3.2bit.txt cat faSize.calJac3.2bit.txt # 2914958544 bases (162452744 N's 2752505800 real 1439244378 upper # 1313261422 lower) in 14205 sequences in 1 files # %45.05 masked total, %47.71 masked real rm /gbdb/calJac3/calJac3.2bit ln -s `pwd`/calJac3.2bit /gbdb/calJac3/calJac3.2bit ########################################################################## ## WINDOWMASKER (DONE- 2013-06-21 - Hiram) mkdir /hive/data/genomes/calJac3/bed/windowMasker cd /hive/data/genomes/calJac3/bed/windowMasker time nice -n +19 doWindowMasker.pl -buildDir=`pwd` -workhorse=hgwdev \ -dbHost=hgwdev calJac3 > do.log 2>&1 & # real 267m27.325s # Masking statistics faSize.calJac3.cleanWMSdust.txt # 2914958544 bases (162452744 N's 2752505800 real 1739687951 # upper 1012817849 lower) in 14205 sequences in 1 files # Total size: mean 205206.5 sd 5236162.4 # min 1031 (chr17_ACFV01172182_random) max 210400635 (chr1) median 4713 # %34.75 masked total, %36.80 masked real # how much does this window masker and repeat masker overlap: featureBits -countGaps calJac3 rmsk windowmaskerSdust \ > fb.calJac3.rmsk.windowmaskerSdust.txt 2>&1 # 752340288 bases of 2914958544 (25.810%) in intersection ######################################################################### # MAKE 11.OOC FILES FOR BLAT (DONE - 2010-02-03 - Hiram) ssh kolossus # numerator is calJac3 gapless bases "real" as reported by faSize # denominator is hg17 gapless bases as reported by featureBits, # 1024 is threshold used for human -repMatch: calc \( 2752505800 / 2897310462 \) \* 1024 # ( 2752505800 / 2897310462 ) * 1024 = 972.821510 # ==> use -repMatch=950 according to size scaled down from 1024 for human. # and rounded down to nearest 50 cd /hive/data/genomes/calJac3 blat calJac3.2bit /dev/null /dev/null -tileSize=11 \ -makeOoc=jkStuff/calJac3.11.ooc -repMatch=950 # Wrote 32908 overused 11-mers to jkStuff/calJac3.11.ooc mkdir /hive/data/staging/data/calJac3 cp -p calJac3.2bit chrom.sizes jkStuff/calJac3.11.ooc \ /hive/data/staging/data/calJac3 gapToLift -bedFile=jkStuff/nonBridgedGaps.bed calJac3 \ jkStuff/calJac3.nonBridged.lft ########################################################################## # BLATSERVERS ENTRY (DONE - 2009-12-23 - Hiram) # After getting a blat server assigned by the Blat Server Gods, ssh hgwdev hgsql -e 'INSERT INTO blatServers (db, host, port, isTrans, canPcr) \ VALUES ("calJac3", "blat12", "17804", "1", "0"); \ INSERT INTO blatServers (db, host, port, isTrans, canPcr) \ VALUES ("calJac3", "blat12", "17805", "0", "1");' \ hgcentraltest # test it with some sequence ############################################################################ # reset position to RHO location as found from blat of hg19 RHO gene hgsql -e \ 'update dbDb set defaultPos="chr15:67439781-67471107" where name="calJac3";' \ hgcentraltest ############################################################################ # genbank run (DONE - 2010-02-04 - Hiram) ssh hgwdev cd $HOME/kent/src/hg/makeDb/genbank # edit etc/genbank.conf to add this section just before calJac1: # Marmoset calJac3.serverGenome = /hive/data/genomes/calJac3/calJac3.2bit calJac3.clusterGenome = /scratch/data/calJac3/calJac3.2bit calJac3.ooc = /scratch/data/calJac3/calJac3.11.ooc calJac3.lift = no calJac3.perChromTables = no calJac3.refseq.mrna.native.pslCDnaFilter = ${ordered.refseq.mrna.native.pslCDnaFilter} calJac3.refseq.mrna.xeno.pslCDnaFilter = ${ordered.refseq.mrna.xeno.pslCDnaFilter} calJac3.genbank.mrna.native.pslCDnaFilter = ${ordered.genbank.mrna.native.pslCDnaFilter} calJac3.genbank.mrna.xeno.pslCDnaFilter = ${ordered.genbank.mrna.xeno.pslCDnaFilter} calJac3.genbank.est.native.pslCDnaFilter = ${ordered.genbank.est.native.pslCDnaFilter} calJac3.genbank.est.xeno.pslCDnaFilter = ${ordered.genbank.est.xeno.pslCDnaFilter} calJac3.downloadDir = calJac3 calJac3.refseq.mrna.native.load = yes calJac3.refseq.mrna.xeno.load = yes calJac3.refseq.mrna.xeno.loadDesc = yes cvs ci -m "adding marmoset calJac3" etc/genbank.conf make etc-update ssh genbank screen # control this business with a screen since it takes a while cd /cluster/data/genbank time nice -n +19 bin/gbAlignStep -initial calJac3 & # var/build/logs/2010.02.10-09:19:20.calJac3.initalign.log # real 238m59.390s ssh hgwdev cd /cluster/data/genbank time ./bin/gbDbLoadStep -drop -initialLoad calJac3 & # logFile: var/dbload/hgwdev/logs/2010.02.11-10:41:22.dbload.log # real 25m55.474s # enable daily alignment and update of hgwdev cd ~/kent/src/hg/makeDb/genbank cvsup # add calJac3 to: etc/align.dbs etc/hgwdev.dbs cvs ci -m "Adding calJac3 - Marmoset - Callithrix jacchus" \ etc/align.dbs etc/hgwdev.dbs make etc-update # done - 2010-02-11 - Hiram ############################################################################ # running cpgIsland business (DONE - 2010-02-11 - Hiram) mkdir /hive/data/genomes/calJac3/bed/cpgIsland cd /hive/data/genomes/calJac3/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 ../../calJac3.2bit:$C stdout \ | maskOutFa stdin hard hardMaskedFa/${C}.fa done ssh swarm cd /hive/data/genomes/calJac3/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: 14205 of 14205 jobs # CPU time in finished jobs: 205s 3.41m 0.06h 0.00d 0.000 y # IO & Wait Time: 38701s 645.02m 10.75h 0.45d 0.001 y # Average job time: 3s 0.05m 0.00h 0.00d # Longest finished job: 23s 0.38m 0.01h 0.00d # Submission to last job: 202s 3.37m 0.06h 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 cd /hive/data/genomes/calJac3/bed/cpgIsland hgLoadBed calJac3 cpgIslandExt -tab \ -sqlTable=$HOME/kent/src/hg/lib/cpgIslandExt.sql cpgIsland.bed # Reading cpgIsland.bed # Loaded 32732 elements of size 10 # Sorted # Saving bed.tab # Loading calJac3 # cleanup rm -fr hardMaskedFa ############################################################################ # LASTZ Human Swap (DONE - 2010-02-11 - Hiram) # original alignment on hg19: cd /hive/data/genomes/hg19/bed/lastzCalJac3.2010-02-11 cat fb.hg19.chainCalJac3Link.txt # 2047068864 bases of 2897316137 (70.654%) in intersection # and for this swap mkdir /hive/data/genomes/calJac3/bed/blastz.hg19.swap cd /hive/data/genomes/calJac3/bed/blastz.hg19.swap time nice -n +19 doBlastzChainNet.pl -verbose=2 \ /hive/data/genomes/hg19/bed/lastzCalJac3.2010-02-11/DEF \ -swap -syntenicNet \ -workhorse=hgwdev -smallClusterHub=memk -bigClusterHub=swarm \ -chainMinScore=5000 -chainLinearGap=medium > swap.log 2>&1 & # real 120m42.991s cat fb.calJac3.chainHg19Link.txt # 2030475813 bases of 2752505800 (73.768%) in intersection ##################################################################### # LASTZ Rhesus Swap (DONE - 2010-02-11 - Hiram) # original alignment to rheMac2 cd /hive/data/genomes/rheMac2/bed/lastzCalJac3.2010-02-11 cat fb.rheMac2.chainCalJac3Link.txt # 1871513554 bases of 2646704109 (70.711%) in intersection # and for this swap mkdir /hive/data/genomes/calJac3/bed/blastz.rheMac2.swap cd /hive/data/genomes/calJac3/bed/blastz.rheMac2.swap time nice -n +19 doBlastzChainNet.pl -verbose=2 \ /hive/data/genomes/rheMac2/bed/lastzCalJac3.2010-02-11/DEF \ -swap -syntenicNet \ -workhorse=hgwdev -smallClusterHub=memk -bigClusterHub=swarm \ -chainMinScore=5000 -chainLinearGap=medium > swap.log 2>&1 & # real 142m34.894s cat fb.calJac3.chainHg19Link.txt # 1916431926 bases of 2752505800 (69.625%) in intersection ############################################################################ # LASTZ Chimp Swap (DONE - 2010-02-11 - Hiram) # original alignment to panTro2 cd /hive/data/genomes/panTro2/bed/lastzCalJac3.2010-02-11 cat fb.panTro2.chainCalJac3Link.txt # 2016331285 bases of 2909485072 (69.302%) in intersection # and this swap run mkdir /hive/data/genomes/calJac3/bed/blastz.panTro2.swap cd /hive/data/genomes/calJac3/bed/blastz.panTro2.swap time nice -n +19 doBlastzChainNet.pl -verbose=2 \ /hive/data/genomes/panTro2/bed/lastzCalJac3.2010-02-11/DEF \ -swap -noLoadChainSplit -syntenicNet \ -workhorse=hgwdev -smallClusterHub=memk -bigClusterHub=swarm \ -chainMinScore=5000 -chainLinearGap=medium > swap.log 2>&1 & # real 118m42.203s cat fb.calJac3.chainHg19Link.txt # 1990168262 bases of 2752505800 (72.304%) in intersection ############################################################################ # LASTZ Orangutan Swap (DONE - 2010-02-11 - Hiram) # original alignment to ponAbe2 cd /hive/data/genomes/ponAbe2/bed/lastzCalJac3.2010-02-11 cat fb.ponAbe2.chainCalJac3Link.txt # 2086557592 bases of 3093572278 (67.448%) in intersection # and this swap run mkdir /hive/data/genomes/calJac3/bed/blastz.ponAbe2.swap cd /hive/data/genomes/calJac3/bed/blastz.ponAbe2.swap time nice -n +19 doBlastzChainNet.pl -verbose=2 \ /hive/data/genomes/ponAbe2/bed/lastzCalJac3.2010-02-11/DEF \ -swap -syntenicNet \ -workhorse=hgwdev -smallClusterHub=memk -bigClusterHub=swarm \ -chainMinScore=5000 -chainLinearGap=medium > swap.log 2>&1 & # real 146m12.301s cat fb.calJac3.chainHg19Link.txt # 1978857628 bases of 2752505800 (71.893%) in intersection ##################################################################### # lastz Dog CanFam2 Swap (DONE - 2010-02-14 - Hiram) # original alignment cd /hive/data/genomes/canFam2/bed/lastzCalJac3.2010-02-12 cat fb.canFam2.chainCalJac3Link.txt # 1363307334 bases of 2384996543 (57.162%) in intersection # and for the swap mkdir /hive/data/genomes/calJac3/bed/blastz.canFam2.swap cd /hive/data/genomes/calJac3/bed/blastz.canFam2.swap time nice -n +19 doBlastzChainNet.pl -verbose=2 \ /hive/data/genomes/canFam2/bed/lastzCalJac3.2010-02-12/DEF \ -swap -syntenicNet \ -workhorse=hgwdev -smallClusterHub=memk -bigClusterHub=swarm \ -chainMinScore=3000 -chainLinearGap=medium > swap.log 2>&1 & # real 129m56.144s cat fb.calJac3.chainHg19Link.txt # 1397333116 bases of 2752505800 (50.766%) in intersection ######################################################################### # lastz Mouse Mm9 Swap (DONE - 2010-02-14 - Hiram) # original alignment to mouse cd /hive/data/genomes/mm9/bed/lastzCalJac3.2010-02-12 cat fb.mm9.chainCalJac3Link.txt # 859869647 bases of 2620346127 (32.815%) in intersection # and for the swap mkdir /hive/data/genomes/calJac3/bed/blastz.mm9.swap cd /hive/data/genomes/calJac3/bed/blastz.mm9.swap time nice -n +19 doBlastzChainNet.pl -verbose=2 \ /hive/data/genomes/mm9/bed/lastzCalJac3.2010-02-12/DEF \ -swap -syntenicNet \ -workhorse=hgwdev -smallClusterHub=memk -bigClusterHub=swarm \ -chainMinScore=3000 -chainLinearGap=medium > swap.log 2>&1 & # real 90m38.739s cat fb.calJac3.chainHg19Link.txt # 861811978 bases of 2752505800 (31.310%) in intersection ##################################################################### # lastz Opossum monDom5 Swap (DONE - 2010-02-14 - Hiram) # original alignment to Opossum cd /hive/data/genomes/monDom5/bed/lastzCalJac3.2010-02-12 cat fb.monDom5.chainCalJac3Link.txt # 216197506 bases of 3501660299 (6.174%) in intersection # and for the swap mkdir /hive/data/genomes/calJac3/bed/blastz.monDom5.swap cd /hive/data/genomes/calJac3/bed/blastz.monDom5.swap time nice -n +19 doBlastzChainNet.pl -verbose=2 \ /hive/data/genomes/monDom5/bed/lastzCalJac3.2010-02-12/DEF \ -swap -noLoadChainSplit -syntenicNet \ -workhorse=hgwdev -smallClusterHub=memk -bigClusterHub=swarm \ -chainMinScore=5000 -chainLinearGap=loose > swap.log 2>&1 & # real 110m13.435s cat fb.calJac3.chainMonDom5Link.txt # 217614612 bases of 2752505800 (7.906%) in intersection ############################################################################## # HUMAN (hg18) PROTEINS TRACK (working braney...) # bash if not using bash shell already cd /cluster/data/calJac3 mkdir /cluster/data/calJac3/blastDb awk '{if ($2 > 1000000) print $1}' chrom.sizes > 1meg.lst twoBitToFa -seqList=1meg.lst calJac3.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 calJac3.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 # 3275 mkdir -p /cluster/data/calJac3/bed/tblastn.hg18KG cd /cluster/data/calJac3/bed/tblastn.hg18KG echo ../../blastDb/*.nsq | xargs ls -S | sed "s/\.nsq//" > query.lst wc -l query.lst # 3275 query.lst # we want around 350000 jobs calc `wc /cluster/data/hg18/bed/blat.hg18KG/hg18KG.psl | awk '{print $1}'`/\(350000/`wc query.lst | awk '{print $1}'`\) # 36727/(350000/3275) = 343.659786 mkdir -p kgfa split -l 344 /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 # 107 107 1391 kg.lst mkdir -p blastOut for i in `cat kg.lst`; do mkdir blastOut/`basename $i .fa`; done tcsh cd /cluster/data/calJac3/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/calJac3/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/calJac3/bed/tblastn.hg18KG gensub2 query.lst kg.lst blastGsub blastSpec para create blastSpec # para try, check, push, check etc. para time # Completed: 100076 of 100076 jobs # CPU time in finished jobs: 2219419s 36990.31m 616.51h 25.69d 0.070 y # IO & Wait Time: 445463s 7424.39m 123.74h 5.16d 0.014 y # Average job time: 27s 0.44m 0.01h 0.00d # Longest finished job: 319s 5.32m 0.09h 0.00d # Submission to last job: 2791s 46.52m 0.78h 0.03d ssh swarm cd /cluster/data/calJac3/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: 254 of 254 jobs # CPU time in finished jobs: 590077s 9834.62m 163.91h 6.83d 0.019 y # IO & Wait Time: 20749s 345.81m 5.76h 0.24d 0.001 y # Average job time: 2405s 40.08m 0.67h 0.03d # Longest finished job: 36592s 609.87m 10.16h 0.42d # Submission to last job: 36604s 610.07m 10.17h 0.42d cd /cluster/data/calJac3/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: 38517 failed: 0 errors: 0 # load table ssh hgwdev cd /cluster/data/calJac3/bed/tblastn.hg18KG hgLoadPsl calJac3 blastHg18KG.psl # check coverage featureBits calJac3 blastHg18KG # 18646943 bases of 332311746 (5.611%) in intersection featureBits calJac3 blastHg18KG ensGene -enrichment # blastHg18KG 5.611%, ensGene 9.520%, both 4.993%, cover 88.97%, enrich 9.35x rm -rf blastOut #end tblastn ############################################################################## # papHam1 Baboon LASTZ/CHAIN/NET (DONE - 2010-02-15 - Hiram) screen # use a screen to manage this multi-day job mkdir /hive/data/genomes/calJac3/bed/lastzPapHam1.2010-02-15 cd /hive/data/genomes/calJac3/bed/lastzPapHam1.2010-02-15 cat << '_EOF_' > DEF # baboon vs. marmoset # same paramters as human hg19 vs marmoset calJac3 BLASTZ=lastz # maximum M allowed with lastz is only 254 BLASTZ_M=254 BLASTZ_Q=/scratch/data/blastz/human_chimp.v2.q BLASTZ_O=600 BLASTZ_E=150 BLASTZ_K=4500 BLASTZ_Y=15000 BLASTZ_T=2 # TARGET: Marmoset (calJac3) SEQ1_DIR=/scratch/data/calJac3/calJac3.2bit SEQ1_LEN=/scratch/data/calJac3/chrom.sizes SEQ1_LIMIT=50 SEQ1_CHUNK=20000000 SEQ1_LAP=10000 # QUERY: Baboon papHam1 SEQ2_DIR=/scratch/data/papHam1/papHam1.2bit SEQ2_LEN=/scratch/data/papHam1/chrom.sizes SEQ2_CHUNK=20000000 SEQ2_LIMIT=300 SEQ2_LAP=0 BASE=/hive/data/genomes/calJac3/bed/lastzPapHam1.2010-02-15 TMPDIR=/scratch/tmp '_EOF_' # << this line keeps emacs coloring happy time nice -n +19 $HOME/kent/src/hg/utils/automation/doBlastzChainNet.pl \ `pwd`/DEF \ -verbose=2 -syntenicNet -chainMinScore=5000 -chainLinearGap=medium \ -workhorse=hgwdev -smallClusterHub=memk -bigClusterHub=swarm \ > do.log 2>&1 & # failed lastz run, finished manually # real 287m24.258s time nice -n +19 doBlastzChainNet.pl `pwd`/DEF \ -continue=cat \ -verbose=2 -syntenicNet -chainMinScore=5000 -chainLinearGap=medium \ -workhorse=hgwdev -smallClusterHub=memk -bigClusterHub=swarm \ > cat.log 2>&1 & # real 158m17.502s cat fb.calJac3.chainPapHam1Link.txt # 1928203329 bases of 2752505800 (70.053%) in intersection time doRecipBest.pl -workhorse=hgwdev -buildDir=`pwd` \ calJac3 papHam1 > rbest.log 2>&1 # real 232m mkdir /hive/data/genomes/papHam1/bed/blastz.calJac3.swap cd /hive/data/genomes/papHam1/bed/blastz.calJac3.swap time nice -n +19 doBlastzChainNet.pl -verbose=2 \ /hive/data/genomes/calJac3/bed/lastzPapHam1.2010-02-15/DEF \ -swap -syntenicNet \ -workhorse=hgwdev -smallClusterHub=pk -bigClusterHub=swarm \ -chainMinScore=5000 -chainLinearGap=medium > swap.log 2>&1 & # real 791m46.765s cat fb.papHam1.chainCalJac3Link.txt # 1908519637 bases of 2741867288 (69.607%) in intersection ############################################################################## # tarSyr1 Tarsier LASTZ/CHAIN/NET (DONE - 2010-02-21 - Hiram) screen # use a screen to manage this multi-day job mkdir /hive/data/genomes/calJac3/bed/lastzTarSyr1.2010-02-21 cd /hive/data/genomes/calJac3/bed/lastzTarSyr1.2010-02-21 cat << '_EOF_' > DEF # tarsier vs. marmoset # same paramters as human hg19 vs tarsier tarSyr1 BLASTZ=lastz # maximum M allowed with lastz is only 254 BLASTZ_M=254 # TARGET: Marmoset (calJac3) SEQ1_DIR=/scratch/data/calJac3/calJac3.2bit SEQ1_LEN=/scratch/data/calJac3/chrom.sizes SEQ1_LIMIT=50 SEQ1_CHUNK=20000000 SEQ1_LAP=10000 # QUERY: Tarsier tarSyr1 SEQ2_DIR=/scratch/data/tarSyr1/tarSyr1.2bit SEQ2_LEN=/scratch/data/tarSyr1/chrom.sizes SEQ2_CHUNK=20000000 SEQ2_LIMIT=300 SEQ2_LAP=0 BASE=/hive/data/genomes/calJac3/bed/lastzTarSyr1.2010-02-21 TMPDIR=/scratch/tmp TMPDIR=/scratch/tmp '_EOF_' # << this line keeps emacs coloring happy time nice -n +19 doBlastzChainNet.pl \ `pwd`/DEF \ -verbose=2 -syntenicNet -chainMinScore=3000 -chainLinearGap=medium \ -workhorse=hgwdev -smallClusterHub=pk -bigClusterHub=swarm \ > do.log 2>&1 & cat fb.calJac3.chainTarSyr1Link.txt # 1286219755 bases of 2752505800 (46.729%) in intersection time doRecipBest.pl -workhorse=hgwdev -buildDir=`pwd` \ calJac3 tarSyr1 > rbest.log 2>&1 & # real 532m mkdir /hive/data/genomes/tarSyr1/bed/blastz.calJac3.swap cd /hive/data/genomes/tarSyr1/bed/blastz.calJac3.swap time nice -n +19 doBlastzChainNet.pl -verbose=2 \ /hive/data/genomes/calJac3/bed/lastzTarSyr1.2010-02-21/DEF \ -swap -syntenicNet \ -workhorse=hgwdev -smallClusterHub=memk -bigClusterHub=swarm \ -chainMinScore=3000 -chainLinearGap=medium > swap.log 2>&1 & # real 120m42.991s cat fb.tarSyr1.chainCalJac3Link.txt # 2030475813 bases of 2752505800 (73.768%) in intersection ##################################################################### # micMur1 Mouse lemur LASTZ/CHAIN/NET (DONE - 2010-02-17,22 - Hiram) # Mouse lemur ( Microcebus murinus) screen # use a screen to manage this multi-day job mkdir /hive/data/genomes/calJac3/bed/lastzMicMur1.2010-02-17 cd /hive/data/genomes/calJac3/bed/lastzMicMur1.2010-02-17 cat << '_EOF_' > DEF # mouse lemur vs. marmoset # same paramters as human hg19 vs Mouse lemur micMur1 BLASTZ=lastz # maximum M allowed with lastz is only 254 BLASTZ_M=254 # TARGET: Marmoset (calJac3) SEQ1_DIR=/scratch/data/calJac3/calJac3.2bit SEQ1_LEN=/scratch/data/calJac3/chrom.sizes SEQ1_LIMIT=5 SEQ1_CHUNK=200000000 SEQ1_LAP=10000 # QUERY: Mouse lemur micMur1 SEQ2_DIR=/hive/data/genomes/micMur1/micMur1.2bit SEQ2_LEN=/hive/data/genomes/micMur1/chrom.sizes SEQ2_CHUNK=20000000 SEQ2_LIMIT=300 SEQ2_LAP=0 BASE=/hive/data/genomes/calJac3/bed/lastzMicMur1.2010-02-17 TMPDIR=/scratch/tmp '_EOF_' # << this line keeps emacs coloring happy time nice -n +19 doBlastzChainNet.pl -verbose=2 \ `pwd`/DEF \ -verbose=2 -syntenicNet -chainMinScore=3000 -chainLinearGap=medium \ -workhorse=hgwdev -smallClusterHub=memk -bigClusterHub=swarm \ > do.log 2>&1 & # real 5502m6.707s # some kluster difficulties, finished cat run manually, then: time nice -n +19 doBlastzChainNet.pl -verbose=2 \ -continue=chainRun `pwd`/DEF \ -verbose=2 -syntenicNet -chainMinScore=3000 -chainLinearGap=medium \ -workhorse=hgwdev -smallClusterHub=pk -bigClusterHub=swarm \ > chainRun.log 2>&1 & # real 374m19.587s calJac3 micMur1 02-17 cat fb.calJac3.chainMicMur1Link.txt # 1258616069 bases of 2752505800 (45.726%) in intersection time doRecipBest.pl -buildDir=`pwd` calJac3 micMur1 > rbest.log 2>&1 # real 235m55.179s calJac3 micMur1 mkdir /hive/data/genomes/micMur1/bed/blastz.calJac3.swap cd /hive/data/genomes/micMur1/bed/blastz.calJac3.swap time nice -n +19 doBlastzChainNet.pl -verbose=2 \ /hive/data/genomes/calJac3/bed/lastzMicMur1.2010-02-17/DEF \ -swap -syntenicNet \ -workhorse=hgwdev -smallClusterHub=memk -bigClusterHub=swarm \ -chainMinScore=3000 -chainLinearGap=medium > swap.log 2>&1 & # real 455m11.215s micMur1 calJac3 swap cat fb.micMur1.chainCalJac3Link.txt # 1243785262 bases of 1852394361 (67.145%) in intersection ##################################################################### # otoGar1 Bushbaby LASTZ/CHAIN/NET (DONE - 2010-02-17,22 - Hiram) screen # use a screen to manage this multi-day job mkdir /hive/data/genomes/calJac3/bed/lastzOtoGar1.2010-02-17 cd /hive/data/genomes/calJac3/bed/lastzOtoGar1.2010-02-17 cat << '_EOF_' > DEF # bushbaby vs. marmoset # same paramters as human hg19 vs Bushbaby otoGar1 BLASTZ=lastz # maximum M allowed with lastz is only 254 BLASTZ_M=254 # TARGET: Marmoset (calJac3) SEQ1_DIR=/scratch/data/calJac3/calJac3.2bit SEQ1_LEN=/scratch/data/calJac3/chrom.sizes SEQ1_LIMIT=5 SEQ1_CHUNK=200000000 SEQ1_LAP=10000 # QUERY: Bushbaby otoGar1 SEQ2_DIR=/scratch/data/otoGar1/otoGar1.rmsk.2bit SEQ2_LEN=/scratch/data/otoGar1/chrom.sizes SEQ2_CHUNK=20000000 SEQ2_LIMIT=300 SEQ2_LAP=0 BASE=/hive/data/genomes/calJac3/bed/lastzOtoGar1.2010-02-17 TMPDIR=/scratch/tmp '_EOF_' # << this line keeps emacs coloring happy time nice -n +19 doBlastzChainNet.pl -verbose=2 \ `pwd`/DEF \ -verbose=2 -syntenicNet -chainMinScore=3000 -chainLinearGap=medium \ -workhorse=hgwdev -smallClusterHub=memk -bigClusterHub=swarm \ > do.log 2>&1 & # real 4722m38.163s # memk failed at the cat run, finish it manually, then, continuing: time nice -n +19 doBlastzChainNet.pl -verbose=2 \ -continue=chainRun `pwd`/DEF \ -verbose=2 -syntenicNet -chainMinScore=3000 -chainLinearGap=medium \ -workhorse=hgwdev -smallClusterHub=pk -bigClusterHub=swarm \ > chainRun.log 2>&1 & # real 285m58.314s cat fb.calJac3.chainOtoGar1Link.txt # 1176505967 bases of 2752505800 (42.743%) in intersection time doRecipBest.pl -workhorse=hgwdev -buildDir=`pwd` \ calJac3 otoGar1 > rbest.log 2>&1 & # real 332m14.375s calJac3 otoGar1 mkdir /hive/data/genomes/otoGar1/bed/blastz.calJac3.swap cd /hive/data/genomes/otoGar1/bed/blastz.calJac3.swap time nice -n +19 doBlastzChainNet.pl -verbose=2 \ /hive/data/genomes/calJac3/bed/lastzOtoGar1.2010-02-17/DEF \ -swap -syntenicNet \ -workhorse=hgwdev -smallClusterHub=pk -bigClusterHub=swarm \ -chainMinScore=3000 -chainLinearGap=medium > swap.log 2>&1 & # real 310m4.077s cat fb.otoGar1.chainCalJac3Link.txt # 1158531484 bases of 1969052059 (58.837%) in intersection ##################################################################### ## 8-Way Multiz (DONE - 2009-06-09,2009-11-10 - Hiram) mkdir /hive/data/genomes/calJac3/bed/multiz8way cd /hive/data/genomes/calJac3/bed/multiz8way /cluster/bin/phast/tree_doctor \ --prune-all-but=calJac1,hg19,panTro2,rheMac2,ponAbe2,mm9,canFam2,monDom5 \ --rename="calJac1 -> calJac3 " \ /hive/data/genomes/hg19/bed/multiz46way/fixedTree/46wayFixed.nh > 8way.nh # *carefully* edit 8way.nh to get calJac3 at the top of this picture # resulting in this tree: (calJac3:0.066389,((rheMac2:0.057695,(ponAbe2:0.018342, (hg19:0.006591,panTro2:0.006639):0.012126):0.014256):0.010000, (mm9:0.352605,(canFam2:0.193569,monDom5:0.581923):0.020666) :0.088210):0.000001); # Use this specification in the phyloGif tool: # http://genome.ucsc.edu/cgi-bin/phyloGif # to obtain a gif image for htdocs/images/phylo/calJac3_8way.gif /cluster/bin/phast/all_dists 8way.nh > 8way.distances.txt # Use this output to create the table below, with this perl script: cat << '_EOF_' > sizeStats.pl #!/usr/bin/env perl use strict; use warnings; open (FH, "grep -y calJac3 8way.distances.txt | sort -k3,3n|") or die "can not read 8way.distances.txt"; my $count = 0; while (my $line = <FH>) { chomp $line; my ($calJac3, $D, $dist) = split('\s+', $line); my $chain = "chain" . ucfirst($D); my $B="/hive/data/genomes/calJac3/bed/lastz.$D/fb.calJac3." . $chain . "Link.txt"; my $chainLinkMeasure = `awk '{print \$5}' ${B} 2> /dev/null | sed -e "s/(//; s/)//"`; chomp $chainLinkMeasure; $chainLinkMeasure = 0.0 if (length($chainLinkMeasure) < 1); $chainLinkMeasure =~ s/\%//; my $swapFile="/hive/data/genomes/${D}/bed/lastz.calJac3/fb.${D}.chainCalJac3Link.txt"; my $swapMeasure = "N/A"; if ( -s $swapFile ) { $swapMeasure = `awk '{print \$5}' ${swapFile} 2> /dev/null | sed -e "s/(//; s/)//"`; chomp $swapMeasure; $swapMeasure = 0.0 if (length($swapMeasure) < 1); $swapMeasure =~ s/\%//; } my $orgName= `hgsql -N -e "select organism from dbDb where name='$D';" hgcentraltest`; chomp $orgName; if (length($orgName) < 1) { $orgName="N/A"; } ++$count; if ($swapMeasure eq "N/A") { printf "# %02d %.4f - %s %s\t(%% %.3f) (%s)\n", $count, $dist, $orgName, $D, $chainLinkMeasure, $swapMeasure } else { printf "# %02d %.4f - %s %s\t(%% %.3f) (%% %.3f)\n", $count, $dist, $orgName, $D, $chainLinkMeasure, $swapMeasure } } close (FH); '_EOF_' # << happy emacs chmod +x ./sizeStats.pl ./sizeStats.pl # # If you can fill in all the numbers in this table, you are ready for # the multiple alignment procedure # # featureBits chainLink measures # chainCalJac3Link # distance on calJac3 on other # 01 0.1090 - Orangutan ponAbe2 (% 71.893) (% 67.448) # 02 0.1094 - Human hg19 (% 73.768) (% 70.654) # 03 0.1094 - Chimp panTro2 (% 72.304) (% 69.302) # 04 0.1341 - Rhesus rheMac2 (% 69.625) (% 70.711) # 05 0.3688 - Dog canFam2 (% 50.766) (% 57.162) # 06 0.5072 - Mouse mm9 (% 31.310) (% 32.815) # 07 0.7572 - Opossum monDom5 (% 7.906) (% 6.174) # create species list and stripped down tree for autoMZ sed 's/[a-z][a-z]*_//g; s/:[0-9\.][0-9\.]*//g; s/;//; /^ *$/d' \ 8way.nh > tmp.nh echo `cat tmp.nh` > tree-commas.nh echo `cat tree-commas.nh` | sed 's/ //g; s/,/ /g' > tree.nh sed 's/[()]//g; s/,/ /g' tree.nh > species.list # bash shell syntax here ... mkdir -p mafLinks H="/hive/data/genomes/calJac3/bed/" for S in `sed -e "s/calJac3 //" species.list` do echo $S ls -og ${H}/lastz.${S}/axtChain/calJac3.${S}.synNet.maf.gz ln -s ${H}/lastz.${S}/axtChain/calJac3.${S}.synNet.maf.gz \ mafLinks/${S}.maf.gz done HERE=`pwd` export HERE PATH=${HERE}/penn:${PATH} export PATH rm -fr tmp mkdir -p tmp cd mafLinks time ../penn/autoMZ + T=${HERE}/tmp \ E=calJac3 "`cat ../tree.nh`" *.sing.maf result.maf # real 3584m8.094s mkdir /gbdb/calJac3/multiz8way ln -s `pwd`/mafLinks/result.maf /gbdb/calJac3/multiz8way/multiz8way.maf cd /scratch/tmp # Loaded 7475045 mafs in 1 files from /gbdb/calJac3/multiz8way time nice -n +19 hgLoadMaf calJac3 multiz8way time nice -n +19 hgLoadMafSummary -minSize=30000 -mergeGap=1500 \ -maxSize=200000 calJac3 multiz8waySummary multiz8way.maf | hgLoadMafSummary calJac1 -minSize=30000 -mergeGap=1500 \ -maxSize=200000 multiz9waySummary stdin # Created 1313222 summary blocks from 34128178 components # and 7475045 mafs from multiz8way.maf # real 8m36.016s ############################################################################## # gorGor2 Bushbaby LASTZ/CHAIN/NET (DONE - 2010-02-22,24 - Hiram) screen # use a screen to manage this multi-day job mkdir /hive/data/genomes/calJac3/bed/lastzOtoGar1.2010-02-22 cd /hive/data/genomes/calJac3/bed/lastzOtoGar1.2010-02-22 cat << '_EOF_' > DEF # Gorilla vs. marmoset # same paramters as human hg19 vs other nearby primates # without all the extra blastz parameters BLASTZ=lastz # maximum M allowed with lastz is only 254 BLASTZ_M=254 # TARGET: Marmoset (calJac3) SEQ1_DIR=/scratch/data/calJac3/calJac3.2bit SEQ1_LEN=/scratch/data/calJac3/chrom.sizes SEQ1_LIMIT=20 SEQ1_CHUNK=200000000 SEQ1_LAP=10000 # QUERY: Gorilla gorGor2 SEQ2_DIR=/scratch/data/gorGor2/gorGor2.2bit SEQ2_LEN=/scratch/data/gorGor2/chrom.sizes SEQ2_CHUNK=12000000 SEQ2_LIMIT=300 SEQ2_LAP=0 BASE=/hive/data/genomes/calJac3/bed/lastzGorGor2.2010-02-22 TMPDIR=/scratch/tmp '_EOF_' # << this line keeps emacs coloring happy time nice -n +19 doBlastzChainNet.pl -verbose=2 \ `pwd`/DEF \ -verbose=2 -syntenicNet -chainMinScore=3000 -chainLinearGap=medium \ -workhorse=hgwdev -smallClusterHub=memk -bigClusterHub=swarm \ > do.log 2>&1 & # Elapsed time: 1956m3.678s cat fb.calJac3.chainGorGor2Link.txt # 2101356280 bases of 2752505800 (76.343%) in intersection time doRecipBest.pl -workhorse=hgwdev -buildDir=`pwd` \ calJac3 gorGor2 > rbest.log 2>&1 # about 4h16m mkdir /hive/data/genomes/gorGor2/bed/blastz.calJac3.swap cd /hive/data/genomes/gorGor2/bed/blastz.calJac3.swap time nice -n +19 doBlastzChainNet.pl -verbose=2 \ /hive/data/genomes/calJac3/bed/lastzGorGor2.2010-02-22/DEF \ -swap -syntenicNet \ -workhorse=hgwdev -smallClusterHub=pk -bigClusterHub=swarm \ -chainMinScore=3000 -chainLinearGap=medium > swap.log 2>&1 & # real 250m57.089s cat fb.gorGor2.chainCalJac3Link.txt # 2135885920 bases of 2829687208 (75.481%) in intersection ##################################################################### ## 13-Way Multiz (DONE - 2010-02-23 - Hiram) ## Re-work 13-Way Multiz: (DONE 2011-04-19 - Chin) ## See redmine track #162 for more detail mv /hive/data/genomes/calJac3/bed/multiz13way \ /hive/data/genomes/calJac3/bed/multiz13way.2010-05-23 mkdir /hive/data/genomes/calJac3/bed/multiz13way cd /hive/data/genomes/calJac3/bed/multiz13way /cluster/bin/phast/tree_doctor \ --prune-all-but=calJac1,hg19,panTro2,rheMac2,ponAbe2,gorGor1,micMur1,otoGar1,papHam1,tarSyr1,mm9,canFam2,monDom5 \ --rename="calJac1 -> calJac3 ; gorGor1 -> gorGor2 " \ /hive/data/genomes/hg19/bed/multiz46way/fixedTree/46wayFixed.nh > 13way.nh # Note: # The rearrangement was not correct, the 13way.nh produced by tree_doctor can be used directly # for the multiz alignment without distance data. The distance will be generated (by 4d) # and used in the phastCon and phyloP steps. # rearrange calJac3 to the top, get some help from tree_doctor: /cluster/bin/phast/tree_doctor --name-ancestors --reroot calJac3 \ --with-branch 13way.nh # edit out the ancestors, and move calJac3 from the bottom to # the top, resulting in this tree: (calJac3:0.066389,(((((hg19:0.006591,panTro2:0.006639):0.002184, gorGor2:0.009411):0.009942,ponAbe2:0.018342):0.014256, (rheMac2:0.036199,papHam1:0.040000):0.021496):0.010000, ((((monDom5:0.581923,canFam2:0.193569):0.020666,mm9:0.352605):0.019992, (micMur1:0.091452,otoGar1:0.128984):0.035463):0.011307, tarSyr1:0.135169):0.056911):0.000001); # more rearranging after seeing what the distance table looks like # below to get them appearing as much as possible in their # distance order top to bottom: (calJac3:0.066389,(((ponAbe2:0.018342, ((hg19:0.006591,panTro2:0.006639):0.002184, gorGor2:0.009411):0.009942):0.014256, (rheMac2:0.036199,papHam1:0.040000):0.021496):0.010000, (tarSyr1:0.135169,((micMur1:0.091452,otoGar1:0.128984):0.035463, (mm9:0.352605, (canFam2:0.193569,monDom5:0.581923):0.020666):0.019992):0.011307):0.056911) :0.000001); # Use 13way.nh directly cat 13way.nh ((((((((((hg19:0.006591,panTro2:0.006639):0.002184,gorGor2:0.009411):0.009942,ponAbe2:0.018342):0.014256,(rheMac2:0.036199,papHam1:0.040000):0.021496):0.010000,calJac3:0.066389):0.056911,tarSyr1:0.135169):0.011307,(micMur1:0.091452,otoGar1:0.128984):0.035463):0.019992,mm9:0.352605):0.020666,canFam2:0.193569):0.156024,monDom5:0.425899); # Use this specification in the phyloGif tool after changing the names: /cluster/bin/phast/tree_doctor \ --rename="calJac3 -> Marmoset ; ponAbe2 -> Orangutan ; hg19 -> Human ; panTro2 -> Chimp ; gorGor2 -> Gorilla ; rheMac2 -> Rhesus ; papHam1 -> Baboon ; tarSyr1 -> Tarsier ; micMur1 -> Mouse_lemur ; otoGar1 -> Bushbaby ; canFam2 -> Dog ; mm9 -> Mouse ; monDom5 -> Opossum " 13way.nh # use the tree provided by Bob K. for display: cat << '_EOF_' > commonNames.13way.nh ((((((Marmoset,((((Human,Chimp),Gorilla),Orangutan), (Baboon,Rhesus))),Tarsier), (MouseLemur, Bushbaby)),Mouse),Dog),Opossum); '_EOF_' # << happy emacs # use http://genome.ucsc.edu/cgi-bin/phyloGif # to generate calJac3_13way.png image and copy over to # htdocs/images/phylo/ cp calJac3_13way.png /usr/local/apache/htdocs/images/phylo # Note the old a gif image for htdocs/images/phylo/calJac3_13way.gif # is wrong, remove it # rm /usr/local/apache/htdocs/images/phylo/calJac3_13way.gif /cluster/bin/phast/all_dists 13way.nh > 13way.distances.tmp # swap the calJac3 from field 2 to field 1 cat 13way.distances.tmp | \ awk '{if ($2 == "calJac3") print $2, " " $1, $3; else print $0;}' | \ sort > 13way.distances.txt # make sure all symlinks lastz.DB -> lastzDb-date # exist here and at the swap locations, the perl script expects this # in order to find featureBits numbers. # This is for our own bookkeeping business so we can # determine what is done. If there is any zero in the table, # then, we miss some steps such as no chain/net for a calJac3 to hg19. # Use 13way.distances.txt to create the table below # with this perl script: cat << '_EOF_' > sizeStats.pl #!/usr/bin/env perl use strict; use warnings; open (FH, "grep -y calJac3 13way.distances.txt | sort -k3,3n|") or die "can not read 13way.distances.txt"; my $count = 0; while (my $line = <FH>) { chomp $line; my ($calJac3, $D, $dist) = split('\s+', $line); my $chain = "chain" . ucfirst($D); my $B="/hive/data/genomes/calJac3/bed/lastz.$D/fb.calJac3." . $chain . "Link.txt"; my $chainLinkMeasure = `awk '{print \$5}' ${B} 2> /dev/null | sed -e "s/(//; s/)//"`; chomp $chainLinkMeasure; $chainLinkMeasure = 0.0 if (length($chainLinkMeasure) < 1); $chainLinkMeasure =~ s/\%//; my $swapFile="/hive/data/genomes/${D}/bed/lastz.calJac3/fb.${D}.chainCalJac3Link.txt"; my $swapMeasure = "N/A"; if ( -s $swapFile ) { $swapMeasure = `awk '{print \$5}' ${swapFile} 2> /dev/null | sed -e "s/(//; s/)//"`; chomp $swapMeasure; $swapMeasure = 0.0 if (length($swapMeasure) < 1); $swapMeasure =~ s/\%//; } my $orgName= `hgsql -N -e "select organism from dbDb where name='$D';" hgcentraltest`; chomp $orgName; if (length($orgName) < 1) { $orgName="N/A"; } ++$count; if ($swapMeasure eq "N/A") { printf "# %02d %.4f - %s %s\t(%% %.3f) (%s)\n", $count, $dist, $orgName, $D, $chainLinkMeasure, $swapMeasure } else { printf "# %02d %.4f - %s %s\t(%% %.3f) (%% %.3f)\n", $count, $dist, $orgName, $D, $chainLinkMeasure, $swapMeasure } } close (FH); '_EOF_' # << happy emacs chmod +x ./sizeStats.pl ./sizeStats.pl # 01 0.1090 - Orangutan ponAbe2 (% 71.893) (% 67.448) # 02 0.1094 - Human hg19 (% 73.768) (% 70.654) # 03 0.1094 - Chimp panTro2 (% 72.304) (% 69.302) # 04 0.1100 - Gorilla gorGor2 (% 76.343) (N/A) # 05 0.1341 - Rhesus rheMac2 (% 69.625) (% 70.711) # 06 0.1379 - Baboon papHam1 (% 70.053) (% 69.607) # 07 0.2585 - Tarsier tarSyr1 (% 46.729) (N/A) # 08 0.2615 - Mouse lemur micMur1 (% 45.726) (% 67.145) # 09 0.2991 - Bushbaby otoGar1 (% 42.743) (% 58.837) # 10 0.3688 - Dog canFam2 (% 50.766) (% 57.162) # 11 0.5072 - Mouse mm9 (% 31.310) (% 32.815) # 12 0.7572 - Opossum monDom5 (% 7.906) (% 6.174) # create species list and stripped down tree for autoMZ # Note: the 13way.nh can not be used directly becasue calJac3 is # not on top. Insted use commonNames.13way.nh provided by Bob # to create species list and stripped down tree. # sed 's/[a-z][a-z]*_//g; s/:[0-9\.][0-9\.]*//g; s/;//; /^ *$/d' \ # 13way.nh > tmp.nh sed -e 's/Marmoset/calJac3/' -e 's/Orangutan/ponAbe2/' \ -e 's/Human/hg19/' -e 's/Chimp/panTro2/' \ -e 's/Gorilla/gorGor2/' -e 's/Rhesus/rheMac2/' \ -e 's/Baboon/papHam1/' -e 's/Tarsier/tarSyr1/' \ -e 's/MouseLemur/micMur1/' -e 's/Bushbaby/otoGar1/' \ -e 's/Dog/canFam2/' -e 's/Mouse/mm9/' \ -e 's/Opossum/monDom5/' \ commonNames.13way.nh | sed 's/ //g; s/,/ /g; s/;//g;' > tree.nh sed 's/[()]//g; s/,/ /g; s/;//g;' tree.nh > species.list # create a tree of species name with comma fo 4d step later sed -e 's/Marmoset/calJac3/' -e 's/Orangutan/ponAbe2/' \ -e 's/Human/hg19/' -e 's/Chimp/panTro2/' \ -e 's/Gorilla/gorGor2/' -e 's/Rhesus/rheMac2/' \ -e 's/Baboon/papHam1/' -e 's/Tarsier/tarSyr1/' \ -e 's/MouseLemur/micMur1/' -e 's/Bushbaby/otoGar1/' \ -e 's/Dog/canFam2/' -e 's/Mouse/mm9/' \ -e 's/Opossum/monDom5/' \ commonNames.13way.nh | sed 's/ //g; s/;//g;' > tree-commas.nh # collect the single whole mafs into one place for splitting: mkdir singleMafs cd singleMafs ln -s ../../lastz.ponAbe2/axtChain/calJac3.ponAbe2.synNet.maf.gz . ln -s ../../lastz.hg19/axtChain/calJac3.hg19.synNet.maf.gz . ln -s ../../lastz.panTro2/axtChain/calJac3.panTro2.synNet.maf.gz . ln -s ../../lastz.gorGor2/axtChain/calJac3.gorGor2.synNet.maf.gz . ln -s ../../lastz.rheMac2/axtChain/calJac3.rheMac2.synNet.maf.gz . ln -s ../../lastz.papHam1/mafRBestNet/calJac3.papHam1.rbest.maf.gz . ln -s ../../lastz.tarSyr1/mafRBestNet/calJac3.tarSyr1.rbest.maf.gz . ln -s ../../lastz.micMur1/mafRBestNet/calJac3.micMur1.rbest.maf.gz . ln -s ../../lastz.otoGar1/mafRBestNet/calJac3.otoGar1.rbest.maf.gz . ln -s ../../lastz.mm9/axtChain/calJac3.mm9.synNet.maf.gz . ln -s ../../lastz.canFam2/axtChain/calJac3.canFam2.synNet.maf.gz . ln -s ../../lastz.monDom5/axtChain/calJac3.monDom5.synNet.maf.gz . mkdir /hive/data/genomes/calJac3/bed/multiz13way/splitMaf cd /hive/data/genomes/calJac3/bed/multiz13way/splitMaf for D in tarSyr1 papHam1 otoGar1 micMur1 do mkdir ${D} mafSplit -useHashedName=8 -byTarget /dev/null ${D}/ \ ../singleMafs/calJac3.${D}.rbest.maf.gz done for D in gorGor2 ponAbe2 hg19 panTro2 rheMac2 mm9 canFam2 monDom5 do mkdir ${D} mafSplit -useHashedName=8 -byTarget /dev/null ${D}/ \ ../singleMafs/calJac3.${D}.synNet.maf.gz done cd /hive/data/genomes/calJac3/bed/multiz13way mkdir run mkdir penn cp -p /cluster/bin/penn/multiz.2008-11-25/multiz penn cp -p /cluster/bin/penn/multiz.2008-11-25/maf_project penn cp -p /cluster/bin/penn/multiz.2008-11-25/autoMZ penn cd /hive/data/genomes/calJac3/bed/multiz13way/run # set the db and pairs directories here cat > autoMultiz.csh << '_EOF_' #!/bin/csh -ef set db = calJac3 set topDir = /hive/data/genomes/$db/bed/multiz13way set c = $1 set result = $2 set pennBin = $topDir/penn set run = `/bin/pwd` set tmp = /scratch/tmp/$db/multiz.$c set pairs = $topDir/splitMaf /bin/rm -fr $tmp /bin/mkdir -p $tmp /bin/cp -p $topDir/tree.nh $topDir/species.list $tmp pushd $tmp > /dev/null foreach s (`/bin/sed -e "s/^$db //" species.list`) set in = $pairs/$s/$c.maf set out = $db.$s.sing.maf if (-e $in.gz) then /bin/zcat $in.gz > $out if (! -s $out) then echo "##maf version=1 scoring=autoMZ" > $out endif else if (-e $in) then /bin/ln -s $in $out else echo "##maf version=1 scoring=autoMZ" > $out endif end set path = ($pennBin $path); rehash $pennBin/autoMZ + T=$tmp E=$db "`cat tree.nh`" $db.*.sing.maf $c.maf \ > /dev/null popd > /dev/null /bin/rm -f $result /bin/cp -p $tmp/$c.maf $result /bin/rm -fr $tmp /bin/rmdir --ignore-fail-on-non-empty /scratch/tmp/$db '_EOF_' # << happy emacs chmod +x autoMultiz.csh cat << '_EOF_' > template #LOOP ./autoMultiz.csh $(root1) {check out line+ /hive/data/genomes/calJac3/bed/multiz13way/run/maf/$(root1).maf} #ENDLOOP '_EOF_' # << happy emacs ssh swarm cd /hive/data/genomes/calJac3/bed/multiz13way/run mkdir maf find ../splitMaf -type f | grep "/[0-9][0-9][0-9].maf" \ | xargs -L 1 basename | sort -u > chr.part.list gensub2 chr.part.list single template jobList para -ram=8g create jobList # para try check time push stuff # put the split mafs back together into a single result # Completed: 256 of 256 jobs # CPU time in finished jobs: 312697s 5211.62m 86.86h 3.62d 0.010 y # IO & Wait Time: 23882s 398.03m 6.63h 0.28d 0.001 y # Average job time: 1315s 21.91m 0.37h 0.02d # Longest finished job: 25624s 427.07m 7.12h 0.30d # Submission to last job: 25750s 429.17m 7.15h 0.30d ssh hgwdev cd /hive/data/genomes/calJac3/bed/multiz13way/run head -q -n 1 maf/000.maf > calJac3.13way.maf for F in maf/*.maf do grep -h -v "^#" ${F} >> calJac3.13way.maf done tail -q -n 1 maf/000.maf >> calJac3.13way.maf # -rw-rw-r-- 1 30553880367 Apr 21 09:40 calJac3.13way.maf #split the /hive/data/genomes/calJac3/bed/multiz13way/run/calJac3.13way.maf # to other folder with fullName for phastCons # load into database ssh hgwdev mkdir -p /gbdb/calJac3/multiz13way/maf cd /hive/data/genomes/calJac3/bed/multiz13way/run ln -s `pwd`/calJac3.13way.maf \ /gbdb/calJac3/multiz13way/maf/multiz13way.maf # this generates an immense multiz13way.tab file in the directory # where it is running. Best to run this over in scratch. cd /scratch/tmp time nice -n +19 hgLoadMaf \ -pathPrefix=/gbdb/calJac3/multiz13way/maf calJac3 multiz13way # Indexing and tabulating /gbdb/calJac3/multiz13way/maf/multiz13way.maf # Loading multiz13way into database # Loaded 13258003 mafs in 1 files from /gbdb/calJac3/multiz13way/maf # real 16m48.119s # load summary table time nice -n +19 cat /gbdb/calJac3/multiz13way/maf/*.maf \ | hgLoadMafSummary calJac3 -minSize=30000 -verbose=2 \ -mergeGap=1500 -maxSize=200000 multiz13waySummary stdin # Indexing and tabulating stdin # Created 2327841 summary blocks from 99934514 components # and 13258003 mafs from stdin # Loading into calJac3 table multiz13waySummary... # Loading complete # real 16m7.795s wc -l multiz13way*.tab # 13258003 multiz13way.tab # 2327841 multiz13waySummary.tab # 15585844 total rm multiz13way*.tab ############################################################################ # GAP ANNOTATE MULTIZ13WAY MAF AND LOAD TABLES (DONE 2011-05-04 - Chin) # mafAddIRows has to be run on single chromosome maf files, it does not # function correctly when more than one reference sequence # are in a single file. # use the /hive/data/genomes/calJac3/bed/multiz13way/run/calJac3.13way.maf # to split cd /hive/data/genomes/calJac3/bed/multiz13way/ mkdir mafSplit cd mafSplit mafSplit -byTarget -useFullSequenceName /dev/null . ../run/calJac3.13way.maf # Splitting 1 files by target sequence -- ignoring first argument # /dev/null ls -lt | wc -l # 10348 # got 10348 mafs named after their chrom/scaff .maf # although there are over 14205 chroms and scaffolds (wc -l # chrom.sizes), # some are too small or have nothing aligning. mkdir /hive/data/genomes/calJac3/bed/multiz13way/anno cd /hive/data/genomes/calJac3/bed/multiz13way/anno # most of these will already exist from previous multiple # alignments # remove the echo from in front of the twoBitInfo command to get # them # to run if this loop appears to be correct for DB in `cat ../species.list` do CDIR="/hive/data/genomes/${DB}" if [ ! -f ${CDIR}/${DB}.N.bed ]; then echo "creating ${DB}.N.bed" echo twoBitInfo -nBed ${CDIR}/${DB}.2bit ${CDIR}/${DB}.N.bed else ls -og ${CDIR}/${DB}.N.bed fi done for DB in `sed -e "s/calJac3 //" ../species.list` do echo "${DB} " ln -s /hive/data/genomes/${DB}/${DB}.N.bed ${DB}.bed echo ${DB}.bed >> nBeds ln -s /hive/data/genomes/${DB}/chrom.sizes ${DB}.len echo ${DB}.len >> sizes done # make sure they all are successful symLinks: ls -ogrtL mkdir mafSplit for F in ../mafSplit/*.maf do B=`basename ${F}` mafAddIRows -nBeds=nBeds ${F} /hive/data/genomes/calJac3/calJac3.2bit \ mafSplit/${B} echo $B done # about 5h 20m # combine into one file head -q -n 1 mafSplit/chr1.maf > calJac3.13way.maf for F in mafSplit/*.maf do grep -h -v "^#" ${F} done >> calJac3.13way.maf # these maf files do not have the end marker, this does nothing: # tail -q -n 1 mafSplit/chr1.maf >> calJac3.13way.maf # How about an official end marker: echo "##eof maf" >> calJac3.13way.maf # Load into database # by loading this into the table multiz13waySummary, it will # replace # the previously loaded table with the unannotated mafs rm /gbdb/calJac3/multiz13way/multiz13way.maf # remove old symlink ln -s `pwd`/calJac3.13way.maf /gbdb/calJac3/multiz13way/multiz13way.maf cd /scratch/tmp time nice -n +19 hgLoadMaf calJac3 multiz13way # Indexing and tabulating /gbdb/calJac3/multiz13way/multiz13way.maf # Loading multiz13way into database # Loaded 14067989 mafs in 1 files from /gbdb/calJac3/multiz13way # real 17m32.053s time nice -n +19 hgLoadMafSummary -verbose=2 -minSize=30000 \ -mergeGap=1500 -maxSize=200000 calJac3 multiz13waySummary \ /gbdb/calJac3/multiz13way/multiz13way.maf # Indexing and tabulating /gbdb/calJac3/multiz13way/multiz13way.maf # Created 2327841 summary blocks from 99934514 components and # 14067989 mafs from /gbdb/calJac3/multiz13way/multiz13way.maf # Loading into calJac3 table multiz13waySummary... # Loading complete # real 21m4.436s # by loading this into the table multiz13waySummary, it will # replace # the previously loaded table with the unannotated mafs # remove the multiz13way*.tab files in this /scratch/tmp directory wc -l multiz13way*.tab # 14067989 multiz13way.tab rm multiz13way*.tab ##################################################################### # all.joiner update, downloads and in pushQ - (DONE - 2010-04-01 - Hiram) cd $HOME/kent/src/hg/makeDb/schema # fixup all.joiner until this is a clean output joinerCheck -database=calJac3 -all all.joiner mkdir /hive/data/genomes/calJac3/goldenPath cd /hive/data/genomes/calJac3/goldenPath time nice -n +19 makeDownloads.pl calJac3 > do.log 2>&1 # real 22m30.329s # now ready for pushQ entry mkdir /hive/data/genomes/calJac3/pushQ cd /hive/data/genomes/calJac3/pushQ time nice -n +19 makePushQSql.pl calJac3 > calJac3.pushQ.sql 2> stderr.out real 2m52.193s # check for errors in stderr.out, some are OK, e.g.: # WARNING: calJac3 does not have seq # WARNING: calJac3 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 calJac3.pushQ.sql hgwbeta:/tmp ssh hgwbeta cd /tmp hgsql qapushq < calJac3.pushQ.sql # in that pushQ entry walk through each entry and see if the # sizes will set properly ############################################################################ # ctgPos2 track - showing clone sequence locations on chromosomes # (DONE - 2010-04-01 - Hiram) mkdir /hive/data/genomes/calJac3/bed/ctgPos2 cd /hive/data/genomes/calJac3/bed/ctgPos2 cat << '_EOF_' > agpToCtgPos2.pl #!/usr/bin/env perl use warnings; use strict; my $argc = scalar(@ARGV); if ($argc != 1) { printf STDERR "usage: zcat your.files.agp.gz | agpToCtgPos2.pl /dev/stdin > ctgPos2.tab\n"; exit 255; } my $agpFile = shift; open (FH, "<$agpFile") or die "can not read $agpFile"; while (my $line = <FH>) { next if ($line =~ m/^#/); chomp $line; my @a = split('\s+', $line); next if ($a[4] =~ m/^N$/); my $chrSize = $a[2]-$a[1]+1; my $ctgSize = $a[7]-$a[6]+1; die "sizes differ $chrSize != $ctgSize\n$line\n" if ($chrSize != $ctgSize); printf "%s\t%d\t%s\t%d\t%d\t%s\n", $a[5], $chrSize, $a[0], $a[1]-1, $a[2], $a[4]; } close (FH); '_EOF_' # << happy emacs export S=../../genbank/Primary_Assembly/assembled_chromosomes/AGP cut -f2 ${S}/chr2acc | while read ACC do C=`grep "${ACC}" ${S}/chr2acc | cut -f1` zcat ${S}/AGP/chr${C}.agp.gz \ | sed -e "s/^${ACC}/chr${C}/" done | ./agpToCtgPos2.pl /dev/stdin > ctgPos2.tab hgLoadSqlTab calJac3 ctgPos2 $HOME/kent/src/hg/lib/ctgPos2.sql ctgPos2.tab ############################################################################ # N-SCAN gene predictions (nscanGene) - (2010-04-26 markd) # obtained NSCAN predictions from michael brent's group # at WUSTL cd /cluster/data/calJac3/bed/nscan wget -nv http://mblab.wustl.edu/predictions/marmoset/calJac3/readme.txt wget -nv http://mblab.wustl.edu/predictions/marmoset/calJac3/calJac3.prot.fa wget -nv http://mblab.wustl.edu/predictions/marmoset/calJac3/calJac3.gtf gzip calJac3.* chmod a-w * # load track gtfToGenePred -genePredExt calJac3.gtf.gz stdout | hgLoadGenePred -genePredExt calJac3 nscanGene stdin hgPepPred calJac3 generic nscanPep calJac3.prot.fa.gz rm *.tab # marmoset/calJac3/trackDb.ra, add: track nscanGene override informant Marmoset N-SCAN used mouse (mm9) as the informant for conservation. PASA clustered human ESTs were mapped to the marmoset genome with the use of the human/marmoset synteny chain track. These remapped clusters were then filtered for splice sites, and used as the EST track on marmoset. # veryify top-level search spec, should produce no results: hgsql -Ne 'select name from nscanGene' calJac3 | egrep -v '^chr[0-9a-zA-Z_]+\.([0-9]+|pasa)((\.[0-9a-z]+)?\.[0-9a-z]+)?$' |head ############################################################################ # Fixing references to old calJac1 date to be to March 2009 - (2010-05-12 mary) # trackDb cd /cluster/home/mary/kent/src/hg/makeDb/trackDb/marmoset find -name "*.html" | xargs grep 2007 # downloads cd /usr/local/apache/htdocs-hgdownload/goldenPath/calJac3 find -name "*" | xargs grep 2007 # hgcentral hgsql hgcentraltest update dbDb set description="March 2009 (WUGSC 3.2/calJac3)" where name="calJac3" ############################################################################ ## Annotate 13-way multiple alignment with gene annotations ## (DONE - 2010-05-17 - Chin) ## use ensGene for calJac3 (DONE - 2010-07-22 - Chin) ## Re-work (DONE 2011-05-05 - Chin) # Gene frames ## survey all genomes to see what type of gene track to use ssh hgwdev mkdir /hive/data/genomes/calJac3/bed/multiz13way/frames cd /hive/data/genomes/calJac3/bed/multiz13way/frames # # survey all the genomes to find out what kinds of gene tracks # they have cat << '_EOF_' > showGenes.csh #!/bin/csh -fe foreach db (`cat ../species.list`) echo -n "${db}: " set tables = `hgsql $db -N -e "show tables like '%Gene%'"` foreach table ($tables) if ($table == "ensGene" || $table == "refGene" || $table == "mgcGenes" || \ $table == "knownGene" || $table == "xenoRefGene" || $table == "nscanGene" ) then set count = `hgsql $db -N -e "select count(*) from $table"` echo -n "${table}: ${count}, " endif end set orgName = `hgsql hgcentraltest -N -e \ "select scientificName from dbDb where name='$db'"` set orgId = `hgsql calJac3 -N -e \ "select id from organism where name='$orgName'"` if ($orgId == "") then echo "Mrnas: 0" else set count = `hgsql calJac3 -N -e "select count(*) from gbCdnaInfo where organism=$orgId"` echo "Mrnas: ${count}" endif end '_EOF_' # << happy emacs chmod +x ./showGenes.csh ./showGenes.csh # calJac3: ensGene: 55137, nscanGene: 28723, refGene: 53, xenoRefGene: 252282, Mrnas: 3789 # hg19: ensGene: 167071, knownGene: 77614, mgcGenes: 31354, nscanGene: 67448, refGene: 38868, xenoRefGene: 136935, Mrnas: 295362 # panTro2: ensGene: 41488, nscanGene: 22251, refGene: 1457, xenoRefGene: 233381, Mrnas: 2257 # gorGor2: xenoRefGene: 240083, Mrnas: 1 # ponAbe2: ensGene: 31566, nscanGene: 23499, refGene: 3549, xenoRefGene: 241018, Mrnas: 0 # papHam1: xenoRefGene: 267869, Mrnas: 71 # rheMac2: ensGene: 42820, nscanGene: 22003, refGene: 2563, xenoRefGene: 303624, Mrnas: 7297 # tarSyr1: ensGene: 45832, Mrnas: 8 # micMur1: ensGene: 37458, Mrnas: 53 # otoGar1: ensGene: 36560, Mrnas: 1 # mm9: ensGene: 93805, knownGene: 55419, mgcGenes: 26707, nscanGene: 21895, refGene: 28429, xenoRefGene: 132694, Mrnas: 261336 # canFam2: ensGene: 30914, nscanGene: 24615, refGene: 1205, xenoRefGene: 225414, Mrnas: 3476 # monDom5: ensGene: 35113, nscanGene: 45946, refGene: 467, xenoRefGene: 228506, Mrnas: 1719 # rearrange that output to create 3 sections: #1. knownGenes: hg19 mm9 #2. ensGene: calJac3,ponAbe2 panTro2 rheMac2 tarSyr1 micMur1 otoGar1 canFam2 monDom5 #3. xenoRefGene: papHam1 gorGor2 mkdir genes # knownGene echo "hg19 mm9" \ | sed -e "s/ */ /g" > knownGene.list for DB in hg19 mm9 do hgsql -N -e "select name,chrom,strand,txStart,txEnd,cdsStart,cdsEnd,exonCount,exonStarts,exonEnds from knownGene" ${DB} \ | genePredSingleCover stdin stdout | gzip -2c \ > /scratch/tmp/${DB}.tmp.gz mv /scratch/tmp/${DB}.tmp.gz genes/$DB.gp.gz echo "${DB} done" done echo "calJac3 ponAbe2 panTro2 rheMac2 tarSyr1 micMur1 otoGar1 canFam2 monDom5" \ | sed -e "s/ */ /g" > ensGene.list # ensGene for DB in calJac3 ponAbe2 panTro2 rheMac2 tarSyr1 micMur1 otoGar1 canFam2 monDom5 do hgsql -N -e "select name,chrom,strand,txStart,txEnd,cdsStart,cdsEnd,exonCount,exonStarts,exonEnds from ensGene" ${DB} \ | genePredSingleCover stdin stdout | gzip -2c \ > /scratch/tmp/${DB}.tmp.gz mv /scratch/tmp/${DB}.tmp.gz genes/$DB.gp.gz echo "${DB} done" done echo "papHam1 gorGor2" > xenoRef.list # xenoRefGene for DB in papHam1 gorGor2 do hgsql -N -e "select name,chrom,strand,txStart,txEnd,cdsStart,cdsEnd,exonCount,exonStarts,exonEnds from xenoRefGene" ${DB} \ | genePredSingleCover stdin stdout | gzip -2c \ > /scratch/tmp/${DB}.tmp.gz mv /scratch/tmp/${DB}.tmp.gz genes/$DB.gp.gz echo "${DB} done" done # verify counts for genes are reasonable: for T in genes/*.gz do echo -n "# $T: " zcat $T | cut -f1 | sort | uniq -c | wc -l done # genes/calJac3.gp.gz: 20843 # genes/canFam2.gp.gz: 19245 # genes/gorGor2.gp.gz: 22444 # genes/hg19.gp.gz: 20597 # genes/micMur1.gp.gz: 16240 # genes/mm9.gp.gz: 20905 # genes/monDom5.gp.gz: 19188 # genes/otoGar1.gp.gz: 15388 # genes/panTro2.gp.gz: 19757 # genes/papHam1.gp.gz: 26534 # genes/ponAbe2.gp.gz: 19895 # genes/rheMac2.gp.gz: 21049 # genes/tarSyr1.gp.gz: 13617 ssh hgwdev cd /hive/data/genomes/calJac3/bed/multiz13way/frames cat ../anno/calJac3.13way.maf \ | nice -n +19 genePredToMafFrames calJac3 stdin stdout \ `sed -e "s#\([a-zA-Z0-9]*\)#\1 genes/\1.gp.gz#g" ../species.list` \ | gzip > multiz13way.mafFrames.gz # verify there are frames on everything: zcat multiz13way.mafFrames.gz | awk '{print $4}' | sort | uniq -c # 200715 calJac3 # 232951 canFam2 # 220393 gorGor2 # 205490 hg19 # 205084 micMur1 # 229058 mm9 # 200391 monDom5 # 212293 otoGar1 # 202695 panTro2 # 212325 papHam1 # 205185 ponAbe2 # 202823 rheMac2 # 157578 tarSyr1 ssh hgwdev cd /hive/data/genomes/calJac3/bed/multiz13way/frames time hgLoadMafFrames calJac3 multiz13wayFrames multiz13way.mafFrames.gz # real 0m28.294s # check it: hgsql -N -e "select src from multiz13wayFrames;" calJac3 \ | sort | uniq -c | grep calJac3 # 200715 calJac3 # enable the trackDb entries: # frames multiz13wayFrames # irows on ############################################################################# ## create upstream refGene maf files mkdir /hive/data/genomes/calJac3/bed/multiz13way/downloads/maf cd /hive/data/genomes/calJac3/bed/multiz13way/downloads/maf # bash script #!/bin/sh for S in 1000 2000 5000 do echo "making upstream${S}.maf" featureBits calJac3 xenoRefGene:upstream:${S} -fa=/dev/null -bed=stdout \ | perl -wpe 's/_up[^\t]+/\t0/' | sort -k1,1 -k2,2n \ | /cluster/bin/$MACHTYPE/mafFrags calJac3 multiz13way \ stdin stdout \ -orgs=/hive/data/genomes/calJac3/bed/multiz13way/species.list \ | gzip -c > upstream${S}.maf.gz echo "done upstream${S}.maf.gz" done # making upstream1000.maf # done upstream1000.maf.gz # making upstream2000.maf # done upstream2000.maf.gz # making upstream5000.maf # done upstream5000.maf.gz # maf files cp -p /hive/data/genomes/calJac3/bed/multiz13way/run/calJac3.13way.maf . gzip calJac3.13way.maf # mv up one level to downloads mv *.* ../ cd .. rm -r maf # download stuff cd /hive/data/genomes/calJac3/bed/multiz13way cp 13way.nh 13way.nh.original cp tree.nh 13way.nh cp 13way.nh downloads/. cp commonNames.13way.nh downloads/. cd downloads # create /hive/data/genomes/calJac3/bed/multiz13way/downloads/README.txt md5sum * > md5sum.txt mkdir -p /usr/local/apache/htdocs-hgdownload/goldenPath/calJac3/multiz13way cd /usr/local/apache/htdocs-hgdownload/goldenPath/calJac3/multiz13way ln -s /hive/data/genomes/calJac3/bed/multiz13way/downloads/*.gz . ln -s /hive/data/genomes/calJac3/bed/multiz13way/downloads/*.nh . ln -s /hive/data/genomes/calJac3/bed/multiz13way/downloads/README.txt . ln -s /hive/data/genomes/calJac3/bed/multiz13way/downloads/md5sum.txt . ######################################################################### # Phylogenetic tree from 13-way (DONE 2010-05-28 - Chin) # Re-work (DONE 2011-04-25 - Chin) # We need one tree for all chroms # use the /hive/data/genomes/calJac3/bed/multiz13way/run/calJac3.13way.maf # to split cd /hive/data/genomes/calJac3/bed/multiz13way/ mkdir mafSplit cd mafSplit mafSplit -byTarget -useFullSequenceName /dev/null . ../run/calJac3.13way.maf # Splitting 1 files by target sequence -- ignoring first argument /dev/null ls -lt | wc -l # 10348 # got 10348 mafs named after their chrom/scaff .maf # although there are over 14205 chroms and scaffolds (wc -l # chrom.sizes), # some are too small or have nothing aligning. mkdir /hive/data/genomes/calJac3/bed/multiz13way/4d cd /hive/data/genomes/calJac3/bed/multiz13way/4d # using ensGene for calJac3, only transcribed genes and nothing # from the randoms and other misc. hgsql calJac3 -Ne \ "select * from ensGene WHERE cdsEnd > cdsStart;" | cut -f 2-20 \ | egrep -E -v "chrM|random|chrUn" \ > ensGene.gp # verify chromosome selection, should just be the ordinary chroms: cut -f2 ensGene.gp | sort | uniq -c wc -l *.gp # 42677 ensGene.gp genePredSingleCover ensGene.gp stdout | sort > ensGeneNR.gp wc -l ensGeneNR.gp # 19912 ensGeneNR.gp ssh memk mkdir /hive/data/genomes/calJac3/bed/multiz13way/4d/run cd /hive/data/genomes/calJac3/bed/multiz13way/4d/run mkdir ../mfa # whole chrom mafs version, using new version of # uses memory-efficient version of phast, from Melissa Hubisz at Cornell # mjhubisz at gmail.com cat << '_EOF_' > 4d.csh #!/bin/csh -fe set PHASTBIN = /cluster/bin/phast.build/cornellCVS/phast.2010-12-30/bin set r = "/hive/data/genomes/calJac3/bed/multiz13way" set c = $1 set infile = $r/mafSplit/$2 set outfile = $3 cd /scratch/tmp # 'clean' maf perl -wpe 's/^s ([^.]+)\.\S+/s $1/' $infile > $c.maf awk -v C=$c '$2 == C {print}' $r/4d/ensGene.gp > $c.gp.tmp sed "s/\t$c\t/\tcalJac3.$c\t/g" $c.gp.tmp > $c.gp set NL=`wc -l $c.gp| gawk '{print $1}'` if ("$NL" != "0") then set PHASTBIN = /cluster/bin/phast.build/cornellCVS/phast.2010-12-30/bin $PHASTBIN/msa_view --4d --features $c.gp -i MAF $c.maf -o SS > $c.ss $PHASTBIN/msa_view -i SS --tuple-size 1 $c.ss > $r/4d/run/$outfile else echo "" > $r/4d/run/$outfile endif rm -f $c.gp.tmp $c.gp $c.maf $c.ss '_EOF_' # << happy emacs chmod +x 4d.csh ls -1S /hive/data/genomes/calJac3/bed/multiz13way/mafSplit/*.maf | \ egrep -E -v "chrM|random|chrUn" \ | sed -e "s#.*multiz13way/mafSplit/##" \ > maf.list cat << '_EOF_' > template #LOOP 4d.csh $(root1) $(path1) {check out line+ ../mfa/$(root1).mfa} #ENDLOOP '_EOF_' # << happy emacs gensub2 maf.list single template stdout | tac > jobList para create jobList para try ... check ... push ... etc para time # Completed: 24 of 24 jobs # CPU time in finished jobs: 2913s 48.55m 0.81h 0.03d 0.000 y # IO & Wait Time: 187s 3.12m 0.05h 0.00d 0.000 y # Average job time: 129s 2.15m 0.04h 0.00d # Longest finished job: 238s 3.97m 0.07h 0.00d # Submission to last job: 421s 7.02m 0.12h 0.00d # combine mfa files ssh hgwdev cd /hive/data/genomes/calJac3/bed/multiz13way/4d # but first clean out junk 1-byte leftovers from above process. # Only 24 (real chrom) out of 2623 files have real data. cd mfa find -type f -size 1c | xargs -iX rm X find -type f -size 0c | xargs -iX rm X cd /hive/data/genomes/calJac3/bed/multiz13way/4d #want comma-less species.list /cluster/bin/phast.build/cornellCVS/phast.2010-12-30/bin/msa_view \ --aggregate "`cat ../species.list`" mfa/*.mfa | sed s/"> "/">"/ \ > 4d.all.mfa # check they are all in there: grep "^>" 4d.all.mfa # >calJac3 # >hg19 # >panTro2 # >gorGor2 # >ponAbe2 # >papHam1 # >rheMac2 # >tarSyr1 # >micMur1 # >otoGar1 # >mm9 # >canFam2 # >monDom5 # cat ../tree-commas.nh # ((((((calJac3,((((hg19,panTro2),gorGor2),ponAbe2),(papHam1,rheMac2))),tarSyr1),(micMur1,otoGar1)),mm9),canFam2),monDom5) # use phyloFit to create tree model (output is phyloFit.mod) /cluster/bin/phast.build/cornellCVS/phast.2010-12-30/bin/phyloFit \ --EM --precision MED --msa-format FASTA --subst-mod REV \ --tree ../tree-commas.nh 4d.all.mfa # Reading alignment from 4d.all.mfa ... # Extracting sufficient statistics ... # Compacting sufficient statistics ... # Fitting tree model to 4d.all.mfa using REV ... # Writing model to phyloFit.mod ... # Done. ############################################################################# # phastCons 13-way (DONE 2010-08-11 - Chin) # Re-work with phast suite 2011-12-30 version (DONE 2011-04-27 - Chin) # split 13way mafs into 10M chunks and generate sufficient # statistics # files for # phastCons ssh swarm mkdir -p /hive/data/genomes/calJac3/bed/multiz13way/cons/msa.split/ss cd /hive/data/genomes/calJac3/bed/multiz13way/cons/msa.split mkdir done cat << '_EOF_' > doSplit.csh #!/bin/csh -ef set c = $1 set MAF = /hive/data/genomes/calJac3/bed/multiz13way/mafSplit/$c.maf set WINDOWS = /hive/data/genomes/calJac3/bed/multiz13way/cons/msa.split/ss/$c set WC = `cat $MAF | wc -l` set NL = `grep "^#" $MAF | wc -l` if ( -s $2 ) then exit 0 endif if ( -s $2.running ) then exit 0 endif date >> $2.running rm -fr $WINDOWS mkdir $WINDOWS pushd $WINDOWS > /dev/null if ( $WC != $NL ) then /cluster/bin/phast.build/cornellCVS/phast.2010-12-30/bin/msa_split \ $MAF -i MAF -o SS -r $WINDOWS/$c -w 10000000,0 -I 1000 -B 5000 endif popd > /dev/null date >> $2 rm -f $2.running '_EOF_' # << happy emacs chmod +x doSplit.csh cat << '_EOF_' > template #LOOP doSplit.csh $(root1) {check out line+ done/$(root1).done} #ENDLOOP '_EOF_' # << happy emacs # do the easy ones first to see some immediate results ls -1S -r ../../mafSplit \ | sed -e "s/.maf//; s/calJac3_//" > maf.list gensub2 maf.list single template jobList para -ram=8g create jobList para try ... check ... etc # Completed: 10347 of 10347 jobs # CPU time in finished jobs: 3233s 53.88m 0.90h 0.04d 0.000 y # IO & Wait Time: 48949s 815.82m 13.60h 0.57d 0.002 y # Average job time: 5s 0.08m 0.00h 0.00d # Longest finished job: 323s 5.38m 0.09h 0.00d # Submission to last job: 574s 9.57m 0.16h 0.01d # some scaffolds were too small to produce output. # expected 10769, found 4749 cd /hive/data/genomes/calJac3/bed/multiz13way/cons/msa.split/ss find -type f | wc -l # 4913 cd ../../done rm *.done # Run phastCons # This job is I/O intensive in its output files, beware where this # takes place or do not run too many at once. ssh swarm mkdir -p /hive/data/genomes/calJac3/bed/multiz13way/cons/run.cons cd /hive/data/genomes/calJac3/bed/multiz13way/cons/run.cons # there are going to be several different phastCons runs using # this same script. They trigger off of the current working # directory # $cwd:t which is the "grp" in this script. It is # all cat << '_EOF_' > doPhast.csh #!/bin/csh -fe set PHASTBIN = /cluster/bin/phast.build/cornellCVS/phast.2010-12-30/bin set c = $1 set cX = $1:r set f = $2 set len = $3 set cov = $4 set rho = $5 set grp = $cwd:t set cons = /hive/data/genomes/calJac3/bed/multiz13way/cons set tmp = $cons/tmp/$f mkdir -p $tmp set ssSrc = $cons set useGrp = "$grp.mod" ln -s $ssSrc/msa.split/ss/$c/$f.ss $tmp ln -s $cons/$grp/$grp.mod $tmp pushd $tmp > /dev/null $PHASTBIN/phastCons $f.ss $useGrp \ --rho $rho --expected-length $len --target-coverage $cov --quiet \ --seqname $c --idpref $c --most-conserved $f.bed --score > $f.pp popd > /dev/null mkdir -p pp/$c bed/$c sleep 4 touch pp/$c bed/$c rm -f pp/$c/$f.pp rm -f bed/$c/$f.bed mv $tmp/$f.pp pp/$c mv $tmp/$f.bed bed/$c rm -fr $tmp '_EOF_' # << happy emacs chmod a+x doPhast.csh # this template will serve for all runs # root1 == chrom name, file1 == ss file name without .ss suffix cat << '_EOF_' > template #LOOP ../run.cons/doPhast.csh $(root1) $(file1) 45 0.3 0.3 {check out line+ pp/$(root1)/$(file1).pp} #ENDLOOP '_EOF_' # << happy emacs ls -1S ../msa.split/ss/*/* | sed -e 's/.ss$//' > ss.list # Create parasol batch and run it # run for all species cd /hive/data/genomes/calJac3/bed/multiz13way/cons mkdir -p all cd all # Using the .mod tree cp -p ../../4d/phyloFit.mod ./all.mod gensub2 ../run.cons/ss.list single ../run.cons/template jobList para -ram=8g create jobList para try ... check ... push ... etc. # Completed: 4913 of 4913 jobs # CPU time in finished jobs: 7560s 126.01m 2.10h 0.09d 0.000 y # IO & Wait Time: 44348s 739.13m 12.32h 0.51d 0.001 y # Average job time: 11s 0.18m 0.00h 0.00d # Longest finished job: 57s 0.95m 0.02h 0.00d # Submission to last job: 295s 4.92m 0.08h 0.00d # create Most Conserved track ssh hgwdev cd /hive/data/genomes/calJac3/bed/multiz13way/cons/all cat bed/*/*.bed | sort -k1,1 -k2,2n \ | awk '{printf "%s\t%d\t%d\tlod=%d\t%s\n", $1, $2, $3, $5, $5;}' \ > tmpMostConserved.bed /cluster/bin/scripts/lodToBedScore tmpMostConserved.bed > mostConserved.bed # load into database time nice -n +19 hgLoadBed calJac3 phastConsElements13way mostConserved.bed # Reading mostConserved.bed # Loaded 1531649 elements of size 5 # Sorted # Creating table definition for phastConsElements13way # Saving bed.tab # Loading calJac3 # # real 0m12.260s # Try for 5% overall cov, and 70% CDS cov featureBits calJac3 -enrichment ensGene:cds phastConsElements13way # ensGene:cds 1.221%, phastConsElements13way 4.704%, #both 0.845%, cover 69.19%, enrich 14.71x # hg19 for comparison # refGene:cds 1.196%, phastConsElements46way 5.065%, # both 0.888%, cover 74.22%, enrich 14.65x # ensGene:cds 1.278%, phastConsElements46way 5.065%, # both 0.910%, cover 71.23%, enrich 14.06x # knownGene:cds 1.252%, phastConsElements46way 5.065%, # both 0.905%, cover 72.29%, enrich 14.27x # Create merged posterier probability file and wiggle track data # files cd /hive/data/genomes/calJac3/bed/multiz13way/cons/all mkdir downloads find ./pp -type f | sed -e "s#^./##; s#\.# d #g; s#-# m #;" \ | sort -k1,1 -k3,3n | sed -e "s# d #.#g; s# m #-#g;" | xargs cat \ | gzip -c > downloads/phastCons13way.wigFix.gz # check integrity of data with wigToBigWig time (zcat downloads/phastCons13way.wigFix.gz \ | wigToBigWig -verbose=2 stdin /hive/data/genomes/calJac3/chrom.sizes \ phastCons13way.bw) > bigWig.log 2>&1 & grep real bigWig.log # real 28m26.948s grep VmPeak bigWig.log # pid=28269: VmPeak: 31261764 kB bigWigInfo phastCons13way.bw # version: 4 # isCompressed: yes # isSwapped: 0 # primaryDataSize: 4,004,563,748 # primaryIndexSize: 84,708,528 # zoomLevels: 10 # chromCount: 4648 # basesCovered: 2,302,681,633 # mean: 0.118453 # min: 0.000000 # max: 1.000000 # std: 0.238431 # encode those files into wiggle data time (zcat downloads/phastCons13way.wigFix.gz \ | wigEncode stdin phastCons13way.wig phastCons13way.wib) # Converted stdin, upper limit 1.00, lower limit 0.00 # real 10m35.994s du -hsc *.wi? # 2.2G phastCons13way.wib # 255M phastCons13way.wig # 2.4G total # Load gbdb and database with wiggle. ln -s `pwd`/phastCons13way.wib /gbdb/calJac3/multiz13way/phastCons13way.wib nice hgLoadWiggle -pathPrefix=/gbdb/calJac3/multiz13way calJac3 \ phastCons13way phastCons13way.wig # Connected to database calJac3 for track phastCons13way # Creating wiggle table definition in calJac3.phastCons13way # Saving wiggle.tab # Loading calJac3 # use to set trackDb.ra entries for wiggle min and max wigTableStats.sh calJac3 phastCons13way ## db.table min max mean count sumData stdDev viewLimits # calJac3.phastCons13way 0 1 0.118453 2302681633 2.7276e+08 0.238431 # viewLimits=0:1 # Create histogram to get an overview of all the data hgWiggle -doHistogram -db=calJac3 \ -hBinSize=0.001 -hBinCount=1000 -hMinVal=0.0 -verbose=2 \ phastCons13way >& histogram.data # create plot of histogram: #orig set terminal png small color x000000 xffffff xc000ff x66ff66 xffff00 x00ffff cat << '_EOF_' | gnuplot > histo.png set terminal png small size 1000,600 x000000 xffffff xc000ff x66ff66 xffff00 x00ffff set size 1.4, 0.8 set key left box set grid noxtics set grid ytics set title " Marmoset calJac3 Histogram phastCons13way track" set xlabel " phastCons13way score" set ylabel " Relative Frequency" set y2label " Cumulative Relative Frequency (CRF)" set y2range [0:1] set y2tics set yrange [0:0.02] plot "histogram.data" using 2:5 title " RelFreq" with impulses, \ "histogram.data" using 2:7 axes x1y2 title " CRF" with lines '_EOF_' # << happy emacs display histo.png & # download stuff mkdir -p /usr/local/apache/htdocs-hgdownload/goldenPath/calJac3/phastCons13way cd /usr/local/apache/htdocs-hgdownload/goldenPath/calJac3/phastCons13way ln -s /hive/data/genomes/calJac3/bed/multiz13way/cons/all/all.mod . ln -s /hive/data/genomes/calJac3/bed/multiz13way/cons/all/histo.png . ln -s /hive/data/genomes/calJac3/bed/multiz13way/cons/all/README.txt . ln -s /hive/data/genomes/calJac3/bed/multiz13way/cons/all/downloads/phastCons13way.wigFix.gz . md5sum phastCons13way.wigFix.gz > md5sum.txt md5sum all.mod >> md5sum.txt md5sum histo.png >> md5sum.txt md5sum README.txt >> md5sum.txt ############################################################################# # Quality track (DONE - 2010-07-14 - Hiram) mkdir -p /hive/data/genomes/calJac3/bed/qual/genbank cd /hive/data/genomes/calJac3/bed/qual/genbank cd /hive/data/genomes/calJac3/bed/qual grep "^chr" ../../calJac3.agp \ | egrep -v "contig|fragment|telo|centrom|clone" \ | awk '{print $6}' | sort -u > contig.name.list cd genbank wget --timestamping \ ftp://ftp.ncbi.nih.gov/genbank/wgs/wgs.ACFV.*.qscore.gz zcat wgs.ACFV.*.qscore.gz | grep "^>" | sed -e "s/^>gb|//; s/|.*//" \ sort -u > contig.list wc -l ../contig.name.list contig.list # 201514 ../contig.name.list # 202483 contig.list # not all the contigs have quality scores here, but almost all: comm -12 ../contig.name.list contig.list | wc -l # 201433 cat << '_EOF_' > extractContigs.pl #!/usr/bin/env perl use strict; use warnings; sub usage() { printf STDERR "usage: ./ctContigs.pl ../contig.name.list > calJac3.qual.fa\n"; printf STDERR "after collecting the contig names from the given list\n"; printf STDERR "this is going to read wgs.ACFV.*.qscore.gz\n"; printf STDERR "and output any fasta elements that match the name list\n"; } my $argc = scalar(@ARGV); if ($argc != 1) { usage; exit 255; } my $file = shift; my %ctgNames; my $readNames = 0; open (FH, "<$file") or die "can not open $file"; while (my $line = <FH>) { chomp $line; $ctgNames{$line} = 1; ++$readNames; } close (FH); printf STDERR "read $readNames from $file\n"; my $ctgCount = 0; foreach my $key (keys %ctgNames) { ++$ctgCount; } printf STDERR "found $ctgCount contig names in $file\n"; open (NF,">not.found.list") or die "can not write to not.found.list"; my $notFound = 0; my $ctgsDone = 0; for (my $i = 1; $i < 15; ++$i) { my $fname = sprintf("wgs.ACFV.%d.qscore.gz", $i); my $goodElement = 0; open (FH, "zcat $fname|") or die "can not read $fname"; while (my $line = <FH>) { if ($line =~ m/^>/) { chomp $line; my $name = $line; $name =~ s/^>gb\|//; $name =~ s/\|.*//; # printf STDERR "name: %s\n", $name; if (exists($ctgNames{$name})) { $goodElement = 1; printf ">%s\n", $name; ++$ctgsDone; } else { printf NF "%s:%s\n", $name, $line; ++$notFound; $goodElement = 0; } } elsif ($goodElement > 0) { $line =~ s/^\s+//; printf "%s", $line; } } close (FH); } close(NF); printf STDERR "found $ctgsDone during processing\n"; printf STDERR "not found $notFound during processing\n"; '_EOF_' # << happy emacs chmod +x ./extractContigs.pl time ./extractContigs.pl ../contig.name.list \ | gzip -c > ../calJac3.contigs.qual.fa.gz # read 201514 from ../contig.name.list # found 201514 contig names in ../contig.name.list # found 201433 during processing # not found 1050 during processing cd .. time qaToQac calJac3.contigs.qual.fa.gz calJac3.contigs.qual.qac # real 1m28.846s time qacAgpLift -mScore=0 ../../calJac3.agp \ calJac3.contigs.qual.qac calJac3.qual.qac # real 1m16.994s time qacToWig -fixed calJac3.qual.qac stdout | gzip -c \ > calJac3.qual.wigFix.gz # real 7m21.203s # avoid memory size problem sizeG=50000000 export sizeG ulimit -d $sizeG ulimit -v $sizeG time wigToBigWig calJac3.qual.wigFix.gz ../../chrom.sizes calJac3.qual.bw # real 29m7.143s bigWigInfo -minMax calJac3.qual.bw # 0.000000 97.000000 bigWigInfo calJac3.qual.bw # version: 4 # isCompressed: yes # isSwapped: 0 # primaryDataSize: 587,271,637 # primaryIndexSize: 91,661,572 # zoomLevels: 10 # chromCount: 14205 # basesCovered: 2,914,958,544 # mean: 87.707709 # min: 0.000000 # max: 97.000000 # std: 25.751106 mkdir /gbdb/calJac3/bbi ln -s `pwd`/calJac3.qual.bw /gbdb/calJac3/bbi hgsql calJac3 -e 'drop table if exists qualityBw; \ create table qualityBw \ (fileName varchar(255) not null); \ insert into qualityBw values ("/gbdb/calJac3/bbi/calJac3.qual.bw");' ############################################################################# # calJac3 - Marmoset - Ensembl Genes version 58 (DONE - 2010-07-20 - hiram) ssh hgwdev cd /hive/data/genomes/calJac3 cat << '_EOF_' > calJac3.ensGene.ra # required db variable db calJac3 nameTranslation "s/^\([0-9XY][0-9]*\)/chr\1/; s/^MT/chrM/; s/^Un/chrUn/" '_EOF_' # << happy emacs doEnsGeneUpdate.pl -ensVersion=58 calJac3.ensGene.ra ssh hgwdev cd /hive/data/genomes/calJac3/bed/ensGene.58 featureBits calJac3 ensGene # 50565085 bases of 2752505800 (1.837%) in intersection ######################################################################### # phyloP conservation for 13-way (Done 2010-08-31 - Chin) # Re-Work ( DONE 2011-04-29 - Chin) # # split SS files into 1M chunks, this business needs smaller files # to complete mkdir /hive/data/genomes/calJac3/bed/multiz13way/consPhyloP cd /hive/data/genomes/calJac3/bed/multiz13way/consPhyloP mkdir run.split cd run.split cat << '_EOF_' > doSplit.csh #!/bin/csh -ef set c = $1 set MAF = /hive/data/genomes/calJac3/bed/multiz13way/mafSplit/$c.maf set WINDOWS = /hive/data/genomes/calJac3/bed/multiz13way/consPhyloP/run.split/ss/$c set WC = `cat $MAF | wc -l` set NL = `grep "^#" $MAF | wc -l` if ( -s $2 ) then exit 0 endif if ( -s $2.running ) then exit 0 endif date >> $2.running rm -fr $WINDOWS mkdir $WINDOWS pushd $WINDOWS > /dev/null if ( $WC != $NL ) then /cluster/bin/phast.build/cornellCVS/phast.2010-12-30/bin/msa_split \ $MAF -i MAF -o SS -r $WINDOWS/$c -w 1000000,0 -I 1000 -B 5000 endif popd > /dev/null date >> $2 rm -f $2.running '_EOF_' # << happy emacs chmod +x doSplit.csh ls -1S -r ../../mafSplit | sed -e "s/.maf//" > maf.list cat << '_EOF_' > template #LOOP doSplit.csh $(path1) {check out exists+ done/$(path1).done} #ENDLOOP '_EOF_' # << happy emacs mkdir ss done ssh swarm cd /hive/data/genomes/calJac3/bed/multiz13way/consPhyloP/run.split gensub2 maf.list single template jobList para -ram=8g create jobList para try para push para check para time # Checking finished jobs # Completed: 10347 of 10347 jobs # CPU time in finished jobs: 3454s 57.57m 0.96h 0.04d 0.000 y # IO & Wait Time: 68218s 1136.97m 18.95h 0.79d 0.002 y # Average job time: 7s 0.12m 0.00h 0.00d # Longest finished job: 354s 5.90m 0.10h 0.00d # Submission to last job: 654s 10.90m 0.18h 0.01d # run phyloP with --method LRT ssh swarm mkdir /hive/data/genomes/calJac3/bed/multiz13way/consPhyloP/run.phyloP cd /hive/data/genomes/calJac3/bed/multiz13way/consPhyloP/run.phyloP # Adjust model file base composition background and rate matrix to # be # representative of the chromosomes in play grep BACKGROUND ../../cons/all/all.mod | awk '{printf "%0.3f\n", $3 + $4}' # 0.548 /cluster/bin/phast.build/cornellCVS/phast.2010-12-30/bin/modFreqs \ ../../cons/all/all.mod 0.548 > all.mod cat << '_EOF_' > doPhyloP.csh #!/bin/csh -fe set PHASTBIN = /cluster/bin/phast.build/cornellCVS/phast.2010-12-30/bin set f = $1 set out = $2 set cName = $f:r:r set chrDir = $f:r set n = $f:r:e set grp = $cwd:t set cons = /hive/data/genomes/calJac3/bed/multiz13way/consPhyloP set tmp = $cons/tmp/$grp/$f rm -fr $tmp mkdir -p $tmp set ssSrc = "$cons/run.split/ss/$chrDir/$f" set useGrp = "$grp.mod" ln -s $cons/run.phyloP/$grp.mod $tmp pushd $tmp > /dev/null $PHASTBIN/phyloP --method LRT --mode CONACC --wig-scores --chrom $cName \ -i SS $useGrp $ssSrc.ss > $f.wigFix popd > /dev/null mkdir -p $out:h sleep 4 mv $tmp/$f.wigFix $out rm -fr $tmp '_EOF_' # << happy emacs chmod +x doPhyloP.csh # Create list of chunks find ../run.split/ss -type f | sed -e "s/.ss$//; s#^../run.split/ss/##" \ > ss.list # Create template file # file1 == $chr/$chunk/file name without .ss suffix cat << '_EOF_' > template #LOOP ../run.phyloP/doPhyloP.csh $(file1) {check out line+ wigFix/$(dir1)/$(file1).wigFix} #ENDLOOP '_EOF_' # << happy emacs ###################### Running all species ####################### # setup run for all species mkdir /hive/data/genomes/calJac3/bed/multiz13way/consPhyloP/all cd /hive/data/genomes/calJac3/bed/multiz13way/consPhyloP/all rm -fr wigFix mkdir wigFix gensub2 ../run.phyloP/ss.list single ../run.phyloP/template jobList para create jobList para try ... check ... push ... etc ... para time # Completed: 7411 of 7411 jobs # CPU time in finished jobs: 129126s 2152.10m 35.87h 1.49d 0.004 y # IO & Wait Time: 189014s 3150.24m 52.50h 2.19d 0.006 y # Average job time: 43s 0.72m 0.01h 0.00d # Longest finished job: 127s 2.12m 0.04h 0.00d # Submission to last job: 454s 7.57m 0.13h 0.01d ssh hgwdev cd /hive/data/genomes/calJac3/bed/multiz13way/consPhyloP/all find ./wigFix -type f \ | sed -e "s#^./##; s/\./ /g; s/-/ - /g" \ | sort -k1,1 -k2,2n | sed -e "s/ - /-/g; s/ /./g" > wigFile.list cat wigFile.list | xargs cat | sed -e "s/__[0-9]//" \ | wigEncode stdin phyloP13way.wig \ phyloP13way.wib > wigEncode.log 2>&1 & cat wigEncode.log # Converted stdin, upper limit 1.28, lower limit -9.33 # good test to make sure no overlapping coordinates, bigWig: # consumes massive amount of memory, in bash raise your memory # limits: ulimit -d 188743680 ulimit -v 188743680 time cat wigFile.list | xargs cat | sed -e "s/__[0-9]//" \ | wigToBigWig stdin ../../../../chrom.sizes phyloP13way.bw & # real 38m1.704s # if you wanted to use the bigWig file, loading bigWig table: ln -s `pwd`/phyloP13way.bw /gbdb/calJac3/bbi hgsql calJac3 -e 'drop table if exists phyloP13wayAll; \ create table phyloP13wayAll \ (fileName varchar(255) not null); \ insert into phyloP13wayAll values ("/gbdb/calJac3/bbi/phyloP13way.bw");' # loading the wiggle table: ln -s `pwd`/phyloP13way.wib /gbdb/calJac3/multiz13way time nice -n +19 hgLoadWiggle -pathPrefix=/gbdb/calJac3/multiz13way calJac3 \ phyloP13wayAll phyloP13way.wig # Connected to database calJac3 for track phyloP13wayAll # Creating wiggle table definition in calJac3.phyloP13wayAll # Saving wiggle.tab # Loading calJac3 # # real 0m27.247s # create download files: cat << '_EOF_' > mkDown.csh #!/bin/csh -fe foreach F (`cat wigFile.list`) set C = $F:h:t:r cat $F | sed -e "s/__[0-9]//" >> downloads/${C}.wigFix end '_EOF_' # << happy emacs chmod +x ./mkDown.csh mkdir downloads time ./mkDown.csh # real 24m44.973s time gzip downloads/chr*.wigFix # real 30m14.304s cd downloads md5sum chr*.wigFix.gz >> md5sum.txt cd .. wigTableStats.sh calJac3 phyloP13wayAll # db.table min max mean count sumData # calJac3.phyloP13wayAll -9.328 1.283 0.0881734 2302681633 2.03035e+08 # stdDev 0.642244 # viewLimits=-3.12305:1.283 # that range is: 9.328 + 1.283 = 10.611 # Create histogram to get an overview of all the data time nice -n +19 hgWiggle -doHistogram \ -hBinSize=0.010611 -hBinCount=1000 -hMinVal=-9.32 -verbose=2 \ -db=calJac3 phyloP13wayAll > histogram.data 2>&1 # real 2m55.034s # create plot of histogram: cat << '_EOF_' | gnuplot > histo.png set terminal png small size 1000,600 x000000 xffffff xc000ff x66ff66 xffff00 x00ffff set size 1.4, 0.8 set key left box set grid noxtics set grid ytics set title " Marmoset calJac3 phyloP13way track" set xlabel " phyloP13way score" set ylabel " Relative Frequency" set y2label " Cumulative Relative Frequency (CRF)" set y2range [0:1] set y2tics set yrange [0:0.04] set xrange [-2:2] plot "histogram.data" using 2:5 title " RelFreq" with impulses, \ "histogram.data" using 2:7 axes x1y2 title " CRF" with lines '_EOF_' # << happy emacs display histo.png & # download stuff mkdir -p /usr/local/apache/htdocs-hgdownload/goldenPath/calJac3/phyloP13way cd /usr/local/apache/htdocs-hgdownload/goldenPath/calJac3/phyloP13way mkdir vertebrate cd vertebrate ln -s /hive/data/genomes/calJac3/bed/multiz13way/consPhyloP/all/downloads/chr*.wigFix.gz . ln -s /hive/data/genomes/calJac3/bed/multiz13way/consPhyloP/all/downloads/md5sum.txt . cd /usr/local/apache/htdocs-hgdownload/goldenPath/calJac3/phyloP13way ln -s /hive/data/genomes/calJac3/bed/multiz13way/consPhyloP/run.phyloP/all.mod . ln -s /hive/data/genomes/calJac3/bed/multiz13way/consPhyloP/README.txt md5sum all.mod > md5sum.txt md5sum README.txt >> md5sum.txt # packs tracks ############################################################################ # Add 13-Way conservation to pushQ /hive/data/genomes/calJac3/bed/multiz13way/quals/calJac3] ln -s /hive/data/genomes/calJac3/bed/qual/calJac3.qual.qac . /hive/data/genomes/calJac3/bed/multiz13way/quals/calJac3] ln -s /hive/data/genomes/calJac3/*.agp . qacAddGapIdx calJac3.agp calJac3.qual.qac calJac3.qac calJac3.qdx /hive/data/genomes/calJac3/bed/multiz13way/quals] time mafAddQRows 10way.list ../run/calJac3.13way.maf calJac3.quals.13way.maf ############################################################################ ## Fix gold/gap tables (DONE - 2010-09-22 - Hiram) # the hgGoldGapGl program has been fixed for chrom names > 16 characters mkdir /hive/data/genomes/calJac3/bed/goldGap cd /hive/data/genomes/calJac3/bed/goldGap # record table status before loading anything: hgsql -e "show table status;" calJac3 > before.status # verify SQL create statements with -noLoad option: hgGoldGapGl -noLoad -verbose=2 -noGl \ calJac3 /hive/data/genomes/calJac3/calJac3.agp # if that looks OK, load the two tables: hgGoldGapGl -verbose=2 -noGl \ calJac3 /hive/data/genomes/calJac3/calJac3.agp # verify only gold and gap tables have changed hgsql -e "show table status;" calJac3 > after.status # it appears only the gold table actually changed diff before.status after.status < gap MyISAM 10 Dynamic 187214 44 8366704 > gap MyISAM 10 Dynamic 187214 44 8366704 < gold MyISAM 10 Dynamic 194974 46 9129232 > gold MyISAM 10 Dynamic 201523 47 9496340 # gold and gap tables together now cover the entire genome: featureBits -countGaps -or calJac3 gold gap # 2914958544 bases of 2914958544 (100.000%) in intersection ############################################################################ # BacEnd track (DONE - 2010-07-12 - Hiram) # data obtained from Pat Minx (pminx at watson dot wustl dot edu) # and Wes Warren (wwarren at watson dot wustl dot edu) mkdir /hive/data/genomes/calJac3/bed/bacEnds cd /hive/data/genomes/calJac3/bed/bacEnds wget --timestamping \ http://genome.wustl.edu/pub/organism/Primates/Callithrix_jacchus/end_sequences/marmoset_bac_ends.3.2.bed.gz gunzip marmoset_bac_ends.3.2.bed.gz egrep -v "^track|^browser" marmoset_bac_ends.3.2.bed \ | hgLoadBed calJac3 bacEnds stdin checkTableCoords -verboseBlocks calJac3 bacEnds wget --timestamping \ http://genome.wustl.edu/pub/organism/Primates/Callithrix_jacchus/end_sequences/marmoset_singletons.3.2.bed.gz gunzip marmoset_singletons.3.2.bed.gz sort -k1,1 -k2,2n marmoset_singletons.3.2.bed \ | hgLoadBed calJac3 bacEndCalJac3Singles stdin checkTableCoords -verboseBlocks calJac3 bacEndCalJac3Singles ############################################################################ # Create liftOver files to and from: calJac3 <-> hg18 (DONE 2011-01-10 - Chin) # original alignment on hg19: cd /hive/data/genomes/hg18/bed/lastzCalJac3.2010-12-20 # and for this swap mkdir /hive/data/genomes/calJac3/bed/blastz.hg18.swap cd /hive/data/genomes/calJac3/bed/blastz.hg18.swap time nice -n +19 doBlastzChainNet.pl -verbose=2 \ /hive/data/genomes/hg18/bed/lastzCalJac3.2010-12-20/DEF \ -swap -syntenicNet \ -stop net \ -workhorse=hgwdev -smallClusterHub=swarm -bigClusterHub=swarm \ -chainMinScore=5000 -chainLinearGap=medium > swap.log 2>&1 & # real 48m54s cd /hive/data/genomes/calJac3/bed/ ln -s blastz.hg18.swap lastz.hg18 cd /hive/data/genomes/calJac3/bed/lastz.hg18/axtChain cp calJac3.hg18.over.chain.gz ../../liftOver/. cd /hive/data/genomes/calJac3/bed/liftOver mv calJac3.hg18.over.chain.gz calJac3ToHg18.over.chain.gz md5sum *.gz > md5sum.txt cd /usr/local/apache/htdocs-hgdownload/goldenPath/calJac3/liftOver ln -s /hive/data/genomes/calJac3/bed/liftOver/calJac3ToHg18.over.chain.gz . ############################################################################ # construct ucscToEnsembl chrom name translation (2011-04-20 - Hiram) # the Ensembl genes v62 have new chrom names, construct translation: mkdir /hive/data/genomes/calJac3/bed/ucscToEnsembl cd /hive/data/genomes/calJac3/bed/ucscToEnsembl zcat ../ensGene.62/download/Callithrix_jacchus.C_jacchus3.2.1.62.gtf.gz \ | awk '{print $1}' | sort -u \ | sed -e 's/^\([0-9XY][0-9]*\)/chr\1/; s/\.1$//;' | sort > t.t cat t.t | while read C do WC=`grep "${C}" ../../chrom.sizes | wc -l` EN=`echo ${C} | sed -e "s/^chr//; s/^\(ACFV.*\)/\1.1/; s/^\(GL.*\)/\1.1/;"` if [ "${WC}" -eq 1 ]; then UCSC=`grep "${C}" ../../chrom.sizes | awk '{print $1}'` echo -e "${UCSC}\t${EN}" else WC=`egrep "^${C} " ../../chrom.sizes | wc -l` if [ "${WC}" -eq 1 ]; then UCSC=`egrep "^${C} " ../../chrom.sizes | awk '{print $1}'` echo -e "${UCSC}\t${EN}" else echo "ERROR: one only not found for: ${C} - ${WC}" fi fi done | sort > ensemblGenes.tab cut -f1 ensemblGenes.tab | sort > ucsc.done.list cut -f1 ../../chrom.sizes | sort > ucsc.all.list comm -13 ucsc.done.list ucsc.all.list | while read UCSC do EN=`echo $UCSC | sed -e 's/^chr[0-9XYU][n0-9]*_//; s/_random//; s/$/.1/'` echo -e "${UCSC}\t${EN}" done > noEnsemblGenes.tab sort ensemblGenes.tab noEnsemblGenes.tab > ucscToEnsembl.tab cat << '_EOF_' > ucscToEnsembl.sql # UCSC to Ensembl chr name translation CREATE TABLE ucscToEnsembl ( ucsc varchar(255) not null, # UCSC chromosome name ensembl varchar(255) not null, # Ensembl chromosome name #Indices PRIMARY KEY(ucsc(21)) ); '_EOF_' hgsql calJac3 < ucscToEnsembl.sql hgsql calJac3 \ -e 'LOAD DATA LOCAL INFILE "ucscToEnsembl.tab" INTO TABLE ucscToEnsembl' ############################################################################ # calJac3 - Marmoset - Ensembl Genes version 62 (DONE - 2011-04-21 - hiram) # need a new lift file to translate names cd /hive/data/genomes/calJac3/jkStuff grep "^chr" ../chrom.sizes | while read L do size=`echo $L | awk '{print $2}'` ucName=`echo $L | awk '{print $1}'` ensName=`hgsql -N -e 'select ensembl from ucscToEnsembl where ucsc="'${ucName}'";' calJac3` echo -e "0\t${ensName}.1\t$size\t$ucName\t$size\n" done > ens.62.lft # place that file name in the calJac3.ensGene.ra file, then: ssh hgwdev cd /hive/data/genomes/calJac3 cat << '_EOF_' > calJac3.ensGene.ra # required db variable db calJac3 # nameTranslation "s/^\([0-9XY][0-9]*\)/chr\1/; s/^MT/chrM/; s/^Un/chrUn/; # s/^GL\([0-9][0-9]*\).1/chrUn_GL\1/;" # name translation in Ensembl v62 liftUp /hive/data/genomes/calJac3/jkStuff/ens.62.lft '_EOF_' # << happy emacs doEnsGeneUpdate.pl -ensVersion=62 calJac3.ensGene.ra ssh hgwdev cd /hive/data/genomes/calJac3/bed/ensGene.62 featureBits calJac3 ensGene # 51944550 bases of 2752505800 (1.887%) in intersection ############################################################################ # hgPal downloads (DONE 2011-05-20 braney) # FASTA from 13way for refGene, ensGene, nscanGene ssh hgwdev screen bash rm -rf /cluster/data/calJac3/bed/multiz13way/pal mkdir /cluster/data/calJac3/bed/multiz13way/pal cd /cluster/data/calJac3/bed/multiz13way/pal for i in `cat ../species.list`; do echo $i; done > order.lst mz=multiz13way gp=refGene db=calJac3 mkdir exonAA exonNuc ppredAA ppredNuc for j in `sort -nk 2 /cluster/data/$db/chrom.sizes | awk '{print $1}'` do echo "date" echo "mafGene -chrom=$j $db $mz $gp order.lst stdout | \ gzip -c > ppredAA/$j.ppredAA.fa.gz" echo "mafGene -chrom=$j -noTrans $db $mz $gp order.lst stdout | \ gzip -c > ppredNuc/$j.ppredNuc.fa.gz" echo "mafGene -chrom=$j -exons -noTrans $db $mz $gp order.lst stdout | \ gzip -c > exonNuc/$j.exonNuc.fa.gz" echo "mafGene -chrom=$j -exons $db $mz $gp order.lst stdout | \ gzip -c > exonAA/$j.exonAA.fa.gz" done > $gp.jobs time sh -x $gp.jobs > $gp.jobs.log 2>&1 & sleep 1 tail -f $gp.jobs.log # real 5m40.876s # user 2m32.229s # sys 2m48.508s mz=multiz13way gp=refGene db=calJac3 zcat exonAA/*.gz | gzip -c > $gp.$mz.exonAA.fa.gz zcat exonNuc/*.gz | gzip -c > $gp.$mz.exonNuc.fa.gz zcat ppredAA/*.gz | gzip -c > $gp.$mz.ppredAA.fa.gz zcat ppredNuc/*.gz | gzip -c > $gp.$mz.ppredNuc.fa.gz rm -rf exonAA exonNuc ppredAA ppredNuc # we're only distributing exons at the moment mz=multiz13way gp=refGene db=calJac3 pd=/usr/local/apache/htdocs/goldenPath/$db/$mz/alignments mkdir -p $pd ln -s `pwd`/$gp.$mz.exonAA.fa.gz $pd/$gp.exonAA.fa.gz ln -s `pwd`/$gp.$mz.exonNuc.fa.gz $pd/$gp.exonNuc.fa.gz mz=multiz13way gp=ensGene db=calJac3 mkdir exonAA exonNuc ppredAA ppredNuc for j in `sort -nk 2 /cluster/data/$db/chrom.sizes | awk '{print $1}'` do echo "date" echo "mafGene -chrom=$j $db $mz $gp order.lst stdout | \ gzip -c > ppredAA/$j.ppredAA.fa.gz" echo "mafGene -chrom=$j -noTrans $db $mz $gp order.lst stdout | \ gzip -c > ppredNuc/$j.ppredNuc.fa.gz" echo "mafGene -chrom=$j -exons -noTrans $db $mz $gp order.lst stdout | \ gzip -c > exonNuc/$j.exonNuc.fa.gz" echo "mafGene -chrom=$j -exons $db $mz $gp order.lst stdout | \ gzip -c > exonAA/$j.exonAA.fa.gz" done > $gp.jobs time sh -x $gp.jobs > $gp.jobs.log 2>&1 & sleep 1 tail -f $gp.jobs.log # real 109m32.085s # user 11m13.926s # sys 7m8.591s mz=multiz13way gp=ensGene db=calJac3 zcat exonAA/*.gz | gzip -c > $gp.$mz.exonAA.fa.gz zcat exonNuc/*.gz | gzip -c > $gp.$mz.exonNuc.fa.gz zcat ppredAA/*.gz | gzip -c > $gp.$mz.ppredAA.fa.gz zcat ppredNuc/*.gz | gzip -c > $gp.$mz.ppredNuc.fa.gz rm -rf exonAA exonNuc ppredAA ppredNuc # we're only distributing exons at the moment mz=multiz13way gp=ensGene db=calJac3 pd=/usr/local/apache/htdocs/goldenPath/$db/$mz/alignments mkdir -p $pd ln -s `pwd`/$gp.$mz.exonAA.fa.gz $pd/$gp.exonAA.fa.gz ln -s `pwd`/$gp.$mz.exonNuc.fa.gz $pd/$gp.exonNuc.fa.gz mz=multiz13way gp=nscanGene db=calJac3 mkdir exonAA exonNuc ppredAA ppredNuc for j in `sort -nk 2 /cluster/data/$db/chrom.sizes | awk '{print $1}'` do echo "date" echo "mafGene -chrom=$j $db $mz $gp order.lst stdout | \ gzip -c > ppredAA/$j.ppredAA.fa.gz" echo "mafGene -chrom=$j -noTrans $db $mz $gp order.lst stdout | \ gzip -c > ppredNuc/$j.ppredNuc.fa.gz" echo "mafGene -chrom=$j -exons -noTrans $db $mz $gp order.lst stdout | \ gzip -c > exonNuc/$j.exonNuc.fa.gz" echo "mafGene -chrom=$j -exons $db $mz $gp order.lst stdout | \ gzip -c > exonAA/$j.exonAA.fa.gz" done > $gp.jobs time sh -x $gp.jobs > $gp.jobs.log 2>&1 & sleep 1 tail -f $gp.jobs.log # real 76m43.791s # user 6m21.521s # sys 6m46.179s mz=multiz13way gp=nscanGene db=calJac3 zcat exonAA/*.gz | gzip -c > $gp.$mz.exonAA.fa.gz zcat exonNuc/*.gz | gzip -c > $gp.$mz.exonNuc.fa.gz zcat ppredAA/*.gz | gzip -c > $gp.$mz.ppredAA.fa.gz zcat ppredNuc/*.gz | gzip -c > $gp.$mz.ppredNuc.fa.gz rm -rf exonAA exonNuc ppredAA ppredNuc # we're only distributing exons at the moment mz=multiz13way gp=nscanGene db=calJac3 pd=/usr/local/apache/htdocs/goldenPath/$db/$mz/alignments mkdir -p $pd ln -s `pwd`/$gp.$mz.exonAA.fa.gz $pd/$gp.exonAA.fa.gz ln -s `pwd`/$gp.$mz.exonNuc.fa.gz $pd/$gp.exonNuc.fa.gz ############################################################################ # LASTZ Gorilla calJac3 (DONE - 2011-10-21,27 - Hiram) mkdir /hive/data/genomes/calJac3/bed/lastzGorGor3.2011-10-21 cd /hive/data/genomes/calJac3/bed/lastzGorGor3.2011-10-21 cat << '_EOF_' > DEF # Gorilla vs Marmoset # maximum M allowed with lastz is only 254 BLASTZ_M=254 BLASTZ_Q=/scratch/data/blastz/human_chimp.v2.q BLASTZ_O=600 BLASTZ_E=150 # other parameters on advice from Webb BLASTZ_K=4500 BLASTZ_Y=15000 BLASTZ_T=2 # TARGET: Marmoset calJac3 SEQ1_DIR=/scratch/data/calJac3/calJac3.2bit SEQ1_LEN=/scratch/data/calJac3/chrom.sizes SEQ1_CHUNK=10000000 SEQ1_LAP=10000 SEQ1_LIMIT=60 # QUERY: Gorilla gorGor3 SEQ2_DIR=/scratch/data/gorGor3/gorGor3.2bit SEQ2_LEN=/scratch/data/gorGor3/chrom.sizes SEQ2_CHUNK=10000000 SEQ2_LAP=0 SEQ2_LIMIT=300 BASE=/hive/data/genomes/calJac3/bed/lastzGorGor3.2011-10-21 TMPDIR=/scratch/tmp '_EOF_' # << happy emacs # establish a screen to control this job screen time nice -n +19 doBlastzChainNet.pl -verbose=2 \ `pwd`/DEF \ -syntenicNet \ -noLoadChainSplit -chainMinScore=5000 -chainLinearGap=medium \ -workhorse=hgwdev -smallClusterHub=encodek -bigClusterHub=swarm \ > do.log 2>&1 & # after recovering manually after a hive filesystem outage time nice -n +19 doBlastzChainNet.pl -verbose=2 \ `pwd`/DEF \ -continue=cat -syntenicNet \ -noLoadChainSplit -chainMinScore=5000 -chainLinearGap=medium \ -workhorse=hgwdev -smallClusterHub=encodek -bigClusterHub=swarm \ > cat.log 2>&1 & # real 147m39.439s cat fb.calJac3.chainGorGor3Link.txt # 1928768328 bases of 2752505800 (70.073%) in intersection cd /hive/data/genomes/calJac3/bed ln -s lastzGorGor3.2011-10-21 lastz.gorGor3 # running the swap mkdir /hive/data/genomes/gorGor3/bed/blastz.calJac3.swap cd /hive/data/genomes/gorGor3/bed/blastz.calJac3.swap time nice -n +19 doBlastzChainNet.pl -verbose=2 \ /hive/data/genomes/calJac3/bed/lastzGorGor3.2011-10-21/DEF \ -swap -syntenicNet \ -noLoadChainSplit -chainMinScore=5000 -chainLinearGap=medium \ -workhorse=hgwdev -smallClusterHub=encodek -bigClusterHub=swarm \ > swap.log 2>&1 & # real 86m1.559s cat fb.gorGor3.chainCalJac3Link.txt # 1933784568 bases of 2822760080 (68.507%) in intersection ############################################################################ # LASTZ Rhesus rheMac3 SWAP (DONE - 2012-07-19 - Chin) cd /hive/data/genomes/rheMac3/bed/lastzCalJac3.2012-07-19 cat fb.rheMac3.chainCalJac3Link.txt # 1893656467 bases of 2639145830 (71.753%) in intersection # running the swap mkdir /hive/data/genomes/calJac3/bed/blastz.rheMac3.swap cd /hive/data/genomes/calJac3/bed/blastz.rheMac3.swap time nice -n +19 doBlastzChainNet.pl -verbose=2 \ /hive/data/genomes/rheMac3/bed/lastzCalJac3.2012-07-19/DEF \ -swap -syntenicNet \ -workhorse=hgwdev -smallClusterHub=memk -bigClusterHub=swarm \ -chainMinScore=5000 -chainLinearGap=medium > swap.log 2>&1 & # real 84m31.344s cat fb.calJac3.chainRheMac3Link.txt # 1931533122 bases of 2752505800 (70.174%) in intersection cd /hive/data/genomes/calJac3/bed ln -s blastz.rheMac3.swap lastz.rheMac3 ############################################################################ # create ucscToINSDC name mapping (DONE - 2013-08-16 - Hiram) mkdir /hive/data/genomes/calJac3/bed/ucscToINSDC cd /hive/data/genomes/calJac3/bed/ucscToINSDC # copying these scripts from the previous load and improving them # with each instance ./translateNames.sh ./verifyAll.sh ./join.sh # verify the track link to INSDC functions ############################################################################## # genscan - (DONE - 2013-11-12 - Hiram) mkdir /hive/data/genomes/calJac3/bed/genscan cd /hive/data/genomes/calJac3/bed/genscan time doGenscan.pl calJac3 > do.log 2>&1 # real 20m7.222s one broken job: ./lastRunGsBig.csh chr13 000 gtf/000/chr13.gtf pep/000/chr13.pep subopt/000/chr13.bed # real 38m21.531s # run 'para time > run.time' on ku, then continuing: time doGenscan.pl -continue=makeBed calJac3 > makeBed.log 2>&1 # real 4m47.592s cat fb.calJac3.genscan.txt # 56633537 bases of 2752505800 (2.058%) in intersection cat fb.calJac3.genscanSubopt.txt # 56321789 bases of 2752505800 (2.046%) in intersection ######################################################################### ############################################################################## # TransMap V3 tracks. see makeDb/doc/transMapTracks.txt (2014-12-21 markd) ############################################################################## # LASTZ chimp panTro4 (DONE - 2015-02-09 - Hiram) mkdir /hive/data/genomes/calJac3/bed/lastzPanTro4.2015-02-09 cd /hive/data/genomes/calJac3/bed/lastzPanTro4.2015-02-09 cat << '_EOF_' > DEF # Marmoset vs chimp # parameters obtained from a tuning run of lastz_D # /hive/users/hiram/multiz/100way/calJac3.panTro4/calJac3.panTro4.tuning.top300.txt BLASTZ=/cluster/bin/penn/lastz-distrib-1.03.54/bin/lastz BLASTZ_T=2 BLASTZ_O=400 BLASTZ_E=30 BLASTZ_M=50 BLASTZ_X=1000 BLASTZ_Y=3400 BLASTZ_Q=/hive/data/genomes/calJac3/bed/lastzPanTro4.2015-02-09/calJac3.panTro4.tuning.Q.txt # A C G T # A 100 -177 -40 -144 # C -177 94 -185 -40 # G -40 -185 94 -177 # T -144 -40 -177 100 # TARGET: Marmoset calJac3 SEQ1_DIR=/hive/data/genomes/calJac3/calJac3.2bit SEQ1_LEN=/hive/data/genomes/calJac3/chrom.sizes SEQ1_CHUNK=20000000 SEQ1_LIMIT=100 SEQ1_LAP=10000 # QUERY: chimp panTro4 SEQ2_DIR=/hive/data/genomes/panTro4/panTro4.2bit SEQ2_LEN=/hive/data/genomes/panTro4/chrom.sizes SEQ2_CHUNK=20000000 SEQ2_LIMIT=150 SEQ2_LAP=0 BASE=/hive/data/genomes/calJac3/bed/lastzPanTro4.2015-02-09 TMPDIR=/dev/shm '_EOF_' # << happy emacs # establish a screen to control this job time (doBlastzChainNet.pl `pwd`/DEF -verbose=2 \ -chainMinScore=3000 -chainLinearGap=medium \ -workhorse=hgwdev -smallClusterHub=ku -bigClusterHub=ku \ -syntenicNet) > do.log 2>&1 # real 323m11.555s cat fb.calJac3.chainPanTro4Link.txt # 2076330774 bases of 2752505800 (75.434%) in intersection time (doRecipBest.pl -buildDir=`pwd` calJac3 panTro4) > rbest.log 2>&1 & # real 55m50.180s # running the swap mkdir /hive/data/genomes/panTro4/bed/blastz.calJac3.swap cd /hive/data/genomes/panTro4/bed/blastz.calJac3.swap time (doBlastzChainNet.pl -verbose=2 \ /hive/data/genomes/calJac3/bed/lastzPanTro4.2015-02-09/DEF \ -swap -syntenicNet -chainMinScore=3000 -chainLinearGap=medium \ -workhorse=hgwdev -smallClusterHub=ku \ -bigClusterHub=ku) > swap.log 2>&1 # real 114m27.647s cat fb.panTro4.chainCalJac3Link.txt # 2113187925 bases of 2902338967 (72.810%) in intersection time (doRecipBest.pl -buildDir=`pwd` panTro4 calJac3) > rbest.log 2>&1 & # real 44m18.149s ############################################################################## # LASTZ chimp panTro5 (DONE - 2017-08-08 - Hiram) mkdir /hive/data/genomes/calJac3/bed/lastzPanTro5.2017-08-08 cd /hive/data/genomes/calJac3/bed/lastzPanTro5.2017-08-08 printf '# Marmoset vs chimp BLASTZ=/cluster/bin/penn/lastz-distrib-1.03.66/bin/lastz BLASTZ_M=254 # TARGET: Marmoset calJac3 SEQ1_DIR=/hive/data/genomes/calJac3/calJac3.2bit SEQ1_LEN=/hive/data/genomes/calJac3/chrom.sizes SEQ1_CHUNK=40000000 SEQ1_LIMIT=100 SEQ1_LAP=10000 # QUERY: chimp panTro5 SEQ2_DIR=/hive/data/genomes/panTro5/panTro5.2bit SEQ2_LEN=/hive/data/genomes/panTro5/chrom.sizes SEQ2_CHUNK=20000000 SEQ2_LIMIT=200 SEQ2_LAP=0 BASE=/hive/data/genomes/calJac3/bed/lastzPanTro5.2017-08-08 TMPDIR=/dev/shm ' > DEF # establish a screen to control this job time (doBlastzChainNet.pl `pwd`/DEF -verbose=2 \ -chainMinScore=3000 -chainLinearGap=medium \ -workhorse=hgwdev -smallClusterHub=ku -bigClusterHub=ku \ -syntenicNet) > do.log 2>&1 # real 676m38.873s cat fb.calJac3.chainPanTro5Link.txt # 2208118247 bases of 2752505800 (80.222%) in intersection time (doRecipBest.pl -workhorse=hgwdev -buildDir=`pwd` calJac3 panTro5) \ > rbest.log 2>&1 & # real 937m46.656s # running the swap mkdir /hive/data/genomes/panTro5/bed/blastz.calJac3.swap cd /hive/data/genomes/panTro5/bed/blastz.calJac3.swap time (doBlastzChainNet.pl -verbose=2 \ /hive/data/genomes/calJac3/bed/lastzPanTro5.2017-08-08/DEF \ -swap -syntenicNet -chainMinScore=3000 -chainLinearGap=medium \ -workhorse=hgwdev -smallClusterHub=ku \ -bigClusterHub=ku) > swap.log 2>&1 # real 140m3.455s cat fb.panTro5.chainCalJac3Link.txt # 2319213349 bases of 3132620660 (74.034%) in intersection time (doRecipBest.pl -workhorse=hgwdev -buildDir=`pwd` panTro5 calJac3) \ > rbest.log 2>&1 & # real 850m40.412s ############################################################################## # LASTZ chimp rheMac8 (DONE - 2017-08-08 - Hiram) mkdir /hive/data/genomes/calJac3/bed/lastzRheMac8.2017-08-08 cd /hive/data/genomes/calJac3/bed/lastzRheMac8.2017-08-08 printf '# Marmoset vs rhesus BLASTZ=/cluster/bin/penn/lastz-distrib-1.03.66/bin/lastz BLASTZ_M=254 # TARGET: Marmoset calJac3 SEQ1_DIR=/hive/data/genomes/calJac3/calJac3.2bit SEQ1_LEN=/hive/data/genomes/calJac3/chrom.sizes SEQ1_CHUNK=40000000 SEQ1_LIMIT=100 SEQ1_LAP=10000 # QUERY: rhesus rheMac8 SEQ2_DIR=/hive/data/genomes/rheMac8/rheMac8.2bit SEQ2_LEN=/hive/data/genomes/rheMac8/chrom.sizes SEQ2_CHUNK=20000000 SEQ2_LIMIT=1000 SEQ2_LAP=0 BASE=/hive/data/genomes/calJac3/bed/lastzRheMac8.2017-08-08 TMPDIR=/dev/shm ' > DEF # establish a screen to control this job time (doBlastzChainNet.pl `pwd`/DEF -verbose=2 \ -chainMinScore=3000 -chainLinearGap=medium \ -workhorse=hgwdev -smallClusterHub=ku -bigClusterHub=ku \ -syntenicNet) > do.log 2>&1 # real 635m57.336s cat fb.calJac3.chainRheMac8Link.txt # 2183155505 bases of 2752505800 (79.315%) in intersection time (doRecipBest.pl -workhorse=hgwdev -buildDir=`pwd` calJac3 rheMac8) \ > rbest.log 2>&1 & # real 1056m51.268s # running the swap mkdir /hive/data/genomes/rheMac8/bed/blastz.calJac3.swap cd /hive/data/genomes/rheMac8/bed/blastz.calJac3.swap time (doBlastzChainNet.pl -verbose=2 \ /hive/data/genomes/calJac3/bed/lastzRheMac8.2017-08-08/DEF \ -swap -syntenicNet -chainMinScore=3000 -chainLinearGap=medium \ -workhorse=hgwdev -smallClusterHub=ku \ -bigClusterHub=ku) > swap.log 2>&1 # real 265m36.306s cat fb.rheMac8.chainCalJac3Link.txt # 2208994431 bases of 3142093174 (70.303%) in intersection time (doRecipBest.pl -workhorse=hgwdev -buildDir=`pwd` rheMac8 calJac3) \ > rbest.log 2>&1 & # real 1183m6.027s -############################################################################ +############################################################################## +# LIFTOVER TO calJac4 (DONE - 2020-09-04 - Hiram) + ssh hgwdev + mkdir /hive/data/genomes/calJac3/bed/blat.calJac4.2020-09-04 + cd /hive/data/genomes/calJac3/bed/blat.calJac4.2020-09-04 + doSameSpeciesLiftOver.pl -verbose=2 \ + -debug -bigClusterHub=ku -dbHost=hgwdev -workhorse=hgwdev \ + -target2Bit=/hive/data/genomes/calJac3/calJac3.2bit \ + -targetSizes=/hive/data/genomes/calJac3/chrom.sizes \ + -ooc=/hive/data/genomes/calJac3/jkStuff/calJac3.11.ooc \ + calJac3 calJac4 + time (doSameSpeciesLiftOver.pl -verbose=2 \ + -bigClusterHub=ku -dbHost=hgwdev -workhorse=hgwdev \ + -target2Bit=/hive/data/genomes/calJac3/calJac3.2bit \ + -targetSizes=/hive/data/genomes/calJac3/chrom.sizes \ + -ooc=/hive/data/genomes/calJac3/jkStuff/calJac3.11.ooc \ + calJac3 calJac4) > doLiftOverToCalJac4.log 2>&1 + # real 257m56.056s + + # see if the liftOver menus function in the browser from calJac3 to calJac4 + +##############################################################################