c5dc7044d76420830c68a063b195599f655e8e5d hiram Fri Jun 12 12:03:39 2020 -0700 running repeat masker refs #25720 diff --git src/hg/makeDb/doc/panPan3/initialBuild.txt src/hg/makeDb/doc/panPan3/initialBuild.txt index 1feb8de..f883b86 100644 --- src/hg/makeDb/doc/panPan3/initialBuild.txt +++ src/hg/makeDb/doc/panPan3/initialBuild.txt @@ -221,167 +221,175 @@ # 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 -XXX - running - Fri Jun 12 10:29:45 PDT 2020 - # real 17m28.042s + # 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 (TBD - 2019-11-19 - Hiram) +# 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 - # 220 between_scaffolds_na - # 639 within_scaffold_map + # 64 between_scaffolds_na + # 619 within_scaffold_map # how much defined by NCBI: awk '{print $3-$2}' *.bed | ave stdin | grep -w total - # total 45844299.000000 + # 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 45844299.000000 + # total 36551040.000000 # equal amounts, no need to adjust the gap table ############################################################################## -# cpgIslands on UNMASKED sequence (TBD - 2019-11-19 - Hiram) +# 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 4m13.285s + # real 4m48.939s cat fb.panPan3.cpgIslandExtUnmasked.txt - # 28001209 bases of 2999027915 (0.934%) in intersection + # 29049864 bases of 3015350297 (0.963%) in intersection ############################################################################# -# cytoBandIdeo - (TBD - 2019-11-19 - Hiram) +# 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 (TBD - 2019-11-19 - Hiram) +# 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 2m48.092s + # real 22m22.873s cat panPan3.keySignature.txt - # 10c42ee6ea4a90775c5da9d8b83854aa + # 25a9ac3ecf4f0b51c5a213bd22ddd73f ############################################################################# -# gapOverlap (TBD - 2019-11-19 - Hiram) +# 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 1m50.075s + # real 1m49.512s # there is only one: wc -l bed.tab # 1 bed.tab cut -f2- bed.tab -chrX 66519024 66677016 chrX:66519025-66677016 1000 + 66519024 66677016 0 2 1000,1000 0,156992 +#chr7 144505525 144507630 chr7:144505526-144507630 803 + +# 144505525 144507630 0 2 803,803 0,1302 cat fb.panPan3.gapOverlap.txt - # 2000 bases of 3044872214 (0.000%) in intersection + # 1606 bases of 3051901337 (0.000%) in intersection ############################################################################# -# tandemDups (TBD - 2019-11-19 - Hiram) +# 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 & +XXX - running - Fri Jun 12 11:33:27 PDT 2020 # real 188m34.598s cat fb.panPan3.tandemDups.txt # 155315479 bases of 3044872214 (5.101%) in intersection bigBedInfo panPan3.tandemDups.bb | sed -e 's/^/# /;' # version: 4 # fieldCount: 13 # hasHeaderExtension: yes # isCompressed: yes # isSwapped: 0 # extraIndexCount: 0 # itemCount: 2,822,307 # primaryDataSize: 72,710,994 # primaryIndexSize: 292,560 # zoomLevels: 9 # chromCount: 5335 # basesCovered: 1,635,503,835 # meanDepth (of bases covered): 14.396921 # minDepth: 1.000000 # maxDepth: 381.000000 # std of depth: 29.341113 ######################################################################### -# ucscToINSDC and ucscToRefSeq table/track (TBD - 2019-11-19 - Hiram) +# 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 2m50.723s + # real 6m11.312s cat refseqPanPan3.keySignature.txt - # 10c42ee6ea4a90775c5da9d8b83854aa + # 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/genbank/vertebrate_mammalian/Gorilla_gorilla/all_assembly_versions/GCA_008122165.1_Kamilah_GGO_v0/GCA_008122165.1_Kamilah_GGO_v0_genomic.fna.gz panPan3.genbank.2bit + 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 & +XXX - running - Fri Jun 12 12:03:06 PDT 2020 # real 3m11.098s cat genbankPanPan3.keySignature.txt # 84734b343949ddf1e28b453d25d3ddf7 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 \ @@ -460,68 +468,69 @@ fi printf "# checking $t: $c0 =? $c1 $ok\n" done # checking refseq: 5486 =? 5486 OK # checking genbank: 5486 =? 5486 OK # verify chrM is here properly: grep chrM panPan3.chromAlias.tab # NC_011120.1 chrM refseq # X93347.1 chrM genbank hgLoadSqlTab panPan3 chromAlias ~/kent/src/hg/lib/chromAlias.sql \ panPan3.chromAlias.tab ######################################################################### -# fixup search rule for assembly track/gold table (TBD - 2019-11-19 - Hiram) - cd ~/kent/src/hg/makeDb/trackDb/gorilla/panPan3 +# 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 - 1 NC_.1 - 6344 SRLZ.1 + | sed -e 's/[0-9][0-9]*//;' | sort | uniq -c | sed -e 's/^/# /;' +# 1 NC_.1 +# 4975 SSBP.1 - # implies a rule: '[NS][CR][L0-9_][Z0-9][0-9]+(\.[0-9]+)?' + # 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 - # 6345 + # 4976 hgsql -N -e "select frag from gold;" panPan3 \ - | egrep -e '[NS][CR][L0-9_][Z0-9][0-9]+(\.[0-9]+)?' | wc -l - # 6345 + | 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][CR][L0-9_][Z0-9][0-9]+(\.[0-9]+)?' | wc -l + | 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][CR][L0-9_][Z0-9][0-9]+(\.[0-9]+)? +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 (TBD - 2019-11-19 - Hiram) +# 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 +XXX - running - Fri Jun 12 11:39:31 PDT 2020 # real 415m10.888s cat faSize.rmsk.txt # 3044872214 bases (45844299 N's 2999027915 real 1492721019 upper # 1506306896 lower) in 5486 sequences in 1 files # Total size: mean 555025.9 sd 8254410.8 min 202 (chrUn_NW_022154607v1) # max 219763114 (chr1) median 29244 # %49.47 masked total, %50.23 masked real egrep -i "versi|relea" do.log # RepeatMasker version development-$Id: RepeatMasker,v 1.332 2017/04/17 19:01:11 rhubley Exp $ # grep version of RepeatMasker$ /hive/data/staging/data/RepeatMasker/RepeatMasker # February 01 2017 (open-4-0-8) 1.332 version of RepeatMasker # grep RELEASE /hive/data/staging/data/RepeatMasker/Libraries/RepeatMaskerLib.embl # CC Dfam_Consensus RELEASE 20181026; * @@ -530,37 +539,38 @@ time featureBits -countGaps panPan3 rmsk # 1506305481 bases of 3044872214 (49.470%) in intersection # real 0m41.857s # 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 1506305481.000000 # real 0m22.271s ########################################################################## -# running simple repeat (TBD - 2019-11-19 - Hiram) +# 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 +XXX - running - Fri Jun 12 11:40:59 PDT 2020 # real 243m33.198s cat fb.simpleRepeat # 260449789 bases of 2999027915 (8.684%) 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