# install required utilities
brew install mafft seqkit brewsci/bio/snp-dists xmlstarlet gawk

# download 335 aligned sequences of SARS-like viruses from GenBank
curl -Lso sarslike.fa 'https://drive.google.com/uc?export=download&id=1j-YFiMYG4DkVKSget2fYW-gaJDy6NCkW'

# make a TSV matrix for number of nucleotide changes without counting insertions or deletions
snp-dists sarslike.fa>sarslike.dist

# download protein sequences and XML metadata (downloading the XML can take a few minutes)
curl 'https://eutils.ncbi.nlm.nih.gov/entrez/eutils/efetch.fcgi?db=nuccore&rettype=fasta_cds_aa&id='$(seqkit seq -ni sarslike.fa|paste -sd, -)>sarslike.aa
curl 'https://eutils.ncbi.nlm.nih.gov/entrez/eutils/efetch.fcgi?db=nuccore&retmode=xml&id='$(seqkit seq -ni sarslike.fa|paste -sd, -)>sarslike.xml

# convert metadata from XML to TSV
xml fo -D sarslike.xml|xml sel -t -m //GBSeq -v GBSeq_accession-version -o $'\t' -v GBSeq_definition -o $'\t' -v GBSeq_create-date -o $'\t' -v './/GBQualifier[GBQualifier_name="collection_date"]/GBQualifier_value' -o $'\t' -v '(.//GBAuthor)[1]' -o ... -v '(.//GBAuthor)[last()]' -o $'\t' -v '(.//GBReference_title[text()!="Direct Submission"])[last()]' -o $'\n'>sarslike.tsv

# `tab \\t` is like `column -ts$'\t'` but it doesn't get thrown off by empty fields
tab(){ awk '{if(NF>m)m=NF;for(i=1;i<=NF;i++){a[NR][i]=$i;l=length($i);if(l>b[i])b[i]=l}}END{for(h in a){for(i=1;i<=m;i++)printf(i==m?"%s\n":"%-"(b[i]+n)"s",a[h][i])}}' "${1+FS=$1}" "n=${2-1}";}

# download RmYN02 spike sequence from GenBank
curl 'https://eutils.ncbi.nlm.nih.gov/entrez/eutils/efetch.fcgi?db=nuccore&rettype=fasta_cds_aa&id=MW201982.1'>RmYN02.spike.fa

# align spike proteins (`--reorder` places similar sequences next to each other)
seqkit grep -nrip 'spike|surface|S protein' sarslike.aa|cat - RmYN02.spike.fa|mafft --reorder ->spike2.aln

# print number of letter changes from first sequence to every other sequence in an aligned set of sequences
dif1(){ awk 'NR==1{split($0,a,"");l=length;next}{split($0,b,"");n=0;for(i=1;i<=l;i++)if(a[i]!=b[i])n++;print n}' "$@";}

# colorize amino acids
cola(){ sed -f <(awk '{print"s/"$1"/\033[38;2;0;0;0m\033[48;2;"$2";"$3";"$4"m"$1"\033[0m/g"}' <(printf %s\\n 'A 242 121 121' 'C 242 182 121' 'D 242 242 121' 'E 182 242 121' 'F 121 242 151' 'G 121 242 222' 'H 121 182 242' 'I 121 121 242' 'K 182 121 242' 'L 242 121 242' 'M 255 191 191' 'N 255 223 191' 'P 255 255 191' 'Q 223 255 191' 'R 191 255 207' 'S 191 255 244' 'T 191 223 255' 'V 191 191 255' 'W 223 191 255' 'Y 255 191 255' '- 140 140 140')) "$@";}

# print region of spike protein around furin cleavage site, number of letter changes in the region from to Wuhan-Hu-1, number of nucleotide changes from Wuhan-Hu-1 in whole genome, and metadata about samples
x=NC_045512.2;sub=$(seqkit subseq -r708:740 spike2.aln);paste <(seqkit seq -ni spike2.aln|sed 's/^lcl|//;s/_prot_.*//') <(dif1 <(seqkit grep -rp "$x"<<<"$sub"|seqkit seq -s;seqkit seq -s<<<"$sub"))|awk -F\\t 'NR==FNR{a[$1]=$2;next}{print a[$1]"\t"$0}' <(awk -F\\t 'NR==1{for(i=2;i<=NF;i++)if($i==x)break;next}{print$1 FS$i}' x=$x sarslike.dist) -|awk 'NR==FNR{a[$1]=$2 FS$3 FS$4 FS$1 FS$5": "$6;next}{print$3,$1,a[$2]}' {,O}FS=\\t <(cat sarslike.tsv <(printf 'MW201982.1\tRmYN02\n')) -|sed 's/, complete genome//'|tab \\t|cut -c-168|paste -d' ' <(seqkit seq -s<<<"$sub"|cola) -