99599c7130790107ff0de9f043930da6aa7fddf1
angie
  Mon Nov 16 16:35:58 2020 -0800
Scripts for automating SARS-CoV-2 Phylogeny tracks (refs #26530): fetching
sequences and metadata from several public sources, mapping GISAID IDs to
public seq IDs, downloading the latest release of the phylogenetic tree from
github.com/roblanf/sarscov2phylo/ , making VCFs from GISAID and public
sequences, and using github.com/yatisht/usher to resolve ambiguous alleles,
make protobuf files for hgPhyloPlace, and add public sequences that have not
been mapped to GISAID sequences to the sarscov2phylo tree for a comprehensive
public tree+VCF.

This is still not fully otto-mated because certain crucial inputs like
GISAID sequences still must be downloaded using a web browser, but the goal
is to automate as much as possible and maybe someday have it fully cron-driven.

There are two main top-level scripts which call other scripts, which may in turn
call scripts, in this hierarchy:

updateIdMapping.sh
getCogUk.sh
getNcbi.sh
searchAllSarsCov2BioSample.sh
bioSampleIdToText.sh
bioSampleTextToTab.pl
gbMetadataAddBioSample.pl
fixNcbiFastaNames.pl

updateSarsCov2Phylo.sh
getRelease.sh
processRelease.sh
cladeLineageColors.pl
mapPublic.sh
extractUnmappedPublic.sh
addUnmappedPublic.sh

many of the above:
util.sh

publicCredits.sh will hopefully be folded into updateSarsCov2Phylo.sh when I
figure out how to automate fetching of author/institution metadata from NCBI
and COG-UK.

diff --git src/hg/utils/otto/sarscov2phylo/mapPublic.sh src/hg/utils/otto/sarscov2phylo/mapPublic.sh
new file mode 100755
index 0000000..299f41f
--- /dev/null
+++ src/hg/utils/otto/sarscov2phylo/mapPublic.sh
@@ -0,0 +1,94 @@
+#!/bin/bash
+set -beEu -x -o pipefail
+
+#	Do not modify this script, modify the source tree copy:
+#	kent/src/hg/utils/otto/sarscov2phylo/mapPublic.sh
+
+usage() {
+    echo "usage: $0 releaseLabel problematic_sites_sarsCov2.vcf epiToPublic"
+}
+
+if [ $# != 3 ]; then
+  usage
+  exit 1
+fi
+
+releaseLabel=$1
+problematicSitesVcf=$2
+epiToPublic=$3
+
+usherDir=~angie/github/usher
+usher=$usherDir/build/usher
+matToVcf=$usherDir/build/matToVcf
+
+echo "releaseLabel=$releaseLabel problematicSitesVcf=$problematicSitesVcf epiToPublic=$epiToPublic"
+
+# Now extract public-sequence-only subtree and VCF.  Start by making new names with shortened year:
+tawk '{if ($4 == "") { print $1, $3 "|" $2;} else { print $1, $3 "|" $2 "|" $4;} }' $epiToPublic \
+| sed -re 's/20([0-9][0-9])(-[0-9-]+)?$/\1\2/; s/, /_/g; s/[, ]/_/g;' \
+    > epiToPublicName
+wc -l epiToPublicName
+
+join -t$'\t' tree.renaming epiToPublicName | cut -f 2,3 > treeToPublic
+wc -l treeToPublic
+
+phyloRenameAndPrune ft_SH.reroot.collapsed.nh treeToPublic ft_SH.reroot.public.nh
+sed -re 's/,/,\n/g' ft_SH.reroot.public.nh | wc -l
+
+# Run usher to collapse (having removed GISAID-only variants) and make protobuf for public set too,
+# on unmasked VCF, and use matToVcf to get ambig-resolved VCF from that for hgTracks
+time vcfRenameAndPrune gisaid-$releaseLabel.unfiltered.vcf.gz treeToPublic stdout \
+| gzip -c \
+    > public-$releaseLabel.unfiltered.vcf.gz
+time $usher -c -u \
+    -v public-$releaseLabel.unfiltered.vcf.gz \
+    -t ft_SH.reroot.public.nh \
+    -o public-$releaseLabel.notMasked.pb
+time $matToVcf -i public-$releaseLabel.notMasked.pb \
+    -v public-$releaseLabel.vcf
+bgzip -f public-$releaseLabel.vcf
+tabix -p vcf public-$releaseLabel.vcf.gz
+
+# Run usher on masked public VCF to get protobuf for usher/hgPhyloPlace:
+time vcfFilter -excludeVcf=mask.vcf public-$releaseLabel.unfiltered.vcf.gz \
+| gzip -c \
+    > public-$releaseLabel.unfiltered.masked.vcf.gz
+time $usher -c -u \
+    -v public-$releaseLabel.unfiltered.masked.vcf.gz \
+    -t ft_SH.reroot.public.nh \
+    -o public-$releaseLabel.masked.pb
+mv uncondensed-final-tree.nh ft_SH.reroot.collapsed.public-$releaseLabel.nwk
+
+#*** Should we get rid of ambigToN... or keep it around as an option?  Maybe inferring missing
+#*** info is too ambitious/optimistic for some purposes?
+
+# Make allele-frequency-filtered versions
+# Disregard pipefail just for this pipe because head prevents the cat from completing:
+set +o pipefail
+sampleCount=$(zcat public-$releaseLabel.vcf.gz | head | grep ^#CHROM | sed -re 's/\t/\n/g' \
+              | tail -n+10 | wc -l)
+set -o pipefail
+minAc001=$(( (($sampleCount + 999) / 1000) ))
+vcfFilter -minAc=$minAc001 -rename public-$releaseLabel.vcf.gz \
+    > public-$releaseLabel.minAf.001.vcf
+wc -l public-$releaseLabel.minAf.001.vcf
+bgzip -f public-$releaseLabel.minAf.001.vcf
+tabix -p vcf public-$releaseLabel.minAf.001.vcf.gz
+
+minAc01=$(( (($sampleCount + 99) / 100) ))
+vcfFilter -minAc=$minAc01 -rename public-$releaseLabel.minAf.001.vcf.gz \
+    > public-$releaseLabel.minAf.01.vcf
+wc -l public-$releaseLabel.minAf.01.vcf
+bgzip -f public-$releaseLabel.minAf.01.vcf
+tabix -p vcf public-$releaseLabel.minAf.01.vcf.gz
+
+# And the lineage colors:
+# TODO: add Nextstrain clade colors
+for source in gisaid lineage nextstrain; do
+    Source=${source^}
+    zcat ${source}Colors.gz \
+    | subColumn -miss=/dev/null 1 stdin treeToPublic stdout \
+    | grep -v EPI_ISL > public${Source}Colors
+    wc -l public${Source}Colors
+    gzip -f public${Source}Colors
+done