Workflow of Thompson et al. (2023)

Individual whole blood from US (RNA)

Author

Harmon Bhasin

Published

September 12, 2024

This is another potential study of this series. In this post, I analyze Thompson 2023, a dataset with 417 samples of whole blood from healthy and SARS-CoV-2 infected individuals in the US.

I’d like to thank Lenni for giving me feedback on the notebook and Will for providing me with a boilerplate rmarkdown file. This notebook uses the MGS Workflow v2.2.1 (note that this is old).

1 The raw data

1.1 About

This dataset is composed of 417 samples of whole blood from the US. In total, there were 353 SARS-CoV-2 infected individuals and 64 healthy individuals that contributed to this study. For each sample, they performed RNA-sequencing, with Illumina NovaSeq 6000, producing 2x100 bp reads.

Code 
low_sample_number <- c('SRR21924256')
# Import libraries and extract metadata from sample names
libraries <- read_csv(libraries_path, show_col_types = FALSE) %>% filter(!sample %in% low_sample_number)
#meta_data <- read_csv('/Users/harmonbhasin/work/securebio/nao-harmon/thijssen2023/preprocessing/SraRunTable.txt') %>%
#  select('Library Name', Run) %>% rename(sample = Run, library = 'Library Name')
#libraries <- libraries %>%
#  left_join(meta_data, by = 'sample') %>%
#  select(sample, library=library.y) %>% 
#  filter(library != 'PCL21') %>%
#  mutate(library = factor(library, levels = c(paste0('PCL', 1:20))))
#libraries

1.1.1 Sample + library preparation

This is paraphrased from the paper:

Patients were seen between April and June 2020 for COVID-19. The whole blood used for RNA-seq was collected in Tempus RNA Blood Tubes (Thermo Fisher Scientific, 4342792). As soon as possible after blood collection, tubes were shaken and stored at −80 °C.

RNA extraction, library preparation and sequencing were performed as described previously 28. In brief, frozen blood samples were thawed, and total RNA was extracted using a modification of the MagMax protocol for Stabilized Blood Tubes RNA Isolation Kit. Samples yielding sufficient RNA (>50 ng) were barcoded and prepared for pooled whole transcriptome sequencing using the TruSeq Stranded Total RNA Library Prep Gold, which is designed to remove ribosomal, globin and mitochondrial RNA. Libraries were amplified with 15 cycles of PCR, pooled and sequenced on a NovaSeq 6000 (Illumina) using Sprime flow cells with 100-bp paired-end reads, targeting a mean of 50 million read pairs per sample.

1.2 Quality control metrics

In total, these 417 samples contained 27.7B read pairs. The samples had 30.7M - 133.9M (mean ~66.4M) read pairs each. The number of reads looks pretty good, although a few samples have a much lower number of reads. The total number of base pairs also looks reasonable, matching the trends of the number of reads. The duplication rate is very high. Adapter content is variable, but most samples have a low adapter content (below 5%). As we’d expect, we see higher quality Q scores near the beginning of the read and a gradual decline towards the end of the read, however all positions seem to have a Q score of 35 which means that our reads are ~99.97% accurate. When looking at the Q score over all sequences, we can see a sharp peak around 35, which corresponds to the previous plot, indicating high read quality.

Code 
# Prepare data
basic_stats_raw_metrics <- basic_stats_raw %>%
  select(library,
         `# Read pairs` = n_read_pairs,
         `Total base pairs\n(approx)` = n_bases_approx,
         `% Duplicates\n(FASTQC)` = percent_duplicates) %>%
  pivot_longer(-library, names_to = "metric", values_to = "value") %>%
  mutate(metric = fct_inorder(metric)) 

# Set up plot templates

