51d2baad200a4ef6f339fc30ec67444117153e1e hiram Thu Jan 23 14:18:21 2020 -0800 now adding ncbiRefSeqSelectCurated track refs #24748 diff --git src/hg/utils/automation/doNcbiRefSeq.pl src/hg/utils/automation/doNcbiRefSeq.pl index 43900b1..67962c7 100755 --- src/hg/utils/automation/doNcbiRefSeq.pl +++ src/hg/utils/automation/doNcbiRefSeq.pl @@ -326,64 +326,79 @@ export annotationRelease=`zcat \$ncbiGffGz | head -100 | grep ^#.annotation-source | sed -e 's/.*annotation-source //'` if [ "\$annotationRelease" == "" ]; then export annotationRelease=\$asmId fi export versionDate=`ls -L --full-time \$ncbiGffGz | awk '{print \$6;}'` echo "\$annotationRelease (\$versionDate)" > ncbiRefSeqVersion.txt # this produces the genePred in NCBI coordinates # 8/23/17: gff3ToGenePred quits over illegal attribute SO_type... make it legal (so_type): zcat \$ncbiGffGz \\ | sed -re 's/([;\\t])SO_type=/\\1so_type=/;' \\ | gff3ToGenePred $warnOnly -refseqHacks -attrsOut=\$asmId.attrs.txt \\ -unprocessedRootsOut=\$asmId.unprocessedRoots.txt stdin \$asmId.gp genePredCheck \$asmId.gp +zcat \$ncbiGffGz \\ + | egrep 'tag=(RefSeq|MANE) Select' | cut -f9- | tr ';' '\\n' \\ + | grep 'Name=' | grep -v NP_ | cut -d= -f2 | sort -u \\ + > \$asmId.refseqSelectTranscripts.txt + # extract labels from semi-structured text in gbff COMMENT/description sections: zcat \$downloadDir/\${asmId}_rna.gbff.gz \\ | (grep ' :: ' || true) \\ | perl -wpe 's/\\s+::.*//; s/^\\s+//;' \\ | sort -u \\ > pragmaLabels.txt # extract cross reference text for refLink \$gff3ToRefLink \$downloadDir/\$asmId.raFile.txt \$ncbiGffGz pragmaLabels.txt 2> \$db.refLink.stderr.txt \\ | sort > \$asmId.refLink.tab # converting the NCBI coordinates to UCSC coordinates liftUp -extGenePred -type=.gp stdout $localLiftFile drop \$asmId.gp \\ | gzip -c > \$asmId.\$db.gp.gz $genePredCheckDb \$asmId.\$db.gp.gz # curated subset of all genes (zegrep "^[NY][MRP]_" \$asmId.\$db.gp.gz || true) > \$db.curated.gp # may not be any curated genes if [ ! -s \$db.curated.gp ]; then rm -f \$db.curated.gp +else + cat \$db.curated.gp | fgrep -f \$asmId.refseqSelectTranscripts.txt - \\ + > \$db.refseqSelect.curated.gp + # may not be any refseqSelect.curated genes + if [ ! -s \$db.refseqSelect.curated.gp ]; then + rm -f \$db.refseqSelect.curated.gp + fi fi # predicted subset of all genes (zegrep "^X[MR]_" \$asmId.\$db.gp.gz || true) > \$db.predicted.gp # not curated or predicted subset of all genes, the left overs (zegrep -v "^[NXY][MRP]_" \$asmId.\$db.gp.gz || true) > \$db.other.gp # curated and predicted without leftovers: (zegrep "^[NXY][MRP]_" \$asmId.\$db.gp.gz || true) > \$db.ncbiRefSeq.gp if [ -s \$db.curated.gp ]; then $genePredCheckDb \$db.curated.gp + if [ -s \$db.refseqSelect.curated.gp ]; then + $genePredCheckDb \$db.refseqSelect.curated.gp + fi fi if [ -s \$db.predicted.gp ]; then $genePredCheckDb \$db.predicted.gp fi if [ -s \$db.other.gp ]; then $genePredCheckDb \$db.other.gp fi # join the refLink metadata with curated+predicted names cut -f1 \$db.ncbiRefSeq.gp | sort -u > \$asmId.\$db.name.list join -t\$'\\t' \$asmId.\$db.name.list \$asmId.refLink.tab > \$asmId.\$db.ncbiRefSeqLink.tab # Make bigBed with attributes in extra columns for ncbiRefSeqOther: twoBitInfo $dbTwoBit stdout | sort -k2,2n > \$db.chrom.sizes genePredToBed \$db.other.gp stdout | sort -k1,1 -k2n,2n > \$db.other.bed @@ -493,30 +508,45 @@ ### curated only if present if [ -s process/\$db.curated.gp ]; then genePredToBigGenePred process/\$db.curated.gp stdout | sort -k1,1 -k2,2n > \$db.ncbiRefSeqCurated.bigGp bedToBigBed -type=bed12+8 -tab -as=bigGenePred.as -extraIndex=name \\ \$db.ncbiRefSeqCurated.bigGp \$db.chrom.sizes \\ \$db.ncbiRefSeqCurated.bb rm -f \$db.ncbiRefSeqCurated.bigGp bigBedInfo \$db.ncbiRefSeqCurated.bb | egrep "^itemCount:|^basesCovered:" \\ | sed -e 's/,//g' > \$db.ncbiRefSeqCurated.stats.txt LC_NUMERIC=en_US /usr/bin/printf "# ncbiRefSeqCurated %s %'d %s %'d\\n" `cat \$db.ncbiRefSeqCurated.stats.txt` | xargs echo ~/kent/src/hg/utils/automation/gpToIx.pl process/\$db.curated.gp \\ | sort -u > \$asmId.ncbiRefSeqCurated.ix.txt ixIxx \$asmId.ncbiRefSeqCurated.ix.txt \$asmId.ncbiRefSeqCurated.ix{,x} rm -f \$asmId.ncbiRefSeqCurated.ix.txt +### and refseqSelect if exists (a subset of curated) + if [ -s process/\$db.refseqSelect.curated.gp ]; then + genePredToBigGenePred process/\$db.refseqSelect.curated.gp stdout | sort -k1,1 -k2,2n > \$db.ncbiRefSeqSelectCurated.bigGp + bedToBigBed -type=bed12+8 -tab -as=bigGenePred.as -extraIndex=name \\ + \$db.ncbiRefSeqSelectCurated.bigGp \$db.chrom.sizes \\ + \$db.ncbiRefSeqSelectCurated.bb + rm -f \$db.ncbiRefSeqSelectCurated.bigGp + bigBedInfo \$db.ncbiRefSeqSelectCurated.bb | egrep "^itemCount:|^basesCovered:" \\ + | sed -e 's/,//g' > \$db.ncbiRefSeqSelectCurated.stats.txt + LC_NUMERIC=en_US /usr/bin/printf "# ncbiRefSeqSelectCurated %s %'d %s %'d\\n" `cat \$db.ncbiRefSeqSelectCurated.stats.txt` | xargs echo + ~/kent/src/hg/utils/automation/gpToIx.pl process/\$db.refseqSelect.curated.gp \\ + | sort -u > \$asmId.ncbiRefSeqSelectCurated.ix.txt + ixIxx \$asmId.ncbiRefSeqSelectCurated.ix.txt \$asmId.ncbiRefSeqSelectCurated.ix{,x} + rm -f \$asmId.ncbiRefSeqSelectCurated.ix.txt + fi fi ### predicted only if present if [ -s process/\$db.predicted.gp ]; then genePredToBigGenePred process/\$db.predicted.gp stdout | sort -k1,1 -k2,2n > \$db.ncbiRefSeqPredicted.bigGp bedToBigBed -type=bed12+8 -tab -as=bigGenePred.as -extraIndex=name \\ \$db.ncbiRefSeqPredicted.bigGp \$db.chrom.sizes \\ \$db.ncbiRefSeqPredicted.bb rm -f \$db.ncbiRefSeqPredicted.bigGp bigBedInfo \$db.ncbiRefSeqPredicted.bb | egrep "^itemCount:|^basesCovered:" \\ | sed -e 's/,//g' > \$db.ncbiRefSeqPredicted.stats.txt LC_NUMERIC=en_US /usr/bin/printf "# ncbiRefSeqPredicted %s %'d %s %'d\\n" `cat \$db.ncbiRefSeqPredicted.stats.txt` | xargs echo ~/kent/src/hg/utils/automation/gpToIx.pl process/\$db.predicted.gp \\ | sort -u > \$asmId.ncbiRefSeqPredicted.ix.txt ixIxx \$asmId.ncbiRefSeqPredicted.ix.txt \$asmId.ncbiRefSeqPredicted.ix{,x} @@ -593,30 +623,34 @@ bedToBigBed -type=bed12+13 -tab -as=bigPsl.as -extraIndex=name \\ \$asmId.bigPsl \$db.chrom.sizes \$asmId.bigPsl.bb rm -f \$asmId.bigPsl _EOF_ ); } else { $bossScript->add(<<_EOF_ # loading the genePred tracks, all genes in one, and subsets hgLoadGenePred -genePredExt \$db ncbiRefSeq process/\$db.ncbiRefSeq.gp $genePredCheckDb ncbiRefSeq if [ -s process/\$db.curated.gp ]; then hgLoadGenePred -genePredExt \$db ncbiRefSeqCurated process/\$db.curated.gp $genePredCheckDb ncbiRefSeqCurated + if [ -s process/\$db.refseqSelect.curated.gp ]; then + hgLoadGenePred -genePredExt \$db ncbiRefSeqSelect process/\$db.refseqSelect.curated.gp + $genePredCheckDb ncbiRefSeqSelect + fi fi if [ -s process/\$db.predicted.gp ]; then hgLoadGenePred -genePredExt \$db ncbiRefSeqPredicted process/\$db.predicted.gp $genePredCheckDb ncbiRefSeqPredicted fi mkdir -p $gbdbDir ln -f -s `pwd`/process/\$db.other.bb $gbdbDir/ncbiRefSeqOther.bb hgBbiDbLink \$db ncbiRefSeqOther $gbdbDir/ncbiRefSeqOther.bb ln -f -s `pwd`/process/ncbiRefSeqOther.ix{,x} $gbdbDir/ ln -f -s `pwd`/process/ncbiRefSeqVersion.txt $gbdbDir/ # select only coding genes to have CDS records