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.

Testing Limits of Self-Identification / Relatedness using Genomic FASTQ Files

Color asked me sign a HIPAA release in order to get access to raw genomic data, which included a FASTQ file with ~15,000 reads.  So, I thought it might be useful to get an idea about how few reads (from random low-coverage Whole Genome Sequencing) can be used to identify myself.

To be clear, I think most rules are meant to take possible future advances into consideration.  So, just because I can't identify myself, doesn't mean somebody else can't identify me with fewer reads (with current of future methods).  Nevertheless, if I can identify myself, then I think that there is a good chance others could probably identify themselves with a similar number of reads and/or variants (and possibly fewer reads/variants).

Down-Sampling 1000 Genomes Omni SNP Chip Data

I compared relatedness estimates for myself with the following genotypes:  1) Veritas WGS, 2) 23andMe SNP chip, 3) Genes for Good SNP chip, and 4) Nebula lcWGS (along with the matching positions from the 1000 Genomes Omni SNP chip).

Perhaps more importantly, I also show kinship/relationship estimates (from plink) for 1000 Genomes samples for parent-to-child relationships as well as more distant relationships:

As you can see, there is a bit more variability in the parent-to-child estimates with a few thousand variants.  The self-identification estimates (among pairs of my 4 samples) were always greater than 0.45, but there is noticeable overlap in the kinship estimates for 1000 Genomes parent-to-child and more distant relatives when you drop down to only using 19 variants.

So, making sure you didn't get false positives for close relationships may be important, particularly with smaller numbers of variants.  If you have SNP chip or regular Whole Genome Sequencing data, then identifying yourself would also be easier than having 2 low-coverage Whole Genome Sequencing datasets.

However, if I can get 1000s (or perhaps even 100s) of variant calls, I am currently most interested in how accurate those calls can be.

Gencove and STITCH Imputed Self-Identification

I have earlier posts showing that the Gencove imputed variants from Nebula were not acceptable for individual variant calls, but I think they provided reasonable broad ancestry and relatedness results.  To be fair, I don't believe Nebula is currently providing low coverage Whole Genome Sequencing results anymore, opting for much higher coverage (like regular Whole Genome Sequencing).  However, Color provided me with considerably fewer lcWGS reads than Nebula (and Color also has a pre-print about lcWGS Polygenic Risk Scores that I was concerned about).

So, I was interested in testing what imputed variants I could get if I uploaded FASTQ files for Gencove analysis myself (as well as an open-source option called STITCH).

There is also more information about running STITCH (as well as more statistics for Gencove variant concordance) within this subfolder (and this subfolder/README) on the human GitHub page.  Essentially, the performance of the human lcWGS looks good at 0.1x (if not better than the earlier Gencove genotypes that were provided to me from Nebula), but there is a drop in performance with the cat lcWGS.

I ran the STITCH analysis on a local computer, so the run-time was longer than Gencove (between 1 day and 1 week, depending upon the number of reference samples - hence, I would start with running STITCH with ~99 reference samples in the future).  However, if you were willing to pay to run analysis on the cloud (or use more local computing power), I think the run-time would be more similar if each chromosome was analyzed in parallel.  Also, STITCH is open-source, and doesn't have any limits on the minimum or maximum number of reads that can be processed.  The performance also looks similar with ~5 million 100 bp paired-end reads, so the window for more accurate results that can be returned from Gencove may be around 2 million reads.  So, I think using STITCH can have advantages in a research setting.

I welcome alternative suggestions of (open-source) methods to try, but would tentatively come up with these suggestions (for 100 bp paired-end reads, with random / even coverage across the genome):

greater than 1 million reads: good chance of self-identification