g_basic <- ggplot(basic_stats_raw_metrics, aes(x = library, y = value)) +
  geom_col(position = "dodge") +
  scale_x_discrete() +
  scale_y_continuous(expand = c(0, 0)) +
  expand_limits(y = c(0, 100)) +
  facet_grid(metric ~ ., scales = "free", space = "free_x", switch = "y") +
  theme_kit + theme(
    axis.title.y = element_blank(),
    strip.text.y = element_text(face = "plain")
  )
g_basic

Code 
# Set up plotting templates
g_qual_raw <- ggplot(mapping=aes(linetype=read_pair, group=interaction(sample,read_pair))) + 
  scale_linetype_discrete(name = "Read Pair") +
  guides(color=guide_legend(nrow=2,byrow=TRUE),
         linetype = guide_legend(nrow=2,byrow=TRUE)) +
  theme_base

# Visualize adapters
g_adapters_raw <- g_qual_raw + 
  geom_line(aes(x=position, y=pc_adapters), data=adapter_stats_raw) +
  scale_y_continuous(name="% Adapters", limits=c(0,NA),
                     breaks = seq(0,10,1), expand=c(0,0)) +
  scale_x_continuous(name="Position", limits=c(0,NA),
                     breaks=seq(0,500,20), expand=c(0,0)) +
  facet_grid(.~adapter)
g_adapters_raw

Code 
# Visualize quality
g_quality_base_raw <- g_qual_raw +
  geom_hline(yintercept=25, linetype="dashed", color="red") +
  geom_hline(yintercept=30, linetype="dashed", color="red") +
  geom_line(aes(x=position, y=mean_phred_score), data=quality_base_stats_raw) +
  scale_y_continuous(name="Mean Phred score", expand=c(0,0), limits=c(10,45)) +
  scale_x_continuous(name="Position", limits=c(0,NA),
                     breaks=seq(0,500,20), expand=c(0,0))
g_quality_base_raw

Code 
g_quality_seq_raw <- g_qual_raw +
  geom_vline(xintercept=25, linetype="dashed", color="red") +
  geom_vline(xintercept=30, linetype="dashed", color="red") +
  geom_line(aes(x=mean_phred_score, y=n_sequences), data=quality_seq_stats_raw) +
  scale_x_continuous(name="Mean Phred score", expand=c(0,0)) +
  scale_y_continuous(name="# Sequences", expand=c(0,0))
g_quality_seq_raw

2 Preprocessing

2.1 Summary

The average fraction of reads at each stage in the preprocessing pipeline is shown in the following table. Adapter trimming & filtering doesn’t lose too many reads. Deduplication loses us about 55% of reads on average, which is good considering the high duplication content seen earlier. Ribodepletion only removes about 2% of reads in total, which makes sense as they did ribosomal depletion during sample preparation.

Code 
# TODO: Group by pool size as well
# Count read losses
n_reads_rel <- basic_stats %>% 
  select(sample, stage, percent_duplicates, n_read_pairs) %>%
  group_by(sample) %>% 
  arrange(sample, stage) %>%
  mutate(p_reads_retained = n_read_pairs / lag(n_read_pairs),
         p_reads_lost = 1 - p_reads_retained,
         p_reads_retained_abs = n_read_pairs / n_read_pairs[1],
         p_reads_lost_abs = 1-p_reads_retained_abs,
         p_reads_lost_abs_marginal = p_reads_lost_abs - lag(p_reads_lost_abs))
n_reads_rel_display <- n_reads_rel %>% 
  rename(Stage=stage) %>% 
  group_by(Stage) %>% 
  summarize(`% Total Reads Lost (Cumulative)` = paste0(round(min(p_reads_lost_abs*100),1), "-", round(max(p_reads_lost_abs*100),1), " (mean ", round(mean(p_reads_lost_abs*100),1), ")"),
            `% Total Reads Lost (Marginal)` = paste0(round(min(p_reads_lost_abs_marginal*100),1), "-", round(max(p_reads_lost_abs_marginal*100),1), " (mean ", round(mean(p_reads_lost_abs_marginal*100),1), ")"), .groups="drop") %>% 
  filter(Stage != "raw_concat") %>%
  mutate(Stage = Stage %>% as.numeric %>% factor(labels=c("Trimming & filtering", "Deduplication", "Initial ribodepletion", "Secondary ribodepletion")))
