Metagenomics Classifications Across Human And Cat Samples

As previously described, part of the reason why I collected a 2nd Whole Genome Sequencing (WGS) sample from Basepaws in order to have raw data to re-analyze for 2 time points (which include Dental Health Test reports).

One of the points raised in that blog post was to improve coverage for germline variants.  I do see some evidence for improvement in Bastu's heterozygous dilute variant.  However, I primarily wanted to better understand and assess the metagenomics classifications (related to the Oral/Dental Health Test), and those are the primary focus for this blog post.

The exact format for the FASTQ files for the Basepaws raw data was different for the 2nd WGS sample than the 1st WGS sample.  The 1st WGS sample had a format more similar to what was provided by other companies, so I wrote a script to confirm the file format for the 2nd WGS samples.  This requires some familiarity with programming (even if mostly to comment or uncomment lines to run), so I am not sure how much it will help with reformatting for all customers.  Nevertheless, you can see the code and details for the file reformatting and the analysis below here.  A summary of data that I am sharing (and I would be interested to see from others) is provided here.

There are different strategies for metagenomics analysis.  For Basepaws, I believe that either regular coverage WGS or low-coverage WGS (lcWGS) is being used for oral microbiome classifications.  There are some human fecal samples where I describe the details of the sequencing design for samples collected from multiple companies here.  The human oral samples are also a mix of different companies and different sequencing design (such as WGS/lcWGS versus 16S Amplicon-Seq).  All cat fecal samples come from PetQCheck.

While I have not looked in as much depth to the sequencing type and re-analysis details for the other raw data, I have at least one oral and at least one fecal sample for myself and my cat Bastu.

The oral samples for myself are from uBiome, Bristle Health, and any WGS data collected for any purpose (Veritas, Sequencing.com, Nebula, etc.).  The oral samples for my cat Bastu only come from Basepaws.

The fecal samples with raw data for myself come from uBiome, Psomagen/Kean, and thryve/Ombre.  The fecal samples for my cat Bastu only come from PetQCheck.

For previous human comparisons, the company / sequencing design had a noticeable effect on the results.

However, if I go back and use a method that is at least capable of providing some classification for all of the sequencing types (Kraken2/Bracken), then I can re-analyze all of those data types together and visualize how they cluster:

I have added clustering for the species (human or cat) and the collection site (oral or fecal).  If Pearson Dissimilarity (1 - Pearson Correlation Coefficient) is used as the distance metric for clustering, then samples cluster primarily by sample site (oral versus fecal).  For example, all fecal samples intended for microbiome analysis do tend to have Bacteroides classification rates.  

Additionally, Veritas only provided reads aligned to human chromosomes in a BAM (Binary Alignment Mapping) alignment file that I converted to paired-end FASTQ files (the reads used for the classifications above), so I expect all accurate assignments should be relatively close to 0% for that sample.  So, this may be one reason to be skeptical about differences in the Basepaws WGS samples with this classification method, which show differences from each other and vary at a level similar to the Veritas samples should really close to 100% true human reads.

I also created another heatmap, only looking at assignments with the Top 20 genera assignments with highest average read fractions (calculated as percentages after removing Eukaryotic reads, which is also how the values for the first heatmap was calculated):

In that view, it looks like there might be additional Prevotella classifications for the cat fecal samples.  However, I have not looked into the details of the PetQCheck sequencing design, to make sure it matches at least one set of human samples.

Additionally, independent of overall clustering, I think it is encouraging that there is greater detection of Faecalibacterium in the fecal samples than the oral samples (for both species).

Also, unlike the human oral sample that should have only contained human (hg19) reads, the human oral samples showed some similarities that were noticeable beyond the differences in sequencing design.

For advanced readers, you can see some additional testing for alternate normalizations here.  However, I think the main conclusions are the same, and (if anything) there is possibly a less clear distinction for the negative control approximation (and the Basepaws samples) with the other normalizations.

Also more for the advanced readers, but you can also see some troubleshooting where I looked for even coverage in representative sequence alignments here.  If a classification is made for a large fraction of reads but coverage is uneven, then I would say those classifications are more likely to be false positives.  For example, in the above images (or here), you can see some higher reads given Staphylococcus and Pseudomonas assignments in the negative control approximation.  Those representative genomes also have uneven coverage in these Bowtie2 alignments or these BWA-MEM alignments.

Finally, as yet another option for more advanced readers, you can see some alternative methods of running Kraken2/Bracken.  Perhaps of note, there is supposed to be an implementation of KrakenUniq within later version of Kracken2 (referred to as "Kraken2Uniq", where I added the parameters "--report-minimizer-data --minimum-hit-groups 3" based upon this protocol):

Without adding the minimizer parameter for the KrakenUnique implementation within Kraken2, increasing the minimum number of hit groups to 10 ("--minimum-hit-groups 10") also reduces similarity between the approximate human negative control and the Basepaws oral Whole Genome Sequencing samples (at least with the latest versions of Kraken2 and Bracken):