0.1 - 1 million reads: intermediate chance of self-identification (perhaps similar to patient's initials, if it narrows down a set of family members?).  Potentially "good" chance of self-identification with other methods and/or future developments.

less than 0.1 million reads: respect general privacy and allow for future improvements, but additional challenges may be countered.  There still may also be sensitive and/or informative rare variants.

I also added the results for GLIMPSE lcWGS imputations (Rubinacci et al. 2020).  These are human results, but the performance was a little lower than STITCH (more similar to the Gencove results for my cat, but lower than the Gencove results for myself).  However, it probably should be noted that I did not specify my ancestry for GLIMPSE but I did specify a limited number of populations for STITCH.  So, if you don't know the  ancestry (or the ancestry used might  confound the results), then that loss of concordance may be  OK.  Also, I think the  GLIMPSE run-time was shorter (within 1 day) and it used the full set of 1000 Genomes samples as the reference set.

Recovery/Observation of Variants in 1 Read

Even if I don't exactly know what is most likely to be able to self-identify myself, I can try to get some idea of the best-case scenario in terms of even having 1x coverage at a potentially informative variant position.

The numbers here are a little different than the 1st section: I am looking for places where my genome varies from the reference genome, and I am considering a larger number of sites.  Nevertheless, I was curious about roughly how many reads it took to recover 500 or 1000 variants from a couple variant lists:

Notice the wide range of any SNPs that can be recovered versus a set of potentially informative SNPs.  However, it looks like you may very roughly notice problems with self-identification (for even future methods) with less than ~250,000 reads (matching the STITCH results above).  This is noticeably more than the ~15,000 of reads that I had to sign a HIPAA release to get from Color, but that lower limit (for any SNPs, called "WGS SNPs" with the gray line) was ~15,000 reads (very similar to what I was provided from Color, albeit single-end instead of paired-end).

In reality, you have sequencing error and a false discovery rate to consider when calling variants from 1 read, the nucleotide distribution at each position is not random between the 4 nucleotides, linkage (non-independence) of variants, and you have 2 copies of chromosomes for each position in the genome reference.  However, if you over-simplified things and asked how many combinations of 4 nucleotides (or even 2 nucleotides) create more unique sequences than the world population, that is noticeably less than the 500 or 1000 thresholds added to the plot above.

So, if you consider rare SNPs (instead of calculating a relatedness estimate, with common SNPs), perhaps you could identify yourself with less than 50,000 reads?  Either way, if you give some rough estimate allowing for future improvements in technology, I would feel safe exercising extra caution with data that has at least 100,000 reads (collected randomly / evenly across the genome).  I also believe that erring on the side of caution for data with fewer reads is probably wise as a preventative measure, but I think the possible applications for that data is lower.

Closing Thoughts

If anybody has knowledge of any other strategies, I am interested in hearing about them.  For example, I think there may be something relevant from Gilly et al. 2019, but I don't currently have a 3rd imputation benchmark set up.  I have also tried to ask a similar question on this Biostars discussion, since it looks like Gencove is no longer freely available (even though I was able to conduct the analysis above using a free trial).  I have all of this data publicly available on my Personal Genome Project page.

I am also not saying that it is not important to consider privacy for samples with less than any of the numbers of reads that I mention above (with or without a way to self-identify myself with current methods).  For example, the FASTQ files have information about the machine, run, and barcode in them (even with only 1 read).  So, if the consumer genomics company had a map between samples and customers, then perhaps that is worth keeping in mind for privacy conversations.  Likewise, if the smaller number of variants includes disease-related variants, perhaps that is also worth considering.

I don't want to cause unnecessary alarm: as mentioned above, I have made my own data public.  However, you do have to take the type of consent into consideration when working with FASTQ files (for data deposit and data sharing).  For example, you currently need "explicit consent" for either public or controlled access of samples collected after January 25th, 2015.

Finally, I would like to thank Robert Davies for the assistance that he provided (in terms of talking about the general idea, as well as attempting to use STITCH for genotype annotations), which you can see from this GitHub discussion.  I would also like to thank several individuals for helping me learn more about the consent requirements for data deposit.

Additional References

I am interested to hear feedback from others, potentially including expansion of this list.

However, if it might help with discussion, here are some possibly useful references:

Selected Genomic Identifiability Studies (or at least relevant publications):

• Sholl et al. 2016 - Supplemental Methods describe using 48 SNPs to "confirm patient identity and eliminate sample mix-up and cross-contamination".
• McGuire et al. 2008 - article describing genomics and privacy with some emphasis on medical records
• Oestreich et al. 2021 - article generally discussing genomics identifiability and privacy
• Ziegenhain and Sandberg 2021 - in theory, considers methodology to provide anonymized processed data.  I have not tried this myself, but this could at best maximize downstream analysis.  That might be useful in that it expands processed data that could be shared with caveats.  However, some things require accurate sequencing reads as unaltered raw data.  Modified sequences should not be represented as such "raw" data.
• Wan et al. 2022 - article generally discussing genomics identifiability and privacy
• Russell et al. 2022 - instead of sequencing coverage (such as in this blog post), this preprint describes the impact of the amount of starting DNA material on microarray genotyping for forensics analysis.
• Kim and Rosenberg 2022 - preprint describing characteristics affecting identifiability for STR (Short Tandem Repeat) analysis
• Popli et al. 2022 - a preprint describing kinship estimates in low coverage sequencing data (and the amount of data for relatedness estimates is the topic for most of the content in this blog post)

While related to the more general topic, I think the goals of Lippert et al. 2017 and Venkatesaramani et al. 2021 and are somewhat different than what I was trying to compare (including image analysis).

I also have some notes in this blog post, some of which are from peer reviewed publications and some from other sources (such as NIH and HHS website).  Again, I would certainly like to learn more.

In addition to the publications for STITCH/GLIMPSE/Gencove, other imputation / low-coverage analysis studies include Martin et al. 2021, Emde et al. 2021 and GLIMPSE2.  Hanks et al. 2022 also compares microarray genotyping and low coverage imputation to Whole Genome sequencing.  If it expected that low coverage analysis includes enough markers to be useful, then I think you are either directly or indirectly saying that level of coverage is sufficient to identify the individual (which I think is a criteria that is likely easier to meet than clinical utility).

I certainly don't want to cause any undue concern.  I think some public data is important for the scientific community, but I think it is appropriate for most individuals to agree to controlled access data sharing.  Nevertheless, I think this is an important topic, which might need additional communication in the scientific community.

Change Log:

3/8/2020 - public post
3/9/2020 - add comment about simplified unique sequence calculation
4/8/2020 - add STITCH results + minor changes
4/9/2020 - minor changes
4/16/2020 - minor change
5/1/2020 - add link for GLIMPSE (before analysis)
7/28/2020 - add GLIMPSE results

1/15/2023 - add additional references for other studies

3/16/2023 - add 48 SNP verification

PRS Results from my Genomics Data (mostly from impute.me)

I haven't had a whole lot of personal experience with Polygenic Risk Score (PRS) estimates, so I thought it was interesting when I found a couple options for re-analysis of my own genomics data (for selected examples):

Association	SNP chip (impute.me) (Folkersen et al. 2020)	Other Re-Analysis Options	23andMe Results
Type 2 Diabetes (No)	(Type 2 Diabetes, 146 variants) Average / Above Average (23andMe-V3, 12/19) Average / Above Average (AncestryDNA, 12/19)	MySeq 1.000 risk ratio [error] (Nebula lcWGS) 0.955 risk ratio (Genos Exome, 3 variants) 1.089 risk ratio (Veritas WGS, 6 variants)	"Typical Risk" of 23% (directly from 23andMe, PRS with 1,244 loci) [actually, slightly lower than normal] Reduces to less than 1% when age, height, weight, fast food consumption, and exercise rate are taken into consideration (also from 23andMe)
Ulcerative Colitis (once, so I think really "no")	(23 variants, and 116 variants) Both Below Average and Above Average Risk, for different PRS (23andMe-V3, 12/19) Both Below Average and Above Average Risk, for different PRS (AncestryDNA, 12/19)
Anxiety Disorder (Yes, but getting better)	(6 variants) Average / Above Average (23andMe-V3, 12/19) Average / Above Average (AncestryDNA, 12/19)
Migraine (Periodic)	(26 variants, and 21 variants) 2 PRS (Average and Above Average) (23andMe-V3, 12/19) 2 PRS (Average and Above Average) (AncestryDNA, 12/19)
Eye Color (Light Brown)		DNA.land Likely to have Brown Eyes (23andMe-V3, 12/19) Likely to have Brown Eyes (23andMe-V3_V5, 12/19) Likely to have Brown Eyes (AncestryDNA, 12/19)	23andMe reports that I am expected to have "brown or hazel eyes" based upon 1 SNP (rs12913832)
Hair Color (Light Brown)	See Below (Roughly 25% Red and 50% Blonde)		For "Light or Dark Hair", 23andMe reports that I have "Likely Dark" Hair (using 42 SNPs) For "Red Hair" 23andMe reports that I am "Unlikely to have red hair" (using 3 MC1R SNPs: rs1805007, rs1805008, and another custom MCR1 probe)
Height (180 cm)	See Below	DNA.land 171 cm: "Likely Taller than Average" (23andMe-V3, 12/19) 171 cm: "Likely Taller than Average" (23andMe-V3_V5, 12/19) 171 cm: "Likely Taller than Average" (AncestryDNA, 12/19)

Individual SNP risks were reported (from impute.me).  While I had a bit of a hard time finding the precise overall risk estimate (without trying to sum / multiply separate risks), this might be OK in terms of getting a sense of whether I was an outlier or not.  For example, being above or below average for "Type 2 Diabetes" seemed to vary (unless you say most people were under something like a null distribution for "average" risk).  In other words, I thought the following plots (which you could see for various traits) were interesting:

impute.me Type 2 Diabetes PRS (23andMe V3)

impute.me Ulcerative Colitis (1st entry, 23andMe V3)

impute.me Ulcerative Colitis (2nd entry, 23andMe V3)

impute.me Anxiety Disorder PRS (23andMe V3)

impute.me Migraine-Broad PRS (23andMe V3)

impute.me Migraine PRS (23andMe V3)

impute.me Hair Color (23andMe V3 + Ancestry DNA, respectively)

impute.me Height (23andMe V3)

I thought the anxiety disorder result was interesting for 2 reasons.  First, I have had issues with anxiety problems (for example, you can click here for notes, even though they are primarily related to PatientsLikeMe).  Second, notice the environmental component is larger than the genetics component.  This matches my concerns that I expressed in this review of "blueprint".  For example, I would say the predictive power from birth has some notable limitations (such as difficulties in the need to take medication at any given point in your life).

While I am not sure if the exact right term was used (since I thought "Ulcerative Colitis" was a condition, rather than a symptom).  However, I was hospitalized for Ulcerative Colitis (even though that was a one time occurrence caused from E. coli with Shiga toxin).

I also get migraines.

I don't have Type 2 Diabetes, but I provided that because I also had other PRS results to compare.  Similarly, if others have suggestions where I can quickly compare to the impute.me PRS results, please let me know and I would be very happy to add them!

For example, I did add DNA.land (and 23andMe) Eye Color and Height based upon a Twitter response.  While I think height is one of the more heritable traits, DNA.land couldn't guess my actual height within a few inches (and there is a noticeable spread of points for the impute.me plot above).  Even though DNA.land gave lower confidence to other predictions, I would say these have been "fair" rather than "high" confidence (and everything else probably should have been "low" confidence).  I am close to the diagonal for the impute.me plot, but I don't know if the scale is 1:1.  For example, my DNA.land height prediction was off by 3-4 inches.  However, to be fair, note that the highest and lowest percentiles for high don't have overlap (there are not any points in the upper-left or bottom-right regions of the scatter plot, even though those make up a smaller fraction of the population).

For comparison, here is the distribution of score for DNA.land (where my true height was greater than anything on the density distribution - perhaps because this was height scaled for female percentiles?):

For impute.me, the predicted hair color shows blondness on the x-axis and redness on the y-axis.  The cyan circle is my actual color (which I filled in), and the while circle is my predicted color.  I think my hair color used to be lighter than it is now (and I think the shade that I reported for myself was a bit too dark), so that is closer to the genetic prediction (perhaps half-way between).

It may be worth noting that 23andMe could predict that I had brown hair and eyes (although I think that covers most people and you need the more rare traits to better calculate accuracy - for example, Francis Collins said that his 23andMe report indicated he had brown eyes when he really had blue eyes, at least 10 years ago).

Again, for comparison, here is the distribution of DNA.land scores for eye color:

I didn't add the AncestryDNA density plots since they looked qualitatively similar to the 23andMe V3 plots (and, on another computer, I had an issue with the percent variance explained appearing in a pie chart that was harder to read).  I also originally intended to test my updated 23andMe genotypes (V3+V5), but I got an error saying that data was already uploaded (from my V3 chip).  However, perhaps I can test those results later, and see if they are still similar.

With a $5 donation, the turn-around time for processing was 1-3 days.

For Genos Exome and Veritas WGS data, I used the BWA-MEM Re-Aligned GATK Variant calls.  However, I think the main conclusion from looking at my diabetes results was that I was of average risk, and I don't believe my own genetic diabetes PRS risk assessment was great without taking additional factors into consideration (for 23andMe, that was a difference between 23% and 1%, after considering BMI, diet, and exercise).

This essentially matches Supplementary Figure S12 for this paper (whose title I respectfully believe can give the reader the wrong impression, and there is at least one objective error that I believe needs to be corrected), where absolute risk explained was usually very low (usually explaining less than 15% of the variation for a trait).  You can also see that the variability explained by "this score" for the impute.me PRS above is estimated to be less than half of the genetic component.

I think the preprint by Brockman et al. 2021 might also have some additional relevant information for this discussion.

In somewhat different contexts, you can also see some notes / concerns about percentiles / indices in the posts on Nebula and basepaws lcWGS results.

Change Log:

12/5/2019 - public post
12/7/2019 - add DNA.land results based upon Twitter reply from Debbie Kennett; revise wording in post
12/8/2019 - mention possible scaling for female height; also fix date for previous log entry.
6/25/2020 - add links to posts with Nebula and basepaws results.  Minor formatting changes.
7/7/2020 - add reference to impute.me paper

4/22/2021 - add reference to another paper

2/4/2024 - change column labels to be more precise

Concerns About Using Low-Coverage Sequencing for Trait or Health Results

This is a subset of my notes from my Nebula lcWGS sequencing on GitHub:

NOTE (2/24/2020): Nebula is currently offering 30x sequencing.  So, my concerns about the low coverage Whole Genome Sequencing (lcWGS) at ~0.5x are probably less relevant for that particular company.  However, if you get lcWGS from another company, then this information is probably relevant.

Concerns about Specific Variants

While I very much support providing FASTQ, BAM and VCF data, one of my concerns about the Nebula results was the use of low-coverage sequencing.

So, one of the first things that I did was visualize the alignments for some of my more confidently understood variants from previous data (using IGV).

For the two alignments below, the Genos Exome is the top alignment, the Nebula low-coverage alignment is in the middle, and the Veritas Whole Genome Sequencing (WGS, regular-coverage) is at the bottom.

My cystic fibrosis variant (rs121908769):

My APOE Alzhiemer's risk variant (rs429358, Nebula alignment in middle, variant is red-blue bar in the right-most exon):

For APOE, I zoomed out from the screenshot so that you could get a better perspective of the error rate per-read at other positions around the gene.

You could see my cystic fibrosis variant in the 1 read covered at that position, but you can't see any reads with the APOE variant.  My concern about the use of low-coverage sequencing is due to imputation (at least for traits).  Even though this APOE variant is somewhat common (I believe ~15% of the population), the imputation failed to identify me as having that variant.  You see that from the .vcf files

My APOE Alzhiemer's risk variant:

19      45411941        rs429358        T       C       .       PASS    .       GT:RC:AC:GP:DS  0/0:0:0:0.923102,0.0768962,1.71523e-06:0.0768996

As described in the gVCF header:

GT = Genotype
RC = Count of Reads with Ref Allele
AC = Count of reads with Alt Allele
GP = Genotype Probability: Pr(0/0), Pr(1/0), Pr(1/1)
DS = Estimated Alternate Allele Dosage

The "0/0" (for genotype/GT in the last column) means that low-coverage imputation couldn't detect my APOE variant.  In other words, I believe Nebula incorrectly estimates my genotype to be 0/0 with a probability of 92.3%, and the probably for the true genotype was 7.7%.  I also see a blog post mentioning that these probabilities are provided to users through the web-interface, although I am having difficulty in finding them without the gVCF (and you won't see them in the PDFs that I have uploaded in this section).