n_reads_rel_display

2.2 Quality control metrics

Cleaning seems to get rid of the library contamination and the polyg contamination. Polya contamination still is present, but I don’t think we’ve found a solution to remove it as yet. Q score remain the same during read cleaning when looking at the positions, with the end of the read actually improving in score. Q scores across all sequences look pretty much the same throughout cleaning.

Code 
g_qual <- ggplot(mapping=aes(linetype=read_pair, group=interaction(sample,read_pair))) + 
  scale_linetype_discrete(name = "Read Pair") +
  guides(color=guide_legend(nrow=2,byrow=TRUE),
         linetype = guide_legend(nrow=2,byrow=TRUE)) +
  theme_base

# Visualize adapters
g_adapters <- g_qual + 
  geom_line(aes(x=position, y=pc_adapters), data=adapter_stats) +
  scale_y_continuous(name="% Adapters", limits=c(0,20),
                     breaks = seq(0,50,10), expand=c(0,0)) +
  scale_x_continuous(name="Position", limits=c(0,NA),
                     breaks=seq(0,140,20), expand=c(0,0)) +
  facet_grid(stage~adapter)
g_adapters

Code 
# Visualize quality
g_quality_base <- g_qual +
  geom_hline(yintercept=25, linetype="dashed", color="red") +
  geom_hline(yintercept=30, linetype="dashed", color="red") +
  geom_line(aes(x=position, y=mean_phred_score), data=quality_base_stats) +
  scale_y_continuous(name="Mean Phred score", expand=c(0,0), limits=c(10,45)) +
  scale_x_continuous(name="Position", limits=c(0,NA),
                     breaks=seq(0,140,20), expand=c(0,0)) +
  facet_grid(stage~.)
g_quality_base

Code 
g_quality_seq <- g_qual +
  geom_vline(xintercept=25, linetype="dashed", color="red") +
  geom_vline(xintercept=30, linetype="dashed", color="red") +
  geom_line(aes(x=mean_phred_score, y=n_sequences), data=quality_seq_stats) +
  scale_x_continuous(name="Mean Phred score", expand=c(0,0)) +
  scale_y_continuous(name="# Sequences", expand=c(0,0)) +
  facet_grid(stage~., scales = "free_y")
g_quality_seq

These plots below show the trends from above in each sample. Deduplication seems to drop by a lot post deduplication,which is a good thing given it was so high earlier. Mean read length remains relatively constant throuhgout. Some samples still seem to have high ribosomal content, however I’m not too worried about this.

Code 
g_stage_trace <- ggplot(basic_stats, aes(x=stage, group=sample)) +
  theme_kit

# Plot reads over preprocessing
g_reads_stages <- g_stage_trace +
  geom_line(aes(y=n_read_pairs)) +
  scale_y_continuous("# Read pairs", expand=c(0,0), limits=c(0,NA))
g_reads_stages

Code 
# Plot relative read losses during preprocessing
g_reads_rel <- ggplot(n_reads_rel, 
                      aes(x=stage, group=sample)) +
  geom_line(aes(y=p_reads_lost_abs_marginal)) +
  scale_y_continuous("% Total Reads Lost", expand=c(0,0), 
                     labels = function(x) x*100) +
  theme_kit
g_reads_rel

Code 
stage_dup <- basic_stats %>% group_by(stage) %>% 
  summarize(dmin = min(percent_duplicates), dmax=max(percent_duplicates),
            dmean=mean(percent_duplicates), .groups = "drop")