Using a larger database can substantially reduce the run-time.  However, there is some additional information uploaded on GitHub for that option, if using a reduced set of total reads for some samples.  You can see a large difference in the Komagataeibacter classification rate for the larger database versus the more stringent requirements for the smaller "MiniKraken2 database.

In other words, I think these are important take-home messages:

• Human and cat fecal samples form a cluster for that collection site, and they separate by species within that cluster (although PetQCheck is the only source for any cat fecal samples).
• The Veritas sample was provided to me as chromosome alignments, rather than reads.  By definition, this means that the reads were aligned to the human genome, so I consider "Vertias_WGS" to be a negative control approximation.
• Both the human fecal and oral samples contain a mix of Whole Genome Sequencing (WGS) and targeted sequencing for the 16S bacterial gene (16S Amplicon-Seq).  So, when these samples cluster together, I consider that a relatively robust result (even if there are also additional factors contributing to the exact abundance values, such as the company and the sequencing design).
• On the contrary, the 2 Basepaws samples are less similar to each other than the different companies / libraries for my human oral samples.  Additionally, if there no no additional troubleshooting, then there is relative clustering between the Basepaws oral cat samples (particularly Basepaws_WGS2) and the human negative control approximation.

In term of trying to understand the Basepaws cat oral samples, there is a noticeable negative shift in the fragment size estimated by Picard (particularly for the first sample):

If the classifications can be reproducibly run in the same way, then I hope sharing reprocessed data might help with helping develop more of an opinion of the cat Oral/Dental Health Test from Basepaws.  I have started a discussion thread here.  In the Basepaws oral health test preprint, I also have additional comments, and I would be interested to hear feedback from others.  My understanding is that Kraken2/Bracken was also used in the preprint, but I think something about the exact way that the program was run was different.  Also, if I understand correctly, then I don't think the program can take the paired-end read information into consideration if the file is in an interleaved format (as provided for my cat's 2nd Basepaws WGS sample); however, I am not sure how important that is.  Within the GitHub discussion, I also provide some information for my human oral microbiome report from Bristle Health.

I would say the Basepaws Oral/Dental Health Test risk assessments were not a great match for my cat's individual dental health assessed by the vet.  However, I understand that the risk assessment is not meant to provide a diagnosis, and I am therefore still not completely sure what I think about the results.

That said, I feel comfortable saying that I wish raw data was provided for all of the Oral/Dental Health Test results from Basepaws (not only for the Whole Genome Test).

So, I hope this contributes some interest and encourages discussion on the topic!  I have some notes for the Basepaws Oral/Dental Health Test here, but I would be interested in seeing re-processed raw data with similar classifications (if possible).

Change Log:

12/20/2022 - public post date.

1/1/2023 - add content related to Kraken2/Bracken troubleshooting + revised take-home messages.

1/4/2022 - add note for Komagataeibacter assignments in Basepaws WGS1.

Disqus / Twitter Follow-Up: Comparing My 23andMe SNP Chip Concordance with Different Veritas WGS files

I am mostly summarizing my points from the following Twitter discussion:

https://twitter.com/carolinefwright/status/1157219572514209792

The original intention was to add this as a comment in the Disqus comment in the original pre-print discussion.  However, I thought this was fairly long, and I wanted to be able to have a little more control over figure formatting.

For reference, I recommended taking a look at Illumina arrays in an earlier comment, and I mentioned there there are datasets with both 23andMe SNP chip data and high-throughput sequencing data (like myself).

As a success story, that comment was followed up on, and new data was added.  In particular, there was a plot of the error/discordance rate between my own 23andMe data and my Veritas WGS data posted on Twitter.

I think this is great, but I think it may be worth emphasizing that re-processing my Veritas WGS data resulted in better concordance with my Exome data.  Additionally, I have an upgraded V5 chip, so there are actually 2 sets of 23andMe Data (although almost all of these results are based upon my later set of V3_V5 genotypes).

Nevertheless, I actually have 2 .vcf files for my Veritas WGS data (the provided .vcf, and the .vcf that I produced from extracted FASTQ files and reprocessed using BWA-MEM and GATK).

I did some analysis with my V3 23andMe genotypes, but I think that was mostly consistent with my V3_V5 genotypes.  Among probes on both arrays, there were only 5 discordant sites (so, I think that is how the previously lower error rate was reported: SNP chip versus same SNP chip, instead of SNP chip versus WGS).

In contrast, even the best-case scenario for my data seems to have concordance of 97.6-99.2% (with MAF > 0.01 variants), and this was slightly lower for the V3_V5 genotypes.  If I only considered my original V3 genotypes, this would have been better than I had previously reported for my Exome versus WGS data (98-99% for BWA-MEM + GATK re-processed variants).  However, either way, there are over a million probes on the V5 23andMe array.  So, I think something about variants being used for “research purposes” may be relevant (although I will show below that certain sets of variants do have higher reproducibility).

