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/updateIdMapping.sh src/hg/utils/otto/sarscov2phylo/updateIdMapping.sh new file mode 100755 index 0000000..be70d07 --- /dev/null +++ src/hg/utils/otto/sarscov2phylo/updateIdMapping.sh @@ -0,0 +1,66 @@ +#!/bin/bash +set -beEu -x -o pipefail + +# Do not modify this script, modify the source tree copy: +# kent/src/hg/utils/otto/sarscov2phylo/updateIdMapping.sh + +usage() { + echo "usage: $0 metadata_date.tsv.gz sequences_date.fasta.gz" +} + +if [ $# != 2 ]; then + usage + exit 1 +fi + +nextmeta=$1 +nextfasta=$2 + +scriptDir=$(dirname "${BASH_SOURCE[0]}") +source $scriptDir/util.sh + +today=$(date +%F) +ottoDir=/hive/data/outside/otto/sarscov2phylo +mapScriptDir=~/chris_ncov +# Should use a better location than this... +installDir=/hive/users/angie/gisaid + +ncbiDir=$ottoDir/ncbi.$today +mkdir -p $ncbiDir +cd $ncbiDir +$scriptDir/getNcbi.sh >& getNcbi.log + +cogUkDir=$ottoDir/cogUk.$today +mkdir -p $cogUkDir +cd $cogUkDir +$scriptDir/getCogUk.sh >& getCogUk.log + +# Last time I checked, CNCB had not updated since September, just keep using what we have +cncbDir=$ottoDir/cncb.latest + +# Set up input files for Chris's scripts to map GISAID <--> public sequences +cd $mapScriptDir +mkdir input/$today +cd input/$today +ln -sf $cncbDir/cncb.nonGenBank.fasta . +ln -sf $ncbiDir/genbank.fa.xz . +ln -sf $cogUkDir/cog_all.fasta.xz . +ln -sf $nextfasta . +xcat $nextmeta | tail -n+2 | cut -f1,3 > seqToEpi + +cd $mapScriptDir +./build.sh -t $today + +cd $installDir + +gbToDate=$ncbiDir/gbToDate +cogUkToDate=$cogUkDir/cogUkToDate +cncbToDate=$cncbDir/cncbToDate + +join -t$'\t' -a 1 -1 2 -o 1.1,1.2,1.3,2.2 \ + <(sort -k2,2 ~/chris_ncov/epiToPublicIdName.$today.txt) \ + <(sort $gbToDate $cncbToDate $cogUkToDate) \ +| sort \ + > epiToPublicAndDate.$today + +ln -sf epiToPublicAndDate.$today epiToPublicAndDate.latest