be4311c07e14feb728abc6425ee606ffaa611a58
markd
  Fri Jan 22 06:46:58 2021 -0800
merge with master

diff --git src/hg/utils/otto/sarscov2phylo/extractUnmappedPublic.sh src/hg/utils/otto/sarscov2phylo/extractUnmappedPublic.sh
index e7a49fb..bace556 100755
--- src/hg/utils/otto/sarscov2phylo/extractUnmappedPublic.sh
+++ src/hg/utils/otto/sarscov2phylo/extractUnmappedPublic.sh
@@ -1,95 +1,95 @@
 #!/bin/bash
 set -beEu -x -o pipefail
 
 #	Do not modify this script, modify the source tree copy:
 #	kent/src/hg/utils/otto/sarscov2phylo/extractUnmappedPublic.sh
 
 usage() {
     echo "usage: $0 epiToPublicNameAndDate genbank.fa cogUk.fa cncb.fa"
 }
 
 if [ $# != 4 ]; then
   usage
   exit 1
 fi
 
 epiToPublic=$1
 genbankFa=$2
 cogUkFa=$3
 cncbFa=$4
 
 minSize=29000
 
 echo "epiToPublic=$epiToPublic"
 echo "genbank=$genbankFa"
 echo "cogUk=$cogUkFa"
 echo "cncb=$cncbFa"
 
 ottoDir=/hive/data/outside/otto/sarscov2phylo
 gbToDate=$ottoDir/ncbi.latest/gbToDate
 cogUkToDate=$ottoDir/cogUk.latest/cogUkToDate
 cncbToDate=$ottoDir/cncb.latest/cncbToDate
 
 # Outputs for next script (addUnmappedPublic.sh):
 alignedFa=unmapped.aligned.fasta
 renaming=unmapped.renaming
 
 scriptDir=$(dirname "${BASH_SOURCE[0]}")
 
 source $scriptDir/util.sh
 
 # faFilter discards fasta header info after the first word, so before running it,
 # make files that associate IDs and strain names.
 fastaNames $genbankFa | cleanGenbank | sort > gbAccName
 # COG-UK has only names, so no need to associate anything.
 # CNCB puts name first, then ID, so swap the order:
 fastaNames $cncbFa | cleanCncb \
 | tawk '{print $2, $1;}' | sort > cncbAccName
 
 # Filter minSize and exclude sequences that were mapped to GISAID IDs
-join -t$'\t' tree.renaming $epiToPublic | cut -f 3 > mappedIds
+join -t$'\t' epi.renaming $epiToPublic | cut -f 3 > mappedIds
 xcat $genbankFa \
 | faFilter -minSize=$minSize stdin stdout \
 | faSomeRecords -exclude stdin mappedIds genbank.unmapped.fa
 fastaSeqCount genbank.unmapped.fa
 #*** TODO: also exclude COG-UK sequences that are in GenBank (some with incomplete names)
 xcat $cogUkFa \
 faFilter -minSize=$minSize stdin stdout \
 | faSomeRecords -exclude stdin mappedIds cogUk.unmapped.fa
 fastaSeqCount cogUk.unmapped.fa
 #*** TODO: also exclude CNCB sequences that are in GenBank (some with incomplete names)
 # Tweak CNCB's fasta headers around to use acc not name:
 xcat $cncbFa \
 | sed -re 's/^>.*\| */>/;' \
 | faFilter -minSize=$minSize stdin stdout \
 | faSomeRecords -exclude stdin mappedIds cncb.unmapped.fa
 fastaSeqCount cncb.unmapped.fa
 
 cat genbank.unmapped.fa cogUk.unmapped.fa cncb.unmapped.fa > unmapped.fa
 fastaSeqCount unmapped.fa
 
 # Make a file for renaming sequences from just the ID to isolateName|ID|date
 join -t$'\t' <(fastaNames genbank.unmapped.fa | sort) gbAccName > gbAccNameUnmapped
 join -t$'\t' <(fastaNames cncb.unmapped.fa | sort) cncbAccName > cncbAccNameUnmapped
 join -t$'\t' <(sort cncbAccNameUnmapped gbAccNameUnmapped) <(sort $cncbToDate $gbToDate) \
 | grep -vF NC_045512.2 \
 | tawk '{ if ($2 == "") { $2 = "?";} if ($3 == "") { $3 = "?"; } print $1, $2 "|" $1 "|" $3; }' \
 | sed -re 's/, /_/g; s/[, ]/_/g;' \
     > $renaming
 join -t$'\t' <(fastaNames cogUk.unmapped.fa | sort) <(sort $cogUkToDate) \
 | tawk '{print $1, $1 "|" $2;}' \
     >> $renaming
 
 #*** TODO: extract submitter credits from metadata, make acknowledgements file
 
 # Use Rob's script that aligns each sequence to NC_045512.2 and concatenates the results
 # as a multifasta alignment (from which we can extract VCF with SNVs):
 #conda install -c bioconda mafft
 rm -f $alignedFa
 export TMPDIR=/dev/shm
 time bash ~angie/github/sarscov2phylo/scripts/global_profile_alignment.sh \
   -i unmapped.fa \
   -o $alignedFa \
   -t 50
 faSize $alignedFa