g_dup_stages <- g_stage_trace +
  geom_line(aes(y=percent_duplicates)) +
  scale_y_continuous("% Duplicates", limits=c(0,NA), expand=c(0,0))
g_dup_stages

Code 
g_readlen_stages <- g_stage_trace + geom_line(aes(y=mean_seq_len)) +
  scale_y_continuous("Mean read length (nt)", expand=c(0,0), limits=c(0,NA))
g_readlen_stages

Code 
# Calculate reads lost during ribodepletion (approximation for % ribosomal reads)
reads_ribo <- n_reads_rel %>% 
  filter(stage %in% c("dedup", "ribo_secondary")) %>% 
  group_by(sample) %>% 
  summarize(p_reads_ribo=1-n_read_pairs[2]/n_read_pairs[1], .groups = "drop") %>%
  inner_join(libraries, by = 'sample')
reads_ribo_summ <- reads_ribo %>%
  group_by(sample) %>%
  summarize(min=min(p_reads_ribo), max=max(p_reads_ribo),
            mean=mean(p_reads_ribo), .groups = "drop") %>%
  inner_join(libraries, by = 'sample')
g_reads_ribo <- ggplot(reads_ribo, 
                       aes(x=library, y=p_reads_ribo)) +
  geom_point() + 
  scale_y_continuous(name="Approx % ribosomal reads", limits=c(0,1),
                     breaks=seq(0,1,0.2), expand=c(0,0), labels = function(y) y*100)+
  theme_kit
g_reads_ribo

3 Taxonomic composition

3.1 High-level composition

To assess the high-level composition of the reads, I ran the ribodepleted files through Kraken2 and summarized the results with Bracken.

The groups listed below were created by Will:

  • Filtered (removed during cleaning)
  • Duplicate (removed during deduplication)
  • Ribosomal (removed during ribodepletion)
  • Unassigned (non-ribosomal reads that were not assigned to any taxon by Kraken/Bracken)
  • Bacterial (non-ribosomal reads assigned to the Bacteria domain by Kraken/Bracken)
  • Archaeal (non-ribosomal reads assigned to the Archaea domain by Kraken/Bracken)
  • Viral (non-ribosomal reads assigned to the Viruses domain by Kraken/Bracken)
  • Human (non-ribosomal reads assigned to the Eukarya domain by Kraken/Bracken)
Code 
# Prepare plotting templates
g_comp_base <- ggplot(mapping=aes(x=library, y=p_reads, fill=classification)) +
  scale_x_discrete(name="Plasma pool") +
  theme_kit + 
  theme(plot.title = element_text(hjust=0, face="plain", size=rel(1.5)))
scale_y_pc_reads <- purrr::partial(scale_y_continuous, name = "% Reads",
                                   expand = c(0,0), labels = function(y) y*100)
geom_comp <- purrr::partial(geom_col, position = "stack", width = 1)

# Define a color palette for the classification
classification_colors <- brewer.pal(8, "Accent")
names(classification_colors) <- c("Filtered", "Duplicate", "Ribosomal", "Unassigned", "Bacterial", "Archaeal", "Viral", "Human")
scale_fill_classification <- function() {
  scale_fill_manual(values = classification_colors, name = "Classification")
}

# Plot overall composition
g_comp <- g_comp_base + geom_comp(data = comp) +
  scale_y_pc_reads(limits = c(0,1.01), breaks = seq(0,1,0.2)) +
  scale_fill_classification() + 
  ggtitle("Read composition (all reads, all groups)")
g_comp

Code 
g_comp_assigned <- g_comp_base + 
  geom_comp(data = comp_assigned) +
  scale_y_pc_reads(limits = c(0,1.01), breaks = seq(0,1,0.2)) +
  scale_fill_classification() + 
  ggtitle("Read composition (assigned reads, all groups)")
g_comp_assigned

3.2 Total viral content

Total viral fraction average \(1.12 \times 10^{-6}\) across samples. As a fraction of assigned (rather than total) reads, this jumped to \(2.96 \times 10^{-6}\):

