3760b502d693a5e3e1540f5eb190c87a6a0e965f angie Tue Feb 22 11:04:44 2022 -0800 Strip quotes from nextclade tsv output, speed up compression by using more threads & speed up sed by usin conditional. diff --git src/hg/utils/otto/sarscov2phylo/gisaidFromChunks.sh src/hg/utils/otto/sarscov2phylo/gisaidFromChunks.sh index 235aa79..418f74a 100755 --- src/hg/utils/otto/sarscov2phylo/gisaidFromChunks.sh +++ src/hg/utils/otto/sarscov2phylo/gisaidFromChunks.sh @@ -15,44 +15,45 @@ today=$(date +%F) # Run pangolin and nextclade on any chunks that need it cd /hive/users/angie/gisaid/chunks make nextclade.tsv make -j10 cd /hive/users/angie/gisaid # Glom all the chunks together. # Remove initial "hCoV-19/" and remove spaces a la nextmeta (e.g. "Hong Kong" -> "HongKong"). # Strip single quotes (e.g. "Cote d'Ivoire" --> "CotedIvoire"). # Also remove a stray comma in a name that caused Newick parsing error ("Hungary/US-32533w,/2020"). # Keep the strain|epiId|date "full names". time xzcat chunks/gisaid_epi_isl_*.fa.xz \ | sed -re 's@^>hCo[Vv]-19/+@>@; s/[ '"'"',()]//g; s/\r$//;' \ -| xz -T 8 \ +| xz -T 20 \ > gisaid_fullNames_$today.fa.xz # Make tmp files with a fullName key and various columns that we'll join together. fastaNames gisaid_fullNames_$today.fa.xz \ | awk -F\| -vOFS="\t" '{print $0, $1, $2, $3;}' \ | sort -u \ > tmp.first3 # Sequence length faSize -detailed <(xzcat gisaid_fullNames_$today.fa.xz) | sort -u > tmp.lengths # Lineage & clade assignments sort -u chunks/pangolin.tsv \ > tmp.lineage -sort -u chunks/nextclade.tsv \ +sed -re 's/"//g;' chunks/nextclade.tsv \ +| sort -u \ > tmp.clade # Countries -- go back to unstripped sequence names: xzcat chunks/gisaid_epi_isl_*.fa.xz \ | grep ^\> \ | sed -re 's@^>hCo[Vv]-19/+@@;' \ | $scriptDir/gisaidNameToCountry.pl \ | sort -u \ > tmp.country # Join locally computed fields and sort by EPI ID for joining with latest real nextmeta join -t$'\t' -a 1 tmp.first3 tmp.lengths \ | join -t$'\t' -a 1 -o 1.1,1.2,1.3,1.4,1.5,2.2 - tmp.clade \ | join -t$'\t' -a 1 -o 1.1,1.2,1.3,1.4,1.5,1.6,2.2 - tmp.lineage \ | join -t$'\t' -a 1 -o 1.1,1.2,1.3,1.4,1.5,1.6,1.7,2.2 - tmp.country \ | tawk '{print $3, $2, $4, $5, $6, $7, $8;}' \ @@ -70,25 +71,25 @@ #19 pangolin_lineage # Fill in other columns from nextmeta when possible (join on EPI ID since names change over time) set +o pipefail zcat $lastRealNextmeta | head -1 \ > metadata_batch_$today.tsv set -o pipefail zcat $lastRealNextmeta \ | tail -n+2 \ | sort -k 3,3 \ | join -t$'\t' -a 1 -2 3 \ -o 1.2,2.2,1.1,2.4,1.3,2.6,1.7,2.8,2.9,2.10,2.11,2.12,2.13,1.4,2.15,2.16,2.17,1.5,1.6,2.20,2.21,2.22,2.23,2.24,2.25,2.26,2.27 \ tmp.epiToLocalMeta - \ | sort \ >> metadata_batch_$today.tsv wc -l metadata_batch_$today.tsv -gzip -f metadata_batch_$today.tsv +pigz -p 8 -f metadata_batch_$today.tsv # Make fasta with strain-name headers a la nextfasta. xzcat gisaid_fullNames_$today.fa.xz \ -| sed -re 's/\|.*//' \ -| xz -T 8 \ +| sed -re '/^>/ s/\|.*//' \ +| xz -T 20 \ > sequences_batch_$today.fa.xz # Clean up rm tmp.*