A Brazilian wastewater sample from November 2019 with SARS2 reads was probably contaminated - sars2.net



The abstract of a Brazilian paper from 2021 said: [https://www.sciencedirect.com/science/article/pii/S0048969721012651]

Human sewage from Florianopolis (Santa Catarina, Brazil) was analyzed for severe acute respiratory syndrome coronavirus-2 (SARS-CoV2) from October 2019 until March 2020. Twenty five ml of sewage samples were clarified and viruses concentrated using a glycine buffer method coupled with polyethylene glycol precipitation, and viral RNA extracted using a commercial kit. SARS-CoV-2 RNA was detected by RT-qPCR using oligonucleotides targeting N1, S and two RdRp regions. The results of all positive samples were further confirmed by a different RT-qPCR system in an independent laboratory. S and RdRp amplicons were sequenced to confirm identity with SARS-CoV-2. Genome sequencing was performed using two strategies; a sequence-independent single-primer amplification (SISPA) approach, and by direct metagenomics using Illumina's NGS. SARS-CoV-2 RNA was detected on 27th November 2019 (5.49 ± 0.02 log10 SARS-CoV-2 genome copies (GC) L−1), detection being confirmed by an independent laboratory and genome sequencing analysis. The samples in the subsequent three events were positive by all RT-qPCR assays; these positive results were also confirmed by an independent laboratory. The average load was 5.83 ± 0.12 log10 SARS-CoV-2 GC L−1, ranging from 5.49 ± 0.02 log10 GC L−1 (27th November 2019) to 6.68 ± 0.02 log10 GC L−1 (4th March 2020). Our findings demonstrate that SARS-CoV-2 was likely circulating undetected in the community in Brazil since November 2019, earlier than the first reported case in the Americas (21st January 2020).

Guy Gadboit downloaded the raw reads from the paper and wrote a Twitter thread about them. [https://x.com/gadboit/status/1793281441322693002]

Coverage and pileup

Gadboit said that he only got 7% coverage for SARS2 even though the paper said that they got 25% coverage. However I got about 15% coverage:

$ enad()(printf %s\\n "${@-$(cat)}"|while IFS= read x;do curl -s "https://www.ebi.ac.uk/ena/portal/api/filereport?accession=$x&result=read_run&fields=fastq_ftp"|sed 1d|cut -f2|tr \; \\n;done|sed s,^,ftp://,|xargs wget -q)
$ x=SRR13152144;enad $x
$ fastp -35 -i $x\_1.fastq.gz -I $x\_2.fastq.gz -o $x\_1.fq.gz -O $x\_2.fq.gz
$ minimap2 -a --sam-hit-only sars2.fa $x\_[12].fq.gz|samtools sort ->$x.bam
$ samtools coverage $x.bam|column -t
#rname      startpos  endpos  numreads  covbases  coverage  meandepth  meanbaseq  meanmapq
MN908947.3  1         29903   92        4348      14.5403   0.293649   36.7       43.5

Gadboit said that he didn't see any other lineage-defining locations except D614G (A23403G) and C8782. However I actually didn't get any read that covered the position of D614G regardless of whether I trimmed reads or not. I got two reads for 8782 and two reads for 28144, but they both had the lineage B allele. But I didn't get any reads that covered the position of the third mutation of WA1 (18060) or the A0 mutation (29095):

$ samtools mpileup -f sars2.fa SRR13152144.bam|ruby -alne'next unless[8782,18060,28144,29095,23403].include?($F[1].to_i);next if$F[3]==0||$F[4]=="*";s=$F[4].upcase.gsub(/[.,]/,$F[2]).gsub(/\^./,"");s2=s.dup;s.enum_for(:scan,/[+-]\d+/).each{m=Regexp.last_match;m.begin(0).upto(m.end(0)+s[m.begin(0)+1,m.end(0)].to_i){|i|s2[i]=" "}};puts $F[1]+" "+$F[2]+" "+"ACGT".chars.map{|x|q=$F[5].chars.values_at(*s2.chars.each_index.select{|i|s2[i]==x}).map{|c|c.ord-33};q.size.to_s+(q.size==0?"":"["+"%.0f"%(q.sum*1.0/q.size)+"]")}*" "'|sed s/^/$x\ /|sort -n|(echo pos ref A C G T;cat)|column -t
pos    ref  A  C      G  T
8782   C    0  2[34]  0  0
28144  T    0  0      0  2[38]

I got 4 different mutations with MAPQ over 10. The MAPQ of each mutation was around 50, because there were two high-quality reads that supported each mutation:

$ x=SRR13152144;samtools index $x.bam;bcftools mpileup -f sars2.fa $x.bam|bcftools call -mv|grep -v ^\#\#|cut -f-6|column -t
MN908947.3  3392   .   G    C    3.22451
MN908947.3  5706   .   A    G    48.4146
MN908947.3  11437  .   TAA  TA   5.7555
MN908947.3  14408  .   C    T    48.4146
MN908947.3  14422  .   C    T    46.4146
MN908947.3  14895  .   T    G    47.4146

I don't understand why G3392C gets a MAPQ of only 3.2, even though it's supported by two reads with high quality like the other mutations. The quality of both bases is H = 72-33 = 39:

$ samtools mpileup -f sars2.fa $x.bam|awk '$2~/^(3392|5706|14408|14422|14895)$/'|cut -f2-|column -t
3392   G  2  Cc  HH
5706   A  2  Gg  HH
14408  C  2  Tt  HH
14422  C  2  Tt  HF
14895  T  2  Gg  HG

Mutations at GISAID

In a nearly complete set of GISAID submissions with a collection date in 2020, I didn't find any sample with all 4 mutations or even 3 out of 4 mutations, but the two earliest samples with 2 out of the 4 mutations were both from Santa Catarina in Brazil (which was the same state where the wastewater samples were collected):

$ curl -Ls sars2.net/f/gisaid2020.tsv.xz|xz -dc>gisaid2020.tsv
$ printf %s\\n A5706G C14408T C14422T T14895G|awk -F\\t 'NR==FNR{a[$0];next}{n=0;split($12,z,",");for(i in z)if(z[i] in a)n++;if(n>1)print n"\t"$0}' - gisaid2020.tsv|sort -k5|head|cut -f1-8,11-14|tr \\t \|
2|EPI_ISL_541370|hCoV-19/Brazil/SC-FIOCRUZ-770/2020|2020-09-22|2020-03-12|B.1.1.33|Brazil|Santa Catarina|Human|10|C241T,C3037T,A5706G,C14408T,A23403G,T27299C,G28881A,G28882A,G28883C,T29148C||
2|EPI_ISL_476278|hCoV-19/Brazil/SC-L16-CD314/2020|2020-06-25|2020-03-24|B.1.1.33|Brazil|Santa Catarina|Human|11|C241T,C3037T,A5706G,C14408T,G18803T,A23403G,T27299C,G28881A,G28882A,G28883C,T29148C||
2|EPI_ISL_541116|hCoV-19/Spain/PV-IBV-98006728/2020|2020-09-22|2020-07-16|B.1|Spain|Basque Country|Human|13|C241T,C1392T,C1593T,T1623C,C3037T,C5512T,A5706G,G10870T,C12815T,C14408T,A20268G,A23403G,A27831G||

The Brazilian samples from Santa Catarina belonged to the clade B.1.1.33. They had D614G and 3 out of the 4 other B.1 mutations (with G25563T missing).

This shows samples whose set of mutation is a subset of the Brazilian sample with the earliest collection date, so that it shows different possible ancestor strains of the Brazilian sample. The first column shows how many times each mutation set was found in each country, and the second column shows how many total mutations the mutation set has relative to Wuhan-Hu-1. The samples with 9 shared mutations are mostly from Brazil, even though there's also samples from USA, Chile, Canada, Portugal, Australia, Uruguay, and UK. However 7 out of 10 mutations are shared by B.1.1 samples from many countries all over the world:

$ printf %s\\n C241T C3037T A5706G C14408T A23403G T27299C G28881A G28882A G28883C T29148C|awk -F\\t 'NR==FNR{a[$0];next}{split($12,b,",");for(i in b)if(!(b[i]in a))next;o[$11 FS$6 FS$5 FS$12]++}END{for(i in o)print o[i]FS i}' - gisaid2020.tsv|sort -rnk2|head -n50|column -ts$'\t'
1    10  Brazil                B.1.1.33    C241T,C3037T,A5706G,C14408T,A23403G,T27299C,G28881A,G28882A,G28883C,T29148C
9    9   Chile                 B.1.1.33    C241T,C3037T,C14408T,A23403G,T27299C,G28881A,G28882A,G28883C,T29148C
7    9   Canada                B.1.1.33    C241T,C3037T,C14408T,A23403G,T27299C,G28881A,G28882A,G28883C,T29148C
5    9   Portugal              B.1.1.33    C241T,C3037T,C14408T,A23403G,T27299C,G28881A,G28882A,G28883C,T29148C
2    9   Australia             B.1.1.33    C241T,C3037T,C14408T,A23403G,T27299C,G28881A,G28882A,G28883C,T29148C
159  9   Brazil                B.1.1.33    C241T,C3037T,C14408T,A23403G,T27299C,G28881A,G28882A,G28883C,T29148C
11   9   USA                   B.1.1.33    C241T,C3037T,C14408T,A23403G,T27299C,G28881A,G28882A,G28883C,T29148C
1    9   Uruguay               B.1.1.33    C241T,C3037T,C14408T,A23403G,T27299C,G28881A,G28882A,G28883C,T29148C
1    9   United Kingdom        B.1.1.33    C241T,C3037T,C14408T,A23403G,T27299C,G28881A,G28882A,G28883C,T29148C
1    9   Brazil                Unassigned  C241T,C3037T,C14408T,A23403G,T27299C,G28881A,G28882A,G28883C,T29148C
6    8   Switzerland           B.1.1       C241T,C3037T,C14408T,A23403G,G28881A,G28882A,G28883C,T29148C
4    8   Brazil                B.1.1.33    C241T,C14408T,A23403G,T27299C,G28881A,G28882A,G28883C,T29148C
3    8   Netherlands           B.1.1       C241T,C3037T,C14408T,A23403G,G28881A,G28882A,G28883C,T29148C
3    8   Australia             B.1.1.33    C3037T,C14408T,A23403G,T27299C,G28881A,G28882A,G28883C,T29148C
2    8   United Kingdom        B.1.1.33    C241T,C3037T,C14408T,A23403G,G28881A,G28882A,G28883C,T29148C
2    8   United Kingdom        B.1.1       C241T,C3037T,C14408T,A23403G,G28881A,G28882A,G28883C,T29148C
2    8   Czech Republic        B.1.1       C241T,C3037T,C14408T,A23403G,G28881A,G28882A,G28883C,T29148C
2    8   Brazil                B.1.1       C241T,C3037T,C14408T,A23403G,T27299C,G28881A,G28882A,G28883C
19   8   Brazil                B.1.1.33    C3037T,C14408T,A23403G,T27299C,G28881A,G28882A,G28883C,T29148C
1    8   USA                   B.1.1.33    C241T,C3037T,C14408T,A23403G,G28881A,G28882A,G28883C,T29148C
1    8   USA                   B.1.1       C241T,C3037T,C14408T,A23403G,T27299C,G28881A,G28882A,G28883C
1    8   Italy                 B.1.1       C241T,C3037T,C14408T,A23403G,G28881A,G28882A,G28883C,T29148C
1    8   Brazil                Unassigned  C241T,C3037T,C14408T,T27299C,G28881A,G28882A,G28883C,T29148C
1    8   Brazil                B.1.1.33    C241T,C3037T,C14408T,A23403G,T27299C,G28881A,G28882A,G28883C
1    8   Brazil                B.1.1       C241T,C3037T,C14408T,A23403G,G28881A,G28882A,G28883C,T29148C
1    8   Australia             B.1.1.33    C241T,C3037T,C14408T,A23403G,T27299C,G28881A,G28882A,G28883C
97   7   Russia                B.1.1       C241T,C3037T,C14408T,A23403G,G28881A,G28882A,G28883C
92   7   Italy                 B.1.1       C241T,C3037T,C14408T,A23403G,G28881A,G28882A,G28883C
91   7   Spain                 B.1.1       C241T,C3037T,C14408T,A23403G,G28881A,G28882A,G28883C
9    7   Turkey                B.1.1       C241T,C3037T,C14408T,A23403G,G28881A,G28882A,G28883C
9    7   Guam                  B.1.1       C241T,C3037T,C14408T,A23403G,G28881A,G28882A,G28883C
9    7   Chile                 B.1.1       C241T,C3037T,C14408T,A23403G,G28881A,G28882A,G28883C
8    7   United Arab Emirates  B.1.1       C241T,C3037T,C14408T,A23403G,G28881A,G28882A,G28883C
8    7   Poland                B.1.1       C241T,C3037T,C14408T,A23403G,G28881A,G28882A,G28883C
8    7   Finland               B.1.1       C241T,C3037T,C14408T,A23403G,G28881A,G28882A,G28883C
76   7   Germany               B.1.1       C241T,C3037T,C14408T,A23403G,G28881A,G28882A,G28883C
72   7   Australia             B.1.1       C241T,C3037T,C14408T,A23403G,G28881A,G28882A,G28883C
7    7   Madagascar            B.1.1       C241T,C3037T,C14408T,A23403G,G28881A,G28882A,G28883C
6    7   Peru                  B.1.1       C241T,C3037T,C14408T,A23403G,G28881A,G28882A,G28883C
6    7   Mexico                B.1.1       C241T,C3037T,C14408T,A23403G,G28881A,G28882A,G28883C
6    7   Colombia              B.1.1       C241T,C3037T,C14408T,A23403G,G28881A,G28882A,G28883C
6    7   Bangladesh            B.1.1       C241T,C3037T,C14408T,A23403G,G28881A,G28882A,G28883C
56   7   Sweden                B.1.1       C241T,C3037T,C14408T,A23403G,G28881A,G28882A,G28883C
53   7   France                B.1.1       C241T,C3037T,C14408T,A23403G,G28881A,G28882A,G28883C
52   7   Canada                B.1.1       C241T,C3037T,C14408T,A23403G,G28881A,G28882A,G28883C
50   7   Cyprus                B.1.1       C241T,C3037T,C14408T,A23403G,G28881A,G28882A,G28883C
5    7   India                 B.1.1       C241T,C3037T,C14408T,A23403G,G28881A,G28882A,G28883C
5    7   Hungary               B.1.1       C241T,C3037T,C14408T,A23403G,G28881A,G28882A,G28883C
5    7   Ecuador               B.1.1       C241T,C3037T,C14408T,A23403G,G28881A,G28882A,G28883C
4    7   South Korea           B.1.1       C241T,C3037T,C14408T,A23403G,G28881A,G28882A,G28883C

Out of the 5 nucleotide changes I got when I ran pileup on the November 2019 sample, only C14408T was part of the 7 mutations in B.1.1. The second mutation which matched the Brazilian samples on GISAID was A5706G, but the oldest two samples with the mutation are from the state of Santa Catarina in Brazil: [https://cov-spectrum.org/explore/World/AllSamples/from%3D2020-01-01%26to%3D2020-03-31/variants?nucMutations=A5706G]

$ grep A5706G gisaid2020.tsv|sort -t$'\t' -k4|head|cut -f1-8,11-14|tr \\t \|
EPI_ISL_541370|hCoV-19/Brazil/SC-FIOCRUZ-770/2020|2020-09-22|2020-03-12|B.1.1.33|Brazil|Santa Catarina|Braco do Norte|10|C241T,C3037T,A5706G,C14408T,A23403G,T27299C,G28881A,G28882A,G28883C,T29148C||
EPI_ISL_476278|hCoV-19/Brazil/SC-L16-CD314/2020|2020-06-25|2020-03-24|B.1.1.33|Brazil|Santa Catarina||11|C241T,C3037T,A5706G,C14408T,G18803T,A23403G,T27299C,G28881A,G28882A,G28883C,T29148C||
EPI_ISL_525869|hCoV-19/USA/OR-OHSU-0648/2020|2020-09-01|2020-06-16|B.1.134|USA|Oregon|Multnomah County|9|C241T,C3037T,A5706G,C8818T,A12579G,C14408T,C15952T,G20980A,A23403G||
EPI_ISL_525840|hCoV-19/USA/OR-OHSU-0546/2020|2020-09-01|2020-06-16|B.1.134|USA|Oregon|Washington County|9|C241T,C3037T,A5706G,C8818T,A12579G,C14408T,C15952T,G20980A,A23403G||
EPI_ISL_525870|hCoV-19/USA/OR-OHSU-0649/2020|2020-09-01|2020-06-16|B.1.134|USA|Oregon|Washington County|9|C241T,C3037T,A5706G,C8818T,A12579G,C14408T,C15952T,G20980A,A23403G||
EPI_ISL_541116|hCoV-19/Spain/PV-IBV-98006728/2020|2020-09-22|2020-07-16|B.1|Spain|Basque Country|San Sebastian|13|C241T,C1392T,C1593T,T1623C,C3037T,C5512T,A5706G,G10870T,C12815T,C14408T,A20268G,A23403G,A27831G||
EPI_ISL_535485|hCoV-19/SouthAfrica/KRISP-K002622/2020|2020-09-15|2020-07-31|B.1.1.412|South Africa|KwaZulu-Natal||12|C241T,C3037T,C3165A,A5706G,G11447A,C14408T,C22675T,A23403G,G28077T,G28881A,G28882A,G28883C||

March 2020 wastewater sample

The Brazilian study also included a second SRA run, which was labeled "Sewage metratranscriptome from Florianopolis, SC, Brazil (Mar20)" (instead of "Nov19" like the other run). [https://trace.ncbi.nlm.nih.gov/Traces?run=SRR13152143] bcftools call gave only 3 mutations for the March 2020 sample but none of them had a MAPQ above 10:

$ enad()(printf %s\\n "${@-$(cat)}"|while IFS= read x;do curl -s "https://www.ebi.ac.uk/ena/portal/api/filereport?accession=$x&result=read_run&fields=fastq_ftp"|sed 1d|cut -f2|tr \; \\n;done|sed s,^,ftp://,|xargs wget -q)
$ x=SRR13152143;enad $x
$ fastp -35 -i $x\_1.fastq.gz -I $x\_2.fastq.gz -o $x\_1.fq.gz -O $x\_2.fq.gz
$ minimap2 -a --sam-hit-only sars2.fa $x\_[12].fastq.gz|samtools sort ->$x.bam
$ samtools index $x.bam;bcftools mpileup -f sars2.fa $x.bam|bcftools call -mv|grep -v ^\#\#|cut -f2,4-6|column -t
6192  C    T    3.7766
8145  A    G    9.88514
8206  A    G    9.88514

However all three mutations were rare in the set of GISAID submissions with a collection date in 2020, and none of them were found in samples from Brazil, and none of the submissions had even 2 out of the 3 mutations.

The March 2020 sample only got about 5% coverage for SARS2, and it didn't have any reads which covered the positions of lineage A (8782 and 28144), WA1 (18060), lineage A0 (29095), or D614G (23403):

$ samtools coverage SRR13152143.bam|column -t
#rname      startpos  endpos  numreads  covbases  coverage  meandepth  meanbaseq  meanmapq
MN908947.3  1         29903   23        1636      5.47102   0.0934689  37.4       57.8
$ samtools depth -a SRR13152143.bam|awk '$2~/^(8782|18060|28144|29095|23403)$/'|cut -f2-
8782    0
18060   0
23403   0
28144   0
29095   0

If I show all mutations instead of using bcftools call, the March 2020 sample only gets one additional mutation (which is found in 1 out of 2 reads and which has quality 18, so it's probably a sequencing error):

$ x=SRR13152143;samtools mpileup -f sars2.fa $x.bam|puptab|awk '$2!=$3&&$4!=0'|(echo pos ref allele count coverage allele% average_quality;cat)|column -t
pos   ref  allele  count  coverage  allele%  average_quality
6137  A    C       1      2         50.000   18.0
6192  C    T       1      1         100.000  31.0
8145  A    G       1      1         100.000  39.0
8206  A    G       1      1         100.000  39.0

In the set of GISAID submissions with a collection date in 2020, there's a total of 183 submissions where the state is listed as Santa Catarina. They include mutations in a total of 618 different positions, but none of hem are the 4 positions shown in the output above.

The March 2020 sample has zero depth for all positions where there's mutations in the two A5706G samples from Santa Catarina:

$ grep A5706G gisaid2020.tsv|grep Santa\ Catarina|cut -f12|tr , \\n|sed 's/.\(.*\)./\1/'|awk 'NR==FNR{a[$0];next}$2 in a' - <(samtools depth -a SRR13152143.bam)|cut -f2-|column -t
241    0
3037   0
5706   0
14408  0
18803  0
23403  0
27299  0
28881  0
28882  0
28883  0
29148  0

STAT results

In the STAT results of the November 2019 sample, the most abundant leaf nodes were mostly human reads, human reads misclassified under monkeys, or bacteria:

$ curl "https://trace.ncbi.nlm.nih.gov/Traces/sra-db-be/run_taxonomy?acc=SRR13152144&cluster_name=public">SRR13152144.stat
$ stal()(jq -r '[.[]|.tax_table[]|.parent]as$par|[.[]|.tax_table[]|select(.tax_id as$x|$par|index($x)|not)]|sort_by(-.total_count)[]|((.total_count|tostring)+";"+.org)' "$@") # STAT leaf nodes
$ stal SRR13152144.stat|head -n24
29976;Homo sapiens
427;Escherichia coli
257;Gorilla gorilla gorilla
241;Klebsiella pneumoniae
185;Clostridioides difficile
171;Dolosigranulum pigrum (bacteria)
124;Pan troglodytes
115;Pongo abelii
97;unclassified Streptococcus
57;Salmonella enterica
53;Aotus nancymaae (monkey)
53;Pseudomonas aeruginosa (bacteria)
52;Bordetella bronchiseptica E013 (bacteria)
44;Nomascus leucogenys (monkey)
39;Streptococcus pneumoniae
35;Staphylococcus aureus
26;Chlorocebus sabaeus (monkey)
20;Staphylococcus epidermidis (bacteria)
15;Veillonella sp. oral taxon 158 str. F0412
14;Severe acute respiratory syndrome coronavirus 2
8;Oncorhynchus kisutch

babarlelephant had found A5706G earlier

After I had written the previous sections in this file, I found that babarlelephant had posted these tweets in 2021: [<https://x.com/babarlelephant/status/1397204321838813199>]

Half of the time when I think I have made some new discovery, it has already been found by Babar or Daoyu or breakfast_dogs.

Other two mutations in November 2020 sample

The other 2 out of the 4 mutations with MAPQ above 10 don't seem to be very informative. The C14422T mutation seems to have emerged independently in a lineage A sample in China, in B.1 samples in USA, in D.2 samples in Australia, in B.1.237 samples in South Africa, and so on:

$ grep C14422T gisaid2020.tsv |awk -F\\t '!a[$12]++'|head|cut -f1,3-8,11-14|tr \\t \|
EPI_ISL_425457|2020-04-14|2020-04-01|B.29|United Kingdom|England||7|G11083T,C14422T,C14805T,C23707T,G26144T,G26779T,G29140C||
EPI_ISL_576534|2020-10-12|2020-06-29|B.1.1|USA|California|Contra Costa County|11|C241T,C3037T,C14408T,C14422T,G21901A,A23403G,C24863T,C27145T,G28881A,G28882A,G28883C|686-694|
EPI_ISL_1040807|2021-02-21|2020-09-07|B.1.237|South Africa|Western Cape||18|G174T,C241T,C1205A,G2164C,C3037T,C6638T,C11455T,C12115T,C12439T,C14408T,C14422T,C15480T,G18651A,A20268G,A23403G,G23587C,C28854T,C29367T||
EPI_ISL_640120|2020-11-17|2020-09-15|B.1.237|South Africa|Western Cape|Western|21|G174T,C241T,C1205A,G2164C,C3037T,C6638T,C11455T,C12115T,C12439T,C14408T,C14422T,C15480T,T16483A,G16487C,G18651A,A20268G,A23403G,G23587C,C28760A,C28854T,C29367T||

There's 15 submissions with T14895G at covSPECTRUM. [<https://cov-spectrum.org/explore/World/AllSamples/from%3D2020-01-01%26to%3D2024-12-31/variants?nucMutations=T14895G>] In the set of GISAID submissions with a collection date in 2020, the T14895G mutation is found in only one sample that seems to have a lot of errors, because it has over 100 nucleotide changes and many long insertions:

$ grep T14895G gisaid2020.tsv|tr \\t \|

I thought that maybe the two mutations could've been sequencing errors if they were close to the ends of reads. But neither of them was close to the ends of reads. T14895G was close to the middle:

$ samtools view SRR13152144.bam|awk '$4+length($10)+20>=x&&$4<=x' x=14895|cut -f1,10|seqkit tab2fx|cat - sars2.fa|mafft --clustalout --globalpair -
SRR13152144.484 ---------------------------------------gactactatcgttataatcta
SRR13152144.484 ------------------------------------------------------------
MN908947.3      ttctttgctcaggatggtaatgctgctatcagcgattatgactactatcgttataatcta

SRR13152144.484 ccaacaatgtgtgatatcagacaactactatttgtagttgaagttgttgataagtacttt
SRR13152144.484 --------------------acaactactatttgtagttgaagttgttgataagtacttt
MN908947.3      ccaacaatgtgtgatatcagacaactactatttgtagttgaagttgttgataagtacttt

SRR13152144.484 gattgttacgatgggggctgtattaatgctaaccaagtcatcgtcaacaacctagacaaa
SRR13152144.484 gattgttacgatgggggctgtattaatgctaaccaagtcatcgtcaacaacctagacaaa
MN908947.3      gattgttacgatggtggctgtattaatgctaaccaagtcatcgtcaacaacctagacaaa
                ************** *********************************************

SRR13152144.484 tcagctggt---------------------------------------------------
SRR13152144.484 tcagctggttttccatttaataaatggggtaaggctagactttattatga----------
MN908947.3      tcagctggttttccatttaataaatggggtaaggctagactttattatgattcaatgagt

Another Braziliian B.1.1.33 mutation found by alchemytoday

alchemytoday wrote: "I see that one T->C is T16251C which, like A5706G, was also found in Brazil in a B.1.1.33 lineage sample in 2020." [https://x.com/alchemytoday/status/1793313402657308946]

However T16251C is only supported by 1 out of 2 reads that cover the position and it has low quality (16):

$ samtools view SRR13152144.bam|awk '$4+length($10)+20>=x&&$4<=x' x=16251|cut -f1,10|seqkit tab2fx|
cat - sars2.fa|mafft --clustalout --globalpair -
SRR13152144.194 -----------------------------------------------------------t
SRR13152144.194 -----------------------------------------------------------t
MN908947.3      gacatgtattctgttatgcttactaatgataacacttcaaggtattgggaacctgagttt

SRR13152144.194 tatgaagctatgtacacaccgcatacagtcttacaggctgttggggcttgcgttctttgc
SRR13152144.194 tatgaggctatgtacacaccgcatacagacttacaggctgttggggcttgtgttctttgc
MN908947.3      tatgaggctatgtacacaccgcatacagtcttacaggctgttggggcttgtgttctttgc
                *****.********************** *********************.*********
                                                                  ^ T16251C
SRR13152144.194 aattcacagacttcattaaga---------------------------------------
SRR13152144.194 aattcacagacttcattaaga---------------------------------------
MN908947.3      aattcacagacttcattaagatgtggtgcttgcatacgtagaccattcttatgttgtaaa

On GISAID there's only two samples from Brazil which have T16251C. [https://cov-spectrum.org/explore/Brazil/AllSamples/AllTimes/variants?nucMutations=T16251C] One sample is from 2021. The other is B.1.1.33 like the samples with A5706G, but its collection date is only in July 2020. A5706G is found in two sets of B.1.1.33 samples from Brazil: there's 3 samples with collection dates in March-April 2020 from Santa Catarina and Para, and there's also 7 samples with a collection date in August to September 2020 from Bahia. However neither set has T16251C.

The November 2019 run has a total of 25 positions with a mutation relative to Wuhan-Hu-1 if I also include mutations that are only supported by a single low-quality read. Out of 20 mutations that are only supported by a single read, all except 2 have a quality of 20 or below:

$ puptab()(ruby -ane's=$F[4].upcase.gsub(/[.,]/,$F[2]).gsub(/\^./,"").gsub("$","");s2=s.dup;s.enum_for(:scan,/[+-]\d+/).each{m=Regexp.last_match;m.begin(0).upto(m.end(0)+s[m.begin(0)+1,m.end(0)].to_i-1){|i|s2[i]=" "}};s2.gsub!(" ","");t=s2.chars.tally;"ACGT".chars.map{|x|q=$F[5].chars.values_at(*s2.chars.each_index.select{|i|s2[i]==x}).map{|c|c.ord-33};puts [$F[1],$F[2],x,t[x]||0,s2.size,"%.3f"%(100.0*(t[x]||0)/s2.size),q.size==0?"NA":"%.1f"%(q.sum.to_f/q.size)]*" "}' "$@")
$ x=SRR13152144;samtools mpileup -f sars2.fa $x.bam|puptab|awk '$2!=$3&&$4!=0'|(echo pos ref allele count coverage allele% average_quality;cat)|column -t
pos    ref  allele  count  coverage  allele%  average_quality
3392   G    C       2      2         100.000  39.0
5706   A    G       2      2         100.000  39.0
5929   T    G       1      1         100.000  24.0
8576   A    C       1      2         50.000   16.0
8601   T    C       1      2         50.000   16.0
8755   A    T       1      2         50.000   18.0
14408  C    T       2      2         100.000  39.0
14422  C    T       2      2         100.000  38.0
14895  T    G       2      2         100.000  38.5
15106  A    C       1      2         50.000   14.0
15125  T    C       1      2         50.000   18.0
15132  T    A       1      2         50.000   20.0
15133  A    C       1      2         50.000   18.0
15168  G    A       1      2         50.000   16.0
15172  T    C       1      2         50.000   18.0
15189  A    C       1      2         50.000   17.0
15195  T    G       1      2         50.000   14.0
16206  G    A       1      2         50.000   20.0
16229  T    A       1      2         50.000   17.0
16251  T    C       1      2         50.000   16.0
17890  A    G       1      4         25.000   18.0
17893  G    A       1      3         33.333   13.0
19218  T    G       1      1         100.000  15.0
19273  C    G       1      1         100.000  16.0
20488  T    G       1      3         33.333   37.0

However the low-quality mutations listed above didn't include any more mutations which were found in Brazilian B.1.1.33 samples:

$ x=SRR13152144;samtools mpileup -f sars2.fa $x.bam|puptab|awk '$2!=$3&&$4!=0'|awk 'NR==FNR{a[$0];next}$2$1$3 in a' <(grep B\.1\.1\.33 gisaid2020.tsv|grep Brazil|cut -f12|tr , \\n) -|(echo pos ref allele count coverage allele% average_quality;cat)|column -t
pos    ref  allele  count  coverage  allele%  average_quality
5706   A    G       2      2         100.000  39.0
14408  C    T       2      2         100.000  39.0
16251  T    C       1      2         50.000   16.0

Coverage of other B.1.1.33 positions

Positions where B.1.1.33 has mutations had zero depth apart from 14408:

$ awk '$5=="B.1.1.33"' gisaid2020.tsv|cut -f12|awk '{split($0,z,",");for(i in z)a[z[i]]++}END{for(i in a)if(a[i]>NR/2)print i}'|sed 's/.\(.*\)./\1/'|awk 'NR==FNR{a[$0];next}$2 in a' - <(samtools depth -a SRR13152144.bam)|cut -f2-|column -t
241    0
3037   0
14408  2
23403  0
27299  0
28881  0
28882  0
28883  0
29148  0

Apart from 5706 and 14408, there was also zero depth for all positions with mutations in the A5706G samples from Santa Catarina:

$ grep A5706G gisaid2020.tsv|grep Santa\ Catarina|cut -f12|tr , \\n|sed 's/.\(.*\)./\1/'|awk 'NR==FNR{a[$0];next}$2 in a' - <(samtools depth -a SRR13152144.bam)|cut -f2-|column -t
241    0
3037   0
5706   2
14408  2
18803  0
23403  0
27299  0
28881  0
28882  0
28883  0
29148  0

Three mutations that Gadboit found in pangolin viruses

Gadboit posted this tweet: [https://x.com/gadboit/status/1794000682669998429]

A8561C is probably a sequencing error, because it's on the sixth-last base of the read and the base quality is 17. That part of the read got trimmed when I used fastp -35, because it removed segments with low-quality bases from the ends of reads. So I got zero coverage for 8651 when I used trimmed reads:

$ puptab()(ruby -ane's=$F[4].upcase.gsub(/[.,]/,$F[2]).gsub(/\^./,"").gsub("$","");s2=s.dup;s.enum_for(:scan,/[+-]\d+/).each{m=Regexp.last_match;m.begin(0).upto(m.end(0)+s[m.begin(0)+1,m.end(0)].to_i-1){|i|s2[i]=" "}};s2.gsub!(" ","");t=s2.chars.tally;"ACGT".chars.map{|x|q=$F[5].chars.values_at(*s2.chars.each_index.select{|i|s2[i]==x}).map{|c|c.ord-33};puts [$F[1],$F[2],x,t[x]||0,s2.size,"%.3f"%(100.0*(t[x]||0)/s2.size),q.size==0?"NA":"%.1f"%(q.sum.to_f/q.size)]*" "}' "$@")
$ samtools mpileup -f sars2.fa SRR13152144.bam|puptab|awk '$1~/^(5929|8651|16206)$/&&$4'|(echo pos ref al count depth al% average_quality;cat)|column -t
pos    ref  al  count  depth  al%      average_quality
5929   T    G   1      1      100.000  24.0
16206  G    A   1      2      50.000   20.0
16206  G    G   1      2      50.000   16.0

When I didn't trim the reads, I got one read that had A8651C, but the quality of the base with the mutation was 17:

$ minimap2 -a --sam-hit-only sars2.fa SRR13152144_[12].fastq.gz|samtools sort ->SRR13152144.untrim.bam
$ samtools mpileup -f sars2.fa SRR13152144.untrim.bam|puptab|awk '$1~/^(5929|8651|16206)$/&&$4'|(echo pos ref al count depth al% average_quality;cat)|column -t
pos    ref  al  count  depth  al%      average_quality
5929   T    G   1      1      100.000  24.0
8651   A    C   1      1      100.000  17.0
16206  G    A   1      2      50.000   20.0
16206  G    G   1      2      50.000   16.0

The output above shows that G16206A is supported by only one read out of two reads that cover the position, and the quality of the base with the mutation is 20. And the 5929 position is covered by only one read, and the quality of the base with the mutation at 5929 is 24.

So all three of Gadboit's pangolin mutations are likely to be sequencing errors.

This shows that the base with a mutation at 8651 is also surrounded by other low-quality bases, so it got removed by fastp -35 which removes segments with low-quality bases from the ends of reads:

$ samtools view SRR13152144.untrim.bam|awk '$4>8500&&$4<=8651'|cut -f-12|tr \\t \|
$ samtools view SRR13152144.untrim.bam|awk '$4>8500&&$4<=8651'|head -n1|cut -f10,11|tr \\t \\n|cut -c-93
                                  the base at position 8651 has quality 2 = 50-33 = 17 ^
                                    other surrounding bases also have low quality ^^^ ^^   ^

Gadboit also posted this tweet: [https://x.com/gadboit/status/1794284842060034514]

I didn't know if he meant that the bat sequences in his screenshot were also submitted in February 2020, because actually their submission date was in 2021. But Gadboit answered that he only meant that the pangolin sequences were submitted in February 2020.

I told him that the first two mutations are also in Rc-os20. And mutations 1 and 3 are in GZ2021C, GZ2021H, and GZ2021I. And mutations 2 and 3 are in Rc-mk2 and in 10 different Vietnamese sequences from Hassanin et al. 2023:

brew install --cask miniconda;conda init ${SHELL##*/}
conda install -c bioconda nextclade
nextclade dataset get --name sars-cov-2 --output-dir sars-cov-2

curl -sL sars2.net/f/sarbe.fa.xz|gzip -dc>sarbe.fa
nextclade run --input-dataset sars-cov-2 --output-all cladeout <(seqkit seq -g sarbe.fa)
awk -F\\t '$26~/T5929G/&&$26~/A8651C/' cladeout/nextclade.tsv|cut -f2

Gadboit answered: "Yes I've also found some of the muts in those Japanese bat viruses Rc-mk2 and Rc-os20. There are actually four of the 17 in each, but a different four. The PCoVs still seem to be the best match." [https://x.com/gadboit/status/1794309583063978444]

My maximum number of shared mutations was only 3/25, which might be because I trimmed the reads, but there were 15 different sequences on shared first place. Here the first column shows the number of mutations shared with the Brazilian reads, and the second column shows total nucleotide changes from Wuhan-Hu-1 with indels not included:

$ samtools mpileup -f sars2.fa SRR13152144.bam|puptab|awk '$2!=$3&&$4{print$2$1$3}'|awk -F\\t 'NR==FNR{a[$0];next}{n=0;for(i in a)if($26~i)n++;print n,$11,$2}' - <(sed 1d cladeout/nextclade.tsv)|sort -rn|sed -E 's/( orf1ab| surface|,).*//'|head -n20
3 6554 OP837780.1 UNVERIFIED: Sarbecovirus sp. isolate RhGB03
3 5957 DQ412042.1 Bat SARS coronavirus Rf1
3 5893 MK211374.1 Coronavirus BtRl-BetaCoV/SC2018
3 5829 GQ153542.1 Bat SARS coronavirus HKU3-7
3 5815 OK017811.1 Sarbecovirus sp. isolate FJ2021M
3 5805 OK017810.1 Sarbecovirus sp. isolate FJ2021E
3 5805 OK017808.1 Sarbecovirus sp. isolate FJ2021A
3 5794 GQ153543.1 Bat SARS coronavirus HKU3-8
3 5781 OK017809.1 Sarbecovirus sp. isolate FJ2021D
3 5778 OK017807.1 Sarbecovirus sp. isolate AH2021A
3 5528 KF294457.1 Bat SARS-like coronavirus isolate Longquan-140
3 5481 LC663958.1 Sarbecovirus sp. Rc-os20 RNA
3 5417 LC663959.1 Sarbecovirus sp. Rc-mk2 RNA
3 3481 MG772934.1 Bat SARS-like coronavirus isolate bat-SL-CoVZXC21
3 3353 OQ297732.1 MAG: Severe acute respiratory syndrome-related coronavirus isolate GD/L105.18/2022
2 5999 KY770860.1 Sarbecovirus sp. isolate Jiyuan-84
2 5932 DQ648856.1 Bat coronavirus (BtCoV/273/2005)
2 5929 KJ473812.1 BtRf-BetaCoV/HeB2013
2 5921 KJ473813.1 BtRf-BetaCoV/SX2013
2 5910 MZ081377.1 Betacoronavirus sp. RmYN07 strain bat/Yunnan/RmYN07/2020

Mutations relative to WA1

The November 2019 reads cover 3 out of 12 positions where B.1.1.33 has a mutation relative to WA1/proCoV2, but all 3 positions have the allele in B.1.1.33 (which is also the allele in B.1), and all 3 positions are covered by two reads with a high base quality:

$ curl -Ls sars2.net/f/gisaid2020.tsv.xz|xz -dc>gisaid2020.tsv
$ gicon()(cut -f12 "$@"|awk '{split($0,z,",");for(i in z)a[z[i]]++}END{for(i in a)if(a[i]>NR/2)print i}')
$ awk -F\\t '$5=="B.1.1.33"' gisaid2020.tsv|gicon # nucleotide changes in more than 50% of B.1.1.33 samples
$ samtools mpileup -f sars2.fa SRR13152144.bam|awk 'NR==FNR{a[$0];next}$2 in a' <(awk -F\\t '$5=="B.1.1.33"' gisaid2020.tsv|gicon|sed 's/.\(.*\)./\1/';printf %s\\n 8782 18060 28144) -|column -t
MN908947.3  8782   C  2  .,  BE
MN908947.3  14408  C  2  Tt  HH
MN908947.3  28144  T  2  .,  HG

Table of mutations in B.1.1.33 and WA1 posted by Gadboit

Gadboit posted this table which shows all positions with mutations in WA1 or the B.1.1.33 consensus: [https://x.com/gadboit/status/1795723769878843738]

I made a similar table, but I also added the other 3/5 mutations that were called in the November 2019 sample, and I added mutations that were not called in the November 2019 but that were supported by 1/1 reads, and I added the 3 mutations that were called in the March 2020 sample, and I added other mutations that were found in early B.1.1.33 samples from Santa Catarina on GISAID:

GISAID submissions with B.1.1.33

In the set of GISAID submissions with a collection date in 2020, B.1.1.33 is actually by far the most common in Brazil:

$ awk -F\\t '$5=="B.1.1.33"{a[$6]++}END{for(i in a)print a[i],i}' gisaid2020.tsv|sort -rn|head
3366 Brazil
162 USA
71 Chile
51 Uruguay
32 Canada
22 United Kingdom
18 Argentina
16 Paraguay
16 Australia
12 French Guiana

The oldest B.1.1.33 sample by collection date is also from Brazil (even though it might have an incorrect collection date):

$ awk -F\\t '$5=="B.1.1.33"{a[$4" "$6]++}END{for(i in a)print a[i],i}' gisaid2020.tsv|sort -k2|head
1 2020-03-01 Brazil
1 2020-03-07 USA
1 2020-03-07 United Kingdom
1 2020-03-09 Brazil
1 2020-03-09 United Kingdom
1 2020-03-11 Brazil
3 2020-03-12 Brazil
1 2020-03-12 United Kingdom
4 2020-03-13 Brazil
1 2020-03-13 Canada

There's only 5 B.1.1.33 submissions that have an older or same collection date than the earliest submission with A5706G, but 3 out of 5 of them are from Brazil. There's 289 B.1.1 samples with an older or same collection date but only 2 of them are from Brazil. B.1.1.33 is defined by the mutations T27299C and T29148C, where the second mutation might to have arisen earlier, because in the output below the samples that are intermediate to B.1.1 and B.1.1.33 only have the second mutation but not the first one, and the intermediate samples are only found in Europe:

$ reline()(gawk -F\\t '{a[$12][$5]++}END{for(i in a){max=0;for(j in a[i])if(a[i][j]>max){max=a[i][j];out=j FS i}print out}}' ${1-~/gisaid2020.tsv}) # show most common lineage designation for each mutation set, or the first listed designation for ties
$ gian()(cut -f12|tr , \\n|awk -F\\t 'NR==FNR{a[$0];next}{split($12,b,",");for(i in b)if(!(b[i]in a))next}1' - ${1-~/gisaid2020.tsv}) # keep only possible ancestor sequences whose mutation set is a subset of the list of mutations from STDIN
$ grep A5706G gisaid2020.tsv|head -n1|gian|awk -F\\t '$11>=5&&$4<="2020-03-12"'|awk -F\\t 'NR==FNR{a[$1]=$2;next}{b[$11 FS a[$12]FS$6 FS$12]++}END{for(i in b)print b[i]FS i}' <(reline) -|sort -rnk2|csvtk -t pretty -s\ |sed 2d
1  10 B.1.1.33 Brazil         C241T,C3037T,A5706G,C14408T,A23403G,T27299C,G28881A,G28882A,G28883C,T29148C
2  9  B.1.1.33 Brazil         C241T,C3037T,C14408T,A23403G,T27299C,G28881A,G28882A,G28883C,T29148C
1  9  B.1.1.33 United Kingdom C241T,C3037T,C14408T,A23403G,T27299C,G28881A,G28882A,G28883C,T29148C
1  9  B.1.1.33 USA            C241T,C3037T,C14408T,A23403G,T27299C,G28881A,G28882A,G28883C,T29148C
5  8  B.1.1    Switzerland    C241T,C3037T,C14408T,A23403G,G28881A,G28882A,G28883C,T29148C
3  8  B.1.1    Netherlands    C241T,C3037T,C14408T,A23403G,G28881A,G28882A,G28883C,T29148C
2  8  B.1.1    United Kingdom C241T,C3037T,C14408T,A23403G,G28881A,G28882A,G28883C,T29148C
1  8  B.1.1    Italy          C241T,C3037T,C14408T,A23403G,G28881A,G28882A,G28883C,T29148C
1  8  B.1.1    Czech Republic C241T,C3037T,C14408T,A23403G,G28881A,G28882A,G28883C,T29148C
9  7  B.1.1    Vietnam        C241T,C3037T,C14408T,A23403G,G28881A,G28882A,G28883C
9  7  B.1.1    Denmark        C241T,C3037T,C14408T,A23403G,G28881A,G28882A,G28883C
86 7  B.1.1    United Kingdom C241T,C3037T,C14408T,A23403G,G28881A,G28882A,G28883C
8  7  B.1.1    Switzerland    C241T,C3037T,C14408T,A23403G,G28881A,G28882A,G28883C
6  7  B.1.1    Netherlands    C241T,C3037T,C14408T,A23403G,G28881A,G28882A,G28883C
6  7  B.1.1    Israel         C241T,C3037T,C14408T,A23403G,G28881A,G28882A,G28883C
5  7  B.1.1    Spain          C241T,C3037T,C14408T,A23403G,G28881A,G28882A,G28883C
5  7  B.1.1    Portugal       C241T,C3037T,C14408T,A23403G,G28881A,G28882A,G28883C
5  7  B.1.1    Japan          C241T,C3037T,C14408T,A23403G,G28881A,G28882A,G28883C
41 7  B.1.1    Italy          C241T,C3037T,C14408T,A23403G,G28881A,G28882A,G28883C
4  7  B.1.1    Germany        C241T,C3037T,C14408T,A23403G,G28881A,G28882A,G28883C
3  7  B.1.1    New Zealand    C241T,C3037T,C14408T,A23403G,G28881A,G28882A,G28883C
2  7  B.1.1    Mexico         C241T,C3037T,C14408T,A23403G,G28881A,G28882A,G28883C
2  7  B.1.1    Ireland        C241T,C3037T,C14408T,A23403G,G28881A,G28882A,G28883C
2  7  B.1.1    France         C241T,C3037T,C14408T,A23403G,G28881A,G28882A,G28883C
2  7  B.1.1    Brazil         C241T,C3037T,C14408T,A23403G,G28881A,G28882A,G28883C
2  7  B.1.1    Austria        C241T,C3037T,C14408T,A23403G,G28881A,G28882A,G28883C
19 7  B.1.1    USA            C241T,C3037T,C14408T,A23403G,G28881A,G28882A,G28883C
16 7  B.1.1    Sweden         C241T,C3037T,C14408T,A23403G,G28881A,G28882A,G28883C
15 7  B.1.1    Greece         C241T,C3037T,C14408T,A23403G,G28881A,G28882A,G28883C
12 7  B.1.1    Belgium        C241T,C3037T,C14408T,A23403G,G28881A,G28882A,G28883C
1  7  B.1.1    Sri Lanka      C241T,C3037T,C14408T,A23403G,G28881A,G28882A,G28883C
1  7  B.1.1    Singapore      C241T,C3037T,C14408T,A23403G,G28881A,G28882A,G28883C
1  7  B.1.1    Saudi Arabia   C241T,C3037T,C14408T,A23403G,G28881A,G28882A,G28883C
1  7  B.1.1    Morocco        C241T,C3037T,C14408T,A23403G,G28881A,G28882A,G28883C
1  7  B.1.1    Lebanon        C241T,C3037T,C14408T,A23403G,G28881A,G28882A,G28883C
1  7  B.1.1    Cyprus         C241T,C3037T,C14408T,A23403G,G28881A,G28882A,G28883C
1  7  B.1.1    Canada         C241T,C3037T,C14408T,A23403G,G28881A,G28882A,G28883C
1  7  B.1.1    Bolivia        C241T,C3037T,C14408T,A23403G,G28881A,G28882A,G28883C
1  7  B.1.1    Australia      C241T,C3037T,C14408T,A23403G,G28881A,G28882A,G28883C
4  6  B.1.1    Australia      C3037T,C14408T,A23403G,G28881A,G28882A,G28883C
1  6  B.1.1    United Kingdom C3037T,C14408T,A23403G,G28881A,G28882A,G28883C
1  6  B.1.1    Russia         C241T,C3037T,A23403G,G28881A,G28882A,G28883C
1  6  B.1.1    Ireland        C3037T,C14408T,A23403G,G28881A,G28882A,G28883C
1  6  B.1.1    Ireland        C241T,C14408T,A23403G,G28881A,G28882A,G28883C
1  5  B.1.1    Ireland        C3037T,A23403G,G28881A,G28882A,G28883C

However when I didn't exclude samples whose collection date was later than the earliest submission with A5706G, then there were also intermediate samples from Brazil that only had the first mutation but not the second one:

$ grep A5706G gisaid2020.tsv|head -n1|gian|awk -F\\t '$11>=5'|awk -F\\t 'NR==FNR{a[$1]=$2;next}{b[$11 FS a[$12]FS$6 FS$12]++}END{for(i in b)print b[i]FS i}' <(reline) -|awk '$2==8'|csvtk -t pretty -s\ |sed 2d
1  8   USA            C241T,C3037T,C14408T,A23403G,T27299C,G28881A,G28882A,G28883C
1  8   Italy          C241T,C3037T,C14408T,A23403G,G28881A,G28882A,G28883C,T29148C
3  8   Netherlands    C241T,C3037T,C14408T,A23403G,G28881A,G28882A,G28883C,T29148C
4  8   United Kingdom C241T,C3037T,C14408T,A23403G,G28881A,G28882A,G28883C,T29148C
2  8   Czech Republic C241T,C3037T,C14408T,A23403G,G28881A,G28882A,G28883C,T29148C
4  8   Brazil         C241T,C14408T,A23403G,T27299C,G28881A,G28882A,G28883C,T29148C
3  8   Australia      C3037T,C14408T,A23403G,T27299C,G28881A,G28882A,G28883C,T29148C
6  8   Switzerland    C241T,C3037T,C14408T,A23403G,G28881A,G28882A,G28883C,T29148C
1  8   Brazil         C241T,C3037T,C14408T,A23403G,G28881A,G28882A,G28883C,T29148C
19 8   Brazil         C3037T,C14408T,A23403G,T27299C,G28881A,G28882A,G28883C,T29148C
1  8   Australia      C241T,C3037T,C14408T,A23403G,T27299C,G28881A,G28882A,G28883C
1  8   Brazil         C241T,C3037T,C14408T,T27299C,G28881A,G28882A,G28883C,T29148C
3  8   Brazil         C241T,C3037T,C14408T,A23403G,T27299C,G28881A,G28882A,G28883C
1  8   USA            C241T,C3037T,C14408T,A23403G,G28881A,G28882A,G28883C,T29148C

The oldest B.1.1.33 submission from Brazil is also unusual because it has a distance of 6 to its closest neighbors on GISAID, because it has 6 mutations that are missing from other B.1.1.33 samples (except there's a couple of Brazilian B.1.1.33 samples that contain 1 out of 6 of the mutations):

$ awk -F\\t '$5=="B.1.1.33"' gisaid2020.tsv|sed 1q|gicon|gidi|awk '$1<8'|sort -n|head|cut -f1,2,5-8,12,13|tr \\t \||pc
0|EPI_ISL_1739298|2020-03-01|B.1.1.33|Brazil|Rio de Janeiro|15|C241T,C3037T,C5183T,C6113T,C8668T,C14408T,G18756T,T21066A,A23403G,C23683T,T27299C,G28881A,G28882A,G28883C,T29148C
6|EPI_ISL_1139068|2020-03-25|B.1.1.33|Brazil|Mato Grosso do Sul|9|C241T,C3037T,C14408T,A23403G,T27299C,G28881A,G28882A,G28883C,T29148C
6|EPI_ISL_1181594|2020-04-20|B.1.1.33|Brazil|Rio de Janeiro|9|C241T,C3037T,C14408T,A23403G,T27299C,G28881A,G28882A,G28883C,T29148C
6|EPI_ISL_1181603|2020-04-13|B.1.1.33|Brazil|Rio de Janeiro|9|C241T,C3037T,C14408T,A23403G,T27299C,G28881A,G28882A,G28883C,T29148C
6|EPI_ISL_1181606|2020-04-05|B.1.1.33|Brazil|Rio de Janeiro|9|C241T,C3037T,C14408T,A23403G,T27299C,G28881A,G28882A,G28883C,T29148C

However the 6 mutations aren't even that rare, because one of them is found in only 5 submissions, but the others are all found in at least a hundred submissions:

$ awk -F\\t '$5=="B.1.1.33"' gisaid2020.tsv|sed 1q|cut -f12|tr , \\n|sort|comm -23 - <(awk -F\\t '$5=="B.1.1.33"' gisaid2020.tsv|gicon|sort)|awk 'NR==FNR{a[$0];next}$0 in a' - <(cut -f12 gisaid2020.tsv|tr , \\n)|sort|uniq -c|sort -n
      5 T21066A
    163 C6113T
    182 C8668T
    886 G18756T
   1321 C23683T
   1662 C5183T

There's one Brazilian submission from July 2020 and another submission from September 2020 which both share 3 out of 6 of the mutations, and there's also a bunch of pangolin viruses which share 3 out of 6 of the mutations (even though 2 mutations are also shared by several RmYN and RsYN bat viruses which are not shown here):

$ gisum()(awk -F\\t 'NR==FNR{a[$0];next}{split($12,z,",");n=0;for(i in z)if(z[i]in a)n++;print n FS$0}' - ${1-~/gisaid2020.tsv}) # sum of shared mutations
$ printf %s\\n T21066A C6113T C8668T G18756T C23683T C5183T|gisum|awk '$1>1'|sort -rn|cut -f-12|head|tr \\t \|
6|EPI_ISL_1739298|hCoV-19/Brazil/RJ-IC0102020_158-20/2020|2021-04-26|2020-03-01|B.1.1.33|Brazil|Rio de Janeiro||29796|Human|15
3|EPI_ISL_5416113|hCoV-19/Brazil/MG-LBIb2964/2020|2021-10-21|2020-07-15|N.4|Brazil|Minas Gerais|Betim|29393|Human|20
3|EPI_ISL_410543|hCoV-19/pangolin/Guangxi/P3B/2017|2020-02-17|2017|B.1|China|Guangxi||29801|Manis javanica|3888
3|EPI_ISL_410542|hCoV-19/pangolin/Guangxi/P2V/2017|2020-02-17|2017|B.1|China|Guangxi||29795|Manis javanica|4249
3|EPI_ISL_410540|hCoV-19/pangolin/Guangxi/P5L/2017|2020-02-17|2017|B.1|China|Guangxi||29806|Manis javanica|4253
3|EPI_ISL_410539|hCoV-19/pangolin/Guangxi/P1E/2017|2020-02-17|2017|B.1|China|Guangxi||29801|Manis javanica|4257
3|EPI_ISL_410538|hCoV-19/pangolin/Guangxi/P4L/2017|2020-02-17|2017|B.1|China|Guangxi||29806|Manis javanica|4255
3|EPI_ISL_2557352|hCoV-19/Brazil/MG-FIOCRUZ-34167/2020|2021-06-16|2020-09-01|B.1.1.33|Brazil|Minas Gerais|Itaobim|29719|Human|20
2|EPI_ISL_848343|hCoV-19/USA/IL-IDPH-I-000137/2020|2021-01-19|2020-10-08|B.1.578|USA|Illinois|Cook County|29782|Human|22
2|EPI_ISL_848343|hCoV-19/USA/IL-IDPH-I-000137/2020|2021-01-19|2020-10-08|B.1.578|USA|Illinois|Cook County|29782|Human|22

B.1.1 is defined by three consecutive nucleotide changes at positions 28881, 28882, and 28883. If samples that likely have an incorrect collection date are excluded, the earliest sample by collection date with all three mutations is from Moscow, the next sample is from China, and the next 3 samples are from UK (so it's difficult to tell where B.1.1 may have originated, or if the samples from China or Moscow have an incorrect collection date). Most of the earliest sequences are from Europe, but there's also two sequences from Brazil with a collection date on February 26th:

$ grep G28881A,G28882A,G28883C gisaid2020.tsv|awk '$5=="B.1.1"&&!$NF'|head|cut -f1,3-8,11,12|tr \\t \|
EPI_ISL_2737704|2021-06-30|2020-02-15|B.1.1|United Kingdom|England||8|C241T,C3037T,C14408T,A23403G,A27938G,G28881A,G28882A,G28883C
EPI_ISL_2737705|2021-06-30|2020-02-16|B.1.1|United Kingdom|England||8|C241T,C3037T,C14408T,A23403G,A27938G,G28881A,G28882A,G28883C
EPI_ISL_466615|2020-06-11|2020-02-16|B.1.1|United Kingdom|England||10|C241T,C3037T,C14408T,A23223G,A23403G,C25587T,T26483C,G28881A,G28882A,G28883C
EPI_ISL_423662|2020-04-12|2020-02-23|B.1.1|United Kingdom|England||7|C241T,C3037T,C14408T,A23403G,G28881A,G28882A,G28883C
EPI_ISL_2737710|2021-06-30|2020-02-24|B.1.1|United Kingdom|England||8|C241T,C3037T,G11842A,C14408T,A23403G,G28881A,G28882A,G28883C
EPI_ISL_464184|2020-06-11|2020-02-25|B.1.1|United Kingdom|England||8|C241T,C3037T,C14408T,A23403G,A27938G,G28881A,G28882A,G28883C
EPI_ISL_2737721|2021-06-30|2020-02-26|B.1.1|United Kingdom|England||7|C241T,C3037T,C14408T,A23403G,G28881A,G28882A,G28883C