3.3 Taxonomic composition of viruses

There seems to be no dominant virus in the samples. Specifically, this seems to be the dataset with the most viral diversity (when compared to the previous plasma datasets). The threshold for the label “other” are the set of families that make up less than 30% composition in all samples.

Code 
major_threshold <- 0.3

# Identify major viral families
viral_families_major_tab <- viral_families %>% 
  group_by(name, taxid) %>%
  summarize(p_reads_viral_max = max(p_reads_viral), .groups="drop") %>%
  filter(p_reads_viral_max >= major_threshold)
viral_families_major_list <- viral_families_major_tab %>% pull(name)
viral_families_major <- viral_families %>% 
  filter(name %in% viral_families_major_list) %>%
  select(name, taxid, sample, p_reads_viral)
viral_families_minor <- viral_families_major %>% 
  group_by(sample) %>%
  summarize(p_reads_viral_major = sum(p_reads_viral), .groups = "drop") %>%
  mutate(name = "Other", taxid=NA, p_reads_viral = 1-p_reads_viral_major) %>%
  select(name, taxid, sample, p_reads_viral)
viral_families_display <- viral_families_major %>% 
  bind_rows(viral_families_minor) %>%
  arrange(desc(p_reads_viral)) %>% 
  mutate(name = factor(name, levels=c(viral_families_major_list, "Other"))) %>%
  rename(p_reads = p_reads_viral, classification=name) %>%
  inner_join(libraries, by='sample')


pal <- c(brewer.pal(8, 'Dark2'),brewer.pal(8, 'Accent'), brewer.pal(8, 'Set1'), brewer.pal(8, 'Pastel1'), brewer.pal(8, 'Pastel2'))

n_classifications <- length(unique(viral_families_display$classification)) - 1

custom_palette_with_black <- c(
  pal[1:n_classifications],
  "black"
)

# Plot
g_families <- g_comp_base + 
  geom_comp(data=viral_families_display) +
  scale_y_continuous(name="% Viral Reads", limits=c(0,1.01), 
                     breaks = seq(0,1,0.2),
                     expand=c(0,0), labels = function(y) y*100) +
  scale_fill_manual(values = custom_palette_with_black)
g_families

4 Human-infecting virus reads

4.1 Overall relative abundance

I calculated the relative abundance of human-infecting viruses in two ways:

  • First, as the total number of deduplicated human-virus reads in each sample, divided by the number of raw reads (“All reads”). This results in a viral fraction of \(1.5 \times 10^{-6}\) across all samples

  • Second, as a fraction of preprocessed (cleaned, deduplicated, computationally ribodepleted) reads (“Preprocessed reads”). This results in a viral fraction of \(3.2 \times 10^{-6}\) across all samples.

4.2 Overall taxonomy and composition

Composition of HV reads was changed from when looking at all viral reads. First, only 15% of the individuals in this dataset had human infecting viruses detectable in their blood. Second, individuals seem to be dominated by a single viral family. The threshold for the label “other” are the set of families that make up less than 5% composition in all samples.

Code 
threshold_major_family <- 0.05

# Count reads for each human-viral family
hv_family_counts <- hv_reads_family %>% 
  group_by(sample, name, taxid) %>%
  count(name = "n_reads_hv") %>%
  group_by(sample) %>%
  mutate(p_reads_hv = n_reads_hv/sum(n_reads_hv))

# Identify high-ranking families and group others
hv_family_major_tab <- hv_family_counts %>% group_by(name) %>% 
  filter(p_reads_hv == max(p_reads_hv)) %>% filter(row_number() == 1) %>%
  arrange(desc(p_reads_hv)) %>% filter(p_reads_hv > threshold_major_family)