Update (8/5): Nebula support got in touch with me and explained that the blog post is in reference to Nebula Research Library (rather than the "Your Traits" section).  I canceled my subscription, I am still able to confirm that I see this under the "Library" section (rather than "Traits," "Ancestry," or "Microbiome").

Likewise, there was no delTT variant in the VCF, so my cystic fibrosis carrier status would also be a false negative (if that was used in the report), even though you could actually see that deletion in the 1 read aligned at that position (because 1 read wasn't sufficient to have confidence in that variant).

Overall Variant Concordance

I can also use my VCF_recovery.pl script to compare recovery of my Veritas WGS variants in my Nebula gVCF.

If you compare SNPs, then the accuracy is noticeably lower than GATK (and even lower than DeepVariant):

3,071,596 / 3,419,611 (89.8%) full SNP recovery
3,184,641 / 3,419,611 (93.1%) partial SNP recovery

The indels are harder to compare (becuase of the freebayes indel format).  So, in the interests of fairness, I am omiting them here (as I did for comparing the provided Genos Exome versus Veritas WGS variants).  However, instead of comparing the provided Veritas WGS .vcf file, I can try comparing the BWA-MEM re-aligned GATK Veritas WGS .vcf (which also had higher concordance between my Exome and WGS datasets):

3,133,635 / 3,419,611 (91.6%) full SNP recovery
3,248,277 / 3,419,611 (95.0%) partial SNP recovery
164,140 / 217,959 (75.3%) full insertion recovery
180,736 / 217,959 (82.9%) partial insertion recovery
190,452 / 266,479 (71.5%) full deletion recovery
213,131 / 266,479 (80.0%) partial deletion recovery

The GATK recovery is a little better.  However, it is very important to emphasize that the gVCF variants do not have 99% accuracy (even for average accuracy, or even for SNPs).  I think whatever benchmark was used for that calculation was probably over-fit on some training data.  To be fair, I think the average SNP chip concordance (with higher coverage WGS data) is also lower than some people might expect, but it is definitely higher than this lcWGS data.

You can also show similar results with precisonFDA (using the BWA-MEM realigned GATK gVCF, which we expect to have better concordance than the provided Veritas gVCF).

For example, the overall file shows noticably low recall when comparing the Nebula imputed gVCF versus the Veritas WGS BWA-MEM re-aligned gVCF:

and, to be more fair for the Exome versus WGS comparison in the blog post, the trend is similar within RefSeq CDS regions:

The screenshots are smaller than in the blog post because there was no precision-recall plot for the Imputed Nebula gVCF comparisons.

So, I  disagree with the use of low-coverage sequencing for traits, and I would respectfully consider removing this section (or only made available to those with higher-coverage sequencing).

When I was trying to upload my raw data to my Personal Genome Project page, I noticed that they had an option called "genetic data - Gencove low pass (e.g. Nebula Genomics)".  This makes me think discouraging low-coverage sequencing is something that needs to be done more broadly (at least for health traits).

Concerns abut Nebula Library Results

My concerns for the previous sections are probably solved when using the higher coverage sequencing data.  So, unless you were an earlier customer and had the lower coverage sequencing data, you probably don't have to be extra careful about possibly overestimated accuracy in your genotype imputations.

However, there is one thing that I think could still be a problem for customers with higher coverage sequencing data (if the reports are the same).  The concept is similar to my concern about the basepaws breed index (described in this blog post) and/or other Polygenic Risk Scores that I have collected for myself, but I think I can explain my concern with the top 3 percentile results that I received from Nebula:

As you can see from this link, the percentile above was calculated using 13 SNPs.  Seven of the thirteen SNPs are on chromosome 6, and 5/7 of those variants didn't have alignments against the main reference chrososome (for hg19) in my higher coverage Veritas Whole Genome Sequencing data.  Nebula predicted 3-4 of those 7 chromosome 6 variants to be homozygous variants, but I am not sure if these are correct or not (and the nucleotide for rs3763312 was different from the variants in dbSNP).  There were a pair of variants on chromosome 10 where I was predicted to be heterozgyous at both sites.  For the remaining 6 non-chr6 variants, I had imputed genotypes for 2 heterozygous variants, 2 homozygous non-reference variants, and 2 homozygous reference variants (and they matched my higher coverage WGS data).

I am only 34, but I definitely don't have hair that looks like the Google images for this disease.  So, I don't know the expected age of onset, but I think I might never get this condition (even though Nebula says that I am at the 100th percentile).

As you can see from this link, the percentile above was calculated using 15 SNPs.  12 of those SNPs had at least 1 variant from the reference genome and 11/12 of those variants matched by Vertias WGS variants.  The discordant variant was rs6910071, which was homozygous for the variant allele on the Nebula lcWGS imputed variants.  So, this could have been consistent with 90% overall accuracy, but I didn't have coverage for either dataset at this position (so, this isn't the same as using a gVCF to make a homozygous reference genotype call).

