27d9d88ca49e60586c2d0b211971c8a360d862ac hiram Thu Sep 8 23:00:48 2022 -0700 sequenceName alias should properly be called assembly no redmine diff --git src/hg/makeDb/doc/panPan3/initialBuild.txt src/hg/makeDb/doc/panPan3/initialBuild.txt index 52840d4..70f90bd 100644 --- src/hg/makeDb/doc/panPan3/initialBuild.txt +++ src/hg/makeDb/doc/panPan3/initialBuild.txt @@ -1,1075 +1,1075 @@ # for emacs: -*- mode: sh; -*- # This file describes browser build for the panPan3 # Can use existing photograph (otherwise find one before starting here) ######################################################################### # Initial steps, reuse existing photograph (DONE - 2020-06-12 - Hiram) # Might be a different individual, requests made # To start this initialBuild.txt document, from a previous assembly document: mkdir ~/kent/src/hg/makeDb/doc/panPan3 cd ~/kent/src/hg/makeDb/doc/panPan3 sed -e 's/gorGor6/panPan3/g; s/GorGor6/PanPan3/g; s/DONE/TBD/g;' \ ../gorGor6/initialBuild.txt > initialBuild.txt mkdir -p /hive/data/genomes/panPan3/refseq cd /hive/data/genomes/panPan3 # Can use existing photograph cp -p ../panPan2/photoReference.txt ./ cat photoReference.txt photoCreditURL http://a-z-animals.com/animals/bonobo/pictures/2955/ photoCreditName Photo courtesy of Kabir Bakie ## download from NCBI cd /hive/data/genomes/panPan3/refseq time rsync -L -a -P --stats \ rsync://ftp.ncbi.nlm.nih.gov/genomes/all/GCF/013/052/645/GCF_013052645.1_Mhudiblu_PPA_v0/ ./ # sent 3,203 bytes received 6,058,365,683 bytes 56,886,092.83 bytes/sec # total size is 6,056,874,807 speedup is 1.00 # real 1m45.870s # this information is from the top of # panPan3/refseq/*_assembly_report.txt # (aka: panPan3/refseq/GCF_013052645.1_Mhudiblu_PPA_v0_assembly_report.txt ) # Assembly name: Mhudiblu_PPA_v0 # Organism name: Pan paniscus (pygmy chimpanzee) # Isolate: Mhudiblu (Carbone #601152) # Sex: female # Taxid: 9597 # BioSample: SAMN11123633 # BioProject: PRJNA526933 # Submitter: University of Washington # Date: 2020-05-15 # Assembly type: haploid # Release type: major # Assembly level: Chromosome # Genome representation: full # WGS project: SSBP01 # Assembly method: FALCON v. (git hash: 91e700c4) Nov 2015 # Expected final version: no # Genome coverage: 85.0x # Sequencing technology: PacBio RSII; Illumina # GenBank assembly accession: GCA_013052645.1 # RefSeq assembly accession: GCF_013052645.1 # RefSeq assembly and GenBank assemblies identical: no # ## Assembly-Units: ## GenBank Unit Accession RefSeq Unit Accession Assembly-Unit name ## GCA_013052655.1 GCF_013052655.1 Primary Assembly ## GCF_000029195.1 non-nuclear # check assembly size for later reference: faSize G*0_genomic.fna.gz # 3051901337 bases (36551040 N's 3015350297 real 1802830734 upper # 1212519563 lower) in 4293 sequences in 1 files # Total size: mean 710901.8 sd 9447605.7 min 209 (NW_023258064.1) # max 224621958 (NC_048239.1) median 32647 # %39.73 masked total, %40.21 masked real ############################################################################# # establish config.ra file (DONE - Hiram - 2020-06-12) cd /hive/data/genomes/panPan3 ~/kent/src/hg/utils/automation/prepConfig.pl panPan3 mammal bonobo \ refseq/*_assembly_report.txt > panPan3.config.ra # verify it is sane: diff panPan3.config.ra ../panPan2/panPan2.config.ra # fixup commonName, from: commonName Pygmy chimpanzee to commonName Bonobo genomeCladePriority 35 to genomeCladePriority 15 orderKey 17442 to orderKey 2623 # verify it really does look sane cat panPan3.config.ra # config parameters for makeGenomeDb.pl: db panPan3 clade mammal genomeCladePriority 15 scientificName Pan paniscus commonName Bonobo assemblyDate May 2020 assemblyLabel University of Washington assemblyShortLabel Mhudiblu_PPA_v0 orderKey 2623 # mitochondrial sequence included in refseq release # mitoAcc NC_001644.1 mitoAcc none fastaFiles /hive/data/genomes/panPan3/ucsc/*.fa.gz agpFiles /hive/data/genomes/panPan3/ucsc/*.agp # qualFiles none dbDbSpeciesDir bonobo photoCreditURL http://a-z-animals.com/animals/bonobo/pictures/2955/ photoCreditName Photo courtesy of Kabir Bakie ncbiGenomeId 10729 ncbiAssemblyId 7133581 ncbiAssemblyName Mhudiblu_PPA_v0 ncbiBioProject 526933 ncbiBioSample SAMN11123633 genBankAccessionID GCF_013052645.1 taxId 9597 ############################################################################# # setup UCSC named files (DONE - 2020-06-12 - Hiram) mkdir /hive/data/genomes/panPan3/ucsc cd /hive/data/genomes/panPan3/ucsc # check for duplicate sequences: time faToTwoBit -noMask ../refseq/G*0_genomic.fna.gz refseq.2bit # real 0m44.404s twoBitDup refseq.2bit # no output is a good result, otherwise, would have to eliminate duplicates # the scripts creating the fasta here will be using this refseq.2bit file # remove it later time ~/kent/src/hg/utils/automation/ucscCompositeAgp.pl \ ../refseq/G*0_genomic.fna.gz \ ../refseq/*_assembly_structure/Primary_Assembly NC_048239.1 chr1 NC_048240.1 chr2A NC_048241.1 chr2B NC_048242.1 chr3 NC_048243.1 chr4 NC_048244.1 chr5 NC_048245.1 chr6 NC_048246.1 chr7 NC_048247.1 chr8 NC_048248.1 chr9 NC_048249.1 chr10 NC_048250.1 chr11 NC_048251.1 chr12 NC_048252.1 chr13 NC_048253.1 chr14 NC_048254.1 chr15 NC_048255.1 chr16 NC_048256.1 chr17 NC_048257.1 chr18 NC_048258.1 chr19 NC_048259.1 chr20 NC_048260.1 chr21 NC_048261.1 chr22 NC_048262.1 chrX # real 11m31.473s time ~/kent/src/hg/utils/automation/unplacedWithChroms.pl \ ../refseq/*_assembly_structure/Primary_Assembly # processed 4268 sequences into chrUn.fa.gz # real 0m50.843s # there are *no* unlocalized sequences time ~/kent/src/hg/utils/automation/unlocalizedWithChroms.pl \ ../refseq/*_assembly_structure/Primary_Assembly # bash syntax here mitoAcc=`grep "^# mitoAcc" ../panPan3.config.ra | awk '{print $NF}'` printf "# mitoAcc %s\n" "$mitoAcc" # mitoAcc NC_001644.1 zcat \ ../refseq/*_assembly_structure/non-nuclear/assem*/AGP/chrMT.comp.agp.gz \ | grep -v "^#" | sed -e "s/^$mitoAcc/chrM/;" > chrM.agp cat chrM.agp # chrM 1 16563 1 O NC_001644.1 1 16563 + printf ">chrM\n" > chrM.fa twoBitToFa -noMask refseq.2bit:$mitoAcc stdout | grep -v "^>" >> chrM.fa gzip chrM.fa faSize chrM.fa.gz # 16563 bases (0 N's 16563 real 16563 upper 0 lower) in 1 sequences in 1 files # verify fasta and AGPs agree time faToTwoBit *.fa.gz test.2bit # real 0m55.327s cat *.agp | checkAgpAndFa stdin test.2bit 2>&1 | tail -4 # All AGP and FASTA entries agree - both files are valid # and no sequence lost from orginal: twoBitToFa test.2bit stdout | faSize stdin # 3051901337 bases (36551040 N's 3015350297 real 3015350297 upper # 0 lower) in 4293 sequences in 1 files # same numbers as above (except for upper/lower masking) # 3051901337 bases (36551040 N's 3015350297 real 1802830734 upper # 1212519563 lower) in 4293 sequences in 1 files # no longer need these temporary 2bit files rm test.2bit refseq.2bit ############################################################################# # Initial database build (DONE - 2020-06-12 - Hiram) # verify sequence and AGP are OK: cd /hive/data/genomes/panPan3 time (makeGenomeDb.pl -workhorse=hgwdev -dbHost=hgwdev -fileServer=hgwdev \ -stop=agp panPan3.config.ra) > agp.log 2>&1 # real 2m32.140s # then finish it off: time (makeGenomeDb.pl -workhorse=hgwdev -dbHost=hgwdev \ -fileServer=hgwdev -continue=db panPan3.config.ra) > db.log 2>&1 # real 16m31.749s # fixing trackDb step: (new files required in trackDb source fetch) time (~/kent/src/hg/utils/automation/makeGenomeDb.pl -workhorse=hgwdev \ -dbHost=hgwdev -continue=trackDb -fileServer=hgwdev \ panPan3.config.ra) > trackDb.log 2>&1 # real 0m6.470s # check in the trackDb files created in TemporaryTrackDbCheckout/ # and add panPan3 to trackDb/makefile # temporary symlink until masked sequence is available cd /hive/data/genomes/panPan3 ln -s `pwd`/panPan3.unmasked.2bit /gbdb/panPan3/panPan3.2bit ############################################################################# # check gap table vs NCBI gap file (DONE - 2020-06-12 - Hiram) mkdir /hive/data/genomes/panPan3/bed/gap cd /hive/data/genomes/panPan3/bed/gap zgrep -v "^#" ../../refseq/G*_gaps.txt.gz \ | awk '{printf "%s\t%d\t%d\t%s_%s\n", $1,$2-1,$3,$5,$6}' \ | sort -k1,1 -k2,2n > refseq.gap.bed # type survey: cut -f4 *.bed | sort | uniq -c # 64 between_scaffolds_na # 619 within_scaffold_map # how much defined by NCBI: awk '{print $3-$2}' *.bed | ave stdin | grep -w total # total 36551040.000000 # how much in the gap table: hgsql -e 'select * from gap;' panPan3 | awk '{print $4-$3}' \ | ave stdin | grep -w total # total 36551040.000000 # equal amounts, no need to adjust the gap table ############################################################################## # cpgIslands on UNMASKED sequence (DONE - 2020-06-12 - Hiram) mkdir /hive/data/genomes/panPan3/bed/cpgIslandsUnmasked cd /hive/data/genomes/panPan3/bed/cpgIslandsUnmasked time (doCpgIslands.pl -dbHost=hgwdev -bigClusterHub=ku -buildDir=`pwd` \ -tableName=cpgIslandExtUnmasked \ -maskedSeq=/hive/data/genomes/panPan3/panPan3.unmasked.2bit \ -workhorse=hgwdev -smallClusterHub=ku panPan3) > do.log 2>&1 # real 4m48.939s cat fb.panPan3.cpgIslandExtUnmasked.txt # 29049864 bases of 3015350297 (0.963%) in intersection ############################################################################# # cytoBandIdeo - (DONE - 2020-06-12 - Hiram) mkdir /hive/data/genomes/panPan3/bed/cytoBand cd /hive/data/genomes/panPan3/bed/cytoBand makeCytoBandIdeo.csh panPan3 ############################################################################# # run up idKeys files for chromAlias/ncbiRefSeq (DONE - 2020-06-12 - Hiram) mkdir /hive/data/genomes/panPan3/bed/idKeys cd /hive/data/genomes/panPan3/bed/idKeys time (doIdKeys.pl \ -twoBit=/hive/data/genomes/panPan3/panPan3.unmasked.2bit \ -buildDir=`pwd` panPan3) > do.log 2>&1 & # real 22m22.873s cat panPan3.keySignature.txt # 25a9ac3ecf4f0b51c5a213bd22ddd73f ############################################################################# # gapOverlap (DONE - 2020-06-12 - Hiram) mkdir /hive/data/genomes/panPan3/bed/gapOverlap cd /hive/data/genomes/panPan3/bed/gapOverlap time (doGapOverlap.pl \ -twoBit=/hive/data/genomes/panPan3/panPan3.unmasked.2bit panPan3 ) \ > do.log 2>&1 & # real 1m49.512s # there is only one: wc -l bed.tab # 1 bed.tab cut -f2- bed.tab #chr7 144505525 144507630 chr7:144505526-144507630 803 + # 144505525 144507630 0 2 803,803 0,1302 cat fb.panPan3.gapOverlap.txt # 1606 bases of 3051901337 (0.000%) in intersection ############################################################################# # tandemDups (DONE - 2020-06-12 - Hiram) mkdir /hive/data/genomes/panPan3/bed/tandemDups cd /hive/data/genomes/panPan3/bed/tandemDups time (~/kent/src/hg/utils/automation/doTandemDup.pl \ -twoBit=/hive/data/genomes/panPan3/panPan3.unmasked.2bit panPan3) \ > do.log 2>&1 & # real 188m31.667s cat fb.panPan3.tandemDups.txt # 158351117 bases of 3051901337 (5.189%) in intersection bigBedInfo panPan3.tandemDups.bb | sed -e 's/^/# /;' # version: 4 # fieldCount: 13 # hasHeaderExtension: yes # isCompressed: yes # isSwapped: 0 # extraIndexCount: 0 # itemCount: 3,019,527 # primaryDataSize: 77,104,242 # primaryIndexSize: 283,772 # zoomLevels: 9 # chromCount: 4151 # basesCovered: 1,653,271,326 # meanDepth (of bases covered): 15.653055 # minDepth: 1.000000 # maxDepth: 425.000000 # std of depth: 34.287645 ######################################################################### # ucscToINSDC and ucscToRefSeq table/track (DONE - 2020-06-12 - Hiram) # construct idKeys for the refseq sequence mkdir /hive/data/genomes/panPan3/refseq/idKeys cd /hive/data/genomes/panPan3/refseq/idKeys faToTwoBit ../GCF_*0_genomic.fna.gz panPan3.refSeq.2bit time (doIdKeys.pl -buildDir=`pwd` \ -twoBit=`pwd`/panPan3.refSeq.2bit refseqPanPan3) > do.log 2>&1 & # real 6m11.312s cat refseqPanPan3.keySignature.txt # 25a9ac3ecf4f0b51c5a213bd22ddd73f # and the genbank sequence needs keys too: mkdir /hive/data/genomes/panPan3/refseq/idKeysGenbank cd /hive/data/genomes/panPan3/refseq/idKeysGenbank faToTwoBit /hive/data/outside/ncbi/genomes/GCA/013/052/645/GCA_013052645.1_Mhudiblu_PPA_v0/GCA_013052645.1_Mhudiblu_PPA_v0_genomic.fna.gz panPan3.genbank.2bit time (doIdKeys.pl -buildDir=`pwd` \ -twoBit=`pwd`/panPan3.genbank.2bit genbankPanPan3) > do.log 2>&1 & # real 12m24.596s cat genbankPanPan3.keySignature.txt # 8fd9bb74aebe493c87e8a961a031ba3f mkdir /hive/data/genomes/panPan3/bed/chromAlias cd /hive/data/genomes/panPan3/bed/chromAlias join -t$'\t' ../idKeys/panPan3.idKeys.txt \ ../../refseq/idKeysGenbank/genbankPanPan3.idKeys.txt | cut -f2- \ | sort -k1,1 | join -t$'\t' <(sort -k1,1 ../../chrom.sizes) - \ | awk '{printf "%s\t0\t%d\t%s\n", $1, $2, $3}' \ | sort -k1,1 -k2,2n > ucscToINSDC.bed join -t$'\t' ../idKeys/panPan3.idKeys.txt \ ../../refseq/idKeys/refseqPanPan3.idKeys.txt | cut -f2- \ | sort -k1,1 | join -t$'\t' <(sort -k1,1 ../../chrom.sizes) - \ | awk '{printf "%s\t0\t%d\t%s\n", $1, $2, $3}' \ | sort -k1,1 -k2,2n > ucscToRefSeq.bed # should be same line counts throughout: wc -l * ../../chrom.sizes # 4292 ucscToINSDC.bed # 4293 ucscToRefSeq.bed # 4293 ../../chrom.sizes # need to find the accession for the INSDC equivalent to chrM: egrep chrM * # refSeq.bed:chrM 0 16563 NC_001644.1 # lookup that accession at NCBI Entrez: D38116.1 # and add to ucscToINSDC.bed: printf "chrM\t0\t16563\tD38116.1\n" >> ucscToINSDC.bed # verify: grep chrM * ucscToINSDC.bed:chrM 0 16563 D38116.1 ucscToRefSeq.bed:chrM 0 16563 NC_001644.1 export chrSize=`cut -f1 ucscToINSDC.bed | awk '{print length($0)}' | sort -n | tail -1` echo $chrSize # 20 # use the $chrSize in this sed sed -e "s/21/$chrSize/" $HOME/kent/src/hg/lib/ucscToINSDC.sql \ | hgLoadSqlTab panPan3 ucscToINSDC stdin ucscToINSDC.bed # should be the same for ucscToRefSeq: export chrSize=`cut -f1 ucscToRefSeq.bed | awk '{print length($0)}' | sort -n | tail -1` echo $chrSize # 26 sed -e "s/21/$chrSize/" $HOME/kent/src/hg/lib/ucscToINSDC.sql \ | sed -e 's/INSDC/RefSeq/g;' \ | hgLoadSqlTab panPan3 ucscToRefSeq stdin ucscToRefSeq.bed # should be quiet for all OK checkTableCoords panPan3 # should cover %100 entirely: featureBits -countGaps panPan3 ucscToINSDC # 3051901337 bases of 3051901337 (100.000%) in intersection featureBits -countGaps panPan3 ucscToRefSeq # 3051901337 bases of 3051901337 (100.000%) in intersection ######################################################################### # add chromAlias table (DONE - 2020-06-12 - Hiram) mkdir /hive/data/genomes/panPan3/bed/chromAlias cd /hive/data/genomes/panPan3/bed/chromAlias hgsql -N -e 'select chrom,name from ucscToRefSeq;' panPan3 \ | sort -k1,1 > ucsc.refseq.tab hgsql -N -e 'select chrom,name from ucscToINSDC;' panPan3 \ | sort -k1,1 > ucsc.genbank.tab grep -v "^#" \ ../../refseq/G*0_assembly_report.txt \ | awk -F$'\t' '{printf "%s\t%s\n", $7,$1}' \ - | sort > refseq.sequenceName.txt + | sort > refseq.assembly.txt # do *not* need the ones that have identical names to UCSC names - join -2 2 -t$'\t' refseq.sequenceName.txt <(sort -k2,2 ucsc.refseq.tab) \ + join -2 2 -t$'\t' refseq.assembly.txt <(sort -k2,2 ucsc.refseq.tab) \ | awk -F$'\t' '{printf "%s\t%s\n", $3, $2}' \ - | awk -F$'\t' '$1 != $2' | sort > ucsc.sequenceName.tab + | awk -F$'\t' '$1 != $2' | sort > ucsc.assembly.tab wc -l *.tab | sed -e 's/^/# /;' # 4293 ucsc.genbank.tab # 4293 ucsc.refseq.tab -# 4271 ucsc.sequenceName.tab +# 4271 ucsc.assembly.tab ~/kent/src/hg/utils/automation/chromAlias.pl ucsc.*.tab \ > panPan3.chromAlias.tab -for t in refseq genbank sequenceName +for t in refseq genbank assembly do c0=`cat ucsc.$t.tab | wc -l` c1=`grep $t panPan3.chromAlias.tab | wc -l` ok="OK" if [ "$c0" -ne "$c1" ]; then ok="ERROR" fi printf "# checking $t: $c0 =? $c1 $ok\n" done # checking refseq: 4293 =? 4293 OK # checking genbank: 4293 =? 4293 OK -# # checking sequenceName: 4271 =? 4271 OK +# # checking assembly: 4271 =? 4271 OK # verify chrM is here properly: grep chrM panPan3.chromAlias.tab | sed -e 's/^/# /;' # D38116.1 chrM genbank -# MT chrM sequenceName +# MT chrM assembly # NC_001644.1 chrM refseq hgLoadSqlTab panPan3 chromAlias ~/kent/src/hg/lib/chromAlias.sql \ panPan3.chromAlias.tab ######################################################################### # fixup search rule for assembly track/gold table (DONE - 2020-06-12 - Hiram) cd ~/kent/src/hg/makeDb/trackDb/bonobo/panPan3 # preview prefixes and suffixes: hgsql -N -e "select frag from gold;" panPan3 \ | sed -e 's/[0-9][0-9]*//;' | sort | uniq -c | sed -e 's/^/# /;' # 1 NC_.1 # 4975 SSBP.1 # implies a rule: '[NS][CS][B0-9_][P0-9][0-9]+(\.[0-9]+)?' # verify this rule will find them all and eliminate them all: hgsql -N -e "select frag from gold;" panPan3 | wc -l # 4976 hgsql -N -e "select frag from gold;" panPan3 \ | egrep -e '[NS][CS][B0-9_][P0-9][0-9]+(\.[0-9]+)?' | wc -l # 4976 hgsql -N -e "select frag from gold;" panPan3 \ | egrep -v -e '[NS][CS][B0-9_][P0-9][0-9]+(\.[0-9]+)?' | wc -l # 0 # hence, add to trackDb/rhesus/panPan3/trackDb.ra searchTable gold shortCircuit 1 termRegex [NS][CS][B0-9_][P0-9][0-9]+(\.[0-9]+)? query select chrom,chromStart,chromEnd,frag from %s where frag like '%s%%' searchPriority 8 # verify searches work in the position box ########################################################################## # running repeat masker (DONE - 2020-06-12 - Hiram) mkdir /hive/data/genomes/panPan3/bed/repeatMasker cd /hive/data/genomes/panPan3/bed/repeatMasker time (doRepeatMasker.pl -buildDir=`pwd` \ -bigClusterHub=ku -dbHost=hgwdev -workhorse=hgwdev \ -smallClusterHub=ku panPan3) > do.log 2>&1 # real 420m57.419s cat faSize.rmsk.txt # 3051901337 bases (36551040 N's 3015350297 real 1503559806 upper # 1511790491 lower) in 4293 sequences in 1 files # Total size: mean 710901.8 sd 9447605.7 min 209 (chrUn_NW_023258064v1) # max 224621958 (chr1) median 32647 # %49.54 masked total, %50.14 masked real # Original masking from NCBI: # %39.73 masked total, %40.21 masked real egrep -i "versi|relea" do.log # RepeatMasker version development-$Id: RepeatMasker,v 1.332 2017/04/17 19:01:11 rhubley Exp $ # CC Dfam_Consensus RELEASE 20181026; * # CC RepBase RELEASE 20181026; # February 01 2017 (open-4-0-8) 1.332 version of RepeatMasker time featureBits -countGaps panPan3 rmsk # 1511789076 bases of 3051901337 (49.536%) in intersection # real 0m51.367s # why is it different than the faSize above ? # because rmsk masks out some N's as well as bases, the faSize count above # separates out the N's from the bases, it doesn't show lower case N's # faster way to get the same result on high contig count assemblies: time hgsql -N -e 'select genoName,genoStart,genoEnd from rmsk;' panPan3 \ | bedSingleCover.pl stdin | ave -col=4 stdin | grep "^total" # total 1511789076.000000 # real 0m22.577s ########################################################################## # running simple repeat (DONE - 2020-06-12 - Hiram) mkdir /hive/data/genomes/panPan3/bed/simpleRepeat cd /hive/data/genomes/panPan3/bed/simpleRepeat time (doSimpleRepeat.pl -buildDir=`pwd` -bigClusterHub=ku \ -dbHost=hgwdev -workhorse=hgwdev -smallClusterHub=ku \ -trf409=6 panPan3) > do.log 2>&1 # real 328m38.113s cat fb.simpleRepeat # 256626302 bases of 3015350297 (8.511%) in intersection cd /hive/data/genomes/panPan3 # using the Window Masker result: cd /hive/data/genomes/panPan3 # twoBitMask bed/windowMasker/panPan3.cleanWMSdust.2bit \ # -add bed/simpleRepeat/trfMask.bed panPan3.2bit # you can safely ignore the warning about fields >= 13 # add to rmsk after it is done: twoBitMask panPan3.rmsk.2bit \ -add bed/simpleRepeat/trfMask.bed panPan3.2bit # you can safely ignore the warning about fields >= 13 twoBitToFa panPan3.2bit stdout | faSize stdin > faSize.panPan3.2bit.txt cat faSize.panPan3.2bit.txt # 3051901337 bases (36551040 N's 3015350297 real 1501255487 upper # 1514094810 lower) in 4293 sequences in 1 files # Total size: mean 710901.8 sd 9447605.7 min 209 (chrUn_NW_023258064v1) # max 224621958 (chr1) median 32647 # %49.61 masked total, %50.21 masked real rm /gbdb/panPan3/panPan3.2bit ln -s `pwd`/panPan3.2bit /gbdb/panPan3/panPan3.2bit ######################################################################### # CREATE MICROSAT TRACK (DONE - 2020-06-13 - Hiram) ssh hgwdev mkdir /cluster/data/panPan3/bed/microsat cd /cluster/data/panPan3/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 panPan3 microsat microsat.bed # Read 27378 elements of size 4 from microsat.bed ########################################################################## ## WINDOWMASKER (DONE - 2020-06-12 - Hiram) mkdir /hive/data/genomes/panPan3/bed/windowMasker cd /hive/data/genomes/panPan3/bed/windowMasker time (doWindowMasker.pl -buildDir=`pwd` -workhorse=hgwdev \ -dbHost=hgwdev panPan3) > do.log 2>&1 # real 131m20.699s # Masking statistics cat faSize.panPan3.cleanWMSdust.txt # 3051901337 bases (36551040 N's 3015350297 real 1785076452 upper # 1230273845 lower) in 4293 sequences in 1 files # Total size: mean 710901.8 sd 9447605.7 min 209 (chrUn_NW_023258064v1) # max 224621958 (chr1) median 32647 # %40.31 masked total, %40.80 masked real cat fb.panPan3.rmsk.windowmaskerSdust.txt # 882436977 bases of 3051901337 (28.914%) in intersection ########################################################################## # cpgIslands - (DONE - 2020-06-15 - Hiram) mkdir /hive/data/genomes/panPan3/bed/cpgIslands cd /hive/data/genomes/panPan3/bed/cpgIslands time (doCpgIslands.pl -dbHost=hgwdev -bigClusterHub=ku \ -workhorse=hgwdev -smallClusterHub=ku panPan3) > do.log 2>&1 # real 4m8.773s cat fb.panPan3.cpgIslandExt.txt # 21946791 bases of 3015350297 (0.728%) in intersection ############################################################################## # genscan - (DONE - 2020-06-15 - Hiram) mkdir /hive/data/genomes/panPan3/bed/genscan cd /hive/data/genomes/panPan3/bed/genscan time (doGenscan.pl -buildDir=`pwd` -workhorse=hgwdev -dbHost=hgwdev \ -bigClusterHub=ku panPan3) > do.log 2>&1 # real 250m19.909s cat fb.panPan3.genscan.txt # 53138526 bases of 3015350297 (1.762%) in intersection cat fb.panPan3.genscanSubopt.txt # 51253170 bases of 3015350297 (1.700%) in intersection ######################################################################### # Create kluster run files (DONE - 2020-06-15 - Hiram) # numerator is panPan3 gapless bases "real" as reported by: featureBits -noRandom -noHap panPan3 gap # 30700245 bases of 2756992444 (1.114%) 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 \( 2756992444 / 2861349177 \) \* 1024 # ( 2756992444 / 2861349177 ) * 1024 = 986.653529 # ==> use -repMatch=950 according to size scaled down from 1024 for human. # and rounded down to nearest 50 cd /hive/data/genomes/panPan3 time blat panPan3.2bit \ /dev/null /dev/null -tileSize=11 -makeOoc=jkStuff/panPan3.11.ooc \ -repMatch=950 # Wrote 40152 overused 11-mers to jkStuff/panPan3.11.ooc # panPan2 repMatch=900: # Wrote 34618 overused 11-mers to jkStuff/panPan2.11.ooc # panPan1 repMatch=950: # Wrote 31225 overused 11-mers to jkStuff/panPan1.11.ooc # check non-bridged gaps to see what the typical size is: hgsql -N \ -e 'select * from gap where bridge="no" order by size;' panPan3 \ | sort -k7,7nr | ave -col=7 stdin # min 100.000000 # max 100.000000 # they are all 100 sized, 64 gaps # minimum gap size is 100 and produces a reasonable number of lifts gapToLift -verbose=2 -minGap=100 panPan3 jkStuff/panPan3.nonBridged.lft \ -bedFile=jkStuff/panPan3.nonBridged.bed wc -l jkStuff/panPan3.nonBri* # 4357 jkStuff/panPan3.nonBridged.bed # 4357 jkStuff/panPan3.nonBridged.lft ############################################################################## # GENBANK AUTO UPDATE (DONE - 2020-06-15 - Hiram) ssh hgwdev cd $HOME/kent/src/hg/makeDb/genbank git pull # /cluster/data/genbank/data/organism.lst shows: # organism mrnaCnt estCnt refSeqCnt # Pan paniscus 557 0 49 # edit etc/genbank.conf to add panPan3 just before panPan2 # panPan3 (bonobo - Pan paniscus - GCF_013052645.1 Mhudiblu_PPA_v0 15 May 2020) panPan3.serverGenome = /hive/data/genomes/panPan3/panPan3.2bit panPan3.ooc = /hive/data/genomes/panPan3/jkStuff/panPan3.11.ooc panPan3.lift = /hive/data/genomes/panPan3/jkStuff/panPan3.nonBridged.lft panPan3.perChromTables = no panPan3.refseq.mrna.native.pslCDnaFilter = ${finished.refseq.mrna.native.pslCDnaFilter} panPan3.refseq.mrna.xeno.pslCDnaFilter = ${finished.refseq.mrna.xeno.pslCDnaFilter} panPan3.genbank.mrna.native.pslCDnaFilter = ${finished.genbank.mrna.native.pslCDnaFilter} panPan3.genbank.mrna.xeno.pslCDnaFilter = ${finished.genbank.mrna.xeno.pslCDnaFilter} panPan3.genbank.est.native.pslCDnaFilter = ${finished.genbank.est.native.pslCDnaFilter} panPan3.genbank.est.xeno.pslCDnaFilter = ${finished.genbank.est.xeno.pslCDnaFilter} panPan3.downloadDir = panPan3 # defaults yes genbank.mrna.native, genbank.est.native.load, # refseq.mrna.native, refseq.mrna.xeno # DO NOT NEED genbank.mrna.xeno except for human, mouse panPan3.genbank.est.native.load = no # verify the files specified exist before checking in the file: grep ^panPan3 etc/genbank.conf | grep hive | awk '{print $NF}' | xargs ls -og # -rw-rw-r-- 1 160616 Jun 15 09:11 /hive/data/genomes/panPan3/jkStuff/panPan3.11.ooc # -rw-rw-r-- 1 254979 Jun 15 09:13 /hive/data/genomes/panPan3/jkStuff/panPan3.nonBridged.lft # -rw-rw-r-- 1 795168110 Jun 13 10:06 /hive/data/genomes/panPan3/panPan3.2bit git commit -m "Added panPan3 bonobo; refs #25720" etc/genbank.conf git push # update /cluster/data/genbank/: make etc-update # enable daily alignment and update of hgwdev cd ~/kent/src/hg/makeDb/genbank git pull # add panPan3 to: # etc/hgwdev.dbs etc/align.dbs git commit -m "Added panPan3 - bonobo refs #25720" etc/hgwdev.dbs etc/align.dbs git push make etc-update # wait a few days for genbank magic to take place, the tracks will # appear ######################################################################## # lastz/chain/net swap human/hg38 (DONE - 2020-06-15 - Hiram) # original alignment cd /hive/data/genomes/hg38/bed/lastzPanPan3.2020-06-15 cat fb.hg38.chainPanPan3Link.txt # 2897225010 bases of 3110768607 (93.135%) in intersection cat fb.hg38.chainSynPanPan3Link.txt # 2880116277 bases of 3110768607 (92.585%) in intersection time (doRecipBest.pl -load -workhorse=hgwdev -buildDir=`pwd` \ hg38 panPan3) > rbest.log 2>&1 & # real 68m45.187s cat fb.hg38.chainRBest.PanPan3.txt # 2727319248 bases of 3110768607 (87.673%) in intersection # and for the swap: mkdir /hive/data/genomes/panPan3/bed/blastz.hg38.swap cd /hive/data/genomes/panPan3/bed/blastz.hg38.swap time (doBlastzChainNet.pl -verbose=2 \ /hive/data/genomes/hg38/bed/lastzPanPan3.2020-06-15/DEF \ -swap -chainMinScore=3000 -chainLinearGap=medium \ -workhorse=hgwdev -smallClusterHub=hgwdev -bigClusterHub=ku \ -syntenicNet) > swap.log 2>&1 # real 73m50.456s cat fb.panPan3.chainHg38Link.txt # 2777869191 bases of 3015350297 (92.124%) in intersection cat fb.panPan3.chainSynHg38Link.txt # 2768018012 bases of 3015350297 (91.798%) in intersection time (doRecipBest.pl -load -workhorse=hgwdev -buildDir=`pwd` \ panPan3 hg38) > rbest.log 2>&1 & # real 66m52.340s cat fb.panPan3.chainRBest.Hg38.txt # 2732415393 bases of 3015350297 (90.617%) in intersection ########################################################################### # lastz/chain/net swap mouse/mm10 (DONE - 2020-06-15 - Hiram) # original alignment cd /hive/data/genomes/mm10/bed/lastzPanPan3.2020-06-15 cat fb.mm10.chainPanPan3Link.txt # 935579510 bases of 2652783500 (35.268%) in intersection cat fb.mm10.chainSynPanPan3Link.txt # 888900388 bases of 2652783500 (33.508%) in intersection time (doRecipBest.pl -load -workhorse=hgwdev -buildDir=`pwd` mm10 panPan3) \ > rbest.log 2>&1 & # real 329m28.051s cat fb.mm10.chainRBest.PanPan3.txt # 890894306 bases of 2652783500 (33.583%) in intersection # and for the swap: mkdir /hive/data/genomes/panPan3/bed/blastz.mm10.swap cd /hive/data/genomes/panPan3/bed/blastz.mm10.swap time (doBlastzChainNet.pl -verbose=2 \ /hive/data/genomes/mm10/bed/lastzPanPan3.2020-06-15/DEF \ -swap -chainMinScore=3000 -chainLinearGap=medium \ -workhorse=hgwdev -smallClusterHub=hgwdev -bigClusterHub=ku \ -syntenicNet) > swap.log 2>&1 # real 55m23.982s cat fb.panPan3.chainMm10Link.txt # 954214151 bases of 3015350297 (31.645%) in intersection cat fb.panPan3.chainSynMm10Link.txt # 887980807 bases of 3015350297 (29.449%) in intersection time (doRecipBest.pl -load -workhorse=hgwdev -buildDir=`pwd` panPan3 mm10) \ > rbest.log 2>&1 # real 305m11.756s cat fb.panPan3.chainRBest.Mm10.txt # 889360051 bases of 3015350297 (29.494%) in intersection ############################################################################## # augustus gene track (DONE - 2020-06-15 - Hiram) mkdir /hive/data/genomes/panPan3/bed/augustus cd /hive/data/genomes/panPan3/bed/augustus time (doAugustus.pl -buildDir=`pwd` -bigClusterHub=ku \ -species=human -dbHost=hgwdev \ -workhorse=hgwdev panPan3) > do.log 2>&1 # real 124m57.738s cat fb.panPan3.augustusGene.txt # 55220084 bases of 3015350297 (1.831%) in intersection ######################################################################### # ncbiRefSeq (DONE - 2020-06-13 - Hiram) mkdir /hive/data/genomes/panPan3/bed/ncbiRefSeq.2020-06-13 cd /hive/data/genomes/panPan3/bed/ncbiRefSeq.2020-06-13 time (~/kent/src/hg/utils/automation/doNcbiRefSeq.pl -buildDir=`pwd` \ -bigClusterHub=ku -dbHost=hgwdev \ -fileServer=hgwdev -smallClusterHub=hgwdev -workhorse=hgwdev \ GCF_013052645.1_Mhudiblu_PPA_v0 panPan3) > do.log 2>&1 & # real 5m43.954s cat fb.ncbiRefSeq.panPan3.txt # 91843495 bases of 3015350297 (3.046%) in intersection # add: include ../../refSeqComposite.ra alpha # to the bonobo/panPan3/trackDb.ra to turn on the track in the browser featureBits -enrichment panPan3 refGene ncbiRefSeq # refGene 0.004%, ncbiRefSeq 3.046%, both 0.004%, cover 99.99%, enrich 32.83x featureBits -enrichment panPan3 ncbiRefSeq refGene # ncbiRefSeq 3.046%, refGene 0.004%, both 0.004%, cover 0.11%, enrich 32.83x featureBits -enrichment panPan3 ncbiRefSeqCurated refGene # ncbiRefSeqCurated 0.004%, refGene 0.004%, both 0.003%, cover 90.24%, enrich 25764.34x featureBits -enrichment panPan3 refGene ncbiRefSeqCurated # refGene 0.004%, ncbiRefSeqCurated 0.004%, both 0.003%, cover 99.48%, enrich 25764.34x ######################################################################### # LIFTOVER TO panPan2 (DONE - 2020-06-15 - Hiram) ssh hgwdev mkdir /hive/data/genomes/panPan3/bed/blat.panPan2.2020-06-15 cd /hive/data/genomes/panPan3/bed/blat.panPan2.2020-06-15 doSameSpeciesLiftOver.pl -verbose=2 \ -debug -bigClusterHub=ku -dbHost=hgwdev -workhorse=hgwdev \ -ooc=/hive/data/genomes/panPan3/jkStuff/panPan3.11.ooc \ panPan3 panPan2 time (doSameSpeciesLiftOver.pl -verbose=2 \ -bigClusterHub=ku -dbHost=hgwdev -workhorse=hgwdev \ -ooc=/hive/data/genomes/panPan3/jkStuff/panPan3.11.ooc \ panPan3 panPan2) > doLiftOverToPanPan2.log 2>&1 # real 273m17.026s # see if the liftOver menus function in the browser from panPan3 to panPan2 ######################################################################### # LIFTOVER TO panPan1 (DONE - 2020-06-15 - Hiram) ssh hgwdev mkdir /hive/data/genomes/panPan3/bed/blat.panPan1.2020-06-15 cd /hive/data/genomes/panPan3/bed/blat.panPan1.2020-06-15 doSameSpeciesLiftOver.pl -verbose=2 \ -debug -bigClusterHub=ku -dbHost=hgwdev -workhorse=hgwdev \ -ooc=/hive/data/genomes/panPan3/jkStuff/panPan3.11.ooc \ panPan3 panPan1 time (doSameSpeciesLiftOver.pl -verbose=2 \ -bigClusterHub=ku -dbHost=hgwdev -workhorse=hgwdev \ -ooc=/hive/data/genomes/panPan3/jkStuff/panPan3.11.ooc \ panPan3 panPan1) > doLiftOverToPanPan1.log 2>&1 # real 275m17.469s # see if the liftOver menus function in the browser from panPan3 to panPan1 ######################################################################### # BLATSERVERS ENTRY (DONE - 2020-06-16 - 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 ("panPan3", "blat1a", "17900", "1", "0"); \ INSERT INTO blatServers (db, host, port, isTrans, canPcr) \ VALUES ("panPan3", "blat1a", "17901", "0", "1");' \ hgcentraltest # test it with some sequence ############################################################################ ## reset default position similar to panPan2 found via blat of FOXP2 protein ## (DONE - 2020-06-16 - Hiram) ssh hgwdev hgsql -e 'update dbDb set defaultPos="chr7:106336524-106733879" where name="panPan3";' hgcentraltest ############################################################################## # crispr whole genome (DONE - 2020-06-15 - Hiram) mkdir /hive/data/genomes/panPan3/bed/crisprAll cd /hive/data/genomes/panPan3/bed/crisprAll # the large shoulder argument will cause the entire genome to be scanned # this takes a while for a new genome to get the bwa indexing done time (~/kent/src/hg/utils/automation/doCrispr.pl -verbose=2 -stop=ranges \ panPan3 ncbiRefSeq -shoulder=250000000 -tableName=crisprAll \ -fileServer=hgwdev -buildDir=`pwd` -smallClusterHub=hgwdev \ -bigClusterHub=ku \ -workhorse=hgwdev) > ranges.log 2>&1 # real 80m34.161s time (~/kent/src/hg/utils/automation/doCrispr.pl -verbose=2 \ -continue=guides -stop=specScores panPan3 ncbiRefSeq \ -shoulder=250000000 -tableName=crisprAll -fileServer=hgwdev \ -buildDir=`pwd` -smallClusterHub=hgwdev -bigClusterHub=ku \ -workhorse=hgwdev) > specScores.log 2>&1 # Number of specScores: 232496133 # real 6369m47.263s cat guides/run.time | sed -e 's/^/# /;' # Completed: 100 of 100 jobs # CPU time in finished jobs: 17860s 297.67m 4.96h 0.21d 0.001 y # IO & Wait Time: 389s 6.48m 0.11h 0.00d 0.000 y # Average job time: 182s 3.04m 0.05h 0.00d # Longest finished job: 575s 9.58m 0.16h 0.01d # Submission to last job: 583s 9.72m 0.16h 0.01d cat specScores/run.time | sed -e 's/^/# /;' # Completed: 3104861 of 3104861 jobs # CPU time in finished jobs: 310874778s 5181246.29m 86354.10h 3598.09d 9.858 y # IO & Wait Time: 8751069s 145851.16m 2430.85h 101.29d 0.277 y # Average job time: 103s 1.72m 0.03h 0.00d # Longest finished job: 414s 6.90m 0.12h 0.00d # Submission to last job: 314519s 5241.98m 87.37h 3.64d ### remember to get back to hgwdev to run this time (~/kent/src/hg/utils/automation/doCrispr.pl -verbose=2 \ -continue=effScores -stop=load panPan3 ncbiRefSeq \ -shoulder=250000000 -tableName=crisprAll -fileServer=hgwdev \ -buildDir=`pwd` -smallClusterHub=hgwdev -bigClusterHub=ku \ -workhorse=hgwdev) > load.log 2>&1 # real 1633m51.680s cat effScores/run.time | sed -e 's/^/# /;' # Completed: 28661 of 28661 jobs # CPU time in finished jobs: 13143908s 219065.14m 3651.09h 152.13d 0.417 y # IO & Wait Time: 68503s 1141.71m 19.03h 0.79d 0.002 y # Average job time: 461s 7.68m 0.13h 0.01d # Longest finished job: 10393s 173.22m 2.89h 0.12d # Submission to last job: 54285s 904.75m 15.08h 0.63d cat offTargets/run.time | sed -e 's/^/# /;' # Completed: 155244 of 155244 jobs # CPU time in finished jobs: 2205628s 36760.46m 612.67h 25.53d 0.070 y # IO & Wait Time: 825117s 13751.95m 229.20h 9.55d 0.026 y # Average job time: 20s 0.33m 0.01h 0.00d # Longest finished job: 157s 2.62m 0.04h 0.00d # Submission to last job: 5940s 99.00m 1.65h 0.07d # cleaning 2021-04-24 - Hiram time (~/kent/src/hg/utils/automation/doCrispr.pl -verbose=2 \ -continue=cleanup panPan3 -tableName=crisprAll -fileServer=hgwdev \ -buildDir=`pwd` -smallClusterHub=hgwdev -bigClusterHub=ku \ -workhorse=hgwdev) > cleanup.log 2>&1 # real 473m37.404s ######################################################################### # all.joiner update, downloads and in pushQ - (WORKING - 2020-06-15 - Hiram) cd $HOME/kent/src/hg/makeDb/schema # verify all the business is done for release ~/kent/src/hg/utils/automation/verifyBrowser.pl panPan3 # 64 tables in database panPan3 - Bonobo, Pan paniscus # verified 62 tables in database panPan3, 2 extra tables, 24 optional tables # NCBI RefSeq genes 10 optional tables # chainNetRBestHg38 3 optional tables # chainNetRBestMm10 3 optional tables # chainNetSynHg38 3 optional tables # chainNetSynMm10 3 optional tables # gapOverlap 1 optional tables # tandemDups 1 optional tables # 1 crisprAllRanges - extra table # 2 crisprAllTargets - extra table # 9 genbank tables found # verified 29 required tables, 0 missing tables # hg38 chainNet to panPan3 found 3 required tables # mm10 chainNet to panPan3 found 3 required tables # hg38 chainNet RBest and syntenic to panPan3 found 6 optional tables # mm10 chainNet RBest and syntenic to panPan3 found 3 optional tables # liftOver to previous versions: 2, from previous versions: 2 # fixup all.joiner until this is a clean output joinerCheck -database=panPan3 -tableCoverage all.joiner joinerCheck -database=panPan3 -times all.joiner joinerCheck -database=panPan3 -keys all.joiner # when clean, check in: git commit -m 'adding rules for panPan3 refs #25720' all.joiner git push # run up a 'make alpha' in hg/hgTables to get this all.joiner file # into the hgwdev/genome-test system cd /hive/data/genomes/panPan3 time (~/kent/src/hg/utils/automation/makeDownloads.pl panPan3) > downloads.log 2>&1 # real 17m56.213s # now ready for pushQ entry mkdir /hive/data/genomes/panPan3/pushQ cd /hive/data/genomes/panPan3/pushQ time ($HOME/kent/src/hg/utils/automation/makePushQSql.pl -redmineList panPan3) > panPan3.pushQ.sql 2> stderr.out # real 13m55.634s # remove the tandemDups and gapOverlap from the file list: sed -i -e "/tandemDups/d" redmine.panPan3.table.list sed -i -e "/Tandem Dups/d" redmine.panPan3.releaseLog.txt sed -i -e "/gapOverlap/d" redmine.panPan3.table.list sed -i -e "/Gap Overlaps/d" redmine.panPan3.releaseLog.txt # check for errors in stderr.out, some are OK, e.g.: # WARNING: panPan3 does not have seq # WARNING: panPan3 does not have extFile # verify the file list does correctly match to files cat redmine.panPan3.file.list | while read L do eval ls $L > /dev/null done # should be silent, missing files will show as errors # verify database tables, how many to expect: wc -l redmine.panPan3.table.list # 53 redmine.panPan3.table.list # how many actual: awk -F'.' '{printf "hgsql -N %s -e '"'"'show table status like \"%s\";'"'"'\n", $1, $2}' redmine.panPan3.table.list | sh | wc -l # 53 # would be a smaller number actual if some were missing # add the path names to the listing files in the redmine issue # in the three appropriate entry boxes: # /hive/data/genomes/panPan3/pushQ/redmine.panPan3.file.list # /hive/data/genomes/panPan3/pushQ/redmine.panPan3.releaseLog.txt # /hive/data/genomes/panPan3/pushQ/redmine.panPan3.table.list #########################################################################