hv_family_counts_major <- hv_family_counts %>%
  mutate(name_display = ifelse(name %in% hv_family_major_tab$name, name, "Other")) %>%
  group_by(sample, name_display) %>%
  summarize(n_reads_hv = sum(n_reads_hv), p_reads_hv = sum(p_reads_hv), 
            .groups="drop") %>%
  mutate(name_display = factor(name_display, 
                               levels = c(hv_family_major_tab$name, "Other")))
hv_family_counts_display <- hv_family_counts_major %>%
  rename(p_reads = p_reads_hv, classification = name_display) %>%
  inner_join(libraries, by = 'sample')

# Plot
g_hv_family <- g_comp_base + 
  geom_col(data=hv_family_counts_display, position = "stack", width=1) +
  scale_y_continuous(name="% HV Reads", limits=c(0,1.01), 
                     breaks = seq(0,1,0.2),
                     expand=c(0,0), labels = function(y) y*100) +
  scale_fill_brewer(palette = 'Accent', name = "Viral family") +
  labs(title="Family composition of human-viral reads") +
  guides(fill=guide_legend(ncol=4)) +
  theme(plot.title = element_text(size=rel(1.4), hjust=0, face="plain"))
g_hv_family

Code 
# Get most prominent families for text
hv_family_collate <- hv_family_counts %>%
  group_by(name, taxid) %>% 
  summarize(n_reads_tot = sum(n_reads_hv),
            p_reads_max = max(p_reads_hv), .groups="drop") %>% 
  arrange(desc(n_reads_tot))
hv_family_collate

4.3 Analyzing specific families

We now investigate the composition of specific families that had more than 5 viral reads. In investigating individual viral families, to avoid distortions from a few rare reads, I restricted myself to samples where that family made up at least 1% of human-viral reads:

4.3.1 Flaviviridae (Number of reads: 25,997)

Code 
# Flaviviridae
plot_viral_family_histogram(taxid_chosen=11050)

Code 
# Flaviviridae
plot_viral_family_composition(taxid_chosen=11050, threshold_major_species = 0.01)

4.3.2 Anelloviridae (Number of reads: 10,103)

Code 
# Anelloviridae
plot_viral_family_histogram(taxid_chosen=687329)

Code 
# Anelloviridae
plot_viral_family_composition(taxid_chosen=687329, threshold_major_species = 0.01)

4.3.3 Coronaviridae (Number of reads: 1,104)

Code 
# Coronaviridae
plot_viral_family_histogram(taxid_chosen=11118)

Code 
# Coronaviridae
plot_viral_family_composition(taxid_chosen=11118, threshold_major_species = 0.01)

4.3.4 Orthoherpesviridae (Number of reads: 534)

Code 
# Orthoherpesviridae
plot_viral_family_histogram(taxid_chosen=3044472)

Code 
# Orthoherpesviridae
plot_viral_family_composition(taxid_chosen=3044472, threshold_major_species = 0.01)

4.3.5 Parvoviridae (Number of reads: 335)

Code 
# Parvoviridae
plot_viral_family_histogram(taxid_chosen=10780)

Code 
# Parvoviridae
plot_viral_family_composition(taxid_chosen=10780, threshold_major_species = 0.01)

4.3.6 Retroviridae (Number of reads: 224)

Code 
# Retroviridae
plot_viral_family_histogram(taxid_chosen=11632)

Code 
# Retroviridae
plot_viral_family_composition(taxid_chosen=11632, threshold_major_species = 0.01)

4.4 Relative abundance of pathogenic viruses of interest

Each dot represents a sample, colored by viral family. The x-axis shows the relative abundance of human-infecting viruses, and the y-axis shows the species.

Code 
species_family <- result %>% select(species, family) %>% rename('name' = 'species')

play <- hv_species_counts %>% 
  ungroup() %>%
  inner_join(libraries, by = 'sample') %>%
  inner_join(species_family, by = 'name') %>%
  mutate(name = ifelse(name == "Severe acute respiratory syndrome-related coronavirus", "SARS-related coronavirus", name))