I have osteoarthritis in my lower back, but I don't believe that I (currently) have rheumatoid arthritis.

While I could believe that I am at increased risk, it is important to note that the summary only describes 4% of variance in disease risk.  I think this should be described for all of the reports, to give a sense of the predictive power (along with other statistics).

I also noticed most of the variants were not present in ClinVar (when I was using dbSNP to check the hg19 genome coordinates and reference allele).

As you can see from this link, the percentile above was calculated using 56 SNPs.

I have blood test results uploaded on my PatientsLikeMe profile, and I thought that I had a normal CRP result.  However, it appears that I might not have remembered that correctly, and I might need to wait until my next checkup to see if I can test my CRP level.  However, this is something where I think it would be relatively easy to show if being at the 99% percentile substantially affects your observed CRP levels (or whether there are limits to what this score represents).

As you can see from this link, the percentile above was calculated using 6 SNPs.  Nebula predicted that I had a homozygous variant for 1 SNP and heterozygous variant for 1 SNP (and reference genotypes for the other 4 variants).  However, all of these variants looked OK in my higher coverage Veritas WGS data.

I believe that I have been previously reported to be at higher risk for restless leg syndrome, but that might have actually been for deep vein thrombosis / venous thromboembolism in the earlier 23andMe reports (before the FDA required approval for a more select set of results).  I do sometimes have difficult sitting perfectly still at night.  However, this does't happen all of the time, and I have never been diagnosed by a doctor for having this condition.  So, I would currently lean towards saying that I don't have restless leg syndrome.