The trends can vary depending upon what I use for the MAF calculation:

1000 Genomes:

gnomAD:

Kaviar:

I am showing 2 plots each because I have 2 .vcf files for my WGS data (the one provided by Veritas, and the BWA-MEM + GATK re-processed variant file).  While the results are above are mostly similar with either Veritas WGS .vcf, there are some noticeable differences in the exact set of rare variants with different population estimates (which could either be because of the composition of individuals, or with the different sample processing strategies for each project).

When I converted my 23andMe to VCF format, I added a “PASS” status (to keep track of variants within repeats, for example).  If there was no other note, the variant had a “PASS” in the FILTER column.  If I only consider “PASS” variants, this is what those plots look like:

1000 Genomes:



  
gnomAD:



  
Kaviar:
 


If I only consider the “PASS” variants, the accuracy also increases a little (from 98.8-99.7% for my initial V3 genotypes, and 98.009-99.997% with the range of MAF for my V3_V5 genotypes), usually in the far-left column.

Finally, if I start from the full set of variants, but I only look at those that were discordant between my WGS .vcf files, I see better concordance for re-processed variants if they are common (but the trend for the smaller variant sets can vary):

1000 Genomes:



  
gnomAD:



  
Kaviar:






For these plots, the largest number of variants are in the far-right column.  Since that far-right column usually has less SNP chip discordance (less red in the barplot), that is consistent with my earlier conclusion that re-processing the WGS data could produce variants with higher overall concordance (in that situation, between Exome and WGS data).  However, this can clearly vary between individual variants/positions (and the best processing strategy may vary depending upon where you need to be calling variants).

I can't really visualize the SNP chip data, and I kind of have to trust the "NC" status for "No Call" positions (and I don't have access to a more raw form of data, like the intensities).

However, as a general rule, I would always recommend checking your alignments for false negatives or false positives (which you can do with a free genome browser like IGV).  I added some ClinVar annotations to try and find some discordant sites to check, and I've listed a couple below:

1) I was surprised that my re-processed was missing my cystic fibrosis variant (which I have a whole other post about).  However, this was purely a formatting issue.

Namely, I threw out most indel positions (indicated by DI) because figuring out exactly what that represents is more difficult than for the SNPs.  However, with my earlier V3 chip analysis, I manually converted some indels in my code (including my cystic fibrosis indel).  However, I converted to match the freebayes indel format in the provided .vcf for the .

In the other blog post, you can clearly see that I am a cystic fibrosis carrier, even with the re-alignment.  So, I looked up the GATK format for that indel and checked the status of that variant.  Indeed, I do have a variant call for chr7 117149181 . CTT C.  So, this is in fact OK (as long as the annotation software can figure out I have the variant, which I think was part of the problem described in the other blog post).

2) While it wasn't a discordant site between .vcf files, there were only a limited number of total ClinVar pathogenic variants.  So, I happened to notice one that indicated I was homozygous for the pathogenic variant in my 23andMe data (1/1) but a variant was not called at that position with either of the Veritas WGS .vcf files (indicated by a 0/0 in the genotype columns).  In other words, the SNP chip data was consistent with a SNP chip replicate, the WGS variant was robust different processing methods, but the result was different for SNP chip versus WGS.

However, I can check the alignments for both the provided and reprocessed WGS data (as well as provided and reprocessed Exome data), which is what I do below:

The plot above shows alignments in the following order (top to bottom): Genos Exome Provided, Genos Exome Reprocessed, Veritas WGS Provided, Veritas WGS Reprocessed

So, from the alignments, I would be inclined to agree with the WGS variant calls.  If there is some isoform of the gene that is so far diverged the reads wouldn't align, that could be an exception.  However, I don't think that situation would be best described as a SNP.  Also, I don't believe any reports indicated that I had predisposition to Neurofibromatosis (and I don't know anybody in my family that was diagnosed with that disease).  This was a custom 23andMe probe (labeled as i5003284, instead of a typical rsID).  However, the larger ANNOVAR annotation file has the ClinVar information, and I can find a dbSNP ID using the UCSC Genome Browser (for hg19, chr17:29541542).  So, in terms of checking the ANNOVAR annotation, if I actually did have two copies of an NF1 pathogenic variant, rs137854557 does have multiple reports indicating it being pathogenic (less of a confident assertion than my cystic fibrosis variant, but more than most of my other "pathogenic" variants).

As a final note, the data and code for analysis of my V3_V5 23andMe genotype data (and 2 Veritas WGS .vcf files) is available here.

Update Log:

8/10/2019 - public post date
1/15/2020 - add tag for "Converted Twitter Response"