adjusted_play <- play %>% 
  group_by(name) %>%
  mutate(virus_prevalence_num = n_distinct(sample)/n_distinct(libraries),
         total_reads_hv = sum(n_reads_hv)) %>%
  ungroup() %>%
  mutate(name = fct_reorder(name, virus_prevalence_num, .desc=TRUE)) %>% 
  select(name, ra_reads_hv, family, virus_prevalence_num, total_reads_hv, n_reads_hv)

pal <- c(brewer.pal(8, 'Dark2'),brewer.pal(8, 'Accent'))

#, labels = label_log(digits=2)

ra_dot <- ggplot(adjusted_play, aes(x = ra_reads_hv, y=name)) +
  geom_quasirandom(orientation = 'y', aes(color = family)) +
  scale_color_manual(values = pal) + 
    scale_x_log10(
    "Relative abundance human virus reads",
    labels=label_log(digits=3)
  ) +
  labs(y =  "",
       color = 'Viral family') + 
  guides(color = guide_legend(ncol=2)) +
  theme_light() + 
    theme(
    axis.text.y = element_text(size = 10),
    axis.text.x = element_text(size = 12),
    axis.ticks.y = element_blank(),axis.line = element_line(colour = "black"),
    axis.title.x = element_text(size = 14),    
    legend.text = element_text(size = 10),
    legend.title = element_text(size = 10, face="bold"),
    #legend.position = 'bottom',
    legend.position = c(1, 1),  # Move legend to top right
    legend.justification = c(1, 1),  # Align legend to top right
    panel.grid.minor = element_blank(),
    panel.border = element_blank(),
    panel.background = element_blank())
ra_dot

I can then take all of the viruses that we found and look up what they’re responsible for and whether they’re dangerous.

Virus Name Common Name Pathogenic Potential
Human immunodeficiency virus 1 HIV-1 Very high, causes AIDS if left untreated, leading to severe immunodeficiency
SARS-related coronavirus SARS-CoV High, causes severe acute respiratory syndrome with potentially fatal outcomes
Hepacivirus hominis Hepatitis C Virus High, causes chronic liver disease and can lead to cirrhosis and liver cancer
Cytomegalovirus humanbeta5 Human Cytomegalovirus Moderate to high, can cause severe disease in immunocompromised individuals and congenital infections
Lymphocryptovirus humangamma4 Epstein-Barr Virus Moderate, causes infectious mononucleosis and associated with certain cancers
Erythroparvovirus primate1 Human Parvovirus B19 Moderate, causes fifth disease and can be serious in certain populations (pregnant women, immunocompromised)
Roseolovirus humanbeta6a Human Herpesvirus 6A Low to moderate, can cause roseola in children and occasionally more severe complications
Roseolovirus humanbeta6b Human Herpesvirus 6B Low to moderate, can cause roseola in children and occasionally more severe complications
Alphapapillomavirus 2 Human Papillomavirus 2 Low to moderate, can cause common warts and rarely progress to cancer
Betapapillomavirus 5 Human Papillomavirus 5 Low to moderate, associated with skin lesions in immunosuppressed individuals
Betapapillomavirus 2 Human Papillomavirus 9 Low to moderate, associated with cutaneous lesions and rarely with non-melanoma skin cancer
Primate T-lymphotropic virus 2 Human T-cell Lymphotropic Virus 2 Low, rarely associated with neurological disorders
Nupapillomavirus 1 Human Papillomavirus 41 Low, can cause benign cutaneous lesions
Alphapolyomavirus quintihominis Human Polyomavirus 5 Low, not known to cause significant disease in healthy individuals
Pegivirus hominis Human Pegivirus Very low, not known to cause disease in humans
Pegivirus columbiaense Pegivirus C Very low, not known to cause disease in humans

