# for emacs: -*- mode: sh; -*- # Xenopus tropicalis # ftp.ncbi.nlm.nih.gov/genbank/genomes/Eukaryotes/vertebrates_other/ # Xenopus_tropicalis/v4.2/ # # http://www.ncbi.nlm.nih.gov/Traces/wgs/?val=AAMC01 # WGS: AAMC01000001:AAMC01190823 ########################################################################## # Download sequence (DONE - 2011-09-07 - Hiram) mkdir -p /hive/data/genomes/xenTro3/genbank cd /hive/data/genomes/xenTro3/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_other/Xenopus_tropicalis/v4.2/*" faSize Primary_Assembly/unplaced_scaffolds/FASTA/unplaced.scaf.fa.gz 1511717716 bases (153400444 N's 1358317272 real 1358317272 upper 0 lower) in 19549 sequences in 1 files ########################################################################### # fixup contig names a little bit (DONE - 2011-09-07 - Hiram) cd /hive/data/genomes/xenTro3/genbank # removing the .1 from all the GL names zcat Primary_Assembly/unplaced_scaffolds/FASTA/unplaced.scaf.fa.gz \ | sed -e "s/^>gi.[0-9][0-9]*.gb.GL/>GL/; s/\.1| .*//;" \ | gzip -c > xenTro3.fa.gz # real 6m24.993s zcat Primary_Assembly/unplaced_scaffolds/AGP/unplaced.scaf.agp.gz \ | sed -e "s/^GL1\([0-9]*\)\.1/GL1\1/" | gzip -c > xenTro3.agp.gz ########################################################################### # Initial database build (DONE - 2011-09-08 - Hiram) cd /hive/data/genomes/xenTro3/genbank cat << '_EOF_' > xenTro3.config.ra # Config parameters for makeGenomeDb.pl: db xenTro3 clade vertebrate genomeCladePriority 80 scientificName Xenopus (Silurana) tropicalis commonName western clawed frog assemblyDate Nov. 2009 assemblyLabel US DOE Joint Genome Institute (JGI-PGF) (GCA_000004195.1) assemblyShortLabel X. tropicalis v4.2 orderKey 444 mitoAcc NC_006839 fastaFiles /cluster/data/xenTro3/genbank/xenTro3.fa.gz agpFiles /cluster/data/xenTro3/genbank/xenTro3.agp.gz dbDbSpeciesDir xenTro taxId 8364 '_EOF_' # << happy emacs # verify sequence and agp are OK makeGenomeDb.pl -stop=agp xenTro3.config.ra > agp.log 2>&1 time makeGenomeDb.pl -continue=db xenTro3.config.ra > db.log 2>&1 # real 12m7.144s ########################################################################## # running repeat masker (DONE - 2011-09-08 - Hiram) mkdir /hive/data/genomes/xenTro3/bed/repeatMasker cd /hive/data/genomes/xenTro3/bed/repeatMasker time doRepeatMasker.pl -buildDir=`pwd` -noSplit \ -species "Xenopus tropicalis" \ -bigClusterHub=swarm -dbHost=hgwdev -workhorse=hgwdev \ -smallClusterHub=memk xenTro3 > do.log 2>&1 & # real 35m31.331s cat faSize.rmsk.txt # 1511735326 bases (153400444 N's 1358334882 real 1025832944 upper # 332501938 lower) in 19550 sequences in 1 files # %21.99 masked total, %24.48 masked real grep -i versi do.log # RepeatMasker version development-$Id: RepeatMasker,v 1.25 2010/09/08 21:32:26 angie Exp $ # June 30 2010 (open-3-2-9) version of RepeatMasker featureBits -countGaps xenTro3 rmsk # 333219656 bases of 1511735326 (22.042%) in intersection # why is it different than the faSize above ? # because rmsk masks out some N's as well as bases, the count above # separates out the N's from the bases, it doesn't show lower case N's ########################################################################## # running simple repeat (DONE - 2011-09-08 - Hiram) mkdir /hive/data/genomes/xenTro3/bed/simpleRepeat cd /hive/data/genomes/xenTro3/bed/simpleRepeat time doSimpleRepeat.pl -buildDir=`pwd` -bigClusterHub=swarm \ -dbHost=hgwdev -workhorse=hgwdev -smallClusterHub=memk \ xenTro3 > do.log 2>&1 & # about 1 hour 3 minutes # real 23m18.960s # real 198m33.953s cat fb.simpleRepeat # 103101528 bases of 1358346944 (7.590%) in intersection # add to rmsk after it is done: cd /hive/data/genomes/xenTro3 twoBitMask xenTro3.rmsk.2bit \ -add bed/simpleRepeat/trfMask.bed xenTro3.2bit # you can safely ignore the warning about fields >= 13 twoBitToFa xenTro3.2bit stdout | faSize stdin > faSize.xenTro3.2bit.txt cat faSize.xenTro3.2bit.txt # 1511735326 bases (153400444 N's 1358334882 real 1024824487 upper # 333510395 lower) in 19550 sequences in 1 files # %22.06 masked total, %24.55 masked real rm /gbdb/xenTro3/xenTro3.2bit ln -s `pwd`/xenTro3.2bit /gbdb/xenTro3/xenTro3.2bit ######################################################################### # Verify all gaps are marked, add any N's not in gap as type 'other' # (DONE - 2011-09-08 - Hiram) mkdir /hive/data/genomes/xenTro3/bed/gap cd /hive/data/genomes/xenTro3/bed/gap time nice -n +19 findMotif -motif=gattaca -verbose=4 \ -strand=+ ../../xenTro3.unmasked.2bit > findMotif.txt 2>&1 # real 0m15.269s grep "^#GAP " findMotif.txt | sed -e "s/^#GAP //" > allGaps.bed featureBits xenTro3 -not gap -bed=notGap.bed # 1358346944 bases of 1358346944 (100.000%) in intersection featureBits xenTro3 allGaps.bed notGap.bed -bed=new.gaps.bed # 12062 bases of 1358346944 (0.001%) in intersection # what is the highest index in the existing gap table: hgsql -N -e "select ix from gap;" xenTro3 | sort -n | tail -1 # 1466 cat << '_EOF_' > mkGap.pl #!/bin/env perl use strict; use warnings; my $ix=`hgsql -N -e "select ix from gap;" xenTro3 | sort -n | tail -1`; chomp $ix; open (FH,") { my ($chrom, $chromStart, $chromEnd, $rest) = split('\s+', $line); ++$ix; printf "%s\t%d\t%d\t%d\tN\t%d\tother\tyes\n", $chrom, $chromStart, $chromEnd, $ix, $chromEnd-$chromStart; } close (FH); '_EOF_' # << happy emacs chmod +x ./mkGap.pl ./mkGap.pl > other.bed featureBits -countGaps xenTro3 other.bed # 12062 bases of 1511735326 (0.001%) in intersection wc -l other.bed # 11242 hgLoadBed -sqlTable=$HOME/kent/src/hg/lib/gap.sql \ -noLoad xenTro3 otherGap other.bed # starting with this many hgsql -e "select count(*) from gap;" xenTro3 # 171274 hgsql xenTro3 -e 'load data local infile "bed.tab" into table gap;' # result count: hgsql -e "select count(*) from gap;" xenTro3 # 182516 # == 171274 + 11242 # verify we aren't adding gaps where gaps already exist # this would output errors if that were true: gapToLift -minGap=1 xenTro3 nonBridged.lift -bedFile=nonBridged.bed # see example in danRer7.txt # there are no non-bridged gaps here: hgsql -N -e "select bridge from gap;" xenTro3 | sort | uniq -c # 182516 yes ########################################################################## ## WINDOWMASKER (DONE - 2011-09-08 - Hiram) mkdir /hive/data/genomes/xenTro3/bed/windowMasker cd /hive/data/genomes/xenTro3/bed/windowMasker time nice -n +19 doWindowMasker.pl -buildDir=`pwd` -workhorse=hgwdev \ -dbHost=hgwdev xenTro3 > do.log 2>&1 & # about 1h 52m # Masking statistics twoBitToFa xenTro3.wmsk.2bit stdout | faSize stdin # 1511735326 bases (153400444 N's 1358334882 real 845383044 upper # 512951838 lower) in 19550 sequences in 1 files # %33.93 masked total, %37.76 masked real twoBitToFa xenTro3.wmsk.sdust.2bit stdout | faSize stdin # 1511735326 bases (153400444 N's 1358334882 real 835449906 upper # 522884976 lower) in 19550 sequences in 1 files # %34.59 masked total, %38.49 masked real hgLoadBed xenTro3 windowmaskerSdust windowmasker.sdust.bed.gz # Loaded 7010382 elements of size 3 featureBits -countGaps xenTro3 windowmaskerSdust # 811865302 bases of 1799143587 (45.125%) in intersection # eliminate the gaps from the masking featureBits xenTro3 -not gap -bed=notGap.bed # 1701353770 bases of 1701353770 (100.000%) in intersection time nice -n +19 featureBits xenTro3 windowmaskerSdust notGap.bed \ -bed=stdout | gzip -c > cleanWMask.bed.gz # 714075485 bases of 1701353770 (41.971%) in intersection # reload track to get it clean hgLoadBed xenTro3 windowmaskerSdust cleanWMask.bed.gz # Loaded 8199653 elements of size 4 featureBits -countGaps xenTro3 windowmaskerSdust # 714075485 bases of 1799143587 (39.690%) in intersection # mask the sequence with this clean mask zcat cleanWMask.bed.gz \ | twoBitMask ../../xenTro3.unmasked.2bit stdin \ -type=.bed xenTro3.cleanWMSdust.2bit twoBitToFa xenTro3.cleanWMSdust.2bit stdout | faSize stdin \ > xenTro3.cleanWMSdust.faSize.txt cat xenTro3.cleanWMSdust.faSize.txt # 1799143587 bases (97789820 N's 1701353767 real 987278282 upper # 714075485 lower) in 6457 sequences in 1 files # %39.69 masked total, %41.97 masked real # how much does this window masker and repeat masker overlap: featureBits -countGaps xenTro3 rmsk windowmaskerSdust # 280229203 bases of 1511735326 (18.537%) in intersection ######################################################################### # MASK SEQUENCE WITH WM+TRF (DONE - 2011-09-09 - Hiram) cd /hive/data/genomes/xenTro3 twoBitMask -add bed/windowMasker/xenTro3.cleanWMSdust.2bit \ bed/simpleRepeat/trfMask.bed xenTro3.2bit # safe to ignore the warnings about BED file with >=13 fields twoBitToFa xenTro3.2bit stdout | faSize stdin > faSize.xenTro3.txt cat faSize.xenTro3.txt # 1799143587 bases (97789820 N's 1701353767 real 987044085 upper # 714309682 lower) in 6457 sequences in 1 files # %39.70 masked total, %41.98 masked real # create symlink to gbdb ssh hgwdev rm /gbdb/xenTro3/xenTro3.2bit ln -s `pwd`/xenTro3.2bit /gbdb/xenTro3/xenTro3.2bit ######################################################################## # MAKE 11.OOC FILE FOR BLAT/GENBANK (DONE - 2011-09-09 - Hiram) # Use -repMatch=650, based on size -- for human we use 1024 # use the "real" number from the faSize measurement, # hg19 is 2897316137, calculate the ratio factor for 1024: calc \( 1358317272 / 2897316137 \) \* 1024 # ( 1358317272 / 2897316137 ) * 1024 = 480.070804 # round up to 500 cd /hive/data/genomes/xenTro3 blat xenTro3.2bit /dev/null /dev/null -tileSize=11 \ -makeOoc=jkStuff/xenTro3.11.ooc -repMatch=500 # Wrote 29991 overused 11-mers to jkStuff/xenTro3.11.ooc # copy all of this stuff to the klusters: # there are no non-bridged gaps hgsql -N -e "select bridge from gap;" xenTro3 | sort | uniq -c # 182516 yes # cd /hive/data/genomes/xenTro3/jkStuff # gapToLift xenTro3 nonBridged.lift -bedFile=nonBridged.bed cd /hive/data/genomes/xenTro3 mkdir /hive/data/staging/data/xenTro3 cp -p jkStuff/xenTro3.11.ooc chrom.sizes \ xenTro3.2bit /hive/data/staging/data/xenTro3 # request rsync copy from cluster admin ######################################################################### # AUTO UPDATE GENBANK (DONE - 2011-09-11 - Hiram) # examine the file: /cluster/data/genbank/data/organism.lst # for your species to see what counts it has for: # organism mrnaCnt estCnt refSeqCnt # Xenopus (Silurana) tropicalis 18786 1271379 8783 # to decide which "native" mrna or ests you want to specify in genbank.conf # this appears that xenTro3 has plenty of native est's ssh hgwdev cd $HOME/kent/src/hg/makeDb/genbank git pull # edit etc/genbank.conf to add xenTro3 before anoCar1 and commit to GIT # xenTro3 (X. tropicalis) 19550 scaffolds xenTro3.serverGenome = /hive/data/genomes/xenTro3/xenTro3.2bit xenTro3.clusterGenome = /scratch/data/xenTro3/xenTro3.2bit xenTro3.ooc = /scratch/data/xenTro3/xenTro3.11.ooc xenTro3.lift = no xenTro3.refseq.mrna.native.pslCDnaFilter = ${lowCover.refseq.mrna.native.pslCDnaFilter} xenTro3.refseq.mrna.xeno.pslCDnaFilter = ${lowCover.refseq.mrna.xeno.pslCDnaFilter} xenTro3.genbank.mrna.native.pslCDnaFilter = ${lowCover.genbank.mrna.native.pslCDnaFilter} xenTro3.genbank.mrna.xeno.pslCDnaFilter = ${lowCover.genbank.mrna.xeno.pslCDnaFilter} xenTro3.genbank.est.native.pslCDnaFilter = ${lowCover.genbank.est.native.pslCDnaFilter} xenTro3.refseq.mrna.native.load = yes xenTro3.genbank.est.native.load = yes xenTro3.refseq.mrna.xeno.load = no xenTro3.genbank.mrna.xeno.load = no xenTro3.downloadDir = xenTro3 xenTro3.perChromTables = no xenTro3.mgc = yes # xenTro3.upstreamGeneTbl = mgcGenes # xenTro3.upstreamMaf = multiz9way # /hive/data/genomes/xenTro3/bed/multiz9way/species.lst # end of section added to etc/genbank.conf git commit -m "adding xenTro3 frog" genbank.conf git push make etc-update # ~/kent/src/hg/makeDb/genbank/src/lib/gbGenome.c already contains # anoCar genome information, if this is a new species, need to add stuff # there ssh hgwdev # used to do this on "genbank" machine screen # long running job managed in screen cd /cluster/data/genbank time nice -n +19 ./bin/gbAlignStep -initial xenTro3 & # var/build/logs/2011.09.12-09:03:09.xenTro3.initalign.log # real 360m20.657s # load database when finished ssh hgwdev cd /cluster/data/genbank time nice -n +19 ./bin/gbDbLoadStep -drop -initialLoad xenTro3 & # logFile: var/dbload/hgwdev/logs/2011.09.12-20:47:05.dbload.log # real 118m27.671s # enable daily alignment and update of hgwdev (DONE - 2011-09-13 - Hiram) cd ~/kent/src/hg/makeDb/genbank git pull # add xenTro3 to: etc/align.dbs etc/hgwdev.dbs git commit -m "Added xenTro3." etc/align.dbs etc/hgwdev.dbs git push make etc-update ######################################################################### # construct ctgPos2 table (DONE - 2011-09-12 - Hiram) mkdir /hive/data/genomes/xenTro3/bed/ctgPos2 cd /hive/data/genomes/xenTro3/bed/ctgPos2 cat << '_EOF_' > chromSizesToCtgPos.pl #!/usr/bin/env perl use strict; use warnings; my %chromSizes; open (FH, "<../../chrom.sizes") or die "can not read ../../chrom.sizes"; while (my $line = ) { chomp $line; my ($contig, $size) = split('\s', $line); $chromSizes{$contig} = $size; } close (FH); open (FH, "<../../genbank/Primary_Assembly/localID2acc") or die "can not read ../../genbank/Primary_Assembly/localID2acc"; while (my $line = ) { chomp $line; my ($scaffold, $contig) = split('\s', $line); $scaffold =~ s/XENTR//; $contig =~ s/.1$//; my $size = $chromSizes{$contig}; printf "%s\t%d\t%s\t0\t%d\tW\n", $scaffold, $size, $contig, $size; } close (FH); '_EOF_' # << happy emacs chmod +x chromSizesToCtgPos.pl ./chromSizesToCtgPos.pl > xenTro3.ctgPos2.tab echo -e "AY789013\t17610\tchrM\t0\t17610\tF" >> xenTro3.ctgPos2.tab # check how big are the names: awk '{print length($1)}' xenTro3.ctgPos2.tab | sort -rn | head -1 # 14 awk '{print length($3)}' xenTro3.ctgPos2.tab | sort -rn | head # 8 sed -e "s/20/14/; s/16/8/" $HOME/kent/src/hg/lib/ctgPos2.sql > ctgPos2.sql hgLoadSqlTab xenTro3 ctgPos2 ctgPos2.sql xenTro3.ctgPos2.tab ############################################################################ # running cpgIsland business (DONE - 2011-09-13 - Hiram) mkdir /hive/data/genomes/xenTro3/bed/cpgIsland cd /hive/data/genomes/xenTro3/bed/cpgIsland # use a previous binary for this program ln -s ../../../xenTro2/bed/cpgIsland/cpglh.exe . mkdir -p hardMaskedFa cut -f1 ../../chrom.sizes | while read C do echo ${C} twoBitToFa ../../xenTro3.2bit:$C stdout \ | maskOutFa stdin hard hardMaskedFa/${C}.fa done ssh swarm cd /hive/data/genomes/xenTro3/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 | egrep -v "^#seq|^total" | 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 chmod +x runOne gensub2 chr.list single template jobList para create jobList para try para check ... etc para time # Completed: 1133 of 1133 jobs # CPU time in finished jobs: 102s 1.70m 0.03h 0.00d 0.000 y # IO & Wait Time: 4469s 74.49m 1.24h 0.05d 0.000 y # Average job time: 4s 0.07m 0.00h 0.00d # Longest finished job: 10s 0.17m 0.00h 0.00d # Submission to last job: 485s 8.08m 0.13h 0.01d # Transform cpglh output to bed + catDir results | awk '{ $2 = $2 - 1; width = $3 - $2; printf("%s\t%d\t%s\t%s %s\t%s\t%s\t%0.0f\t%0.1f\t%s\t%s\n", $1, $2, $3, $5,$6, width, $6, width*$7*0.01, 100.0*2*$6/width, $7, $9); }' > cpgIsland.bed # verify longest unique chrom name: cut -f1 cpgIsland.bed | awk '{print length($0)}' | sort -rn | head -1 # 8 # update the length 14 in the template to be 16: sed -e "s/14/8/" $HOME/kent/src/hg/lib/cpgIslandExt.sql > cpgIslandExt.sql cd /hive/data/genomes/xenTro3/bed/cpgIsland hgLoadBed xenTro3 cpgIslandExt -tab -sqlTable=cpgIslandExt.sql cpgIsland.bed # Loaded 15341 elements of size 10 featureBits xenTro3 cpgIslandExt # 9355659 bases of 1358334882 (0.689%) in intersection # there should be no output from checkTableCoords: checkTableCoords -verboseBlocks -table=cpgIslandExt xenTro3 # cleanup rm -fr hardMaskedFa ######################################################################### # GENSCAN GENE PREDICTIONS (DONE - 2011-09-13 - Hiram) mkdir /hive/data/genomes/xenTro3/bed/genscan cd /hive/data/genomes/xenTro3/bed/genscan # use a previously existing genscan binary ln -s ../../../xenTro2/bed/genscan/hg3rdParty . # create hard masked .fa files mkdir -p hardMaskedFa cut -f1 ../../chrom.sizes | while read C do echo ${C} twoBitToFa ../../xenTro3.2bit:$C stdout \ | maskOutFa stdin hard hardMaskedFa/${C}.fa done # Generate a list file, genome.list, of all the hard-masked contig chunks: find ./hardMaskedFa/ -type f | sed -e 's#^./##' > genome.list wc -l genome.list # 19550 genome.list # Run on small cluster (more mem than big cluster). ssh swarm cd /hive/data/genomes/xenTro3/bed/genscan # Make 3 subdirectories for genscan to put their output files in mkdir gtf pep subopt # Create template file, template, for gensub2. For example (3-line file): cat << '_EOF_' > template #LOOP /cluster/bin/x86_64/gsBig {check in exists+ $(path1)} {check out exists gtf/$(root1).gtf} -trans={check out exists pep/$(root1).pep} -subopt={check out exists subopt/$(root1).bed} -exe=hg3rdParty/genscanlinux/genscan -par=hg3rdParty/genscanlinux/HumanIso.smat -tmp=/tmp -window=2400000 #ENDLOOP '_EOF_' # << emacs gensub2 genome.list single template jobList para create jobList para try para check ... etc... para time # Completed: 19550 of 19550 jobs # CPU time in finished jobs: 25171s 419.52m 6.99h 0.29d 0.001 y # IO & Wait Time: 51770s 862.83m 14.38h 0.60d 0.002 y # Average job time: 4s 0.07m 0.00h 0.00d # Longest finished job: 176s 2.93m 0.05h 0.00d # Submission to last job: 336s 5.60m 0.09h 0.00d # this did not work, runs out of file handles ? find ./gtf -type f | xargs -n 256 endsInLf -zeroOk # Concatenate results: cd /hive/data/genomes/xenTro3/bed/genscan find ./gtf -type f | xargs cat > genscan.gtf find ./pep -type f | xargs cat > genscan.pep find ./subopt -type f | xargs cat > genscanSubopt.bed # Load into the database (without -genePredExt because no frame info): # Don't load the Pep anymore -- redundant since it's from genomic. ssh hgwdev cd /hive/data/genomes/xenTro3/bed/genscan # to construct a local file with the genePred business: gtfToGenePred genscan.gtf genscan.gp # this produces exactly the same thing and loads the table: ldHgGene -gtf xenTro3 genscan genscan.gtf # Read 47412 transcripts in 293762 lines in 1 files # 47412 groups 9702 seqs 1 sources 1 feature types # 47412 gene predictions hgLoadBed xenTro3 genscanSubopt genscanSubopt.bed # Loaded 369728 elements of size 6 featureBits xenTro3 genscan # 62125158 bases of 1358334882 (4.574%) in intersection # previously: featureBits xenTro2 genscan # 62149060 bases of 1359412157 (4.572%) in intersection featureBits xenTro1 genscan # 63766483 bases of 1381238994 (4.617%) in intersection ######################################################################### # CREATE MICROSAT TRACK (DONE - 2011-09-15 - Hiram ssh hgwdev mkdir /cluster/data/xenTro3/bed/microsat cd /cluster/data/xenTro3/bed/microsat awk '($5==2 || $5==3) && $6 >= 15 && $8 == 100 && $9 == 0 {printf("%s\t%s\t%s\t%dx%s\n", $1, $2, $3, $6, $16);}' \ ../simpleRepeat/simpleRepeat.bed > microsat.bed hgLoadBed xenTro3 microsat microsat.bed # Loaded 13365 elements of size 4 ######################################################################### # create ucscToEnsembl scaffold name mapping (DONE - 2011-09-15 - Hiram) # this allows the "ensembl" blue bar button to appear mkdir /hive/data/genomes/xenTro3/bed/ucscToEnsembl cd /hive/data/genomes/xenTro3/bed/ucscToEnsembl awk '{printf "%s\t%s.1\n", $1, $1}' ../../chrom.sizes > 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(8)) ); '_EOF_' hgsql xenTro3 < ucscToEnsembl.sql hgsql xenTro3 \ -e 'LOAD DATA LOCAL INFILE "ucscToEnsembl.tab" INTO TABLE ucscToEnsembl' # verify the blue bar "ensembl" link is now available ######################################################################### # constructing downloads (DONE - 2011-09-16 - Hiram) cd /hive/data/genomes/xenTro3 time makeDownloads.pl -workhorse=hgwdev xenTro3 # real 14m44.729s # edit the README files in goldenPath/*/README.txt ######################################################################### # create pushQ entry (DONE - 2011-09-16 - Hiram) # first make sure all.joiner is up to date and has this new organism # a keys check should be clean: cd ~/kent/src/hg/makeDb/schema joinerCheck -database=xenTro3 -keys all.joiner mkdir /hive/data/genomes/xenTro3/pushQ cd /hive/data/genomes/xenTro3/pushQ makePushQSql.pl xenTro3 > xenTro3.sql 2> stderr.out # check stderr.out for no significant problems, it is common to see: # WARNING: hgwdev does not have /gbdb/xenTro3/wib/gc5Base.wib # WARNING: hgwdev does not have /gbdb/xenTro3/wib/quality.wib # WARNING: hgwdev does not have /gbdb/xenTro3/bbi/quality.bw # WARNING: xenTro3 does not have seq # WARNING: xenTro3 does not have extFile # which are no real problem # if some tables are not identified: # WARNING: Could not tell (from trackDb, all.joiner and hardcoded lists of # supporting and genbank tables) which tracks to assign these tables to: # ensGtp # ensPep # ensemblSource # ensemblToGeneName # ucscToEnsembl # put them in manually after loading the pushQ entry scp -p xenTro3.sql hgwbeta:/tmp ssh hgwbeta cd /tmp hgsql qapushq < xenTro3.sql ######################################################################### # HUMAN (hg18) PROTEINS TRACK (DONE 2011-09-17 braney ) # bash if not using bash shell already cd /cluster/data/xenTro3 mkdir /cluster/data/xenTro3/blastDb awk '{if ($2 > 1000000) print $1}' chrom.sizes > 1meg.lst twoBitToFa -seqList=1meg.lst xenTro3.2bit temp.fa faSplit gap temp.fa 1000000 blastDb/x -lift=blastDb.lft # 1353 pieces of 1353 written rm temp.fa 1meg.lst awk '{if ($2 <= 1000000) print $1}' chrom.sizes > less1meg.lst twoBitToFa -seqList=less1meg.lst xenTro3.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 # 1808 mkdir -p /cluster/data/xenTro3/bed/tblastn.hg18KG cd /cluster/data/xenTro3/bed/tblastn.hg18KG echo ../../blastDb/*.nsq | xargs ls -S | sed "s/\.nsq//" > query.lst wc -l query.lst # 1808 query.lst # we want around 250000 jobs calc `wc /cluster/data/hg18/bed/blat.hg18KG/hg18KG.psl | awk '{print $1}'`/\(250000/`wc query.lst | awk '{print $1}'`\) # 36727/(250000/1808) = 265.609664 mkdir -p kgfa split -l 266 /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 # 139 139 1807 kg.lst mkdir -p blastOut for i in `cat kg.lst`; do mkdir blastOut/`basename $i .fa`; done tcsh cd /cluster/data/xenTro3/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/xenTro3/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/xenTro3/bed/tblastn.hg18KG gensub2 query.lst kg.lst blastGsub blastSpec para create blastSpec # para try, check, push, check etc. para time # Completed: 251312 of 251312 jobs # CPU time in finished jobs: 8416814s 140280.23m 2338.00h 97.42d 0.267 y # IO & Wait Time: 1310570s 21842.83m 364.05h 15.17d 0.042 y # Average job time: 39s 0.65m 0.01h 0.00d # Longest finished job: 237s 3.95m 0.07h 0.00d # Submission to last job: 19145s 319.08m 5.32h 0.22d cd /cluster/data/xenTro3/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: 139 of 139 jobs # CPU time in finished jobs: 562s 9.37m 0.16h 0.01d 0.000 y # IO & Wait Time: 30974s 516.23m 8.60h 0.36d 0.001 y # Average job time: 227s 3.78m 0.06h 0.00d # Longest finished job: 299s 4.98m 0.08h 0.00d # Submission to last job: 310s 5.17m 0.09h 0.00d cd /cluster/data/xenTro3/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: 41882 failed: 0 errors: 0 # load table ssh hgwdev cd /cluster/data/xenTro3/bed/tblastn.hg18KG hgLoadPsl xenTro3 blastHg18KG.psl # check coverage featureBits xenTro3 blastHg18KG # 20459423 bases of 1358334882 (1.506%) in intersection featureBits xenTro3 blastHg18KG ensGene -enrichment # blastHg18KG 1.506%, ensGene 2.717%, both 1.337%, cover 88.77%, enrich 32.67x rm -rf blastOut #end tblastn ########################################################################## # BLATSERVERS ENTRY (DONE - 2011-09-20 - 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 ("xenTro3", "blat12", "17828", "1", "0"); \ INSERT INTO blatServers (db, host, port, isTrans, canPcr) \ VALUES ("xenTro3", "blat12", "17829", "0", "1");' \ hgcentraltest # test it with some sequence ############################################################################ # set default position the same as was xenTro2 via blat # (DONE - 2011-09-20 - Hiram) hgsql -e \ 'update dbDb set defaultPos="GL172663:565941-578374" where name="xenTro3";' \ hgcentraltest ############################################################################ # Swap lastz galGal3 Chicken (DONE - 2011-09-20 - Hiram) # original alignment cd /hive/data/genomes/galGal3/bed/lastzXenTro3.2011-09-20 cat fb.galGal3.chainXenTro3Link.txt # 54978665 bases of 1042591351 (5.273%) in intersection # and this swap mkdir /hive/data/genomes/xenTro3/bed/blastz.galGal3.swap cd /hive/data/genomes/xenTro3/bed/blastz.galGal3.swap time nice -n +19 doBlastzChainNet.pl -verbose=2 \ /hive/data/genomes/galGal3/bed/lastzXenTro3.2011-09-20/DEF \ -noLoadChainSplit -chainMinScore=5000 -chainLinearGap=loose \ -swap -workhorse=hgwdev -smallClusterHub=memk -bigClusterHub=swarm \ > swap.log 2>&1 & # real 13m36.708s cat fb.xenTro3.chainGalGal3Link.txt # 64250673 bases of 1358334882 (4.730%) in intersection ######################################################################### # Swap lastz anoCar2 Lizard (DONE - 2011-09-21 - Hiram) # original alignment cd /cluster/data/anoCar2/bed/blastz.xenTro3.2011-09-19 cat fb.anoCar2.chainXenTro3Link.txt # 86932478 bases of 1701353770 (5.110%) in intersection # and this swap mkdir /hive/data/genomes/xenTro3/bed/blastz.anoCar2.swap cd /hive/data/genomes/xenTro3/bed/blastz.anoCar2.swap time doBlastzChainNet.pl -verbose=2 \ /cluster/data/anoCar2/bed/blastz.xenTro3.2011-09-19/DEF \ -bigClusterHub=swarm -workhorse=hgwdev -smallClusterHub=memk \ -swap -chainMinScore=5000 \ -chainLinearGap=loose > swap.log 2>&1 & # Elapsed time: 88m21s cat fb.xenTro3.chainAnoCar2Link.txt # 92033566 bases of 1358334882 (6.775%) in intersection cd /hive/data/genomes/xenTro3/bed ln -s blastz.anoCar2.swap lastz.anoCar2 ########################################################################## # Swap lastz hg19 Human (DONE - 2011-09-21 - Hiram) # original alignment cd /hive/data/genomes/hg19/bed/lastzXenTro3.2011-09-20 cat fb.hg19.chainXenTro3Link.txt # 87928753 bases of 2897316137 (3.035%) in intersection # running the swap - DONE - 2011-09-21 mkdir /hive/data/genomes/xenTro3/bed/blastz.hg19.swap cd /hive/data/genomes/xenTro3/bed/blastz.hg19.swap time nice -n +19 doBlastzChainNet.pl -verbose=2 \ /hive/data/genomes/hg19/bed/lastzXenTro3.2011-09-20/DEF \ -chainMinScore=5000 -chainLinearGap=loose \ -workhorse=hgwdev -smallClusterHub=memk -bigClusterHub=swarm \ -swap > swap.log 2>&1 & # real 37m7.001s cat fb.xenTro3.chainHg19Link.txt # 90929066 bases of 1358334882 (6.694%) in intersection ############################################################################## # Swap lastz danRer7 Zebrafish (DONE - 2011-09-21 - Hiram) # original alignment cd /hive/data/genomes/danRer7/bed/lastzXenTro3.2011-09-20 cat fb.danRer7.chainXenTro3Link.txt # 91771881 bases of 1409770109 (6.510%) in intersection cd /hive/data/genomes/danRer7/bed ln -s lastzXenTro3.2011-09-20 lastz.xenTro3 # running the swap mkdir /hive/data/genomes/xenTro3/bed/blastz.danRer7.swap cd /hive/data/genomes/xenTro3/bed/blastz.danRer7.swap time nice -n +19 doBlastzChainNet.pl -verbose=2 \ /hive/data/genomes/danRer7/bed/lastzXenTro3.2011-09-20/DEF \ -chainMinScore=5000 -chainLinearGap=loose \ -workhorse=hgwdev -smallClusterHub=memk -bigClusterHub=swarm \ -swap > swap.log 2>&1 & # real 67m46.506s cat fb.xenTro3.chainDanRer7Link.txt # 96539932 bases of 1358334882 (7.107%) in intersection ######################################################################### # Swap lastz mm9 Mouse (DONE - 2011-09-21 - Hiram) # original alignment cd /hive/data/genomes/mm9/bed/lastzXenTro3.2011-09-20 cat fb.mm9.chainXenTro3Link.txt # 81920795 bases of 2620346127 (3.126%) in intersection # running the swap - DONE - 2011-09-21 mkdir /hive/data/genomes/xenTro3/bed/blastz.mm9.swap cd /hive/data/genomes/xenTro3/bed/blastz.mm9.swap time nice -n +19 doBlastzChainNet.pl -verbose=2 \ /hive/data/genomes/mm9/bed/lastzXenTro3.2011-09-20/DEF \ -chainMinScore=5000 -chainLinearGap=loose \ -workhorse=hgwdev -smallClusterHub=memk -bigClusterHub=swarm \ -swap > swap.log 2>&1 & # real 57m36.529s cat fb.xenTro3.chainMm9Link.txt # 89770014 bases of 1358334882 (6.609%) in intersection ############################################################################## # Swap lastz melGal1 Turkey (DONE - 2011-09-23 - Hiram) # repeated a second time with better lastz parameters, see below # original alignment cd /hive/data/genomes/melGal1/bed/lastzXenTro3.2011-09-22 cat fb.melGal1.chainXenTro3Link.txt # 35993896 bases of 935922386 (3.846%) in intersection # running the swap mkdir /hive/data/genomes/xenTro3/bed/blastz.melGal1.swap cd /hive/data/genomes/xenTro3/bed/blastz.melGal1.swap time nice -n +19 doBlastzChainNet.pl -verbose=2 \ /hive/data/genomes/melGal1/bed/lastzXenTro3.2011-09-22/DEF \ -swap \ -workhorse=hgwdev -smallClusterHub=memk -bigClusterHub=swarm \ -chainMinScore=5000 -chainLinearGap=loose > swap.log 2>&1 # real 3m29.407s cat fb.xenTro3.chainMelGal1Link.txt # 42384698 bases of 1358334882 (3.120%) in intersection cd /hive/data/genomes/xenTro3/bed ln -s blastz.melGal1.swap lastz.melGal1 ####################################################################### # Swap lastz monDom5 Opossum (DONE - 2011-09-23 - Hiram) # the original alignment cd /hive/data/genomes/monDom5/bed/lastzXenTro3.2011-09-22 cat fb.monDom5.chainXenTro3Link.txt # 75079446 bases of 3501660299 (2.144%) in intersection # and the swap mkdir /hive/data/genomes/xenTro3/bed/blastz.monDom5.swap cd /hive/data/genomes/xenTro3/bed/blastz.monDom5.swap time nice -n +19 doBlastzChainNet.pl -verbose=2 \ /hive/data/genomes/monDom5/bed/lastzXenTro3.2011-09-22/DEF \ -swap -syntenicNet \ -workhorse=hgwdev -smallClusterHub=memk -bigClusterHub=swarm \ -chainMinScore=5000 -chainLinearGap=loose > swap.log 2>&1 & # real 52m10.229s cat fb.xenTro3.chainMonDom5Link.txt # 74593612 bases of 1358334882 (5.492%) in intersection cd /hive/data/genomes/xenTro3/bed ln -s blastz.monDom5.swap lastz.monDom5 ######################################################################### # Swap lastz rn4 Rat (DONE - 2011-09-23 - Hiram) # original alignment cd /hive/data/genomes/rn4/bed/lastzXenTro3.2011-09-22 cat fb.rn4.chainXenTro3Link.txt # 66202042 bases of 2571531505 (2.574%) in intersection # and the swap mkdir /hive/data/genomes/xenTro3/bed/blastz.rn4.swap cd /hive/data/genomes/xenTro3/bed/blastz.rn4.swap time nice -n +19 doBlastzChainNet.pl -verbose=2 \ /hive/data/genomes/rn4/bed/lastzXenTro3.2011-09-22/DEF \ -swap -syntenicNet \ -workhorse=hgwdev -smallClusterHub=memk -bigClusterHub=swarm \ -chainMinScore=5000 -chainLinearGap=loose > swap.log 2>&1 & # Elapsed time: 20m33s cat fb.xenTro3.chainRn4Link.txt # 74593612 bases of 1358334882 (5.492%) in intersection cd /hive/data/genomes/xenTro3/bed ln -s blastz.rn4.swap lastz.rn4 ############################################################################# # lastz Turkey melGal1 again (DONE - 2011-09-28 - Hiram) # original alignment cat fb.melGal1.chainXenTro3Link.txt # 49081729 bases of 935922386 (5.244%) in intersection # running the swap mkdir /hive/data/genomes/xenTro3/bed/blastz.melGal1.swap cd /hive/data/genomes/xenTro3/bed/blastz.melGal1.swap time nice -n +19 doBlastzChainNet.pl -verbose=2 \ /hive/data/genomes/melGal1/bed/lastzXenTro3.2011-09-27/DEF \ -swap \ -workhorse=hgwdev -smallClusterHub=memk -bigClusterHub=swarm \ -chainMinScore=5000 -chainLinearGap=loose > swap.log 2>&1 # real 4m9.023s cat fb.xenTro3.chainMelGal1Link.txt # 57087653 bases of 1358334882 (4.203%) in intersection cd /hive/data/genomes/xenTro3/bed ln -s blastz.melGal1.swap lastz.melGal1 ####################################################################### ## 9-Way Multiz (DONE - 2011-09-28 - Hiram) ssh hgwdev mkdir /hive/data/genomes/xenTro3/bed/multiz9way cd /hive/data/genomes/xenTro3/bed/multiz9way # working on: xenTro3 anoCar2 danRer7 galGal3 hg19 melGal1 mm9 monDom5 rn4 kent/src/hg/utils/phyloTrees/50way.nh xenTro3,xenTro2,galGal3,mm9,hg19,gasAcu1,ornAna1 # All distances remain as specified in the 46way.nh /cluster/bin/phast/tree_doctor --prune-all-but \ anoCar1,xenTro2,danRer6,galGal3,melGal1,mm9,hg19,monDom5,rn4 \ $HOME/kent/src/hg/utils/phyloTrees/50way.nh \ | sed -e "s/anoCar1/anoCar2/; s/xenTro2/xenTro3/; s/danRer6/danRer7/" \ > 9way.nh # what that looks like: cat 9way.nh # (((((hg19:0.144018,(mm9:0.084509,rn4:0.091589):0.271974):0.278985, # monDom5:0.340786):0.181168,((melGal1:0.100000,galGal3:0.165536):0.299223, # anoCar2:0.489241):0.105143):0.172371,xenTro3:0.855573):0.311354, # danRer7:0.878701); # rearrange to get xenTro3 on top: cat << '_EOF_' > xenTro3.9way.nh ((xenTro3:0.855573,(((hg19:0.144018,(mm9:0.084509, rn4:0.091589):0.271974):0.278985,monDom5:0.340786):0.181168, ((melGal1:0.100000,galGal3:0.165536):0.299223, anoCar2:0.489241):0.105143):0.172371):0.311354,danRer7:0.878701); '_EOF_' # << happy emacs # convert to species names /cluster/bin/phast/tree_doctor --rename \ "xenTro3->X_tropicalis; hg19->Human; mm9->Mouse; rn4->Rat; monDom5->Opossum; melGal1->Turkey; galGal3->Chicken; anoCar2->Lizard; danRer7->Zebrafish" \ xenTro3.9way.nh > xenTro3.commonNames.9way.nh # ((X_tropicalis:0.855573,(((Human:0.144018, # (Mouse:0.084509,Rat:0.091589):0.271974):0.278985, # Opossum:0.340786):0.181168,((Turkey:0.100000, # Chicken:0.165536):0.299223,Lizard:0.489241):0.105143):0.172371):0.311354, # Zebrafish:0.878701); # Use this specification in the phyloGif tool: # http://genome.ucsc.edu/cgi-bin/phyloGif # to obtain a png image for src/hg/htdocs/images/phylo/xenTro3_9way.png /cluster/bin/phast/all_dists xenTro3.9way.nh > 9way.distances.txt # Use this output to create the table below grep -i xenTro3 9way.distances.txt | sort -k3,3n # xenTro3 melGal1 1.532310 # xenTro3 monDom5 1.549898 # xenTro3 galGal3 1.597846 # xenTro3 anoCar2 1.622328 # xenTro3 hg19 1.632115 # xenTro3 mm9 1.844580 # xenTro3 rn4 1.851660 # xenTro3 danRer7 2.045628 # If you can fill in all the numbers in this table, you are ready for # the multiple alignment procedure # featureBits chainLink measures # chainAnoCar2Link chain linearGap # distance on xenTro3 on other minScore # 1 1.532 - turkey melGal1 (% 4.20) (% 5.24) 5000 loose # 2 1.550 - opossum monDom5 (% 5.49) (% 2.14) 5000 loose # 3 1.598 - chicken galGal3 (% 4.73) (% 5.27) 5000 loose # 4 1.622 - lizard anoCar2 (% 6.78) (% 5.11) 5000 loose # 5 1.632 - human hg19 (% 6.69) (% 3.04) 5000 loose # 6 1.845 - mouse mm9 (% 6.61) (% 3.13) 5000 loose # 7 1.852 - rat rn4 (% 5.49) (% 2.57) 5000 loose # 8 2.046 - zebrafish danRer7 (% 7.11) (% 6.51) 5000 loose # None of this concern for distances matters in building the first step, the # maf files. # create species list and stripped down tree for autoMZ sed 's/[a-z][a-z]*_//g; s/:[0-9\.][0-9\.]*//g; s/;//; /^ *$/d' \ xenTro3.9way.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 # split the maf files into a set of hashed named files # this hash named split keeps the same chr/contig names in the same # named hash file. mkdir /hive/data/genomes/xenTro3/bed/multiz9way/mafSplit cd /hive/data/genomes/xenTro3/bed/multiz9way/mafSplit for D in `sed -e "s/xenTro3 //" ../species.list` do echo "${D}" mkdir $D cd $D echo "mafSplit -byTarget -useHashedName=10 /dev/null . ../../../lastz.${D}/axtChain/xenTro3.${D}.maf.gz" mafSplit -byTarget -useHashedName=8 /dev/null . \ ../../../lastz.${D}/mafNet/xenTro3.${D}.net.maf.gz cd .. done # construct a list of all possible maf file names. # they do not all exist in each of the species directories find . -type f | wc -l # 2048 find . -type f | grep ".maf$" | xargs -L 1 basename | sort -u > maf.list wc -l maf.list # 256 maf.list mkdir /hive/data/genomes/xenTro3/bed/multiz9way/splitRun cd /hive/data/genomes/xenTro3/bed/multiz9way/splitRun mkdir maf run cd run mkdir penn cp -p /cluster/bin/penn/multiz.2009-01-21/multiz penn cp -p /cluster/bin/penn/multiz.2009-01-21/maf_project penn cp -p /cluster/bin/penn/multiz.2009-01-21/autoMZ penn # set the db and pairs directories here cat > autoMultiz.csh << '_EOF_' #!/bin/csh -ef set db = xenTro3 set c = $1 set result = $2 set run = `/bin/pwd` set tmp = /scratch/tmp/$db/multiz.$c set pairs = /hive/data/genomes/xenTro3/bed/multiz9way/mafSplit /bin/rm -fr $tmp /bin/mkdir -p $tmp /bin/cp -p ../../tree.nh ../../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 = ($run/penn $path); rehash $run/penn/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/xenTro3/bed/multiz9way/splitRun/maf/$(root1).maf} #ENDLOOP '_EOF_' # << happy emacs ln -s ../../mafSplit/maf.list maf.list ssh swarm cd /hive/data/genomes/xenTro3/bed/multiz9way/splitRun/run # the tac reverses the list to get the small jobs first gensub2 maf.list single template stdout | tac > jobList para -ram=8g create jobsList # Completed: 256 of 256 jobs # CPU time in finished jobs: 14427s 240.45m 4.01h 0.17d 0.000 y # IO & Wait Time: 1927s 32.12m 0.54h 0.02d 0.000 y # Average job time: 64s 1.06m 0.02h 0.00d # Longest finished job: 280s 4.67m 0.08h 0.00d # Submission to last job: 1536s 25.60m 0.43h 0.02d # assemble into a single maf file cd /hive/data/genomes/xenTro3/bed/multiz9way head -1 splitRun/maf/001.maf > multiz9way.maf for F in splitRun/maf/*.maf do egrep -v "^#" ${F} done >> multiz9way.maf tail -1 splitRun/maf/001.maf >> multiz9way.maf # -rw-rw-r-- 1 1524521286 Sep 28 13:28 multiz9way.maf # Load into database ssh hgwdev cd /hive/data/genomes/xenTro3/bed/multiz9way mkdir /gbdb/xenTro3/multiz9way ln -s `pwd`/multiz9way.maf /gbdb/xenTro3/multiz9way cd /scratch/tmp time nice -n +19 hgLoadMaf xenTro3 multiz9way # Indexing and tabulating /gbdb/xenTro3/multiz9way/multiz9way.maf # Loaded 2591342 mafs in 1 files from /gbdb/xenTro3/multiz9way # real 1m4.383s time nice -n +19 hgLoadMafSummary -verbose=2 -minSize=30000 \ -mergeGap=1500 -maxSize=200000 xenTro3 multiz9waySummary \ /gbdb/xenTro3/multiz9way/multiz9way.maf # Created 516059 summary blocks from 5782695 components and 2591342 mafs # from /gbdb/xenTro3/multiz9way/multiz9way.maf # real 1m2.197s wc -l multiz9way*.tab # 2591342 multiz9way.tab # 516059 multiz9waySummary.tab # 3107401 total rm multiz9way*.tab ####################################################################### # GAP ANNOTATE MULTIZ9WAY MAF AND LOAD TABLES (DONE - 2011-09-28 - Hiram) # 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. mkdir -p /hive/data/genomes/xenTro3/bed/multiz9way/anno/mafSplit cd /hive/data/genomes/xenTro3/bed/multiz9way/anno/mafSplit time mafSplit -byTarget -useFullSequenceName \ /dev/null . ../../multiz9way.maf # real 0m53.756s ls | wc # 9202 9202 119622 cd /hive/data/genomes/xenTro3/bed/multiz9way/anno # a couple of these do not yet have an N.bed file for DB in xenTro3 melGal1 anoCar2 do cd /hive/data/genomes/${DB} twoBitInfo -nBed ${DB}.2bit ${DB}.N.bed done for DB in `cat ../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 ssh swarm cd /hive/data/genomes/xenTro3/bed/multiz9way/anno mkdir result # NEXT TIME: this template should have a check out exists+ statement cat << '_EOF_' > template #LOOP mafAddIRows -nBeds=nBeds $(path1) /hive/data/genomes/xenTro3/xenTro3.2bit result/$(file1) #ENDLOOP '_EOF_' # << happy emacs ls mafSplit/*.maf > maf.list gensub2 maf.list single template jobList # limit jobs to one per node with the ram=8g requirement para -ram=8g create jobList para try ... check ... push ... # Completed: 9202 of 9202 jobs # CPU time in finished jobs: 4291s 71.52m 1.19h 0.05d 0.000 y # IO & Wait Time: 24581s 409.68m 6.83h 0.28d 0.001 y # Average job time: 3s 0.05m 0.00h 0.00d # Longest finished job: 10s 0.17m 0.00h 0.00d # Submission to last job: 146s 2.43m 0.04h 0.00d # verify all result files have some content, look for 0 size files: find . -type f -size 0 # should see none # combine into one file head -q -n 1 result/GL172637.maf > xenTro3.9way.maf for F in result/*.maf do grep -h -v "^#" ${F} done >> xenTro3.9way.maf # these maf files do not have the end marker, this does nothing: # tail -q -n 1 result/GL172637.maf >> xenTro3.9way.maf # How about an official end marker: echo "##eof maf" >> xenTro3.9way.maf # Load into database rm /gbdb/xenTro3/multiz9way/multiz9way.maf # remove old symlink ln -s `pwd`/xenTro3.9way.maf /gbdb/xenTro3/multiz9way/multiz9way.maf cd /scratch/tmp time nice -n +19 hgLoadMaf xenTro3 multiz9way # Loaded 2929343 mafs in 1 files from /gbdb/xenTro3/multiz9way # real 1m21.568s time nice -n +19 hgLoadMafSummary -verbose=2 -minSize=30000 \ -mergeGap=1500 -maxSize=200000 xenTro3 multiz9waySummary \ /gbdb/xenTro3/multiz9way/multiz9way.maf # Created 516059 summary blocks from 5782695 components and 2929343 mafs # from /gbdb/xenTro3/multiz9way/multiz9way.maf # real 1m6.705s wc -l multiz9way*.tab # 2929343 multiz9way.tab # 516059 multiz9waySummary.tab # 3445402 total rm multiz9way*.tab ######################################################################### # MULTIZ9WAY MAF FRAMES (DONE - 2011-05-16 - Hiram) ssh hgwdev mkdir /hive/data/genomes/xenTro3/bed/multiz9way/frames cd /hive/data/genomes/xenTro3/bed/multiz9way/frames mkdir genes # we get more annotations using the ensGene on mm9 and hg19 # instead of knownGene, so, using ensGene for hg19 and mm9 xenTro3 ensGene hg19 knownGene mm9 knownGene rn4 rgdGene2 monDom5 ensGene melGal1 ensGene galGal3 ensGene anoCar2 ensGene danRer7 ensGene #------------------------------------------------------------------------ # get the genes for all genomes # mRNAs with CDS. single select to get cds+psl, then split that up and # create genePred # using ensGene for: for qDB in xenTro3 monDom5 melGal1 galGal3 anoCar2 danRer7 do echo hgsql -N -e \"'select * from 'ensGene\;\" ${qDB} hgsql -N -e "select * from ensGene" ${qDB} | cut -f 2-11 \ | genePredSingleCover stdin stdout | gzip -2c > genes/$qDB.gp.gz done # if you were using knownGene for mm9 hg19 # genePreds; (must keep only the first 10 columns for knownGene) for qDB in mm9 hg19 do echo hgsql -N -e \"'select * from 'knownGene\;\" ${qDB} hgsql -N -e "select * from knownGene" ${qDB} | cut -f 1-10 \ | genePredSingleCover stdin stdout | gzip -2c > genes/$qDB.gp.gz done hgsql -N -e "select * from rgdGene2" rn4 | cut -f 1-10 \ | genePredSingleCover stdin stdout | gzip -2c > genes/rn4.gp.gz # 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/anoCar2.gp.gz: 17757 # genes/danRer7.gp.gz: 25874 # genes/galGal3.gp.gz: 16491 # genes/hg19.gp.gz: 20709 # genes/melGal1.gp.gz: 14040 # genes/mm9.gp.gz: 20905 # genes/monDom5.gp.gz: 19188 # genes/rn4.gp.gz: 18093 # genes/xenTro3.gp.gz: 18416 ssh hgwdev cd /hive/data/genomes/xenTro3/bed/multiz9way/frames time (cat ../anno/xenTro3.9way.maf \ | nice -n +19 genePredToMafFrames xenTro3 stdin stdout \ `sed -e "s#\([a-zA-Z0-9]*\)#\1 genes/\1.gp.gz#g" ../species.list` \ | gzip > multiz9way.mafFrames.gz) # real 2m22.546s # verify there are frames on everything: zcat multiz9way.mafFrames.gz | awk '{print $4}' | sort | uniq -c # 279231 anoCar2 # 372399 danRer7 # 258019 galGal3 # 341339 hg19 # 249007 melGal1 # 322663 mm9 # 257623 monDom5 # 214301 rn4 # 191614 xenTro3 ssh hgwdev cd /hive/data/genomes/xenTro3/bed/multiz9way/frames time hgLoadMafFrames xenTro3 multiz9wayFrames multiz9way.mafFrames.gz # real 0m21.425s ############################################################################ # creating upstream mafs (DONE - 2011-10-04 - Hiram) ssh hgwdev mkdir /hive/data/genomes/xenTro3/goldenPath/multiz9way cd /hive/data/genomes/xenTro3/goldenPath/multiz9way # run this bash script cat << '_EOF_' > mkUpstream.sh DB=xenTro3 GENE=ensGene NWAY=multiz9way export DB GENE for S in 1000 2000 5000 do echo "making upstream${S}.${GENE}.maf" echo "featureBits ${DB} ${GENE}:upstream:${S} -fa=/dev/null -bed=stdout" featureBits ${DB} ${GENE}:upstreamAll:${S} -fa=/dev/null -bed=stdout \ | perl -wpe 's/_up[^\t]+/\t0/' | sort -k1,1 -k2,2n \ | mafFrags ${DB} ${NWAY} stdin stdout \ -orgs=/hive/data/genomes/${DB}/bed/${NWAY}/species.list \ | gzip -c > upstream${S}.${GENE}.maf.gz echo "done upstream${S}.${GENE}.maf.gz" done '_EOF_' # << happy emacs chmod +x mkUpstream.sh time ./mkUpstream.sh # real 17m4.037s # FIXUP/VERIFY the genbank.conf file to indicate this gene choice for the # upstream automatic generation ######################################################################### # Phylogenetic tree from 9-way (DONE - 2011-09-28 - Hiram) # this is a chancy operation with this business. These alignments are # too distant and not much of the gene sequence is matching. mkdir /hive/data/genomes/xenTro3/bed/multiz9way/4d cd /hive/data/genomes/xenTro3/bed/multiz9way/4d # the split annotated maf's are in: ../anno/result/*.maf cd /hive/data/genomes/xenTro3/bed/multiz9way/4d # using ensGene for xenTro3, only transcribed genes and nothing # from the randoms and other misc. hgsql xenTro3 -Ne \ "select * from ensGene WHERE cdsEnd > cdsStart;" | cut -f 2-20 \ > ensGene.gp # 22705 ensGene.gp genePredSingleCover ensGene.gp stdout | sort > ensGeneNR.gp wc -l ensGeneNR.gp # 18416 ensGeneNR.gp ssh encodek mkdir /hive/data/genomes/xenTro3/bed/multiz9way/4d/run cd /hive/data/genomes/xenTro3/bed/multiz9way/4d/run mkdir ../mfa # newer versions of msa_view have a slightly different operation # the sed of the gp file inserts the reference species in the chr name cat << '_EOF_' > 4d.csh #!/bin/csh -fe set PHASTBIN = /cluster/bin/phast.build/cornellCVS/phast.2010-12-30/bin set r = "/hive/data/genomes/xenTro3/bed/multiz9way" set c = $1 set infile = $r/anno/result/$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/ensGeneNR.gp | sed -e "s/\t$c\t/\txenTro3.$c\t/" > $c.gp set NL=`wc -l $c.gp| gawk '{print $1}'` if ("$NL" != "0") then $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 $c.maf $c.ss '_EOF_' # << happy emacs chmod +x 4d.csh ls -1S /hive/data/genomes/xenTro3/bed/multiz9way/anno/result/*.maf \ | sed -e "s#.*multiz9way/anno/result/##" \ > maf.list cat << '_EOF_' > template #LOOP 4d.csh $(root1) $(path1) {check out line+ ../mfa/$(root1).mfa} #ENDLOOP '_EOF_' # << happy emacs # the tac puts the quick jobs at the front gensub2 maf.list single template stdout | tac > jobList para create jobList para try ... check para -maxJob=5 push para time # there are numerous problems in some of the result. Not a whole lot # of sequence is matching in these distant organisms # Completed: 9143 of 9202 jobs # Crashed: 30 jobs # Other count: 29 jobs # CPU time in finished jobs: 711s 11.85m 0.20h 0.01d 0.000 y # IO & Wait Time: 23273s 387.89m 6.46h 0.27d 0.001 y # Average job time: 3s 0.04m 0.00h 0.00d # Longest finished job: 7s 0.12m 0.00h 0.00d # Submission to last job: 2270s 37.83m 0.63h 0.03d # combine mfa files ssh hgwdev cd /hive/data/genomes/xenTro3/bed/multiz9way/4d # remove the broken empty files, size 0 and size 1: find ./mfa -type f -size 0 | xargs rm -f # most interesting, this did not identify files of size 1: # find ./mfa -type f -size 1 ls -og mfa | awk '$3 == 1' > empty.list sed -e "s#^#mfa/##" empty.list | xargs rm -f #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 # >xenTro3 # >hg19 # >mm9 # >rn4 # >monDom5 # >melGal1 # >galGal3 # >anoCar2 # >danRer7 # tree-commas.nh: #((xenTro3,(((hg19,(mm9, rn4)),monDom5), ((melGal1,galGal3), anoCar2))),danRer7) # use phyloFit to create tree model (output is phyloFit.mod) time nice -n +19 \ /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 # real 2m10.469s mv phyloFit.mod all.mod grep TREE all.mod # TREE: ((xenTro3:0.678022,(((hg19:0.160837,(mm9:0.0938355,rn4:0.110316):0.205024):0.213345,monDom5:0.346247):0.158444,((melGal1:0.0616202,galGal3:0.0595626):0.332209,anoCar2:0.500881):0.098368):0.212664):0.492003,danRer7:0.492003); # it actually looks pretty good as a tree estimate ... ######################################################################### # phastCons 9-way (DONE - 2011-09-29 - Hiram) # split 9way mafs into 10M chunks and generate sufficient statistics # files for # phastCons ssh encodek mkdir -p /hive/data/genomes/xenTro3/bed/multiz9way/cons/SS cd /hive/data/genomes/xenTro3/bed/multiz9way/cons/SS mkdir result done cat << '_EOF_' > mkSS.csh #!/bin/csh -ef set c = $1 set MAF = /hive/data/genomes/xenTro3/bed/multiz9way/anno/result/$c.maf set WINDOWS = /hive/data/genomes/xenTro3/bed/multiz9way/cons/SS/result/$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 mkSS.csh cat << '_EOF_' > template #LOOP mkSS.csh $(root1) {check out line+ done/$(root1)} #ENDLOOP '_EOF_' # << happy emacs # do the easy ones first to see some immediate results ls -1S -r ../../anno/result | sed -e "s/.maf//" > maf.list gensub2 maf.list single template jobList para create jobList para try ... check ... etc # Completed: 9202 of 9202 jobs # CPU time in finished jobs: 816s 13.60m 0.23h 0.01d 0.000 y # IO & Wait Time: 24389s 406.49m 6.77h 0.28d 0.001 y # Average job time: 3s 0.05m 0.00h 0.00d # Longest finished job: 10s 0.17m 0.00h 0.00d # Submission to last job: 1276s 21.27m 0.35h 0.01d # not all of them produce results # there isn't much aligned in these MAF files find ./result -type f | wc # 5690 5690 214806 # 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 encodek mkdir -p /hive/data/genomes/xenTro3/bed/multiz9way/cons/run.cons cd /hive/data/genomes/xenTro3/bed/multiz9way/cons/run.cons # there are going to be only one phastCons run using # this same script. It triggers 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/xenTro3/bed/multiz9way/cons set tmp = $cons/tmp/$f mkdir -p $tmp set ssSrc = $cons set useGrp = "$grp.mod" ln -s $ssSrc/SS/result/$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 +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 find ../SS -type f | grep ".ss$" | sed -e 's/.ss$//' > ss.list wc -l ss.list # 3291 ss.list # Create parasol batch and run it # run for all species cd /hive/data/genomes/xenTro3/bed/multiz9way/cons mkdir -p all cd all # Using the .mod tree cp -p ../../4d/all.mod ./all.mod gensub2 ../run.cons/ss.list single ../run.cons/template jobList para create jobList para try ... check ... para push # Completed: 5690 of 5690 jobs # CPU time in finished jobs: 2411s 40.18m 0.67h 0.03d 0.000 y # IO & Wait Time: 37365s 622.76m 10.38h 0.43d 0.001 y # Average job time: 7s 0.12m 0.00h 0.00d # Longest finished job: 20s 0.33m 0.01h 0.00d # Submission to last job: 2030s 33.83m 0.56h 0.02d cd /hive/data/genomes/xenTro3/bed/multiz9way/cons/all find ./bed -type f | grep ".bed$" | xargs cat | 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 xenTro3 phastConsElements9way mostConserved.bed # Loaded 827762 elements of size 5 # real 0m4.880s # on human we often try for 5% overall cov, and 70% CDS cov # most bets are off here for that goal, these alignments are too few # and too far between featureBits xenTro3 -enrichment ensGene:cds phastConsElements9way # --rho 0.3 --expected-length 45 --target-coverage 0.3 # ensGene:cds 2.244%, phastConsElements9way 10.055%, both 1.669%, cover 74.39%, enrich 7.40x featureBits xenTro3 -enrichment mgcGenes:cds phastConsElements9way # mgcGenes:cds 0.741%, phastConsElements9way 10.055%, both 0.619%, cover 83.53%, enrich 8.31x # Create merged posterier probability file and wiggle track data files cd /hive/data/genomes/xenTro3/bed/multiz9way/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/phastCons9way.wigFix.gz # check integrity of data with wigToBigWig time (zcat downloads/phastCons9way.wigFix.gz \ | wigToBigWig -verbose=2 stdin /hive/data/genomes/xenTro3/chrom.sizes \ phastCons9way.bw) > bigWig.log 2>&1 & tail bigWig.log # pid=32277: VmPeak: 2219472 kB # real 2m45.147s bigWigInfo phastCons9way.bw # version: 4 # isCompressed: yes # isSwapped: 0 # primaryDataSize: 489,270,288 # primaryIndexSize: 14,137,652 # zoomLevels: 10 # chromCount: 5690 # basesCovered: 202,032,632 # mean: 0.630799 # min: 0.000000 # max: 1.000000 # std: 0.378048 # encode those files into wiggle data time (zcat downloads/phastCons9way.wigFix.gz \ | wigEncode stdin phastCons9way.wig phastCons9way.wib) # Converted stdin, upper limit 1.00, lower limit 0.00 # real 1m2.832s du -hsc *.wi? # 193M phastCons9way.wib # 43M phastCons9way.wig # 236M total # Load gbdb and database with wiggle. ln -s `pwd`/phastCons9way.wib /gbdb/xenTro3/multiz9way/phastCons9way.wib time nice -n +19 hgLoadWiggle -pathPrefix=/gbdb/xenTro3/multiz9way \ xenTro3 phastCons9way phastCons9way.wig # real 0m4.207s # use to set trackDb.ra entries for wiggle min and max # and verify table is loaded correctly wigTableStats.sh xenTro3 phastCons9way # db.table min max mean count sumData stdDev viewLimits #xenTro3.phastCons9way 0 1 0.630799 202032632 1.27442e+08 0.378048 viewLimits=0:1 # Create histogram to get an overview of all the data time nice -n +19 hgWiggle -doHistogram -db=xenTro3 \ -hBinSize=0.001 -hBinCount=1000 -hMinVal=0.0 -verbose=2 \ phastCons9way > histogram.data 2>&1 # real 0m42.544s # create plot of histogram: cat << '_EOF_' | gnuplot > histo.png set terminal png small x000000 xffffff xc000ff x66ff66 xffff00 x00ffff set size 1.4, 0.8 set key left box set grid noxtics set grid ytics set title " X. tropicalis xenTro3 Histogram phastCons9way track" set xlabel " phastCons9way 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 & ############################################################################# # phyloP for 7-way (DONE - 2011-09-29 - Hiram) # run phyloP with score=LRT ssh encodek mkdir /cluster/data/xenTro3/bed/multiz9way/consPhyloP cd /cluster/data/xenTro3/bed/multiz9way/consPhyloP mkdir run.phyloP cd 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.495 /cluster/bin/phast.build/cornellCVS/phast.2010-12-30/bin/modFreqs \ ../../cons/all/all.mod 0.495 > all.mod # verify, the BACKGROUND should now be paired up: grep BACK all.mod # BACKGROUND: 0.252500 0.247500 0.247500 0.252500 # following the pattern from panTro3 with grp: "all" (no other grp) cat << '_EOF_' > doPhyloP.csh #!/bin/csh -fe set PHASTBIN = /cluster/bin/phast.build/cornellCVS/phast.2010-12-30/bin set f = $1 set file1 = $f:t set out = $2 set cName = $f:t:r set n = $f:r:e set grp = $cwd:t set cons = /hive/data/genomes/xenTro3/bed/multiz9way/consPhyloP set tmp = $cons/tmp/$grp/$f rm -fr $tmp mkdir -p $tmp set ssSrc = "/hive/data/genomes/xenTro3/bed/multiz9way/cons/SS/result/$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 > $file1.wigFix popd > /dev/null mkdir -p $out:h sleep 4 mv $tmp/$file1.wigFix $out rm -fr $tmp rmdir --ignore-fail-on-non-empty $cons/tmp/$grp/$cName rmdir --ignore-fail-on-non-empty $cons/tmp/$grp rmdir --ignore-fail-on-non-empty $cons/tmp '_EOF_' # << happy emacs # Create list of chunks find ../../cons/SS/result -type f | grep ".ss$" \ | sed -e "s/.ss$//; s#^../../cons/SS/result/##" > ss.list # make sure the list looks good wc -l ss.list # 5690 ss.list # Create template file # file1 == $chr/$chunk/file name without .ss suffix cat << '_EOF_' > template #LOOP ../run.phyloP/doPhyloP.csh $(path1) {check out line+ wigFix/$(dir1)/$(file1).wigFix} #ENDLOOP '_EOF_' # << happy emacs ###################### Running all species ####################### # setup run for all species mkdir /hive/data/genomes/xenTro3/bed/multiz9way/consPhyloP/all cd /hive/data/genomes/xenTro3/bed/multiz9way/consPhyloP/all rm -fr wigFix mkdir wigFix gensub2 ../run.phyloP/ss.list single ../run.phyloP/template jobList # the -ram=8g will allow only one job per node to slow this down since # it would run too fast otherwise. Either run on one of the small # klusters or use the -ram=8g on the para create para -ram=8g create jobList para try ... check ... etc ... para -maxJob=64 push para time > run.time # Completed: 5690 of 5690 jobs # CPU time in finished jobs: 9670s 161.17m 2.69h 0.11d 0.000 y # IO & Wait Time: 38146s 635.76m 10.60h 0.44d 0.001 y # Average job time: 8s 0.14m 0.00h 0.00d # Longest finished job: 63s 1.05m 0.02h 0.00d # Submission to last job: 2448s 40.80m 0.68h 0.03d # make downloads mkdir downloads time (find ./wigFix -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/phyloP9way.wigFix.gz) & # real 5m28.973s # check integrity of data with wigToBigWig time (zcat downloads/phyloP9way.wigFix.gz \ | wigToBigWig -verbose=2 stdin /hive/data/genomes/xenTro3/chrom.sizes \ phyloP9way.bw) > bigWig.log 2>&1 & egrep "real|VmPeak" bigWig.log # pid=17458: VmPeak: 2219472 kB # real 2m53.180s bigWigInfo phyloP9way.bw # version: 4 # isCompressed: yes # isSwapped: 0 # primaryDataSize: 366,163,491 # primaryIndexSize: 14,137,652 # zoomLevels: 10 # chromCount: 5690 # basesCovered: 202,032,632 # mean: 0.447745 # min: -3.316000 # max: 2.298000 # std: 0.772169 # encode those files into wiggle data time (zcat downloads/phyloP9way.wigFix.gz \ | wigEncode stdin phyloP9way.wig phyloP9way.wib) & # Converted stdin, upper limit 2.30, lower limit -3.32 # real 1m3.694s du -hsc *.wi? # 193M phyloP9way.wib # 43M phyloP9way.wig # 236M total # Load gbdb and database with wiggle. ln -s `pwd`/phyloP9way.wib /gbdb/xenTro3/multiz9way/phyloP9way.wib nice hgLoadWiggle -pathPrefix=/gbdb/xenTro3/multiz9way xenTro3 \ phyloP9way phyloP9way.wig # use to set trackDb.ra entries for wiggle min and max # and verify table is loaded correctly wigTableStats.sh xenTro3 phyloP9way # db.table min max mean count sumData stdDev viewLimits # xenTro3.phyloP9way -3.316 2.298 0.447745 202032632 9.04591e+07 0.772169 # viewLimits=-3.316:2.298 # Create histogram to get an overview of all the data time nice -n +19 hgWiggle -doHistogram -db=xenTro3 \ -hBinSize=0.001 -hBinCount=1000 -hMinVal=0.0 -verbose=2 \ phyloP9way > histogram.data 2>&1 # real 0m19.396s # find out the range for the 2:5 graph grep -v chrom histogram.data | grep "^[0-9]" | ave -col=5 stdin # Q1 0.000240 # median 0.000383 # Q3 0.000726 # average 0.001000 # min 0.000107 # max 0.041706 # count 1000 # total 1.000011 # standard deviation 0.002765 # create plot of histogram: cat << '_EOF_' | gnuplot > histo.png set terminal png small x000000 xffffff xc000ff x66ff66 xffff00 x00ffff set size 1.4, 0.8 set key left box set grid noxtics set grid ytics set title " X. tropicalis xenTro3 Histogram phyloP9way track" set xlabel " phyloP9way score" set ylabel " Relative Frequency" set y2label " Cumulative Relative Frequency (CRF)" set y2range [0:1] set y2tics set yrange [0:0.007] 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 data for 7-way (DONE - 2011-09-29 - Hiram) mkdir /usr/local/apache/htdocs-hgdownload/goldenPath/xenTro3/phastCons9way cd /usr/local/apache/htdocs-hgdownload/goldenPath/xenTro3/phastCons9way ln -s \ /hive/data/genomes/xenTro3/bed/multiz9way/cons/all/downloads/phastCons9way.wigFix.gz \ ./phastCons9way.wigFix.gz ln -s /hive/data/genomes/xenTro3/bed/multiz9way/cons/all/phastCons9way.bw \ ./phastCons9way.bw ln -s /hive/data/genomes/xenTro3/bed/multiz9way/cons/all/all.mod \ ./vertebrate.mod # use a README from a recent multiz like this, this one is from anoCar2 cd /hive/data/genomes/xenTro3/bed/multiz9way/cons cp -p /hive/data/genomes/anoCar2/bed/multiz7way/cons/README.txt . # need to edit this to finish it off ln -s /hive/data/genomes/xenTro3/bed/multiz9way/cons/README.txt . md5sum * > /hive/data/genomes/xenTro3/bed/multiz9way/cons/md5sum.txt ln -s /hive/data/genomes/xenTro3/bed/multiz9way/cons/md5sum.txt . mkdir /usr/local/apache/htdocs-hgdownload/goldenPath/xenTro3/phyloP9way cd /usr/local/apache/htdocs-hgdownload/goldenPath/xenTro3/phyloP9way ln -s \ /hive/data/genomes/xenTro3/bed/multiz9way/consPhyloP/all/downloads/phyloP9way.wigFix.gz \ ./phyloP9way.wigFix.gz ln -s \ /hive/data/genomes/xenTro3/bed/multiz9way/consPhyloP/all/phyloP9way.bw \ ./phyloP9way.bw ln -s \ /hive/data/genomes/xenTro3/bed/multiz9way/consPhyloP/run.phyloP/all.mod \ ./vertebrate.mod # use a README from a recent multiz like this, this one is from anoCar2 /hive/data/genomes/xenTro3/bed/multiz9way/consPhyloP cp /hive/data/genomes/anoCar2/bed/multiz7way/consPhyloP/README.txt . # need to edit this to finish it off ln -s /hive/data/genomes/xenTro3/bed/multiz9way/consPhyloP/README.txt . md5sum * > /hive/data/genomes/xenTro3/bed/multiz9way/consPhyloP/md5sum.txt ln -s /hive/data/genomes/xenTro3/bed/multiz9way/consPhyloP/md5sum.txt . mkdir /hive/data/genomes/xenTro3/bed/multiz9way/downloads cd /hive/data/genomes/xenTro3/bed/multiz9way/downloads # this was already done elsewhere: grep TREE ../cons/all/all.mod | sed -e "s/TREE: //" > 9way.nh /cluster/bin/phast/tree_doctor --rename \ "xenTro3->X_tropicalis; hg19->Human; mm9->Mouse; rn4->Rat; monDom5->Opossum; melGal1->Turkey; galGal3->Chicken; anoCar2->Lizard; danRer7->Zebrafish" \ 9way.nh > xenTro3.commonNames.9way.nh time cp -p ../anno/xenTro3.9way.maf ./multiz9way.maf # real 1m13.654s time gzip multiz9way.maf # real 7m12.108s mkdir /usr/local/apache/htdocs-hgdownload/goldenPath/xenTro3/multiz9way cd /usr/local/apache/htdocs-hgdownload/goldenPath/xenTro3/multiz9way ln -s /hive/data/genomes/xenTro3/bed/multiz9way/downloads/9way.nh . ln -s \ /hive/data/genomes/xenTro3/bed/multiz9way/downloads/xenTro3.commonNames.9way.nh . ln -s \ /hive/data/genomes/xenTro3/bed/multiz9way/downloads/multiz9way.maf.gz . # use a README from a recent multiz like this, this one is from danRer7 cd /hive/data/genomes/xenTro3/bed/multiz9way/downloads cp -p /hive/data/genomes/anoCar2/bed/multiz7way/downloads/README.txt . # need to edit this to finish it off cd /usr/local/apache/htdocs-hgdownload/goldenPath/xenTro3/multiz9way ln -s /hive/data/genomes/xenTro3/bed/multiz9way/downloads/README.txt . md5sum *.nh *.gz *.txt \ > /hive/data/genomes/xenTro3/bed/multiz9way/downloads/md5sum.txt ln -s /hive/data/genomes/xenTro3/bed/multiz9way/downloads/md5sum.txt . ############################################################################# # hgPal downloads (DONE braney 2011-09-30) # FASTA from 9way for ensGene ssh hgwdev screen bash rm -rf /cluster/data/xenTro3/bed/multiz9way/pal mkdir /cluster/data/xenTro3/bed/multiz9way/pal cd /cluster/data/xenTro3/bed/multiz9way/pal for i in `cat ../species.list`; do echo $i; done > order.lst mz=multiz9way gp=ensGene db=xenTro3 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 29m58.498s # user 8m8.312s # sys 6m7.591s mz=multiz9way gp=ensGene db=xenTro3 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=multiz9way gp=ensGene db=xenTro3 pd=/usr/local/apache/htdocs-hgdownload/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 rat rn5 (DONE - 2012-11-14 - Hiram) # the original alignment: cd /hive/data/genomes/rn5/bed/lastzXenTro3.2012-11-13 cat fb.rn5.chainXenTro3Link.txt # 82011038 bases of 2572853723 (3.188%) in intersection # and this swap: mkdir /hive/data/genomes/xenTro3/bed/blastz.rn5.swap cd /hive/data/genomes/xenTro3/bed/blastz.rn5.swap time nice -n +19 doBlastzChainNet.pl -verbose=2 \ /hive/data/genomes/rn5/bed/lastzXenTro3.2012-11-13/DEF \ -workhorse=hgwdev -smallClusterHub=encodek -bigClusterHub=swarm \ -swap -chainMinScore=5000 -chainLinearGap=loose > swap.log 2>&1 & # real 25m43.688s cat fb.xenTro3.chainRn5Link.txt # 87276772 bases of 1358334882 (6.425%) in intersection # set sym link to indicate this is the lastz for this genome: cd /hive/data/genomes/xenTro3/bed ln -s blastz.rn5.swap lastz.rn5 #########################################################################