# for emacs: -*- mode: sh; -*- # $Id: ailMel1.txt,v 1.15 2010/06/11 17:12:56 chinhli Exp $ # Panda sequence: http://panda.genomics.org.cn/page/panda/download.jsp # Ailuropoda melanoleuca # http://www.ncbi.nlm.nih.gov/genome/2702 # http://www.ncbi.nlm.nih.gov/bioproject/38683 # http://www.ncbi.nlm.nih.gov/Traces/wgs/?val=ACTA00 ########################################################################## # Download sequence (DONE - 2009-12-17 - Hiram) mkdir /hive/data/genomes/ailMel1 cd /hive/data/genomes/ailMel1 mkdir genbank cd genbank wget --timestamping -r --cut-dirs=6 --level=0 -nH -x \ --no-remove-listing -np \ "ftp.ncbi.nlm.nih.gov:genbank/genomes/Eukaryotes/vertebrates_mammals/Ailuropoda_melanoleuca/AilMel_1.0/*" # first pass attempt at sequence from public.genomics.org.cn # 2009-12-17 # NCBI has cleaned the sequence and removed a couple of scaffolds that # are cloning vector sequence mkdir download cd download wget --timestamping \ 'ftp://public.genomics.org.cn/BGI/panda/genome_sequence/panda.scafSeq.gapFilled.noMito.GeneScaffold.fa.gz' wget --timestamping \ 'ftp://public.genomics.org.cn/BGI/panda/genome_sequence/panda_contig100.mapping_scaffolds_nucmer.list.clean.tab.sort.gff.gz' wget --timestamping \ 'ftp://panda.genomics.org.cn/pub/panda/panda.scafSeq.gapFilled.noMito.GeneScaffold.fa.gz' faCount panda.scafSeq.gapFilled.noMito.GeneScaffold.fa.gz > faCount.txt grep "^scaff" faCount.txt | awk '{print $1}' > scaffold.list faSomeRecords panda.scafSeq.gapFilled.noMito.GeneScaffold.fa.gz \ scaffold.list stdout | gzip > scaffolds.fa.gz hgFakeAgp scaffolds.fa.gz stdout | sed -e "s/D/W/" > fake.agp ########################################################################## # Initial genome build (DONE - 2009-12-17 - Hiram) cd /hive/data/genomes/ailMel1 cat << '_EOF_' > ailMel1.config.ra # Config parameters for makeGenomeDb.pl: db ailMel1 clade mammal genomeCladePriority 19 scientificName Ailuropoda melanoleuca commonName Panda assemblyDate Dec. 2009 assemblyLabel BGI-Shenzhen AilMel 1.0 Dec. 2009 orderKey 220 mitoAcc NC_009492 fastaFiles /hive/data/genomes/ailMel1/genbank/Primary_Assembly/unplaced_scaffolds/FASTA/unplaced.scaf.fa.gz agpFiles /hive/data/genomes/ailMel1/genbank/Primary_Assembly/unplaced_scaffolds/AGP/unplaced.scaf.agp.gz # qualFiles none dbDbSpeciesDir panda taxId 9646 makeGenomeDb.pl -workhorse=hgwdev -stop=seq ailMel1.config.ra > seq.out 2>&1 # real 4m4.783s makeGenomeDb.pl -start=agp -stop=agp ailMel1.config.ra > agp.out 2>&1 # real 0m20.968s makeGenomeDb.pl -start=db -stop=db ailMel1.config.ra > db.out 2>&1 # real 6m28.168s makeGenomeDb.pl -start=dbDb -stop=dbDb ailMel1.config.ra > dbDb.out 2>&1 makeGenomeDb.pl -start=trackDb -stop=trackDb ailMel1.config.ra > trackDb.out 2>&1 ########################################################################## # running repeat masker (DONE - 2010-02-02 - Hiram) mkdir /hive/data/genomes/ailMel1/bed/repeatMasker cd /hive/data/genomes/ailMel1/bed/repeatMasker doRepeatMasker.pl -buildDir=`pwd` -noSplit -bigClusterHub=swarm \ -workhorse=hgwdev ailMel1 > do.log 2>&1 cat faSize.rmsk.txt # 2299509015 bases (54196184 N's 2245312831 real 1371210941 upper 874101890 # lower) in 81467 sequences in 1 files # %38.01 masked total, %38.93 masked real ########################################################################## # running simple repeat (DONE - 2010-02-02 - Hiram) mkdir /hive/data/genomes/ailMel1/bed/simpleRepeat cd /hive/data/genomes/ailMel1/bed/simpleRepeat time doSimpleRepeat.pl -buildDir=`pwd` -smallClusterHub=swarm \ -workhorse=hgwdev ailMel1 > do.log 2>&1 & # real 103m30.880s cat fb.simpleRepeat # 25064758 bases of 2245312831 (1.116%) in intersection cd /hive/data/genomes/ailMel1 twoBitMask ailMel1.rmsk.2bit \ -add bed/simpleRepeat/trfMask.bed ailMel1.2bit # you can safely ignore the warning about fields >= 13 twoBitToFa ailMel1.2bit stdout | faSize stdin > faSize.ailMel1.2bit.txt cat faSize.ailMel1.2bit.txt # 2299509015 bases (54196184 N's 2245312831 real 1370617891 upper # 874694940 lower) in 81467 sequences in 1 files # %38.04 masked total, %38.96 masked real rm /gbdb/ailMel1/ailMel1.2bit ln -s `pwd`/ailMel1.2bit /gbdb/ailMel1/ailMel1.2bit ########################################################################## # 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 ("ailMel1", "blat12", "17802", "1", "0"); \ INSERT INTO blatServers (db, host, port, isTrans, canPcr) \ VALUES ("ailMel1", "blat12", "17803", "0", "1");' \ hgcentraltest # test it with some sequence ######################################################################### # MAKE 11.OOC FILES FOR BLAT (DONE - 2010-02-03 - Hiram) ssh kolossus # numerator is ailMel1 gapless bases as reported by faSize # denominator is hg17 gapless bases as reported by featureBits, # 1024 is threshold used for human -repMatch: calc \( 2245312831 / 2897310462 \) \* 1024 # ( 2245312831 / 2897310462 ) * 1024 = 793.563675 # ==> use -repMatch=750 according to size scaled down from 1024 for human. # and rounded down to nearest 50 cd /hive/data/genomes/ailMel1 blat ailMel1.2bit /dev/null /dev/null -tileSize=11 \ -makeOoc=jkStuff/ailMel1.11.ooc -repMatch=750 # Wrote 25467 overused 11-mers to jkStuff/ailMel1.11.ooc mkdir /hive/data/staging/data/ailMel cp -p ailMel1.2bit chrom.sizes jkStuff/ailMel1.11.ooc \ /hive/data/staging/data/ailMel1 ######################################################################### # 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 after canFam2: # ailMel1 (panda) ailMel1.serverGenome = /hive/data/genomes/ailMel1/ailMel1.2bit ailMel1.clusterGenome = /scratch/data/ailMel1/ailMel1.2bit ailMel1.lift = no ailMel1.ooc = /scratch/data/ailMel1/ailMel1.11.ooc ailMel1.refseq.mrna.native.pslCDnaFilter = ${ordered.refseq.mrna.native.pslCDnaFilter} ailMel1.refseq.mrna.xeno.pslCDnaFilter = ${ordered.refseq.mrna.xeno.pslCDnaFilter} ailMel1.genbank.mrna.native.pslCDnaFilter = ${ordered.genbank.mrna.native.pslCDnaFilter} ailMel1.genbank.mrna.xeno.pslCDnaFilter = ${ordered.genbank.mrna.xeno.pslCDnaFilter} ailMel1.genbank.est.native.pslCDnaFilter = ${ordered.genbank.est.native.pslCDnaFilter} ailMel1.refseq.mrna.native.load = yes ailMel1.refseq.mrna.xeno.load = yes ailMel1.genbank.mrna.xeno.load = yes ailMel1.downloadDir = ailMel1 ailMel1.perChromTables = no # Edit src/lib/gbGenome.c to add new species. With these four lines: # static char *ailMelNames[] = {"Ailuropoda melanoleuca", "Canis familiaris", # "Canis sp.", "Canis lupus familiaris", # "Canis lupus", NULL}; # ... later ... # {"ailMel", ailMelNames}, cvs ci -m "adding ailMel - Panda with Dog specifications" src/lib/gbGenome.c make install-server 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 ailMel1 & # real 652m27.871s # var/build/logs/2010.04.09-15:37:11.ailMel1.initalign.log ssh hgwdev cd /cluster/data/genbank time ./bin/gbDbLoadStep -drop -initialLoad ailMel1 & # logFile: var/dbload/hgwdev/logs/2010.04.12-10:13:40.dbload.log # real 30m57.886s # enable daily alignment and update of hgwdev cd ~/kent/src/hg/makeDb/genbank cvsup # add oryCun2 to: etc/align.dbs etc/hgwdev.dbs cvs ci -m "Adding ailMel1 - Panda - Ailuropoda melanoleuca" \ etc/align.dbs etc/hgwdev.dbs make etc-update # done - 2010-04-12 - Hiram ############################################################################ # reset position to RHO location as found from blat of hg19 RHO gene hgsql -e \ 'update dbDb set defaultPos="GL192818.1:558576-566855" where name="ailMel1";' \ hgcentraltest ############################################################################ # LASTZ Opossum monDom5 (DONE - 2010-02-05 - Hiram) # original alignment cd /hive/data/genomes/monDom5/bed/lastzAilMel1.2010-02-04 cat fb.monDom5.chainAilMel1Link.txt # 223510659 bases of 3501660299 (6.383%) in intersection # and for the swap mkdir /hive/data/genomes/ailMel1/bed/blastz.monDom5.swap cd /hive/data/genomes/ailMel1/bed/blastz.monDom5.swap time nice -n +19 doBlastzChainNet.pl -verbose=2 \ /hive/data/genomes/monDom5/bed/lastzAilMel1.2010-02-04/DEF \ -swap -noLoadChainSplit -syntenicNet \ -workhorse=hgwdev -smallClusterHub=memk -bigClusterHub=swarm \ -chainMinScore=5000 -chainLinearGap=loose > swap.log 2>&1 & # real 69m35.464s cat fb.ailMel1.chainMonDom5Link.txt # 211209682 bases of 2245312831 (9.407%) in intersection ######################################################################### # LASTZ Mouse Mm9 (DONE - 2010-02-05 - Hiram) # original alignment cd /hive/data/genomes/mm9/bed/lastzAilMel1.2010-02-04 cat fb.mm9.chainAilMel1Link.txt # 749595031 bases of 2620346127 (28.607%) in intersection # and for the swap mkdir /hive/data/genomes/ailMel1/bed/blastz.mm9.swap cd /hive/data/genomes/ailMel1/bed/blastz.mm9.swap time doBlastzChainNet.pl -verbose=2 \ /hive/data/genomes/mm9/bed/lastzAilMel1.2010-02-04/DEF \ -swap -noLoadChainSplit -bigClusterHub=swarm -smallClusterHub=memk \ -workhorse=hgwdev \ -chainMinScore=3000 -chainLinearGap=medium > swap.log 2>&1 & # real 54m57.140s cat fb.ailMel1.chainMm9Link.txt # 739076250 bases of 2245312831 (32.916%) in intersection ####################################################################### # felCat4 Cat BLASTZ/CHAIN/NET (DONE - 2010-06-07 - Chin) screen # use a screen to manage this multi-day job mkdir /hive/data/genomes/ailMel1/bed/lastzFelCat4.2010-06-07 cd /hive/data/genomes/ailMel1/bed/lastzFelCat4.2010-06-07 cat << '_EOF_' > DEF # panda vs. cat # maximum M allowed with lastz is only 254 BLASTZ_M=254 # QUERY: Panda ailMel1 SEQ2_DIR=/scratch/data/ailMel1/ailMel1.2bit SEQ2_LEN=/scratch/data/ailMel1/chrom.sizes SEQ2_CHUNK=10000000 SEQ2_LAP=10000 SEQ2_LIMIT=300 # TARGET: Cat (felCat4) SEQ1_DIR=/scratch/data/felCat4/felCat4.2bit SEQ1_LEN=/scratch/data/felCat4/chrom.sizes SEQ1_LIMIT=50 SEQ1_CHUNK=20000000 SEQ1_LAP=10000 BASE=/hive/data/genomes/ailMel1/bed/lastzFelCat4.2010-06-07 TMPDIR=/scratch/tmp '_EOF_' # << this line keeps emacs coloring happy time nice -n +19 doBlastzChainNet.pl -verbose=2 \ `pwd`/DEF \ -syntenicNet -noDbNameCheck \ -chainMinScore=3000 -chainLinearGap=medium \ -workhorse=hgwdev -smallClusterHub=memk -bigClusterHub=swarm \ > do.log 2>&1 & # real 2761m37.696s time nice -n +19 doBlastzChainNet.pl -verbose=2 \ `pwd`/DEF \ -continue blastz \ -syntenicNet -noDbNameCheck \ -chainMinScore=3000 -chainLinearGap=medium \ -workhorse=hgwdev -smallClusterHub=memk -bigClusterHub=swarm \ > blastz.log 2>&1 & # real 2909m58.289s # *** All done ! Elapsed time: 2909m58s # *** Make sure that goldenPath/felCat4/vsAilMel1/README.txt is accurate. # *** Add {chain,net}AilMel1 tracks to trackDb.ra if necessary. cat fb.felCat4.chainAilMel1Link.txt # 1503647735 bases of 1990635005 (75.536%) in intersection # run rbest 2010-06-28 (screen 10498 Job 10560) time doRecipBest.pl -workhorse=hgwdev -buildDir=`pwd` \ felCat4 ailMel1 > rbest.log 2>&1 & # real 463m40.359s #swap (Copied from felCat4.txt) mkdir /hive/data/genomes/felCat4/bed/blastz.ailMel1.swap cd /hive/data/genomes/felCat4/bed/blastz.ailMel1.swap time nice -n +19 doBlastzChainNet.pl -verbose=2 \ /hive/data/genomes/ailMel1/bed/lastzFelCat4.2010-06-07/DEF \ -swap -syntenicNet -noDbNameCheck \ -workhorse=hgwdev -smallClusterHub=encodek -bigClusterHub=swarm \ -chainMinScore=3000 -chainLinearGap=medium > swap.log 2>&1 & # real 213m28.468s # the swap run failed due to # mod of # /usr/local/apache/htdocs-hgdownload/goldenPath/ailMel1/liftOver/md5sum.txt # was set to # -r--r--r-- 1 382 Apr 26 17:17 md5sum.txt # correct this to 664 by cluster-admin # and continue from download step time nice -n +19 doBlastzChainNet.pl -verbose=2 \ /hive/data/genomes/ailMel1/bed/lastzFelCat4.2010-06-07/DEF \ -continue=download \ -swap -syntenicNet -noDbNameCheck \ -workhorse=hgwdev -smallClusterHub=encodek -bigClusterHub=swarm \ -chainMinScore=3000 -chainLinearGap=medium > download.log 2>&1 & # real 32m9.204s # *** All done ! Elapsed time: 32m9s # *** Make sure that goldenPath/ailMel1/vsFelCat4/README.txt is accurate. # *** Add {chain,net}FelCat4 tracks to trackDb.ra if necessary. cd /hive/data/genomes/ailMel1/bed/blastz.felCat4.swap cat fb.ailMel1.chainFelCat4Link.txt # 1507273252 bases of 2245312831 (67.130%) in intersection # Due to the SEQ1 and SEQ2 spec'ed in # /hive/data/genomes/ailMel1/bed/lastzFelCat4.2010-06-07/DEF # need to copy the dDirectory from ailMel1 to felCatv17e # without this, 6way step will have problem # cd /hive/data/genomes/ailMel1/bed/blastz.felCat4.swap] # cp -pr * /hive/data/genomes/felCat4/bed/blastz.ailMel1.swap/ # cd /hive/data/genomes/felCat4/bed # ln -s lastzAilMel1.2010-06-07 lastz.ailMel1 # mkdir /hive/data/genomes/felCat4/bed/lastzAilMel1.2010-06-07 # cd /hive/data/genomes/ailMel1/bed/lastzFelCat4.2010-06-07 # cp -rp * /hive/data/genomes/felCat4/bed/lastzAilMel1.2010-06-07/ # cd /hive/data/genomes/felCat4/bed # ln -s /hive/data/genomes/felCat4/bed/lastzAilMel1.2010-06-07 # lastz.ailMel1 # cd /hive/data/genomes/ailMel1/bed # ln -s blastz.felCat4.swap lastz.felCat4 ######################################################################### # 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=ailMel1 -all all.joiner mkdir /hive/data/genomes/ailMel1/goldenPath cd /hive/data/genomes/ailMel1/goldenPath makeDownloads.pl ailMel1 > do.log 2>&1 # now ready for pushQ entry mkdir /hive/data/genomes/ailMel1/pushQ cd /hive/data/genomes/ailMel1/pushQ makePushQSql.pl ailMel1 > ailMel1.pushQ.sql 2> stderr.out # check for errors in stderr.out, some are OK, e.g.: # WARNING: ailMel1 does not have seq # WARNING: ailMel1 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 ailMel1.pushQ.sql hgwbeta:/tmp ssh hgwbeta cd /tmp hgsql qapushq < ailMel1.pushQ.sql # in that pushQ entry walk through each entry and see if the # sizes will set properly ############################################################################ # Marking *all* gaps - they are not all in the AGP file # (DONE - 2010-03-25 - Hiram) mkdir /hive/data/genomes/ailMel1/bed/allGaps cd /hive/data/genomes/ailMel1/bed/allGaps time nice -n +19 findMotif -motif=gattaca -verbose=4 \ -strand=+ ../../ailMel1.2bit > findMotif.txt 2>&1 # real 1m12.153s grep "^#GAP " findMotif.txt | sed -e "s/^#GAP //" > allGaps.bed featureBits ailMel1 -not gap -bed=notGap.bed # 2245312831 bases of 2245312831 (100.000%) in intersection featureBits ailMel1 allGaps.bed notGap.bed -bed=new.gaps.bed # 0 bases of 2245312831 (0.000%) in intersection # real 84m58.382s # no new gaps ######################################################################### # lastz Human hg19 swap (DONE - 2010-02-05 - Hiram) # original alignment cd /hive/data/genomes/hg19/bed/lastzAilMel1.2010-02-04 cat fb.hg19.chainAilMel1Link.txt # 1453400264 bases of 2897316137 (50.164%) in intersection # and the swap mkdir /hive/data/genomes/ailMel1/bed/blastz.hg19.swap cd /hive/data/genomes/ailMel1/bed/blastz.hg19.swap time nice -n +19 doBlastzChainNet.pl -verbose=2 \ /hive/data/genomes/hg19/bed/lastzAilMel1.2010-02-04/DEF \ -swap -noLoadChainSplit -syntenicNet \ -workhorse=hgwdev -smallClusterHub=memk -bigClusterHub=swarm \ -chainMinScore=3000 -chainLinearGap=medium > swap.log 2>&1 & real 124m14.393s cat fb.ailMel1.chainHg19Link.txt # 1411953704 bases of 2245312831 (62.884%) in intersection ######################################################################### # lastz Dog canFam2 swap (DONE - 2010-02-05 - Hiram) # original alignment cd /hive/data/genomes/canFam2/bed/lastzAilMel1.2010-02-04 cat fb.canFam2.chainAilMel1Link.txt # 1791212709 bases of 2384996543 (75.103%) in intersection # and the swap mkdir /hive/data/genomes/ailMel1/bed/blastz.canFam2.swap cd /hive/data/genomes/ailMel1/bed/blastz.canFam2.swap time doBlastzChainNet.pl -verbose=2 \ /hive/data/genomes/canFam2/bed/lastzAilMel1.2010-02-04/DEF \ -swap -noLoadChainSplit -workhorse=hgwdev -bigClusterHub=pk \ -smallClusterHub=memk -chainMinScore=3000 -chainLinearGap=medium \ > swap.log 2>&1 & # real 128m41.005s cat fb.ailMel1.chainCanFam2Link.txt # 1788107935 bases of 2245312831 (79.637%) in intersection ######################################################################### # lastz Rat rn4 swap (DONE - 2010-02-18 - Hiram) # original alignment cd /hive/data/genomes/rn4/bed/lastzAilMel1.2010-02-04 cat fb.rn4.chainAilMel1Link.txt # 708197812 bases of 2571531505 (27.540%) in intersection # and the swap mkdir /hive/data/genomes/ailMel1/bed/blastz.rn4.swap cd /hive/data/genomes/ailMel1/bed/blastz.rn4.swap time nice -n +19 doBlastzChainNet.pl -verbose=2 \ /hive/data/genomes/rn4/bed/lastzAilMel1.2010-02-04/DEF \ -swap -syntenicNet \ -workhorse=hgwdev -smallClusterHub=memk -bigClusterHub=swarm \ -chainMinScore=3000 -chainLinearGap=medium > swap.log 2>&1 & # real 73m37.751s cat fb.ailMel1.chainRn4Link.txt # 695366144 bases of 2245312831 (30.970%) in intersection ######################################################################### # HUMAN (hg18) PROTEINS TRACK (DONE braney 2010-04-19) # bash if not using bash shell already cd /cluster/data/ailMel1 mkdir /cluster/data/ailMel1/blastDb awk '{if ($2 > 1000000) print $1}' chrom.sizes > 1meg.lst twoBitToFa -seqList=1meg.lst ailMel1.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 ailMel1.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 # 2757 mkdir -p /cluster/data/ailMel1/bed/tblastn.hg18KG cd /cluster/data/ailMel1/bed/tblastn.hg18KG echo ../../blastDb/*.nsq | xargs ls -S | sed "s/\.nsq//" > query.lst wc -l query.lst # 2757 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/2757) = 289.303826 mkdir -p kgfa split -l 289 /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 # 128 128 1664 kg.lst mkdir -p blastOut for i in `cat kg.lst`; do mkdir blastOut/`basename $i .fa`; done tcsh cd /cluster/data/ailMel1/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/ailMel1/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/ailMel1/bed/tblastn.hg18KG gensub2 query.lst kg.lst blastGsub blastSpec para create blastSpec # para try, check, push, check etc. para time # Completed: 352896 of 352896 jobs # CPU time in finished jobs: 13166354s 219439.24m 3657.32h 152.39d 0.418 y # IO & Wait Time: 1995282s 33254.69m 554.24h 23.09d 0.063 y # Average job time: 43s 0.72m 0.01h 0.00d # Longest finished job: 207s 3.45m 0.06h 0.00d # Submission to last job: 15726s 262.10m 4.37h 0.18d ssh swarm cd /cluster/data/ailMel1/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=150000 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: 128 of 128 jobs #CPU time in finished jobs: 1401s 23.35m 0.39h 0.02d 0.000 y #IO & Wait Time: 47407s 790.12m 13.17h 0.55d 0.002 y #Average job time: 381s 6.36m 0.11h 0.00d #Longest finished job: 466s 7.77m 0.13h 0.01d #Submission to last job: 470s 7.83m 0.13h 0.01d cd /cluster/data/ailMel1/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: 64037 failed: 0 errors: 0 # load table ssh hgwdev cd /cluster/data/ailMel1/bed/tblastn.hg18KG hgLoadPsl ailMel1 blastHg18KG.psl # check coverage featureBits ailMel1 blastHg18KG # 32308509 bases of 2245312831 (1.439%) in intersection featureBits ailMel1 blastHg18KG xenoRefGene -enrichment # blastHg18KG 1.439%, xenoRefGene 2.186%, both 1.137%, cover 79.05%, enrich 36.17x rm -rf blastOut #end tblastn