For the 2 sets of SNPs that I checked (for alopecia areata and restless leg syndrome), I also visualized my Genos Exome alignment.  However, most of the variants were not covered by sequencing of coding regions.

In general, the journals where these results are published may make some readers think the results are useful.  However, being able to publish a result in a prestigious journal doesn't mean the associations are predictive enough to be clinically meaningful.  Also, even with a more subtle association, being in a prestigious journal doesn't necessarily mean the result can be reproduced.  For example, there are retractions in high impact journals (you can see some in this blog post), and there are objectively wrong conclusions in papers that haven't been retracted (and the science-wide error rate is mentioned in this blog post).  I don't want to cause unnecessarily alarm, but I think it is important to emphasize that time and large sample sizes (and independent validation) are needed to become comfortable with using genetic results to guide your medical treatment.

Nebula does provide a warning: "Disclaimer: Nebula Library is for research, information, and educational use only. This information is not medical advice, nor is it intended to be used for any diagnostic purpose. Please seek the assistance of a health care provider with any questions regarding your health."  However, I think this is easy to miss and the importance may not be fully understood among all customers.

Additional Note #1: I tried to provide a review for Nebula on Trustpilot, but I have encountered some difficulties.  You can see a screenshot of the current review (that was not accepted) here.  While I still haven't gotten a response for the last attempt to submit a review and get an explanation of what I need to change.  While I am not certain if this is the cause, the link that I received to submit a review initially created a review under another name.  To be fair, Nebula did pay me the $10 Amazon gift card (even when I provided a screenshot of a 2-star review, when most are 4- or 5-star reviews), and I hope that this review can eventually be posted (in which case, I will provide a link to that review, instead of this longer explanation).

Additional Note #2: You can see my report to FDA MedWatch (MW5093887) in MAUDE here.  I received an acknowledgement via mail for another report, but I just looked for this report after waiting a while.

Update Log:

8/4/2019 - public post date
8/5/2019 - add update about Nebula research library
8/6/2019 - minor changes
8/10/2019 - add link for 23andMe SNP chip versus WGS concordance
8/14/2019 - minor changes
8/15/2019 - minor changes; add box around APOE variant
8/16/2019 - minor changes
8/17/2019 - revise title (to better emphasize importance, but also presentation of just my own data)
8/16/2019 - minor changes
11/26/2019 - add arrows to Nebula samples in IGV screenshots
1/26/2020 - add screenshot of Trustpilot review
2/24/2020 - mention that Nebula is currently offering higher coverage sequencing
3/19/2020 - add concerns about Nebula Library (to match what I described in my MedWatch report)
3/20/2020 - add links to SNP details for selected library results
3/22/2020 - add link to other PRS post
4/18/2020 - add notes about checking RA post
4/24/2020 - add notes about FDA MedWatch submission
4/26/2020 - minor changes
7/6/2020 - minor changes "Polygenic Risk Score" label, for the last part of the blog post.