This was interesting to see. This is the first study where we’ve seen signficant numbers of human-infecting viruses in whole blood.

4.5 Relative abundance assuming perfect human read removal

Assuming we’re able to perfectly remove all human reads, the average relative abundance of known human infecting virus is \(1.91 \times 10^{-5}\).

5 Conclusion

There were some interesting takeways from this analysis:

  1. Whole blood metagenomics was able to detect HIV-1, SARS-CoV-2, HCV, CMV, EBV, HTLV-2, and other herpesviruses. This is important as some of these are latent viruses.
  2. Far fewer individuals tend to have detectable levels of human infecting viruses, when compared to the plasma datasets.
  3. Important viruses can be detected from RNA-seq data, even if the study wasn’t explictly designed to look at them.

6 Appendix

6.1 Human-infecting virus families, genera, and species

To get a good overview of families, genera, and species, we can look at a Sankey plot where the magnitude of relative abundance, averaged over all samples, is shown in parentheses.

Code 
# Function to create links
create_links <- function(data) {
  family_to_genus <- data %>%
    filter(!is.na(genus)) %>%
    group_by(family, genus) %>%
    summarise(value = n(), .groups = "drop") %>%
    mutate(source = family, target = genus)
  
  genus_to_species <- data %>%
    group_by(genus, species) %>%
    summarise(value = n(), .groups = "drop") %>%
    mutate(source = genus, target = species)

  family_to_species <- data %>%
    filter(is.na(genus)) %>%
    group_by(family, species) %>%
    summarise(value = n(), .groups = "drop") %>%
    mutate(source = family, target = species)

  bind_rows(family_to_genus, genus_to_species, family_to_species) %>%
    filter(!is.na(source))
}

# Function to create nodes
create_nodes <- function(links) {
  data.frame(
    name = c(links$source, links$target) %>% unique()
  )
}

# Function to prepare data for Sankey diagram
prepare_sankey_data <- function(links, nodes) {
  links$IDsource <- match(links$source, nodes$name) - 1
  links$IDtarget <- match(links$target, nodes$name) - 1
  list(links = links, nodes = nodes)
}

# Function to create Sankey plot
create_sankey_plot <- function(sankey_data) {
  sankeyNetwork(
    Links = sankey_data$links, 
    Nodes = sankey_data$nodes,
    Source = "IDsource", 
    Target = "IDtarget",
    Value = "value", 
    NodeID = "name",
    sinksRight = TRUE,
    nodeWidth = 25,
    fontSize = 14,
  )
}

save_sankey_as_png <- function(sankey_plot, width = 1000, height = 800) {
  # Save the plot as an HTML file
  saveWidget(sankey_plot, sprintf('%s/sankey.html',data_dir))
}

format_scientific <- function(x, digits=2) {
  sapply(x, function(val) {
    if (is.na(val) || abs(val) < 1e-15) {
      return("0")
    } else {
      exponent <- floor(log10(abs(val)))
      coef <- val / 10^exponent
      #return(sprintf("%.1f × 10^%d", round(coef, digits), exponent))
      # May or may not be smart, just keeping magnitude
      return(sprintf("10^%d", exponent))
    }
  })
}

data <- result %>% 
  mutate(across(c(genus_n_reads_tot, genus_ra_reads_tot), ~replace_na(., 0)),
         genus = ifelse(is.na(genus), "Unknown Genus", genus)) %>%
  mutate(
  species = paste0(species, sprintf(' (%s)', format_scientific(species_ra_reads_tot))),
  genus = paste0(genus, sprintf(' (%s)', format_scientific(genus_ra_reads_tot))),
  family = paste0(family, sprintf(' (%s)', format_scientific(family_ra_reads_tot)))
)
links <- as.data.frame(create_links(data))
nodes <- create_nodes(links)
sankey_data <- prepare_sankey_data(links, nodes)
sankey <- create_sankey_plot(sankey_data)

sankey