# for emacs: -*- mode: sh; -*- # $Id: bosTau6.txt,v 1.6 2010/05/06 16:27:44 chinhli Exp $ # Bos taurus (bovine) -- Bos_taurus_UMD_3.1 (2009-11-19) # file template copied from oviAri1.txt # Bos taurus (NCBI Project ID: 10709, Accession: GCA_000003055.3) # by International Cow Genomics Consortium (ISGC) # assembly] sequence: # ftp://ftp.ncbi.nlm.nih.gov/genbank/genomes/Eukaryotes/vertebrates_mammals/Bos_taurus/Bos_taurus_UMD_3.1/ ########################################################################## # Download sequence (DONE - 2011-05-02 Chin) mkdir /hive/data/genomes/bosTau6 cd /hive/data/genomes/bosTau6 mkdir genbank cd genbank wget --timestamping -r --cut-dirs=6 --level=0 -nH -x \ --no-remove-listing -np \ "ftp://ftp.ncbi.nlm.nih.gov/genbank/genomes/Eukaryotes/vertebrates_mammals/Bos_taurus/Bos_taurus_UMD_3.1/*" # FINISHED --2011-05-02 15:20:18-- # Downloaded: 229 files, 1.9G in 11m 53s (2.73 MB/s) # Read ASSEMBLY_INFO mkdir ucscChr # stay at genbank directory # fixup the accession names to become UCSC chrom names export S=Primary_Assembly/assembled_chromosomes cut -f2 ${S}/chr2acc | while read ACC do C=`grep "${ACC}" ${S}/chr2acc | cut -f1` echo "${ACC} -> chr${C}" zcat ${S}/AGP/chr${C}.comp.agp.gz \ | sed -e "s/^${ACC}/chr${C}/" | gzip > ucscChr/chr${C}.agp.gz done export S=Primary_Assembly/assembled_chromosomes cut -f2 ${S}/chr2acc | while read ACC do C=`grep "${ACC}" ${S}/chr2acc | cut -f1` echo "${ACC} -> chr${C}" echo ">chr${C}" > ucscChr/chr${C}.fa zcat ${S}/FASTA/chr${C}.fa.gz | grep -v "^>" >> ucscChr/chr${C}.fa gzip ucscChr/chr${C}.fa & done # Check them with faSize faSize Primary_Assembly/assembled_chromosomes/FASTA/chr*.fa.gz # 2660906405 bases (20740230 N's 2640166175 real 2640166175 upper # 0 lower) in 30 sequences in 30 files faSize ucscChr/chr*.fa.gz # 2660906405 bases (20740230 N's 2640166175 real 2640166175 upper # 0 lower) in 30 sequences in 30 files # For unplaced scalfolds, named them as chrUn_xxxxxxxx # where xxxxxx is the original access id as: chrUn_GL340781.1 # The ".1" at the end need to be filter out since # MySQL does not allow "." as part of the table name and # will casue problems in genbank task step later export S=Primary_Assembly/unplaced_scaffolds zcat ${S}/AGP/unplaced.scaf.agp.gz | grep "^#" > ucscChr/chrUn.agp # append the gap records zcat ${S}/AGP/unplaced.scaf.agp.gz | grep -v "^#" \ | sed -e "s/^/chrUn_/" -e "s/\.1//" >> ucscChr/chrUn.agp gzip ucscChr/chrUn.agp & zcat ${S}/FASTA/unplaced.scaf.fa.gz \ | sed -e "s#^>.*|tpg|#>chrUn_#; s#|.*##" -e "s/\.1//" \ | gzip > ucscChr/chrUn.fa.gz # about 3286 sequences in the unplaced zcat ucscChr/chrUn.fa.gz | grep "^>" | wc # 3286 3286 13144 # Check them with faSize faSize Primary_Assembly/unplaced_scaffolds/FASTA/unplaced.scaf.fa.gz # 9499556 bases (40 N's 9499516 real 9499516 upper 0 # lower) in 3286 sequences in 1 files faSize ucscChr/chrUn.fa.gz # 9499556 bases (40 N's 9499516 real 9499516 upper 0 # lower) in 3286 sequences in 1 files ######################################################################### # Initial database build (DONE - 2011-05-09 - Chin) cd /hive/data/genomes/bosTau6 cat << '_EOF_' > bosTau6.config.ra # Config parameters for makeGenomeDb.pl: db bosTau6 clade mammal genomeCladePriority 31 scientificName Bos taurus commonName Cow assemblyDate Nov. 2009 assemblyLabel UMD_3.1 (NCBI project 10708, accession GCA_000003055.3) assemblyShortLabel Bos_taurus_UMD_3.1 orderKey 236 mitoAcc NC_006853 fastaFiles /hive/data/genomes/bosTau6/genbank/ucscChr/chr*.fa.gz agpFiles /hive/data/genomes/bosTau6/genbank/ucscChr/chr*.agp.gz # qualFiles none dbDbSpeciesDir cow taxId 9913 '_EOF_' # << happy emacs time makeGenomeDb.pl -noGoldGapSplit -workhorse=hgwdev bosTau6.config.ra \ > makeGenomeDb.log 2>&1 & # real 21m31.993s # real real 12m42.419s # add the trackDb entries to the source tree, and the 2bit link: ln -s `pwd`/bosTau6.unmasked.2bit /gbdb/bosTau6/bosTau6.2bit # Per instructions in makeGenomeDb.log: #Search for '***' notes in each file in and make corrections (sometimes the #files used for a previous assembly might make a better template): # description.html /cluster/data/bosTau6/html/{trackDb.ra,gap.html,gold.html} # #Then copy these files to your ~/kent/src/hg/makeDb/trackDb/cow/bosTau6 # - cd ~/kent/src/hg/makeDb/trackDb # - edit makefile to add bosTau6 to DBS. # - git add cow/bosTau6/*.{ra,html} # - git commit -m "Added bosTau6 to DBS." makefile # - git commit -m "Initial descriptions for bosTau6." cow/bosTau6/*.{ra,html} # - git pull; git push # - Run make update DBS=bosTau6 and make alpha when done. # - (optional) Clean up /cluster/data/bosTau6/TemporaryTrackDbCheckout ######################################################################### # RepeatMasker (DONE 2011-05-09 - Chin) mkdir /hive/data/genomes/bosTau6/bed/repeatMasker cd /hive/data/genomes/bosTau6/bed/repeatMasker time nice -n +19 doRepeatMasker.pl -buildDir=`pwd` \ -workhorse=hgwdev -bigClusterHub=swarm -noSplit bosTau6 > do.log 2>&1 & # real 355m43.294s cat faSize.rmsk.txt # 2670422299 bases (20740270 N's 2649682029 real 1370807360 upper # 1278874669 lower) in 3317 sequences in 1 files # %47.89 masked total, %48.27 masked real ######################################################################### # simpleRepeats (DONE - 2011-05-09 - Chin) mkdir /hive/data/genomes/bosTau6/bed/simpleRepeat cd /hive/data/genomes/bosTau6/bed/simpleRepeat time nice -n +19 doSimpleRepeat.pl -buildDir=`pwd` -workhorse=hgwdev \ -dbHost=hgwdev -bigClusterHub=swarm -smallClusterHub=memk \ bosTau6 > do.log 2>&1 & # real 13m47.346s # Batch failed after 4 tries on # ./TrfRun.csh /hive/data/genomes/bosTau6/TrfPart/071/071.lst.bed # Command failed: # ssh -x -o 'StrictHostKeyChecking = no' -o 'BatchMode = yes' memk # nice /hive/data/genomes/bosTau6/bed/simpleRepeat/run.cluster/doTrf.csh # 071 failed due to no input TrfPart/071/071.lst.bed # which is chrM and can be ignored. # to continue, create an empty TrfPart/071/071.lst.bed by: touch /hive/data/genomes/bosTau6/TrfPart/071/071.lst.bed # run a para time > run.time on memk to get that file to exist: # Checking finished jobs # Completed: 71 of 72 jobs # Crashed: 1 jobs # CPU time in finished jobs: 12340s 205.67m 3.43h 0.14d 0.000 y # IO & Wait Time: 236s 3.93m 0.07h 0.00d 0.000 y # Average job time: 177s 2.95m 0.05h 0.00d # Longest finished job: 606s 10.10m 0.17h 0.01d # Submission to last job: 984s 16.40m 0.27h 0.01d # continue time doSimpleRepeat.pl -buildDir=`pwd` -bigClusterHub=swarm \ -dbHost=hgwdev -workhorse=hgwdev -smallClusterHub=memk \ -continue=filter bosTau6 > filter.log 2>&1 & # real 0m20.430s cat fb.simpleRepeat # 32390355 bases of 2649685036 (1.222%) in intersection # add to the repeatMasker cd /hive/data/genomes/bosTau6 twoBitMask bosTau6.rmsk.2bit -add bed/simpleRepeat/trfMask.bed bosTau6.2bit # safe to ignore warnings about >=13 fields twoBitToFa bosTau6.2bit stdout | faSize stdin > bosTau6.2bit.faSize.txt cat bosTau6.2bit.faSize.txt # 2670422299 bases (20740270 N's 2649682029 real 1370313596 upper # 1279368433 lower) in 3317 sequences in 1 files # %47.91 masked total, %48.28 masked real # double check with featureBits featureBits -countGaps bosTau6 gap # 20740270 bases of 2670422299 (0.777%) in intersection rm /gbdb/bosTau6/bosTau6.2bit ln -s `pwd`/bosTau6.2bit /gbdb/bosTau6/bosTau6.2bit ######################################################################### # Marking *all* gaps - they are not all in the AGP file # (DONE - 2011-05-10 - Chin) mkdir /hive/data/genomes/bosTau6/bed/allGaps cd /hive/data/genomes/bosTau6/bed/allGaps time nice -n +19 findMotif -motif=gattaca -verbose=4 \ -strand=+ ../../bosTau6.unmasked.2bit > findMotif.txt 2>&1 # real 1m10.942s grep "^#GAP " findMotif.txt | sed -e "s/^#GAP //" > allGaps.bed featureBits bosTau6 -not gap -bed=notGap.bed # 2649685036 bases of 2649685036 (100.000%) in intersection featureBits bosTau6 allGaps.bed notGap.bed -bed=new.gaps.bed # 3007 bases of 2649685036 (0.000%) in intersection # what is the highest index in the existing gap table: hgsql -N -e "select ix from gap;" bosTau6 | sort -n | tail -1 # 17354 # use tcsh and ctrl-c to create the here doc cat << '_EOF_' > mkGap.pl #!/usr/bin/env perl use strict; use warnings; my $ix=`hgsql -N -e "select ix from gap;" bosTau6 | 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 bosTau6 other.bed # 3007 bases of 2649685036 (0.000%) in intersection hgLoadBed -sqlTable=$HOME/kent/src/hg/lib/gap.sql \ -noLoad bosTau6 otherGap other.bed # Loaded 2011 elements of size 8 # adding this many: wc -l bed.tab # 2011 bed.tab # starting with this many hgsql -e "select count(*) from gap;" bosTau6 # 72454 hgsql bosTau6 -e 'load data local infile "bed.tab" into table gap;' # result count: hgsql -e "select count(*) from gap;" bosTau6 # 74465 # == 72454 + 2011 ######################################################################## # Create kluster run files (DONE - 2011-05-10 - Chin) # numerator is bosTau6 gapless bases "real" as reported by: featureBits -noRandom -noHap bosTau6 gap # 20740230 bases of 2640182513 (0.786%) in intersection # denominator is hg19 gapless bases as reported by: # featureBits -noRandom -noHap hg19 gap # 234344806 bases of 2861349177 (8.190%) in intersection # 1024 is threshold used for human -repMatch: calc \( 2640182513 / 2861349177 \) \* 1024 ( 2640182513 / 2861349177 ) * 1024 = 944.850393 # ==> use -repMatch=900 according to size scaled down from 1024 for human. # and rounded down to nearest 50 cd /hive/data/genomes/bosTau6 blat bosTau6.2bit \ /dev/null /dev/null -tileSize=11 -makeOoc=jkStuff/bosTau6.11.ooc \ -repMatch=900 & # Wrote 33615 overused 11-mers to jkStuff/bosTau6.11.ooc mkdir /hive/data/staging/data/bosTau6 cp -p bosTau6.2bit jkStuff/bosTau6.11.ooc /hive/data/staging/data/bosTau6 cp -p chrom.sizes /hive/data/staging/data/bosTau6 # check non-bridged gaps to see what the typical size is: hgsql -N \ -e 'select * from gap where bridge="no" order by size;' bosTau6 \ | sort -k7,7nr # most non-bridges gaps have size = 100, follow pig's example (5000) # decide on a minimum gap for this break gapToLift -verbose=2 -minGap=100 bosTau6 jkStuff/nonBridged.lft \ -bedFile=jkStuff/nonBridged.bed cp -p jkStuff/nonBridged.lft \ /hive/data/staging/data/bosTau6/bosTau6.nonBridged.lft # ask cluster-admin to copy (evry time if any file chsnged) # /hive/data/staging/data/bosTau6 directory to # /scratch/data/bosTau6 on cluster nodes # before porceed to genbank step ######################################################################## # GENBANK AUTO UPDATE (DONE - 2011-05-10 - Chin) ssh hgwdev cd $HOME/kent/src/hg/makeDb/genbank git pull # edit etc/genbank.conf to add bosTau6 just before susScr2 # bosTau6 (Cow) bosTau6.serverGenome = /hive/data/genomes/bosTau6/bosTau6.2bit bosTau6.clusterGenome = /scratch/data/bosTau6/bosTau6.2bit bosTau6.ooc = /scratch/data/bosTau6/bosTau6.11.ooc bosTau6.lift = no bosTau6.perChromTables = no bosTau6.refseq.mrna.native.pslCDnaFilter = ${ordered.refseq.mrna.native.pslCDnaFilter} bosTau6.refseq.mrna.xeno.pslCDnaFilter = ${ordered.refseq.mrna.xeno.pslCDnaFilter} bosTau6.genbank.mrna.native.pslCDnaFilter = ${ordered.genbank.mrna.native.pslCDnaFilter} bosTau6.genbank.mrna.xeno.pslCDnaFilter = ${ordered.genbank.mrna.xeno.pslCDnaFilter} bosTau6.genbank.est.native.pslCDnaFilter = ${ordered.genbank.est.native.pslCDnaFilter} bosTau6.genbank.est.xeno.pslCDnaFilter = ${ordered.genbank.est.xeno.pslCDnaFilter} bosTau6.downloadDir = bosTau6 bosTau6.refseq.mrna.native.load = yes bosTau6.refseq.mrna.xeno.load = yes bosTau6.refseq.mrna.xeno.loadDesc = yes git add etc/genbank.conf git commit -m "Added bosTau6" etc/genbank.conf git push # update /cluster/data/genbank/: make etc-update # Edit src/lib/gbGenome.c to add new species. Skipped screen # control this business with a screen since it takes a while cd /cluster/data/genbank time nice -n +19 ./bin/gbAlignStep -initial bosTau6 & # logFile: var/build/logs/2011.05.10-15:10:42.bosTau6.initalign.log # real 736m26.569s # To re-do, rm the dir first: # /cluster/data/genbank/data/aligned/genbank.183.0/bosTau6 # load database when finished ssh hgwdev cd /cluster/data/genbank time nice -n +19 ./bin/gbDbLoadStep -drop -initialLoad bosTau6 & # logFile: var/dbload/hgwdev/logs/2011.05.11-08:24:53.dbload.log # real 32m5.964s # enable daily alignment and update of hgwdev cd ~/kent/src/hg/makeDb/genbank gitpull # add bosTau6 to: # etc/align.dbs # etc/hgwdev.dbs git add etc/align.dbs git add etc/hgwdev.dbs git commit -m "Added bosTau6 - Cow" etc/align.dbs etc/hgwdev.dbs git push make etc-update ######################################################################### # reset default position (DONE - 2011-05-11 - Chin) # Reset default position to an area of chr3 with a number of genes hgsql -e \ 'update dbDb set defaultPos="chr3:21,592,157-21,596,916" where name="bosTau6";' \ hgcentraltest ############################################################################ # ctgPos2 track - showing clone sequence locations on chromosomes # (DONE 2011-05-11 - Chin) # NOTE - create bosTau6 entry in all.joiner since this is a new species mkdir /hive/data/genomes/bosTau6/bed/ctgPos2 cd /hive/data/genomes/bosTau6/bed/ctgPos2 cat << '_EOF_' > agpToCtgPos2.pl #!/usr/bin/env perl use warnings; use strict; my $argc = scalar(@ARGV); if ($argc != 1) { printf STDERR "usage: zcat your.files.agp.gz | agpToCtgPos2.pl /dev/stdin > ctgPos2.tab\n"; exit 255; } my $agpFile = shift; open (FH, "<$agpFile") or die "can not read $agpFile"; while (my $line = ) { next if ($line =~ m/^#/); chomp $line; my @a = split('\s+', $line); next if ($a[4] =~ m/^U$/); my $chrSize = $a[2]-$a[1]+1; my $ctgSize = $a[7]-$a[6]+1; die "sizes differ $chrSize != $ctgSize\n$line\n" if ($chrSize != $ctgSize); printf "%s\t%d\t%s\t%d\t%d\t%s\n", $a[5], $chrSize, $a[0], $a[1]-1, $a[2], $a[4]; } close (FH); '_EOF_' # << happy emacs chmod 775 agpToCtgPos2.pl export S=../../genbank/Primary_Assembly/assembled_chromosomes cut -f2 ${S}/chr2acc | while read ACC do C=`grep "${ACC}" ${S}/chr2acc | cut -f1` zcat ${S}/AGP/chr${C}.agp.gz \ | sed -e "s/^${ACC}/chr${C}/" done | ./agpToCtgPos2.pl /dev/stdin > ctgPos2.tab hgLoadSqlTab bosTau6 ctgPos2 $HOME/kent/src/hg/lib/ctgPos2.sql ctgPos2.tab # add the track ctgPos2 to src/hg/makeDb/trackDb/cow/bosTau6/trackDb.ra # at src/makeDb/trackdb do "make update DBS=bosTau6" or/and "make alpha" # based on result of # hgsql -N -e "select type from ctgPos2;" bosTau6 | sort | uniq # W # prepare the src/hg/makeDb/trackDb/cow/bosTau6/ctgPos2.html ############################################################################ # running cpgIsland business (DONE -2011-05-16 - Chin) mkdir /hive/data/genomes/bosTau6/bed/cpgIsland cd /hive/data/genomes/bosTau6/bed/cpgIsland cvs -d /projects/compbio/cvsroot checkout -P hg3rdParty/cpgIslands cd hg3rdParty/cpgIslands # needed to fixup this source, adding include to readseq.c: #include "string.h" # and to cpg_lh.c: #include "unistd.h" #include "stdlib.h" # and fixing a declaration in cpg_lh.c sed -e "s#\(extern char\* malloc\)#// \1#" cpg_lh.c > tmp.c mv tmp.c cpg_lh.c make cd ../../ ln -s hg3rdParty/cpgIslands/cpglh.exe mkdir -p hardMaskedFa cut -f1 ../../chrom.sizes | while read C do echo ${C} twoBitToFa ../../bosTau6.2bit:$C stdout \ | maskOutFa stdin hard hardMaskedFa/${C}.fa done ssh swarm cd /hive/data/genomes/bosTau6/bed/cpgIsland mkdir results cut -f1 ../../chrom.sizes > chr.list cat << '_EOF_' > template #LOOP ./runOne $(root1) {check out exists results/$(root1).cpg} #ENDLOOP '_EOF_' # << happy emacs # the faCount business is to make sure there is enough sequence to # work with in the fasta. cpglh.exe does not like files with too many # N's - it gets stuck cat << '_EOF_' > runOne #!/bin/csh -fe set C = `faCount hardMaskedFa/$1.fa | grep ^chr | awk '{print $2 - $7 }'` if ( $C > 200 ) then ./cpglh.exe hardMaskedFa/$1.fa > /scratch/tmp/$1.$$ mv /scratch/tmp/$1.$$ $2 else touch $2 endif '_EOF_' # << happy emacs chmod 775 runOne gensub2 chr.list single template jobList para create jobList para try para check ... etc para time para problems para status # then, kick it with para push # check it with plb # when all are done, para time shows: # Checking finished jobs # Completed: 3317 of 3317 jobs # CPU time in finished jobs: 197s 3.29m 0.05h 0.00d 0.000 y # IO & Wait Time: 9373s 156.21m 2.60h 0.11d 0.000 y # Average job time: 3s 0.05m 0.00h 0.00d # Longest finished job: 18s 0.30m 0.01h 0.00d # Submission to last job: 677s 11.28m 0.19h 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 ssh hgwdev cd /hive/data/genomes/bosTau6/bed/cpgIsland hgLoadBed bosTau6 cpgIslandExt -tab \ -sqlTable=$HOME/kent/src/hg/lib/cpgIslandExt.sql cpgIsland.bed # Loaded 37638 elements of size 10 # Sorted # Creating table definition for cpgIslandExt # Saving bed.tab # Loading bosTau6 # cleanup rm -fr hardMaskedFa ########################################################################### # ornAna1 Platypus LASTZ/CHAIN/NET (DONE - 2011-05-20 - Chin) screen # use screen to control this job mkdir /cluster/data/bosTau6/bed/blastzOrnAna1.2011-05-20 cd /cluster/data/bosTau6/bed/blastzOrnAna1.2011-05-20 cat << '_EOF_' > DEF # Cow vs. Platypus BLASTZ_Y=3400 BLASTZ_L=6000 BLASTZ_K=2200 BLASTZ_Q=/cluster/data/blastz/HoxD55.q BLASTZ_M=50 # TARGET: Cow bosTau6 SEQ1_DIR=/scratch/data/bosTau6/bosTau6.2bit SEQ1_LEN=/cluster/data/bosTau6/chrom.sizes # Maximum number of scaffolds that can be lumped together SEQ1_LIMIT=300 SEQ1_CHUNK=20000000 SEQ1_LAP=0 # QUERY: Platypus ornAna1 SEQ2_DIR=/scratch/data/ornAna1/ornAna1.2bit SEQ2_LEN=/cluster/data/ornAna1/chrom.sizes SEQ2_CHUNK=20000000 SEQ2_LIMIT=300 SEQ2_LAP=0 BASE=/cluster/data/bosTau6/bed/blastzOrnAna1.2011-05-20 TMPDIR=/scratch/tmp '_EOF_' # << happy emacs time nice -n +19 doBlastzChainNet.pl `pwd`/DEF \ -syntenicNet -noDbNameCheck \ -workhorse=hgwdev -smallClusterHub=memk -bigClusterHub=swarm \ -chainMinScore=5000 -verbose=2 -chainLinearGap=loose > do.log 2>&1 & # real 329m59.325s # filter with doRecipBest.pl time doRecipBest.pl -workhorse=hgwdev -buildDir=`pwd` \ bosTau6 ornAna1 > rbest.log 2>&1 # real 15m15.120s cat fb.bosTau6.chainOrnAna1Link.txt # 195910907 bases of 2649682029 (7.394%) in intersection cd /hive/data/genomes/bosTau6/bed ln -s blastzOrnAna1.2011-05-20 lastz.ornAna1 mkdir /cluster/data/ornAna1/bed/blastz.bosTau6.swap cd /cluster/data/ornAna1/bed/blastz.bosTau6.swap time nice -n +19 doBlastzChainNet.pl \ /cluster/data/bosTau6/bed/blastzOrnAna1.2011-05-20/DEF \ -chainMinScore=5000 -verbose=2 -smallClusterHub=memk \ -swap -syntenicNet \ -chainLinearGap=loose -bigClusterHub=swarm > swap.log 2>&1 & # real 59m45.336 cat fb.ornAna1.chainBosTau6Link.txt # 190401000 bases of 1842236818 (10.335%) in intersection # TBD for all chain/net: *** All done ! Elapsed time: 59m18s *** Make sure that goldenPath/ornAna1/vsBosTau6/README.txt is accurate. *** Add {chain,net}BosTau6 tracks to trackDb.ra if necessary. cd /cluster/data/ornAna1/bed ln -s blastz.bosTau6.swap lastz.bosTau6 ##################################################################### # susScr2 Pig BLASTZ/CHAIN/NET (DONE - 2011-05-24 - Chin) screen # use a screen to manage this multi-day job mkdir /hive/data/genomes/bosTau6/bed/lastzSusScr2.2011-05-24 cd /hive/data/genomes/bosTau6/bed/lastzSusScr2.2011-05-24 cat << '_EOF_' > DEF # Pig vs. Cow BLASTZ_M=50 # TARGET: Cow BosTau6 SEQ1_DIR=/scratch/data/bosTau6/bosTau6.2bit SEQ1_LEN=/scratch/data/bosTau6/chrom.sizes SEQ1_CHUNK=10000000 SEQ1_LAP=10000 # QUERY: Pig SusScr2 SEQ2_DIR=/scratch/data/susScr2/susScr2.2bit SEQ2_LEN=/scratch/data/susScr2/chrom.sizes SEQ2_CHUNK=10000000 SEQ2_LAP=0 BASE=/hive/data/genomes/bosTau6/bed/lastzSusScr2.2011-05-24 TMPDIR=/scratch/tmp '_EOF_' # << this line keeps emacs coloring happy time nice -n +19 doBlastzChainNet.pl -verbose=2 \ `pwd`/DEF \ -noLoadChainSplit -syntenicNet \ -workhorse=hgwdev -smallClusterHub=memk -bigClusterHub=swarm \ -chainMinScore=3000 -chainLinearGap=medium > do.log 2>&1 & # real 2198m44.401s cat fb.bosTau6.chainSusScr2Link.txt # 1306738899 bases of 2649682029 (49.317%) in intersection cd /hive/data/genomes/bosTau6/bed ln -s lastzSusScr2.2011-05-24 lastz.susScr2 # then swap mkdir /hive/data/genomes/susScr2/bed/blastz.bosTau6.swap cd /hive/data/genomes/susScr2/bed/blastz.bosTau6.swap time nice -n +19 doBlastzChainNet.pl -verbose=2 \ /hive/data/genomes/bosTau6/bed/lastzSusScr2.2011-05-24/DEF \ -swap -noLoadChainSplit -syntenicNet \ -workhorse=hgwdev -smallClusterHub=memk -bigClusterHub=pk \ -chainMinScore=3000 -chainLinearGap=medium > swap.log 2>&1 & # real 380m53.612s cat fb.susScr2.chainBosTau6Link.txt # 1475323050 bases of 2231298548 (66.119%) in intersection ######################################################################### # SWAP mm9 LASTZ (DONE - 2010-05-17 - Chin) # original alignment done at mm9.txt cd /hive/data/genomes/mm9/bed/lastzBosTau6.2011-05-17 cat fb.mm9.chainBosTau6Link.txt # 699351036 bases of 2620346127 (26.689%) in intersection # Create link cd /hive/data/genomes/mm9/bed ln -s lastzBosTau6.2011-05-17 lastz.bosTau6 # and the swap mkdir /hive/data/genomes/bosTau6/bed/blastz.mm9.swap cd /hive/data/genomes/bosTau6/bed/blastz.mm9.swap time nice -n +19 doBlastzChainNet.pl -verbose=2 \ /hive/data/genomes/mm9/bed/lastzBosTau6.2011-05-17/DEF \ -swap -syntenicNet \ -noLoadChainSplit \ -workhorse=hgwdev -smallClusterHub=memk -bigClusterHub=swarm \ -chainMinScore=3000 -chainLinearGap=medium > swap.log 2>&1 # real 53m5.237s cat fb.bosTau6.chainMm9Link.txt # 688894115 bases of 2649682029 (25.999%) in intersection cd /hive/data/genomes/bosTau6/bed ln -s blastz.mm9.swap lastz.mm9 ############################################################################ # SWAP hg19 LASTZ (DONE 2011-05-16 - Chin) # original alignment cd /hive/data/genomes/hg19/bed/lastzBosTau6.2011-05-16 cat fb.hg19.chainBosTau6Link.txt # 1370696434 bases of 2897316137 (47.309%) in intersection # and the swap mkdir /hive/data/genomes/bosTau6/bed/blastz.hg19.swap cd /hive/data/genomes/bosTau6/bed/blastz.hg19.swap time nice -n +19 doBlastzChainNet.pl -verbose=2 \ /hive/data/genomes/hg19/bed/lastzBosTau6.2011-05-16/DEF \ -swap -syntenicNet \ -noLoadChainSplit \ -workhorse=hgwdev -smallClusterHub=memk -bigClusterHub=swarm \ -chainMinScore=3000 -chainLinearGap=medium > swap.log 2>&1 # real 72m57.889s cat fb.bosTau6.chainHg19Link.txt # 1336966333 bases of 2649682029 (50.458%) in intersection cd /hive/data/genomes/bosTau6/bed ln -s blastz.hg19.swap lastz.hg19 ############################################################################ # SWAP canFam2 LASTZ (DONE 2011-05-25 - Chin) # original alignment cd /hive/data/genomes/canFam2/bed/lastzBosTau6.2011-05-25 cat fb.canFam2.chainBosTau6Link.txt # 1389712912 bases of 2384996543 (58.269%) in intersection # and the swap mkdir /hive/data/genomes/bosTau6/bed/blastz.canFam2.swap cd /hive/data/genomes/bosTau6/bed/blastz.canFam2.swap time nice -n +19 doBlastzChainNet.pl -verbose=2 \ /hive/data/genomes/canFam2/bed/lastzBosTau6.2011-05-25/DEF \ -swap -syntenicNet \ -noLoadChainSplit \ -workhorse=hgwdev -smallClusterHub=memk -bigClusterHub=swarm \ -chainMinScore=3000 -chainLinearGap=medium > swap.log 2>&1 # real 88m38.555s cat fb.bosTau6.chainCanFam2Link.txt # 1404315725 bases of 2649682029 (52.999%) in intersection cd /hive/data/genomes/bosTau6/bed ln -s blastz.canFam2.swap lastz.canFam2 ############################################################################ # SWAP rn4 LASTZ (DONE 2011-05-31 - Chin) cd /hive/data/genomes/rn4/bed/blastzBosTau6.2011-05-31 cat fb.rn4.chainBosTau6Link.txt # 660647959 bases of 2571531505 (25.691%) in intersection # and the swap mkdir /hive/data/genomes/bosTau6/bed/blastz.rn4.swap cd /hive/data/genomes/bosTau6/bed/blastz.rn4.swap time nice -n +19 doBlastzChainNet.pl -verbose=2 \ /hive/data/genomes/rn4/bed/blastzBosTau6.2011-05-31/DEF \ -swap -noLoadChainSplit -syntenicNet \ -workhorse=hgwdev -smallClusterHub=memk -bigClusterHub=swarm \ -chainMinScore=3000 -chainLinearGap=medium > swap.log 2>&1 & # real 49m43.881 cat fb.bosTau6.chainRn4Link.txt # 648232024 bases of 2649682029 (24.465%) in intersection cd /hive/data/genomes/bosTau6/bed ln -s blastz.rn4.swap lastz.rn4 ############################################################################ # SWAP monDom5 LASTZ (DONE 2011-05-31 - Chin) cd /hive/data/genomes/monDom5/bed/lastzBosTau6.2011-05-31 cat fb.monDom5.chainBosTau6Link.txt # 207703303 bases of 3501660299 (5.932%) in intersection # and the swap mkdir /hive/data/genomes/bosTau6/bed/blastz.monDom5.swap cd /hive/data/genomes/bosTau6/bed/blastz.monDom5.swap time nice -n +19 doBlastzChainNet.pl -verbose=2 \ /hive/data/genomes/monDom5/bed/lastzBosTau6.2011-05-31/DEF \ -swap -noLoadChainSplit -syntenicNet \ -workhorse=hgwdev -smallClusterHub=memk -bigClusterHub=swarm \ -chainMinScore=5000 -chainLinearGap=loose > swap.log 2>&1 & # real 52m10.229s cat fb.bosTau6.chainMonDom5Link.txt # 199980646 bases of 2649682029 (7.547%) in intersection cd /hive/data/genomes/bosTau6/bed ln -s blastz.monDom5.swap lastz.monDom5 ######################################################################### # all.joiner update, downloads and in pushQ - (working 2011-06-14 - Chin) cd $HOME/kent/src/hg/makeDb/schema # fixup all.joiner until this is a clean output joinerCheck -database=bosTau6 -all all.joiner mkdir /hive/data/genomes/bosTau6/goldenPath cd /hive/data/genomes/bosTau6/goldenPath makeDownloads.pl bosTau6 > do.log 2>&1 # now ready for pushQ entry mkdir /hive/data/genomes/bosTau6/pushQ cd /hive/data/genomes/bosTau6/pushQ makePushQSql.pl bosTau6 > bosTau6.pushQ.sql 2> stderr.out # check for errors in stderr.out, some are OK, e.g.: # WARNING: hgwdev does not have /gbdb/bosTau6/wib/gc5Base.wib # WARNING: hgwdev does not have /gbdb/bosTau6/wib/quality.wib # WARNING: hgwdev does not have /gbdb/bosTau6/bbi/quality.bw # WARNING: bosTau6 does not have seq # WARNING: bosTau6 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: # tableList # copy it to hgwbeta scp -p bosTau6.pushQ.sql hgwbeta:/tmp ssh hgwbeta cd /tmp hgsql qapushq < bosTau6.pushQ.sql # in that pushQ entry walk through each entry and see if the # sizes will set properly ############################################################################ # bosTau6 - Cow - Ensembl Genes version 64 (DONE - 2011-10-12 - hiram) ssh hgwdev cd /hive/data/genomes/bosTau6 cat << '_EOF_' > bosTau6.ensGene.ra # required db variable db bosTau6 # optional nameTranslation, the sed command that will transform # Ensemble names to UCSC names. With quotes just to make sure. nameTranslation "s/^\([0-9X][0-9]*\)/chr\1/; s/^MT/chrM/; s/^GJ\([0-9]*\).1/chrUn_GJ\1/" '_EOF_' # << happy emacs doEnsGeneUpdate.pl -ensVersion=64 bosTau6.ensGene.ra ssh hgwdev cd /hive/data/genomes/bosTau6/bed/ensGene.64 featureBits bosTau6 ensGene # 42264149 bases of 2649682029 (1.595%) in intersection ############################################################################ # ucscToEnsembl name translation table (DONE - 2011-10-12 - Hiram) mkdir /hive/data/genomes/bosTau6/bed/ucscToEnsembl cd /hive/data/genomes/bosTau6/bed/ucscToEnsembl cut -f1 ../../chrom.sizes | while read C do E=`echo ${C} | sed -e "s/^chrM/MT/; s/chrUn_GJ\([0-9]*\)/GJ\1.1/; s/^chr//;"` echo -e "${C}\t${E}" done > ucscToEnsembl.tab # determine longest UCSC name: cut -f1 ucscToEnsembl.tab | awk '{print length($0)}' | sort -nr | head -1 # 14 # use that size in this SQL definition index: 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(14)) ); '_EOF_' hgLoadSqlTab bosTau6 ucscToEnsembl ucscToEnsembl.sql ucscToEnsembl.tab ############################################################################# # LIFTOVER TO bosTau4 (DONE - 2011-10-20 - Chin) sh hgwdev cd /hive/data/genomes/bosTau6 ln -s jkStuff/bosTau6.11.ooc 11.ooc time nice -n +19 doSameSpeciesLiftOver.pl -verbose=2 \ -bigClusterHub=swarm -dbHost=hgwdev -workhorse=hgwdev \ bosTau6 bosTau4 > doLiftOverToBosTau4.log 2>&1 & # real 110m26.710s pushd /usr/local/apache/htdocs-hgdownload/goldenPath/bosTau6/liftOver/ md5sum bosTau6ToBosTau4.over.chain.gz >> md5sum.txt popd ############################################################################# # LIFTOVER TO bosTau4 ( redo DONE - 2012-05-11 - Chin) sh hgwdev cd /hive/data/genomes/bosTau6 ln -s jkStuff/bosTau6.11.ooc 11.ooc time nice -n +19 doSameSpeciesLiftOver.pl -verbose=2 \ -bigClusterHub=swarm -dbHost=hgwdev -workhorse=hgwdev \ bosTau6 bosTau4 > doLiftOverToBosTau4.log 2>&1 & # real 1590m52.986s pushd /usr/local/apache/htdocs-hgdownload/goldenPath/bosTau6/liftOver/ md5sum bosTau6ToBosTau4.over.chain.gz >> md5sum.txt popd ############################################################################# # LIFTOVER TO bosTau5 (DONE - 2012-05-11 - Chin) sh hgwdev cd /hive/data/genomes/bosTau6 ln -s jkStuff/bosTau6.11.ooc 11.ooc time nice -n +19 doSameSpeciesLiftOver.pl -verbose=2 \ -bigClusterHub=swarm -dbHost=hgwdev -workhorse=hgwdev \ bosTau6 bosTau5 > doLiftOverToBosTau5.log 2>&1 & # real 1567m22.050s pushd /usr/local/apache/htdocs-hgdownload/goldenPath/bosTau6/liftOver/ md5sum bosTau6ToBosTau5.over.chain.gz >> md5sum.txt popd ############################################################################# # LIFTOVER TO bosTau7 (DONE - 2012-05-11 - Chin) sh hgwdev cd /hive/data/genomes/bosTau6 ln -s jkStuff/bosTau6.11.ooc 11.ooc time nice -n +19 doSameSpeciesLiftOver.pl -verbose=2 \ -bigClusterHub=swarm -dbHost=hgwdev -workhorse=hgwdev \ bosTau6 bosTau7 > doLiftOverToBosTau7.log 2>&1 & # real 1594m27.003s pushd /usr/local/apache/htdocs-hgdownload/goldenPath/bosTau6/liftOver/ md5sum bosTau6ToBosTau7.over.chain.gz >> md5sum.txt popd ############################################################################## # SWAP mm10 LASTZ (DONE - 2012-06-19 - Chin) # original alignment done at mm10.txt cd /hive/data/genomes/mm10/bed/lastzBosTau6.2012-06-19 cat fb.mm10.chainBosTau6Link.txt # 700039696 bases of 2652783500 (26.389%) in intersection # set sym link to indicate this is the lastz for this genome: cd /hive/data/genomes/mm10/bed ln -s lastzBosTau6.2012-06-19 lastz.bosTau6 # swap mkdir /hive/data/genomes/bosTau6/bed/blastz.mm10.swap cd /hive/data/genomes/bosTau6/bed/blastz.mm10.swap time nice -n +19 doBlastzChainNet.pl -verbose=2 \ /hive/data/genomes/mm10/bed/lastzBosTau6.2012-06-19/DEF \ -swap -syntenicNet \ -workhorse=hgwdev -smallClusterHub=encodek -bigClusterHub=swarm \ -chainMinScore=3000 -chainLinearGap=medium > swap.log 2>&1 & # real 72m13.925s cat fb.bosTau6.chainMm10Link.txt # 688651806 bases of 2649682029 (25.990%) in intersection # set sym link to indicate this is the lastz for this genome: cd /hive/data/genomes/bosTau6/bed ln -s blastz.mm10.swap lastz.mm10 ######################################################################### # SWAP dog canFam3 LASTZ (DONE - 2012-06-24 - Chin) mkdir /hive/data/genomes/canFam3/bed/lastzBosTau6.2012-06-24 cd /hive/data/genomes/canFam3/bed/lastzBosTau6.2012-06-24 cat fb.canFam3.chainBosTau6Link.txt # 1387159926 bases of 2392715236 (57.974%) in intersection # Create link cd /hive/data/genomes/canFam3/bed ln -s lastzBosTau6.2012-06-24 lastz.bosTau6 # and the swap mkdir /hive/data/genomes/bosTau6/bed/blastz.canFam3.swap cd /hive/data/genomes/bosTau6/bed/blastz.canFam3.swap time nice -n +19 doBlastzChainNet.pl -verbose=2 \ /hive/data/genomes/canFam3/bed/lastzBosTau6.2012-06-24/DEF \ -swap -syntenicNet \ -noLoadChainSplit \ -workhorse=hgwdev -smallClusterHub=memk -bigClusterHub=swarm \ -chainMinScore=3000 -chainLinearGap=medium > swap.log 2>&1 & # real 119m54.003s cat fb.bosTau6.chainCanFam3Link.txt # 1399687351 bases of 2649682029 (52.825%) in intersection cd /hive/data/genomes/bosTau6/bed ln -s blastz.canFam3.swap lastz.canFam3 ############################################################################ # NUMTS TRACK (working 2013-02-05 - Chin) # copy and unzip # http://193.204.182.50/files/tracks.zip and # http://193.204.182.50/files/other%20NumtS%20tracks.zip # to /hive/data/outside/NumtS/2012/tracks mkdir /hive/data/genomes/bosTau6/bed/NumtS cd /hive/data/genomes/bosTau6/bed/NumtS cp /hive/data/outside/NumtS/2012/tracks/bosTau6/* . cat Tracks_nuclear_NumtS_BosTau6 | grep -v "^track" > bosTau6NumtS.bed cat tracks_NumtS_BT6_assembled | grep -v "^track" > bosTau6NumtSAssembled.bed cat Tracks_mit_NumtS_BosTau6 | grep -v "^track" > bosTau6NumtSMitochondrion.bed # load the 3 bed files to bosTau6 hgLoadBed bosTau6 numtSAssembled bosTau6NumtSAssembled.bed hgLoadBed bosTau6 numtS bosTau6NumtS.bed hgLoadBed bosTau6 numtSMitochondrion bosTau6NumtSMitochondrion.bed # Make /gbdb/ links and load bam mkdir /gbdb/bosTau6/NumtS ln -s `pwd`/BT-6_NumtS_seqs.sorted.bam{,.bai} /gbdb/bosTau6/NumtS/ hgBbiDbLink bosTau6 bamAllNumtSSorted /gbdb/bosTau6/NumtS/BT-6_NumtS_seqs.sorted.bam # setup trackDb for bosTau6 ############################################################################