Viral proteins in bulk RNA seq data

Open in Google Colab

In this tutorial, we will align RNA sequencing data collected from SARS-CoV2 infected human iPSC derived cardiomyocytes to viral RdRP protein sequences. Let’s see if we can detect SARS-CoV2 RdRP-like sequences as expected. This is a SMART-seq dataset, but the workflow is the same for bulk RNA sequencing data.

Workflow reference:

Laura Luebbert, Delaney K Sullivan, Maria Carilli, Kristján Eldjárn Hjörleifsson, Alexander Viloria Winnett, Tara Chari, Lior Pachter (2023). Efficient and accurate detection of viral sequences at single-cell resolution reveals novel viruses perturbing host gene expression. bioRxiv 2023.12.11.571168; doi: https://doi.org/10.1101/2023.12.11.571168

Data references:

RNA seq data: https://www.cell.com/cell-reports-medicine/pdf/S2666-3791(20)30068-9.pdf

PalmDB viral protein (RdRP) reference database: https://www.nature.com/articles/s41586-021-04332-2

PalmDB is an uncurated database of viral RNA-dependent RNA polymerase (RdRP) sequences, primarily derived from metagenomic sources. Consequently, some entries may originate from non-viral sources or represent host-derived sequences. We provide example code below demonstrating one method for the masking of host sequences and the subsequent extraction and BLAST analysis of the identified reads.

# SRR numbers of the sequencing dataset(s) to analyze
srr_numbers = ["SRR11777734", "SRR11777735", "SRR11777736", "SRR11777737", "SRR11777738", "SRR11777739"]

# Genus and species of the host (separate by "_")
host_species = "homo_sapiens"

# Here, we will only align a subset of [n_seq_to_keep] sequences from each of these SRRs
# to decrease the runtime of this example notebook
# Set n_seq_to_keep=False to align the complete dataset
n_seq_to_keep = 10000000

# Number of threads to use during alignments
threads = 2

# k-mer length: Increase k for increased specificity (at a potential trade-off with sensitivity)
# k should be an odd integer and should be >=31 for translated alignment
# Since k relates to nucleotides and not amino acids, k corresponds to matching ⌊k/3⌋ amino acids (for k=31, this corresponds to ⌊31/3⌋=10 amino acids)
k = 31

Install software

!pip install -q ffq gget kb_python==0.29.1 anndata==0.10.9

     ━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━ 49.6/49.6 kB 3.0 MB/s eta 0:00:00
  Preparing metadata (setup.py) ... done
  Preparing metadata (setup.py) ... done
   ━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━ 36.5/36.5 MB 29.4 MB/s eta 0:00:00
   ━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━ 129.0/129.0 kB 10.3 MB/s eta 0:00:00
   ━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━ 43.1/43.1 MB 12.2 MB/s eta 0:00:00
   ━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━ 50.4/50.4 kB 3.5 MB/s eta 0:00:00
   ━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━ 3.2/3.2 MB 82.8 MB/s eta 0:00:00
   ━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━ 34.4/34.4 MB 45.3 MB/s eta 0:00:00
   ━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━ 50.9/50.9 MB 8.8 MB/s eta 0:00:00
   ━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━ 2.1/2.1 MB 74.8 MB/s eta 0:00:00
   ━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━ 1.6/1.6 MB 62.9 MB/s eta 0:00:00
   ━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━ 56.9/56.9 kB 4.6 MB/s eta 0:00:00
   ━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━ 22.0/22.0 MB 76.2 MB/s eta 0:00:00
   ━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━ 88.8/88.8 kB 7.2 MB/s eta 0:00:00
   ━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━ 40.6/40.6 kB 2.9 MB/s eta 0:00:00
   ━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━ 83.6/83.6 kB 6.8 MB/s eta 0:00:00
  Building wheel for loompy (setup.py) ... done
  Building wheel for session-info (setup.py) ... done

import anndata
import gget
import glob
import gzip
import json
import matplotlib.colors
import matplotlib.pyplot as plt
import numpy as np
import pandas as pd
import platform
from Bio import SeqIO

%config InlineBackend.figure_format='retina'

# This step will not be necessary after the next kb release
if k<32:
  # Install kallisto v0.50.0 from source
  # This version performs translated search faster while providing the same results as later versions, but it does not support k>31
  # (the features to increase runtime and support k>31 will be incorporated into the next version of kb, so this step will not be necessary after the next release)
  !git clone https://github.com/pachterlab/kallisto.git --branch v0.50.0
  !cd kallisto && mkdir build && cd build && cmake .. && make
  kallisto = "kallisto/build/src/kallisto"
else:
  # Find the kallisto binary installed automatically with kb (or simply remove the --kallisto arguments below)
  kallisto = glob.glob(f"/usr/local/lib/python*/dist-packages/kb_python/bins/{platform.system().lower()}/kallisto/kallisto_k64")[0]

Download the RNA sequencing data

Use ffq to download the sequencing data from the SRR numbers:

# Get ftp download links for raw data with ffq and store results in json file
srr_numbers = " ".join(srr_numbers)
!ffq $srr_numbers \
    --ftp \
    -o ffq.json

# Load ffq output
f = open("ffq.json")
data_json = json.load(f)
f.close()

# Download raw data using FTP links fetched by ffq
for dataset in data_json:
    url = dataset["url"]
    !curl -O $url

[2025-01-11 14:04:53,185]    INFO Parsing run SRR11777734
[2025-01-11 14:04:55,640]    INFO Parsing run SRR11777735
[2025-01-11 14:05:00,165]    INFO Parsing run SRR11777736
[2025-01-11 14:05:02,303]    INFO Parsing run SRR11777737
[2025-01-11 14:05:04,384]    INFO Parsing run SRR11777738
[2025-01-11 14:05:06,424]    INFO Parsing run SRR11777739
  % Total    % Received % Xferd  Average Speed   Time    Time     Time  Current
                                 Dload  Upload   Total   Spent    Left  Speed
100 2962M  100 2962M    0     0  29.6M      0  0:01:39  0:01:39 --:--:-- 30.5M
  % Total    % Received % Xferd  Average Speed   Time    Time     Time  Current
                                 Dload  Upload   Total   Spent    Left  Speed
100 2634M  100 2634M    0     0  13.5M      0  0:03:13  0:03:13 --:--:-- 12.3M
  % Total    % Received % Xferd  Average Speed   Time    Time     Time  Current
                                 Dload  Upload   Total   Spent    Left  Speed
100 2484M  100 2484M    0     0  12.6M      0  0:03:16  0:03:16 --:--:-- 15.3M
  % Total    % Received % Xferd  Average Speed   Time    Time     Time  Current
                                 Dload  Upload   Total   Spent    Left  Speed
100 3018M  100 3018M    0     0  28.6M      0  0:01:45  0:01:45 --:--:-- 30.1M
  % Total    % Received % Xferd  Average Speed   Time    Time     Time  Current
                                 Dload  Upload   Total   Spent    Left  Speed
100 3737M  100 3737M    0     0  29.9M      0  0:02:04  0:02:04 --:--:-- 30.5M
  % Total    % Received % Xferd  Average Speed   Time    Time     Time  Current
                                 Dload  Upload   Total   Spent    Left  Speed
100 3152M  100 3152M    0     0  29.6M      0  0:01:46  0:01:46 --:--:-- 30.5M

Download the viral protein reference sequences

In this case, we are using a version of the PalmDB database v1 that was optimized for viral sequence detection with kallisto as described in this manuscript. These files are stored in the GitHub repository accompanying the manuscript.

# Download the virus ID to taxonomy mapping
!wget https://raw.githubusercontent.com/pachterlab/LSCHWCP_2023/main/PalmDB/ID_to_taxonomy_mapping.csv
# Download the customized transcripts to gene mapping
!wget https://raw.githubusercontent.com/pachterlab/LSCHWCP_2023/main/PalmDB/palmdb_clustered_t2g.txt
# Download the RdRP amino acid sequences
!wget https://raw.githubusercontent.com/pachterlab/LSCHWCP_2023/main/PalmDB/palmdb_rdrp_seqs.fa

--2025-01-11 14:18:55--  https://raw.githubusercontent.com/pachterlab/LSCHWCP_2023/main/PalmDB/ID_to_taxonomy_mapping.csv
Resolving raw.githubusercontent.com (raw.githubusercontent.com)... 185.199.108.133, 185.199.109.133, 185.199.110.133, ...
Connecting to raw.githubusercontent.com (raw.githubusercontent.com)|185.199.108.133|:443... connected.
HTTP request sent, awaiting response... 200 OK
Length: 19705497 (19M) [text/plain]
Saving to: ‘ID_to_taxonomy_mapping.csv’

ID_to_taxonomy_mapp 100%[===================>]  18.79M  --.-KB/s    in 0.1s

2025-01-11 14:18:56 (157 MB/s) - ‘ID_to_taxonomy_mapping.csv’ saved [19705497/19705497]

--2025-01-11 14:18:56--  https://raw.githubusercontent.com/pachterlab/LSCHWCP_2023/main/PalmDB/palmdb_clustered_t2g.txt
Resolving raw.githubusercontent.com (raw.githubusercontent.com)... 185.199.110.133, 185.199.109.133, 185.199.108.133, ...
Connecting to raw.githubusercontent.com (raw.githubusercontent.com)|185.199.110.133|:443... connected.
HTTP request sent, awaiting response... 200 OK
Length: 4561689 (4.3M) [text/plain]
Saving to: ‘palmdb_clustered_t2g.txt’

palmdb_clustered_t2 100%[===================>]   4.35M  --.-KB/s    in 0.07s

2025-01-11 14:18:57 (61.2 MB/s) - ‘palmdb_clustered_t2g.txt’ saved [4561689/4561689]

--2025-01-11 14:18:57--  https://raw.githubusercontent.com/pachterlab/LSCHWCP_2023/main/PalmDB/palmdb_rdrp_seqs.fa
Resolving raw.githubusercontent.com (raw.githubusercontent.com)... 185.199.108.133, 185.199.109.133, 185.199.110.133, ...
Connecting to raw.githubusercontent.com (raw.githubusercontent.com)|185.199.108.133|:443... connected.
HTTP request sent, awaiting response... 200 OK
Length: 35361991 (34M) [text/plain]
Saving to: ‘palmdb_rdrp_seqs.fa’

palmdb_rdrp_seqs.fa 100%[===================>]  33.72M  --.-KB/s    in 0.1s

2025-01-11 14:18:58 (230 MB/s) - ‘palmdb_rdrp_seqs.fa’ saved [35361991/35361991]

Build a reference index from the viral protein sequences and mask host (here, human) sequences

In addition to building a reference index from the protein sequences, we also want to mask host (in this case, human) sequences. For a more conservative host masking workflow that incorporates an additional step for the removal of host sequences, seethis notebook.

The `--aa argument <https://kallisto.readthedocs.io/en/latest/translated/pseudoalignment.html>`__ tells kb that this is an amino acid reference.

The `--d-list argument <https://kallisto.readthedocs.io/en/latest/index/index_generation.html#the-d-list>`__ is the path to the host genome/transcriptome. These sequences will be masked in the index. Here, we are using `gget ref <https://pachterlab.github.io/gget/en/ref.html>`__ to fetch the latest human genome and transcriptome from Ensembl.

We are using --workflow custom here since we do not have a .gtf file for the PalmDB fasta file.

--kallisto provides the path to the kallisto binary v0.50.0 (if k ≤ 31), which we installed above. This argument can be omitted to use the kallisto binary included in kb.

Building the index can take some time, depending on the number of threads used and the size of the D-list (here, ~20 min), but it only needs to be done once.

# Download the host reference genome and transcriptome from Ensembl
!gget ref -w cdna,dna -d $host_species

# Concatenate host genome and transcriptome into one file
host_genome = glob.glob("*.dna*.fa.gz")[0]
host_transcriptome = glob.glob("*.cdna.all.fa.gz")[0]
host_genome_transciptome = f"{host_species}.cdna_dna.fa.gz"
!cat $host_genome $host_transcriptome > $host_genome_transciptome

14:19:06 - INFO - Fetching reference information for homo_sapiens from Ensembl release: 113.
{
    "homo_sapiens": {
        "transcriptome_cdna": {
            "ftp": "http://ftp.ensembl.org/pub/release-113/fasta/homo_sapiens/cdna/Homo_sapiens.GRCh38.cdna.all.fa.gz",
            "ensembl_release": 113,
            "release_date": "2024-08-15",
            "release_time": "14:44",
            "bytes": "76M"
        },
        "genome_dna": {
            "ftp": "http://ftp.ensembl.org/pub/release-113/fasta/homo_sapiens/dna/Homo_sapiens.GRCh38.dna.primary_assembly.fa.gz",
            "ensembl_release": 113,
            "release_date": "2024-08-15",
            "release_time": "20:18",
            "bytes": "841M"
        }
    }
}
  % Total    % Received % Xferd  Average Speed   Time    Time     Time  Current
                                 Dload  Upload   Total   Spent    Left  Speed
100 75.9M  100 75.9M    0     0  19.3M      0  0:00:03  0:00:03 --:--:-- 19.3M
  % Total    % Received % Xferd  Average Speed   Time    Time     Time  Current
                                 Dload  Upload   Total   Spent    Left  Speed
100  841M  100  841M    0     0  29.6M      0  0:00:28  0:00:28 --:--:-- 30.9M

%%time
# Generate the viral protein reference index with masked host sequences using kallisto
reference_index = f"palmdb_{host_species}_dlist_cdna_dna.idx"

!kb ref \
    --kallisto $kallisto \
    --verbose \
    -t $threads \
    --aa \
    -k $k \
    --workflow custom \
    --d-list $host_genome_transciptome \
    -i $reference_index \
    palmdb_rdrp_seqs.fa

[2025-01-11 14:19:49,218]   DEBUG [main] Printing verbose output
[2025-01-11 14:19:51,424]   DEBUG [main] kallisto binary located at /content/kallisto/build/src/kallisto
[2025-01-11 14:19:51,424]   DEBUG [main] bustools binary located at /usr/local/lib/python3.10/dist-packages/kb_python/bins/linux/bustools/bustools
[2025-01-11 14:19:51,424]   DEBUG [main] Creating `tmp` directory
[2025-01-11 14:19:51,424]   DEBUG [main] Namespace(list=False, command='ref', tmp=None, keep_tmp=False, verbose=True, i='palmdb_homo_sapiens_dlist_cdna_dna.idx', g=None, f1=None, include_attribute=None, exclude_attribute=None, f2=None, c1=None, c2=None, d=None, k=31, t=2, d_list='homo_sapiens.cdna_dna.fa.gz', d_list_overhang=1, aa=True, workflow='custom', distinguish=False, make_unique=False, overwrite=False, kallisto='kallisto/build/src/kallisto', bustools='/usr/local/lib/python3.10/dist-packages/kb_python/bins/linux/bustools/bustools', opt_off=False, fasta='palmdb_rdrp_seqs.fa', gtf=None, feature=None, no_mismatches=False, ec_max_size=None, flank=None)
[2025-01-11 14:19:51,425]    INFO [ref_custom] Indexing palmdb_rdrp_seqs.fa to palmdb_homo_sapiens_dlist_cdna_dna.idx
[2025-01-11 14:19:51,425]   DEBUG [ref_custom] kallisto index -i palmdb_homo_sapiens_dlist_cdna_dna.idx -k 31 -t 2 -d homo_sapiens.cdna_dna.fa.gz --aa palmdb_rdrp_seqs.fa
[2025-01-11 14:19:51,526]   DEBUG [ref_custom]
[2025-01-11 14:19:51,526]   DEBUG [ref_custom] [build] loading fasta file palmdb_rdrp_seqs.fa
[2025-01-11 14:19:51,526]   DEBUG [ref_custom] [build] k-mer length: 31
[2025-01-11 14:19:53,831]   DEBUG [ref_custom] KmerStream::KmerStream(): Start computing k-mer cardinality estimations (1/2)
[2025-01-11 14:19:53,831]   DEBUG [ref_custom] KmerStream::KmerStream(): Start computing k-mer cardinality estimations (1/2)
[2025-01-11 14:20:02,359]   DEBUG [ref_custom] KmerStream::KmerStream(): Finished
[2025-01-11 14:20:02,961]   DEBUG [ref_custom] CompactedDBG::build(): Estimated number of k-mers occurring at least once: 37641510
[2025-01-11 14:20:02,962]   DEBUG [ref_custom] CompactedDBG::build(): Estimated number of minimizer occurring at least once: 7877811
[2025-01-11 14:20:15,821]   DEBUG [ref_custom] CompactedDBG::filter(): Processed 87630084 k-mers in 296561 reads
[2025-01-11 14:20:15,821]   DEBUG [ref_custom] CompactedDBG::filter(): Found 37508942 unique k-mers
[2025-01-11 14:20:15,822]   DEBUG [ref_custom] CompactedDBG::filter(): Number of blocks in Bloom filter is 257317
[2025-01-11 14:20:16,328]   DEBUG [ref_custom] CompactedDBG::construct(): Extract approximate unitigs (1/2)
[2025-01-11 14:21:41,073]   DEBUG [ref_custom] CompactedDBG::construct(): Extract approximate unitigs (2/2)
[2025-01-11 14:22:06,663]   DEBUG [ref_custom] CompactedDBG::construct(): Closed all input files
[2025-01-11 14:22:06,663]   DEBUG [ref_custom]
[2025-01-11 14:22:06,663]   DEBUG [ref_custom] CompactedDBG::construct(): Splitting unitigs (1/2)
[2025-01-11 14:22:06,763]   DEBUG [ref_custom]
[2025-01-11 14:22:06,764]   DEBUG [ref_custom] CompactedDBG::construct(): Splitting unitigs (2/2)
[2025-01-11 14:22:07,565]   DEBUG [ref_custom] CompactedDBG::construct(): Before split: 2040769 unitigs
[2025-01-11 14:22:07,566]   DEBUG [ref_custom] CompactedDBG::construct(): After split (1/1): 2040769 unitigs
[2025-01-11 14:22:07,566]   DEBUG [ref_custom] CompactedDBG::construct(): Unitigs split: 1559
[2025-01-11 14:22:07,566]   DEBUG [ref_custom] CompactedDBG::construct(): Unitigs deleted: 0
[2025-01-11 14:22:07,566]   DEBUG [ref_custom]
[2025-01-11 14:22:07,566]   DEBUG [ref_custom] CompactedDBG::construct(): Joining unitigs
[2025-01-11 14:22:18,006]   DEBUG [ref_custom] CompactedDBG::construct(): After join: 2021701 unitigs
[2025-01-11 14:22:18,006]   DEBUG [ref_custom] CompactedDBG::construct(): Joined 19439 unitigs
[2025-01-11 14:22:18,908]   DEBUG [ref_custom] [build] extracting distinguishing flanking k-mers from "homo_sapiens.cdna_dna.fa.gz"
[2025-01-11 14:52:44,153]   DEBUG [ref_custom] [build] identified 130 distinguishing flanking k-mers
[2025-01-11 14:52:54,811]   DEBUG [ref_custom] [build] building MPHF
[2025-01-11 14:53:05,472]   DEBUG [ref_custom] [build] creating equivalence classes ...
[2025-01-11 14:54:13,754]   DEBUG [ref_custom] [build] target de Bruijn graph has k-mer length 31 and minimizer length 23
[2025-01-11 14:54:13,855]   DEBUG [ref_custom] [build] target de Bruijn graph has 2021821 contigs and contains 37541835 k-mers
[2025-01-11 14:54:29,704]   DEBUG [ref_custom]
[2025-01-11 14:54:31,007]    INFO [ref_custom] Finished creating custom index
[2025-01-11 14:54:31,007]   DEBUG [main] Removing `tmp` directory
CPU times: user 14 s, sys: 1.85 s, total: 15.8 s
Wall time: 34min 50s

Alternatively, you can download a precomputed PalmDB reference index for use with kb in which the human (or mouse) genome and transcriptome were masked. You can find links to all available precomputed reference indices here.

Note: These indices were computed with k=31.

# # Download precomputed index with masked human genome and transcriptome
# k=31
# !wget https://data.caltech.edu/records/sh33z-hrx98/files/palmdb_human_dlist_cdna_dna.idx?download=1
# !mv palmdb_human_dlist_cdna_dna.idx?download=1 palmdb_human_dlist_cdna_dna.idx

Align data using kallisto translated search

Create a batch file so we can run all fastq files simultaneously (to learn more about batch files, see Box 7 in the Protocols paper):

if n_seq_to_keep:
    # To decrease the runtime of this example tutorial,
    # we will only align the top [n_seq_to_keep] sequences in each fastq file
    print(f"Number of reads kept per fastq file: {n_seq_to_keep:,}")

with open("batch.txt", "w") as batch_file:
  for filename in np.sort(glob.glob("*.fastq.gz")):

    if n_seq_to_keep:
        # Shorten fastq files (skip this step during a real analysis)
        n_rows_to_keep = n_seq_to_keep * 4
        filename_short = filename.split(".fastq.gz")[0] + "_short.fastq"
        !zcat $filename | head -$n_rows_to_keep > $filename_short
        # Delete original file to save disk space
        !rm $filename
    else:
      filename_short = filename

    # Write batch file in the following format:
    # sample_name \t fastq_filepath
    sample_name = filename.split("/")[-1].split(".")[0]
    batch_file.write(sample_name + "\t" +  filename_short + "\n")

Number of reads kept per fastq file: 10,000,000

The --aa argument tells kb that this is an amino acid reference.

The `-x argument <https://kallisto.readthedocs.io/en/latest/sc/technologies.html>`__ tells kb where to find the barcodes and UMIs in the data.

The --h5ad flag generates a .h5ad file from the count matrix.

NOTE for paired read data: kallisto translated search does not yet support the simultaneous alignment of paired reads, so parity should always be set to ‘single’ (--parity single) and the R1 and R2 files obtained using paired read sequencing should be aligned separately. The resulting matrices from the R1 and R2 reads should provide similar viral counts.

%%time
out_folder = "kb_output"

!kb count \
    --kallisto $kallisto \
    --verbose \
    -t $threads \
    --aa \
    -k $k \
    -x bulk \
    --parity single \
    -i $reference_index \
    -g palmdb_clustered_t2g.txt \
    --h5ad \
    -o $out_folder \
    batch.txt

[2025-01-11 14:57:23,533]   DEBUG [main] Printing verbose output
[2025-01-11 14:57:25,749]   DEBUG [main] kallisto binary located at /content/kallisto/build/src/kallisto
[2025-01-11 14:57:25,749]   DEBUG [main] bustools binary located at /usr/local/lib/python3.10/dist-packages/kb_python/bins/linux/bustools/bustools
[2025-01-11 14:57:25,749]   DEBUG [main] Creating `kb_output/tmp` directory
[2025-01-11 14:57:25,749]   DEBUG [main] Namespace(list=False, command='count', tmp=None, keep_tmp=False, verbose=True, i='palmdb_homo_sapiens_dlist_cdna_dna.idx', g='palmdb_clustered_t2g.txt', x='bulk', o='kb_output', num=False, w=None, r=None, t=2, m='2G', strand=None, inleaved=False, genomebam=False, aa=True, gtf=None, chromosomes=None, workflow='standard', em=False, mm=False, tcc=False, filter=None, filter_threshold=None, c1=None, c2=None, overwrite=False, dry_run=False, batch_barcodes=False, loom=False, h5ad=True, loom_names='barcode,target_name', sum='none', cellranger=False, gene_names=False, N=None, report=False, no_inspect=False, long=False, threshold=0.8, error_rate=None, platform='ONT', kallisto='kallisto/build/src/kallisto', bustools='/usr/local/lib/python3.10/dist-packages/kb_python/bins/linux/bustools/bustools', opt_off=False, k=31, no_validate=False, no_fragment=False, union=False, no_jump=False, quant_umis=False, keep_flags=False, parity='single', fragment_l=None, fragment_s=None, bootstraps=None, matrix_to_files=False, matrix_to_directories=False, fastqs=['batch.txt'])
[2025-01-11 14:57:29,555]    INFO [count] Using index palmdb_homo_sapiens_dlist_cdna_dna.idx to generate BUS file to kb_output from
[2025-01-11 14:57:29,555]    INFO [count]         /content/kb_output/tmp/tmpo_jkrwvs
[2025-01-11 14:57:29,555]   DEBUG [count] kallisto bus -i palmdb_homo_sapiens_dlist_cdna_dna.idx -o kb_output -x BULK -t 2 --aa --batch /content/kb_output/tmp/tmpo_jkrwvs
[2025-01-11 14:57:29,656]   DEBUG [count]
[2025-01-11 14:57:29,656]   DEBUG [count] [bus] will try running read files supplied in batch file
[2025-01-11 14:57:39,802]   DEBUG [count] [index] k-mer length: 31
[2025-01-11 14:58:12,736]   DEBUG [count] [index] number of targets: 296,561
[2025-01-11 14:58:12,737]   DEBUG [count] [index] number of k-mers: 37,541,835
[2025-01-11 14:58:12,737]   DEBUG [count] [index] number of distinguishing flanking k-mers: 130
[2025-01-11 14:58:12,737]   DEBUG [count] [quant] running in single-end mode
[2025-01-11 14:58:12,737]   DEBUG [count] [quant] will process file 1: SRR11777734_short.fastq
[2025-01-11 14:58:12,737]   DEBUG [count] [quant] will process file 1: SRR11777735_short.fastq
[2025-01-11 14:58:12,737]   DEBUG [count] [quant] will process file 1: SRR11777736_short.fastq
[2025-01-11 14:58:12,737]   DEBUG [count] [quant] will process file 1: SRR11777737_short.fastq
[2025-01-11 14:58:12,737]   DEBUG [count] [quant] will process file 1: SRR11777738_short.fastq
[2025-01-11 14:58:12,737]   DEBUG [count] [quant] will process file 1: SRR11777739_short.fastq
[2025-01-11 15:00:21,012]   DEBUG [count] [quant] finding pseudoalignments for all files ...
[2025-01-11 15:02:19,618]   DEBUG [count] [progress] 1M reads processed (0.0% mapped)
[2025-01-11 15:04:15,214]   DEBUG [count] [progress] 2M reads processed (0.0% mapped)
[2025-01-11 15:06:08,798]   DEBUG [count] [progress] 3M reads processed (0.0% mapped)
[2025-01-11 15:08:02,103]   DEBUG [count] [progress] 4M reads processed (0.0% mapped)
[2025-01-11 15:10:22,726]   DEBUG [count] [progress] 5M reads processed (0.0% mapped)
[2025-01-11 15:12:40,731]   DEBUG [count] [progress] 6M reads processed (0.0% mapped)
[2025-01-11 15:14:55,730]   DEBUG [count] [progress] 7M reads processed (0.0% mapped)
[2025-01-11 15:17:12,941]   DEBUG [count] [progress] 8M reads processed (0.0% mapped)
[2025-01-11 15:19:32,636]   DEBUG [count] [progress] 9M reads processed (0.0% mapped)
[2025-01-11 15:21:36,944]   DEBUG [count] [progress] 11M reads processed (0.0% mapped)
[2025-01-11 15:23:48,579]   DEBUG [count] [progress] 12M reads processed (0.0% mapped)
[2025-01-11 15:25:57,169]   DEBUG [count] [progress] 13M reads processed (0.0% mapped)
[2025-01-11 15:27:59,126]   DEBUG [count] [progress] 14M reads processed (0.0% mapped)
[2025-01-11 15:30:11,821]   DEBUG [count] [progress] 15M reads processed (0.0% mapped)
[2025-01-11 15:32:10,526]   DEBUG [count] [progress] 16M reads processed (0.0% mapped)
[2025-01-11 15:34:02,451]   DEBUG [count] [progress] 17M reads processed (0.0% mapped)
[2025-01-11 15:36:01,092]   DEBUG [count] [progress] 18M reads processed (0.0% mapped)
[2025-01-11 15:37:51,728]   DEBUG [count] [progress] 19M reads processed (0.0% mapped)
[2025-01-11 15:39:57,462]   DEBUG [count] [progress] 20M reads processed (0.0% mapped)
[2025-01-11 15:41:49,835]   DEBUG [count] [progress] 21M reads processed (0.0% mapped)
[2025-01-11 15:43:41,889]   DEBUG [count] [progress] 23M reads processed (0.0% mapped)
[2025-01-11 15:45:38,517]   DEBUG [count] [progress] 24M reads processed (0.0% mapped)
[2025-01-11 15:47:33,131]   DEBUG [count] [progress] 25M reads processed (0.0% mapped)
[2025-01-11 15:49:54,888]   DEBUG [count] [progress] 26M reads processed (0.0% mapped)
[2025-01-11 15:51:46,296]   DEBUG [count] [progress] 27M reads processed (0.0% mapped)
[2025-01-11 15:53:50,729]   DEBUG [count] [progress] 28M reads processed (0.0% mapped)
[2025-01-11 15:55:57,382]   DEBUG [count] [progress] 29M reads processed (0.0% mapped)
[2025-01-11 15:57:56,238]   DEBUG [count] [progress] 30M reads processed (0.0% mapped)
[2025-01-11 15:59:34,752]   DEBUG [count] [progress] 31M reads processed (0.0% mapped)
[2025-01-11 16:01:40,671]   DEBUG [count] [progress] 32M reads processed (0.0% mapped)
[2025-01-11 16:03:43,186]   DEBUG [count] [progress] 33M reads processed (0.0% mapped)
[2025-01-11 16:05:41,756]   DEBUG [count] [progress] 35M reads processed (0.0% mapped)
[2025-01-11 16:07:34,144]   DEBUG [count] [progress] 36M reads processed (0.0% mapped)
[2025-01-11 16:09:41,273]   DEBUG [count] [progress] 37M reads processed (0.0% mapped)
[2025-01-11 16:11:44,253]   DEBUG [count] [progress] 38M reads processed (0.0% mapped)
[2025-01-11 16:13:37,678]   DEBUG [count] [progress] 39M reads processed (0.0% mapped)
[2025-01-11 16:15:20,938]   DEBUG [count] [progress] 40M reads processed (0.0% mapped)
[2025-01-11 16:17:30,574]   DEBUG [count] [progress] 41M reads processed (0.0% mapped)
[2025-01-11 16:19:11,738]   DEBUG [count] [progress] 42M reads processed (0.0% mapped)
[2025-01-11 16:21:26,263]   DEBUG [count] [progress] 43M reads processed (0.0% mapped)
[2025-01-11 16:23:10,200]   DEBUG [count] [progress] 44M reads processed (0.0% mapped)
[2025-01-11 16:25:00,529]   DEBUG [count] [progress] 45M reads processed (0.0% mapped)
[2025-01-11 16:26:43,524]   DEBUG [count] [progress] 46M reads processed (0.0% mapped)
[2025-01-11 16:28:33,196]   DEBUG [count] [progress] 47M reads processed (0.0% mapped)
[2025-01-11 16:30:25,325]   DEBUG [count] [progress] 48M reads processed (0.0% mapped)
[2025-01-11 16:32:25,779]   DEBUG [count] [progress] 49M reads processed (0.0% mapped)
[2025-01-11 16:34:15,891]   DEBUG [count] [progress] 50M reads processed (0.0% mapped)
[2025-01-11 16:36:08,055]   DEBUG [count] [progress] 51M reads processed (0.0% mapped)
[2025-01-11 16:37:48,494]   DEBUG [count] [progress] 52M reads processed (0.0% mapped)
[2025-01-11 16:39:31,229]   DEBUG [count] [progress] 53M reads processed (0.0% mapped)
[2025-01-11 16:41:47,782]   DEBUG [count] [progress] 54M reads processed (0.0% mapped)
[2025-01-11 16:43:50,284]   DEBUG [count] [progress] 55M reads processed (0.0% mapped)
[2025-01-11 16:46:10,261]   DEBUG [count] [progress] 56M reads processed (0.0% mapped)
[2025-01-11 16:48:01,779]   DEBUG [count] [progress] 57M reads processed (0.0% mapped)
[2025-01-11 16:50:22,105]   DEBUG [count] [progress] 58M reads processed (0.0% mapped)
[2025-01-11 16:50:50,525]   DEBUG [count] [progress] 59M reads processed (0.0% mapped)              done
[2025-01-11 16:50:50,526]   DEBUG [count] [quant] processed 60,000,000 reads, 18,380 reads pseudoaligned
[2025-01-11 16:50:50,826]   DEBUG [count]
[2025-01-11 16:50:54,552]    INFO [count] Sorting BUS file kb_output/output.bus to kb_output/tmp/output.s.bus
[2025-01-11 16:50:54,552]   DEBUG [count] bustools sort -o kb_output/tmp/output.s.bus -T kb_output/tmp -t 2 -m 2G kb_output/output.bus
[2025-01-11 16:50:56,700]   DEBUG [count] all fits in buffer
[2025-01-11 16:50:58,002]   DEBUG [count] Read in 18380 BUS records
[2025-01-11 16:50:58,003]   DEBUG [count] reading time 0.000307s
[2025-01-11 16:50:58,003]   DEBUG [count] sorting time 0.002648s
[2025-01-11 16:50:58,003]   DEBUG [count] writing time 0.00043s
[2025-01-11 16:50:58,003]    INFO [count] Inspecting BUS file kb_output/tmp/output.s.bus
[2025-01-11 16:50:58,003]   DEBUG [count] bustools inspect -o kb_output/inspect.json kb_output/tmp/output.s.bus
[2025-01-11 16:50:59,106]    INFO [count] Generating count matrix kb_output/counts_unfiltered/cells_x_genes from BUS file kb_output/tmp/output.s.bus
[2025-01-11 16:50:59,107]   DEBUG [count] bustools count -o kb_output/counts_unfiltered/cells_x_genes -g palmdb_clustered_t2g.txt -e kb_output/matrix.ec -t kb_output/transcripts.txt --genecounts --cm kb_output/tmp/output.s.bus
[2025-01-11 16:51:02,377]   DEBUG [count] kb_output/counts_unfiltered/cells_x_genes.mtx passed validation
[2025-01-11 16:51:02,379]    INFO [count] Writing gene names to file kb_output/counts_unfiltered/cells_x_genes.genes.names.txt
[2025-01-11 16:51:02,945] WARNING [count] 99228 gene IDs do not have corresponding valid gene names. These genes will use their gene IDs instead.
[2025-01-11 16:51:02,994]    INFO [count] Reading matrix kb_output/counts_unfiltered/cells_x_genes.mtx
[2025-01-11 16:51:03,502]    INFO [count] Writing matrix to h5ad kb_output/counts_unfiltered/adata.h5ad
[2025-01-11 16:51:03,854]   DEBUG [main] Removing `kb_output/tmp` directory
CPU times: user 50.9 s, sys: 6.98 s, total: 57.8 s
Wall time: 1h 54min 1s

Load generated count matrix

# Open count matrix generated by kb (AnnData object in h5ad format)
adata = anndata.read_h5ad(f"{out_folder}/counts_unfiltered/adata.h5ad")
adata

AnnData object with n_obs × n_vars = 6 × 99228

Add SRR and sample labels to adata.obs based on the batch barcode assigned by kb count:

def read_list_from_txt(filepath):
  with open(filepath, "r") as file:
      my_list = [line.strip() for line in file]
  return my_list

batch_barcode_df = pd.DataFrame()
batch_barcode_df["batch_name"] = read_list_from_txt(f"{out_folder}/matrix.cells")
batch_barcode_df["batch_barcode"] = read_list_from_txt(f"{out_folder}/matrix.sample.barcodes")
batch_barcode_df

	batch_name	batch_barcode
0	SRR11777734	AAAAAAAAAAAAAAAA
1	SRR11777735	AAAAAAAAAAAAAAAC
2	SRR11777736	AAAAAAAAAAAAAAAG
3	SRR11777737	AAAAAAAAAAAAAAAT
4	SRR11777738	AAAAAAAAAAAAAACA
5	SRR11777739	AAAAAAAAAAAAAACC

adata.obs = adata.obs.merge(batch_barcode_df, left_on="barcode", right_on="batch_barcode").set_index("batch_barcode")
# Add metadata for each SRR file
adata.obs["sample"] = ['Infected 1', 'Infected 2', 'Infected 3', 'Control 1', 'Control 2', 'Control 3']
adata.obs

	batch_name	sample
batch_barcode
AAAAAAAAAAAAAAAA	SRR11777734	Infected 1
AAAAAAAAAAAAAAAC	SRR11777735	Infected 2
AAAAAAAAAAAAAAAG	SRR11777736	Infected 3
AAAAAAAAAAAAAAAT	SRR11777737	Control 1
AAAAAAAAAAAAAACA	SRR11777738	Control 2
AAAAAAAAAAAAAACC	SRR11777739	Control 3

Plot counts for SARS-CoV RdRP sequences in each sample

Find the virus IDs for SARS-CoV in the PalmDB viral protein reference:

# Load the PalmDB virus ID to virus taxonomy mapping
# Note: Only the rep_ID will occur in the reference index
u_tax_csv = "ID_to_taxonomy_mapping.csv"
tax_df = pd.read_csv(u_tax_csv)

# Find the virus IDs for SARS-CoV proteins
tax_df[tax_df["species"].str.contains("Severe acute respiratory syndrome")]

	ID	rep_ID	phylum	class	order	family	genus	species	strandedness
0	u1	u1	Pisuviricota	Pisoniviricetes	Nidovirales	Coronaviridae	Betacoronavirus	Severe acute respiratory syndrome-related coro...	+ssRNA
1	u10066	u1	Pisuviricota	Pisoniviricetes	Nidovirales	Coronaviridae	Betacoronavirus	Severe acute respiratory syndrome-related coro...	+ssRNA
2	u10121	u1	Pisuviricota	Pisoniviricetes	Nidovirales	Coronaviridae	Betacoronavirus	Severe acute respiratory syndrome-related coro...	+ssRNA
3	u10124	u1	Pisuviricota	Pisoniviricetes	Nidovirales	Coronaviridae	Betacoronavirus	Severe acute respiratory syndrome-related coro...	+ssRNA
4	u102609	u1	Pisuviricota	Pisoniviricetes	Nidovirales	Coronaviridae	Betacoronavirus	Severe acute respiratory syndrome-related coro...	+ssRNA
...	...	...	...	...	...	...	...	...	...
404	u889	u1	Pisuviricota	Pisoniviricetes	Nidovirales	Coronaviridae	Betacoronavirus	Severe acute respiratory syndrome-related coro...	+ssRNA
405	u9015	u1	Pisuviricota	Pisoniviricetes	Nidovirales	Coronaviridae	Betacoronavirus	Severe acute respiratory syndrome-related coro...	+ssRNA
406	u944	u1	Pisuviricota	Pisoniviricetes	Nidovirales	Coronaviridae	Betacoronavirus	Severe acute respiratory syndrome-related coro...	+ssRNA
407	u9771	u1	Pisuviricota	Pisoniviricetes	Nidovirales	Coronaviridae	Betacoronavirus	Severe acute respiratory syndrome-related coro...	+ssRNA
408	u992	u1	Pisuviricota	Pisoniviricetes	Nidovirales	Coronaviridae	Betacoronavirus	Severe acute respiratory syndrome-related coro...	+ssRNA

409 rows × 9 columns

# Target virus IDs to plot (here, SARS-CoV)
target_ids = tax_df[tax_df["species"].str.contains("Severe acute respiratory syndrome-related coronavirus")]["rep_ID"].unique()
target_ids

array(['u1'], dtype=object)

Create bar plot showing raw counts for SARS-CoV:

fig, ax = plt.subplots(figsize=(5, 5))
fontsize = 14
width = 0.75

samples = adata.obs.index.values
x_labels = adata.obs["sample"].values

# Plot total raw counts
counts = []
labels = samples
for sample in samples:
    counts.append(adata.X[adata.obs.index == sample, adata.var.index.isin(target_ids)].sum())

x = np.arange(len(labels))

ax.bar(x, counts, width=width, color="#003049", edgecolor="black")

# Adjust figure labels and axes
ax.set_yscale("symlog")
ax.set_ylabel(f"Total counts for {target_ids[0]} (SARS-CoV))", fontsize=fontsize)
# ax.set_xlabel("Sample", fontsize=fontsize)
ax.set_xticks(x, x_labels, rotation=45, ha="right")
ax.tick_params(axis="both", labelsize=fontsize-2)
ax.set_title(f"SARS-CoV-2 infected human\niPSC-derived cardiomyocytes", fontsize=fontsize+2)
ax.grid(True, which="both", color="lightgray", ls="--", lw=1)
ax.set_axisbelow(True)

# Save figure
fig.savefig("sars_cov_rdrp_count.png", dpi=300, bbox_inches="tight")

fig.show()

counts

[611.0, 945.0, 734.0, 0.0, 0.0, 0.0]

Visualize other virus-like sequences observed in the dataset

When interpreting the presence of other virus- or RdRP-like sequences (each identified by its own virus ID), keep in mind that there will likely be many RdRP-like sequences introduced by contamination of laboratory reagents. A (non-comprehensive) list of virus IDs observed in blank sequencing data is available here and we will download it below for comparison to the virus IDs observed in our dataset. Another indication that a sequence may have originated from sample contamination is its uniform presence across all sequenced samples.

# Download list of virus IDs seen in blank sequencing reagents
!wget https://raw.githubusercontent.com/pachterlab/LSCHWCP_2023/refs/heads/main/viruses_in_blank_reagents/total_raw_count_per_virus_id_in_laboratory_reagents.csv
blank_df = pd.read_csv("total_raw_count_per_virus_id_in_laboratory_reagents.csv")
blank_df

--2025-01-11 16:51:18--  https://raw.githubusercontent.com/pachterlab/LSCHWCP_2023/refs/heads/main/viruses_in_blank_reagents/total_raw_count_per_virus_id_in_laboratory_reagents.csv
Resolving raw.githubusercontent.com (raw.githubusercontent.com)... 185.199.108.133, 185.199.109.133, 185.199.110.133, ...
Connecting to raw.githubusercontent.com (raw.githubusercontent.com)|185.199.108.133|:443... connected.
HTTP request sent, awaiting response... 200 OK
Length: 54642 (53K) [text/plain]
Saving to: ‘total_raw_count_per_virus_id_in_laboratory_reagents.csv’

total_raw_count_per 100%[===================>]  53.36K  --.-KB/s    in 0.009s

2025-01-11 16:51:18 (5.76 MB/s) - ‘total_raw_count_per_virus_id_in_laboratory_reagents.csv’ saved [54642/54642]

	count	virus
0	5374210.0	u172514
1	4301784.0	u226460
2	752771.0	u237705
3	456424.0	u202260
4	359783.0	u223701
...	...	...
4573	1.0	u169611
4574	1.0	u41840
4575	1.0	u169761
4576	1.0	u169999
4577	1.0	u158491

4578 rows × 2 columns

Plot the top 20 virus IDs observed in this dataset (raw counts are normalized to the total number of viral reads in each sample, and virus IDs also observed in blank sequencing reagents are marked in red):

def plot_top_virs_stacked_bar(adata, top_n=20):
    # Convert sparse matrix to dense if needed
    counts_matrix = adata.X
    if not isinstance(counts_matrix, np.ndarray):
        counts_matrix = counts_matrix.toarray()

    # Calculate total reads per sample
    sample_totals = counts_matrix.sum(axis=1, keepdims=True)

    # Normalize counts by total reads per sample to get fractions
    normalized_counts = counts_matrix / sample_totals

    # Sum normalized counts across all samples for each virus
    vir_totals = normalized_counts.sum(axis=0)

    # Get the top N viruses across all samples
    top_vir_indices = np.argsort(vir_totals)[-top_n:][::-1]
    top_virs = adata.var_names[top_vir_indices]
    top_counts = normalized_counts[:, top_vir_indices]

    sample_ids = adata.obs["sample"].values

    # Create a stacked bar plot
    fig, ax = plt.subplots(figsize=(7, 5))
    colors = ['#e34a33', '#fdbb84', '#fee8c8', '#2b8cbe', '#a6bddb', '#ece7f2']
    fontsize = 14

    bar_bottom = np.zeros(top_n)
    for i, sample_id in enumerate(sample_ids):
        sample_fractions = top_counts[i, :]
        ax.bar(
            top_virs,
            sample_fractions,
            bottom=bar_bottom,
            color=colors[i % len(colors)],
            label=sample_id
        )
        bar_bottom += sample_fractions

    # Adjust figure labels and axes
    ax.set_title(f"Top {top_n} Virus IDs", fontsize=fontsize+2)
    ax.set_xlabel("Virus ID", fontsize=fontsize)
    ax.set_ylabel("Fraction of Total Reads in Sample", fontsize=fontsize)
    ax.set_xticklabels(top_virs, rotation=45, ha='right')
    # Mark virus IDs that have been observed in blank sequencing reagents in red
    for i, label in enumerate(top_virs):
        if label in blank_df["virus"].values:
            ax.get_xticklabels()[i].set_color('red')
    ax.legend(title="Sample", loc='upper right', fontsize=fontsize-2, title_fontsize=fontsize)
    ax.tick_params(axis='both', which='major', labelsize=fontsize-2)
    ax.grid(axis='y', linestyle='--', color="lightgrey", alpha=1)
    ax.set_axisbelow(True)
    ax.margins(x=0.01)
    ax.set_ylim(0, np.max(bar_bottom) * 1.05)

    # Save figure
    fig.savefig("top_virus_IDs_per_sample.png", dpi=300, bbox_inches="tight")

    fig.show()

    return top_virs

top_virs = plot_top_virs_stacked_bar(adata, top_n=20)

Show the predicted taxonomies for each virus ID shown in the plot above (these are the species-like operational taxonomic units (OTUs) provided by PalmDB):

top_vir_tax = tax_df[tax_df["rep_ID"].isin(top_virs)][tax_df.columns[1:]].drop_duplicates()
top_vir_tax['rep_ID'] = pd.Categorical(top_vir_tax['rep_ID'], categories=top_virs, ordered=True)
top_vir_tax = top_vir_tax.sort_values('rep_ID')
top_vir_tax

	rep_ID	phylum	class	order	family	genus	species	strandedness
0	u1	Pisuviricota	Pisoniviricetes	Nidovirales	Coronaviridae	Betacoronavirus	Severe acute respiratory syndrome-related coro...	+ssRNA
35782	u34159	.	.	.	.	.	.	unknown
150602	u100074	Lenarviricota	Howeltoviricetes	Cryppavirales	.	.	.	+ssRNA
34276	u27694	Peploviricota	Herviviricetes	Herpesvirales	Herpesviridae	Varicellovirus	Bubaline alphaherpesvirus 1	dsDNA
66423	u100002	Lenarviricota	Allassoviricetes	Levivirales	.	.	.	+ssRNA
37641	u39566	.	.	.	.	.	.	unknown
36132	u35404	.	.	.	.	.	.	unknown
207936	u100644	Lenarviricota	Amabiliviricetes	Wolframvirales	.	.	.	+ssRNA
106844	u100017	Lenarviricota	Allassoviricetes	Levivirales	Leviviridae	.	.	+ssRNA
37033	u37761	.	.	.	.	.	.	unknown
35795	u34182	.	.	.	.	.	.	unknown
38543	u41991	.	.	.	.	.	.	unknown
263358	u202260	.	.	.	.	.	.	unknown
35684	u33987	.	.	.	.	.	.	unknown
38443	u41838	.	.	.	.	.	.	unknown
148078	u100049	Lenarviricota	Miaviricetes	.	.	.	.	+ssRNA
95765	u100012	Lenarviricota	Allassoviricetes	.	.	.	.	+ssRNA
35685	u33988	.	.	.	.	.	.	unknown
142323	u100041	Kitrinoviricota	.	.	.	.	.	+ssRNA
36058	u34860	.	.	.	.	.	.	unknown

To confirm the presence of a known virus for which a reference genome is available, you can align the data to the viral (nucleotide) reference genome using the standard kb workflow (simply omit the --aa flag in the kb ref and kb count commands and supply the viral reference genome to kb ref instead of the PalmDB protein reference).

Extract reads that aligned to specific virus IDs for further analysis

If this approach fails or no reference genome is available, use kb extract to extract raw reads that aligned to a given virus ID for further analysis (such as BLAST alignment):

# Sequences that matched to [vir_ids_to_extract] will be extracted
# Here, we will extract reads that mapped to the virus ID u100074
# u100074 is interesting because it was only observed in infected samples
vir_ids_to_extract = ["u100074"]

# Provide a single fastq file to extract reads from (here, we are using "Infected 1")
fastq_file = "/content/SRR11777734_short.fastq"

%%time
# Extract reads using kb extract
# kb extract will output a new fastq file containing the extracted reads for each target virus ID
outfolder_extract = "extracted_seqs"
targets = " ".join(vir_ids_to_extract)

!kb extract \
    --kallisto $kallisto \
    --verbose \
    -t $threads \
    --aa \
    -g palmdb_clustered_t2g.txt \
    -i $reference_index \
    -ts $targets \
    -o $outfolder_extract \
    $fastq_file

[2025-01-11 16:51:42,733]   DEBUG [main] Printing verbose output
[2025-01-11 16:51:44,941]   DEBUG [main] kallisto binary located at /content/kallisto/build/src/kallisto
[2025-01-11 16:51:44,941]   DEBUG [main] bustools binary located at /usr/local/lib/python3.10/dist-packages/kb_python/bins/linux/bustools/bustools
[2025-01-11 16:51:44,941]   DEBUG [main] Creating `extracted_seqs/tmp` directory
[2025-01-11 16:51:44,941]   DEBUG [main] Namespace(list=False, command='extract', tmp=None, keep_tmp=False, verbose=True, fastq='/content/SRR11777734_short.fastq', i='palmdb_homo_sapiens_dlist_cdna_dna.idx', targets=['u100074'], target_type='gene', extract_all=False, extract_all_fast=False, extract_all_unmapped=False, mm=False, g='palmdb_clustered_t2g.txt', o='extracted_seqs', t=2, strand=None, aa=True, N=None, kallisto='kallisto/build/src/kallisto', bustools='/usr/local/lib/python3.10/dist-packages/kb_python/bins/linux/bustools/bustools', opt_off=False, k=31)
[2025-01-11 16:51:47,837]    INFO [extract] Performing alignment using kallisto...
[2025-01-11 16:51:47,837]    INFO [extract] Using index palmdb_homo_sapiens_dlist_cdna_dna.idx to generate BUS file to extracted_seqs/tmp from
[2025-01-11 16:51:47,837]    INFO [extract]         /content/SRR11777734_short.fastq
[2025-01-11 16:51:47,837]   DEBUG [extract] kallisto bus -i palmdb_homo_sapiens_dlist_cdna_dna.idx -o extracted_seqs/tmp -x bulk -t 2 --num --aa /content/SRR11777734_short.fastq
[2025-01-11 16:51:47,939]   DEBUG [extract]
[2025-01-11 16:51:58,604]   DEBUG [extract] [index] k-mer length: 31
[2025-01-11 16:52:23,202]   DEBUG [extract] [index] number of targets: 296,561
[2025-01-11 16:52:23,202]   DEBUG [extract] [index] number of k-mers: 37,541,835
[2025-01-11 16:52:23,202]   DEBUG [extract] [index] number of distinguishing flanking k-mers: 130
[2025-01-11 16:52:23,202]   DEBUG [extract] [quant] running in single-end mode
[2025-01-11 16:52:23,202]   DEBUG [extract] [quant] will process file 1: /content/SRR11777734_short.fastq
[2025-01-11 16:54:19,189]   DEBUG [extract] [quant] finding pseudoalignments for all files ...
[2025-01-11 16:56:19,089]   DEBUG [extract] [progress] 1M reads processed (0.0% mapped)
[2025-01-11 16:58:47,495]   DEBUG [extract] [progress] 2M reads processed (0.0% mapped)
[2025-01-11 17:01:01,756]   DEBUG [extract] [progress] 3M reads processed (0.0% mapped)
[2025-01-11 17:03:02,161]   DEBUG [extract] [progress] 4M reads processed (0.0% mapped)
[2025-01-11 17:04:59,634]   DEBUG [extract] [progress] 5M reads processed (0.0% mapped)
[2025-01-11 17:07:01,514]   DEBUG [extract] [progress] 6M reads processed (0.0% mapped)
[2025-01-11 17:09:00,622]   DEBUG [extract] [progress] 7M reads processed (0.0% mapped)
[2025-01-11 17:11:14,126]   DEBUG [extract] [progress] 8M reads processed (0.0% mapped)
[2025-01-11 17:11:21,856]   DEBUG [extract] [progress] 9M reads processed (0.0% mapped)              done
[2025-01-11 17:11:21,856]   DEBUG [extract] [quant] processed 10,000,000 reads, 3,086 reads pseudoaligned
[2025-01-11 17:11:21,956]   DEBUG [extract]
[2025-01-11 17:11:24,972]    INFO [extract] Alignment complete. Beginning extraction of reads using bustools...
[2025-01-11 17:11:24,972]   DEBUG [extract] Replacing transcript entries with -1 for equivalence classes that map to multiple genes from extracted_seqs/tmp/matrix.ec
[2025-01-11 17:11:24,973]   DEBUG [extract] Fetching equivalence classes that map to multiple genes...
[2025-01-11 17:12:14,721]   DEBUG [extract] Found the following equivalence classes that map to multiple genes: 2 ,4 ,5 ,7 ,9 ,13 ,15 ,16 ,17 ,18 ,20 ,22 ,23 ,24 ,25 ,32 ,34 ,35 ,36 ,37 ,38 ,39 ,42 ,48 ,49 ,50 ,53 ,54 ,55 ,56 ,57 ,58 ,59 ,60 ,61 ,65 ,66 ,69 ,70 ,71 ,72 ,74 ,75 ,78 ,80 ,82 ,83 ,87 ,90 ,94 ,96 ,97 ,98 ,100 ,101 ,102 ,103 ,104 ,105 ,107 ,108 ,110 ,114 ,115 ,119 ,120 ,121 ,123 ,128 ,129 ,132 ,133 ,135 ,136 ,138 ,145 ,146 ,147 ,150 ,151 ,152 ,153 ,157 ,158 ,159 ,160 ,161 ,162 ,164 ,166 ,167 ,168 ,170 ,171 ,172 ,173 ,177 ,179 ,182 ,183 ,186 ,187 ,188 ,190 ,191 ,192 ,193 ,194 ,196 ,197 ,200 ,202 ,204 ,205 ,207 ,208 ,210 ,211 ,213 ,215 ,216 ,217 ,219 ,221 ,225 ,226 ,229 ,231 ,233 ,235 ,236 ,237 ,239 ,242 ,245 ,246 ,247 ,248 ,249 ,255 ,256 ,258 ,259 ,263 ,264 ,265 ,267 ,268 ,269 ,271 ,274 ,277 ,279 ,281 ,283 ,284 ,285 ,290 ,291 ,292 ,293 ,296 ,297 ,298 ,299 ,300 ,303 ,304 ,305 ,306 ,307 ,312 ,314 ,315 ,316 ,320 ,323 ,324 ,325 ,328 ,329 ,337 ,338 ,339 ,346 ,349 ,354 ,355 ,357 ,358 ,359 ,360 ,362 ,363 ,365 ,367 ,369 ,370 ,371 ,373 ,374 ,376 ,377 ,381 ,382 ,383 ,384 ,385 ,386 ,387 ,389 ,390 ,391 ,396 ,399 ,400 ,401 ,402 ,403 ,407 ,408 ,413 ,415 ,416 ,418 ,420 ,421 ,423 ,428 ,429 ,430 ,432 ,433 ,435 ,439 ,441 ,444 ,445 ,447 ,448 ,449 ,453 ,455 ,458 ,459 ,461 ,464 ,465 ,467 ,472 ,475 ,478 ,481 ,482 ,485 ,486 ,490 ,492 ,495 ,496 ,497 ,498 ,502 ,503 ,504 ,505 ,507 ,508 ,509 ,513 ,515 ,518 ,519 ,522 ,523 ,525 ,532
[2025-01-11 17:12:14,843]   DEBUG [extract] matrix.ec file where transcript entries were replaced with -1 for equivalence classes that map to multiple genes saved at extracted_seqs/tmp/matrix_no_mm.ec
[2025-01-11 17:12:15,266]    INFO [extract] Extracting reads for following transcripts for gene ID u100074: u100074, u100077, u100080, u100082, u100106, u100192, u100299, u100338, u100387, u100423, u100460, u100513, u100518, u100623, u100634, u100649, u100684, u100772, u100780, u100804, u100895, u100924, u100931, u101016, u101036, u101069, u101090, u101092, u101093, u101105, u101184, u101231, u101235, u101236, u101251, u101300, u101343, u101352, u101374, u101387, u101409, u101411, u101420, u101425, u101426, u101436, u101438, u101442, u101489, u101490, u101586, u101592, u101599, u101626, u101648, u101660, u101719, u101749, u101779, u101815, u101829, u101840, u101881, u101920, u101977, u102031, u102069, u102077, u102112, u102140, u102169, u102232, u102236, u102261, u102294, u102304, u102487, u102494, u102575, u102656, u102764, u102778, u102803, u102952, u103011, u103033, u103040, u103115, u103121, u103176, u103185, u103189, u103369, u103370, u103375, u103455, u103478, u103555, u103685, u103693, u103727, u103854, u104111, u104124, u104213, u104349, u104381, u104625, u104641, u104724, u104760, u104885, u104975, u105128, u105131, u105147, u105154, u105175, u105195, u105197, u105251, u105252, u105254, u105287, u105345, u105353, u105392, u105438, u105486, u105512, u105521, u105570, u105593, u105614, u105629, u105823, u105828, u105844, u105852, u105868, u105873, u105898, u106143, u106176, u106306, u106410, u106450, u106560, u106674, u106680, u106699, u106718, u106719, u106733, u106734, u106756, u106769, u106819, u106824, u106944, u106950, u106956, u106986, u107010, u107033, u107055, u107071, u107110, u107132, u107161, u107326, u107336, u107337, u107351, u107398, u107480, u107558, u107600, u107705, u107722, u107730, u107739, u107767, u107829, u107913, u107919, u107935, u108011, u108046, u108126, u108143, u108178, u108215, u108223, u108297, u108324, u108358, u108423, u108424, u108463, u108465, u108511, u108594, u108670, u108673, u108682, u108689, u108698, u108716, u108758, u108768, u108784, u108802, u108814, u108819, u108859, u108878, u108930, u108939, u108941, u109006, u109012, u109022, u109112, u109135, u109207, u109213, u109252, u109268, u109279, u109341, u109378, u109381, u109393, u109469, u109561, u109626, u109657, u109714, u109773, u109838, u109857, u109858, u109893, u109968, u110056, u110095, u110100, u110150, u110188, u110214, u110217, u110261, u110313, u110323, u110462, u110474, u110536, u110550, u110569, u110581, u110583, u110684, u110695, u110709, u110716, u110768, u110776, u110850, u110870, u110892, u110951, u110959, u110967, u111005, u111081, u111093, u111104, u111152, u111262, u111264, u111280, u111301, u111317, u111398, u111401, u111403, u111405, u111411, u111457, u111481, u111799, u111843, u112320, u112335, u112354, u112361, u112486, u11253, u112577, u112591, u112606, u112643, u112980, u113051, u113091, u113099, u113106, u113115, u113157, u113164, u113210, u113268, u113270, u113303, u113358, u113392, u113468, u113554, u113705, u113747, u113762, u113868, u113875, u113882, u113938, u113972, u113984, u113991, u114025, u114029, u114051, u114171, u114178, u114203, u114268, u114272, u114296, u114309, u114311, u114320, u114404, u114437, u114551, u114591, u114644, u114655, u114735, u114737, u114756, u114768, u114774, u114786, u114824, u114830, u114837, u114849, u114883, u114885, u114903, u115018, u115085, u115115, u115140, u115173, u115204, u115205, u115245, u115279, u115293, u115299, u115305, u115337, u115368, u115391, u115444, u115461, u115463, u115472, u115497, u115543, u11555, u115565, u115577, u115580, u115594, u115609, u115622, u115626, u115643, u115645, u115664, u115686, u115690, u115724, u115739, u115787, u115799, u115899, u115904, u115932, u116004, u116033, u116039, u116052, u116062, u116083, u116084, u116131, u116136, u116146, u116151, u116154, u116155, u116190, u116218, u116271, u116359, u116407, u116462, u116475, u116515, u116624, u116632, u116642, u116658, u116663, u116710, u116721, u116731, u116735, u116805, u116821, u116832, u116868, u116914, u116955, u116978, u117050, u117133, u117188, u117279, u117280, u117294, u117425, u117535, u117626, u117845, u117904, u118083, u118113, u118235, u118288, u118303, u118305, u118356, u118419, u118427, u118429, u118430, u118482, u118487, u118513, u118539, u118579, u118613, u118615, u118618, u118692, u118745, u118754, u118824, u118861, u118872, u118898, u118931, u119026, u119030, u119141, u119182, u119372, u119399, u119401, u119409, u119423, u119449, u119487, u119509, u119535, u119557, u119570, u119610, u119655, u119775, u119785, u119808, u119818, u119849, u119859, u119898, u119939, u119946, u11997, u119996, u120005, u120027, u120057, u120061, u120075, u120077, u120078, u120081, u120082, u120090, u120140, u120143, u120221, u120252, u120254, u120257, u120265, u120270, u120328, u120335, u120350, u120364, u120373, u120409, u120440, u120442, u120469, u120557, u120586, u120599, u120609, u120619, u120637, u120677, u120683, u120767, u120787, u120835, u120836, u120844, u120873, u120875, u120892, u120911, u120912, u120939, u120943, u120945, u120947, u121033, u121097, u121128, u121150, u121200, u121299, u121300, u121310, u121447, u121514, u121545, u121546, u121565, u121597, u121623, u121643, u121653, u121698, u121709, u121777, u121831, u121858, u121867, u121871, u121887, u121908, u121974, u122078, u122086, u122096, u122124, u122130, u122142, u122145, u122190, u122203, u122259, u122264, u122270, u122286, u122390, u122391, u122398, u122422, u122437, u122468, u122483, u122498, u122521, u122532, u122676, u122684, u122737, u122783, u122854, u122968, u123039, u123045, u123052, u123084, u123143, u123165, u123173, u123179, u123183, u123247, u123273, u123277, u123359, u123402, u123455, u123485, u123512, u123536, u123574, u123603, u123653, u123695, u123734, u123793, u123815, u123829, u123836, u123861, u123866, u123930, u123957, u124032, u124089, u124241, u124266, u124278, u124292, u124320, u124385, u124446, u124531, u124532, u124629, u124631, u124717, u124740, u124743, u124771, u124850, u124852, u124857, u124862, u124867, u124878, u124939, u124943, u124981, u125059, u125147, u125194, u125243, u125244, u125444, u125462, u125591, u125647, u125682, u125692, u125706, u125755, u125756, u125771, u125794, u125796, u125831, u125834, u125841, u125856, u125858, u125866, u125872, u125890, u125928, u125985, u126158, u126164, u126169, u126175, u126195, u126197, u126286, u126318, u126371, u126372, u126376, u126508, u126524, u126528, u126534, u126567, u126639, u126655, u126675, u126748, u126861, u126887, u127013, u127241, u127295, u127340, u127430, u127480, u127545, u127575, u127659, u127708, u127709, u127730, u127734, u127803, u127807, u127842, u127866, u127944, u128020, u128028, u128038, u128041, u128053, u128083, u128098, u128128, u128159, u128255, u128262, u128306, u128315, u128329, u128348, u128362, u128416, u128431, u128576, u128590, u128656, u128689, u128742, u128769, u128782, u128825, u128829, u128906, u128949, u128991, u129104, u129123, u129145, u129166, u129174, u129180, u129197, u129213, u129235, u129343, u129432, u129750, u129797, u129842, u129844, u129911, u129943, u129944, u130095, u130228, u130292, u130311, u130400, u130557, u130560, u130674, u130675, u130683, u130865, u130889, u130907, u130911, u130918, u130943, u130946, u131046, u131062, u131097, u131105, u131106, u131146, u131185, u131196, u131240, u131278, u131314, u131333, u131353, u131377, u131501, u131508, u131601, u131672, u131721, u131738, u131754, u131782, u131792, u131841, u131863, u131870, u131886, u131904, u131913, u131927, u132014, u132117, u132124, u132140, u132147, u132158, u132168, u132184, u132256, u132269, u132326, u132332, u132357, u132372, u132376, u132377, u132387, u132393, u132462, u132491, u132511, u132517, u132561, u132584, u132625, u132645, u132662, u132668, u132696, u132729, u132733, u132748, u132749, u132757, u132759, u132776, u132830, u132877, u132891, u132906, u132916, u132919, u132933, u132951, u132972, u133047, u133109, u133145, u133191, u133249, u133346, u133362, u133539, u133787, u133850, u133902, u133905, u133927, u133965, u134049, u134059, u134076, u134101, u134138, u134179, u134191, u134228, u134245, u134246, u134275, u134279, u134288, u134312, u134400, u134446, u134513, u134647, u134673, u134683, u134769, u134801, u134802, u134807, u134819, u134862, u134866, u134868, u134869, u134929, u135038, u135338, u135417, u135445, u135461, u135479, u135537, u135547, u135551, u135587, u135683, u135698, u135867, u136025, u136142, u136161, u136169, u136236, u136277, u136288, u136438, u136464, u136469, u136470, u136584, u136586, u136649, u136654, u136656, u136753, u136776, u136809, u136837, u136927, u136960, u136973, u137075, u137128, u137209, u137224, u137229, u137304, u137339, u137416, u137449, u137515, u137553, u137555, u137594, u137612, u137650, u137699, u137751, u137782, u137795, u137798, u137830, u137838, u137865, u137968, u138025, u138028, u138040, u138044, u138066, u138218, u138229, u138255, u138311, u138364, u138456, u138479, u138504, u138628, u138689, u138736, u138795, u138894, u138895, u139005, u139150, u139202, u139292, u139328, u139430, u139533, u139558, u139567, u139613, u139624, u139852, u139894, u139924, u139925, u139970, u139986, u140051, u140102, u140224, u140245, u140315, u140316, u140464, u140511, u140512, u140548, u140574, u140617, u140667, u140752, u140772, u140791, u140800, u140822, u140828, u140834, u140955, u140962, u140991, u141020, u141057, u141104, u141111, u141146, u141149, u141150, u141171, u141173, u141248, u141253, u141278, u141329, u141332, u141388, u141427, u141459, u141487, u141488, u141497, u141561, u141608, u141610, u141617, u141642, u141652, u141735, u141806, u141829, u141874, u141927, u141929, u141933, u141935, u141943, u141949, u141979, u141981, u142008, u142009, u142056, u142062, u142075, u142141, u142146, u142167, u142177, u142199, u142257, u142311, u142331, u142342, u142355, u142363, u142365, u142377, u142430, u142482, u142486, u142509, u142611, u142686, u142705, u142714, u142745, u142749, u142777, u142794, u142809, u142812, u142813, u142866, u142876, u142888, u142911, u142936, u142972, u142980, u143014, u143026, u143082, u143121, u143125, u143212, u143216, u143371, u143441, u143444, u143534, u143577, u143594, u143626, u143654, u143671, u143672, u143674, u143696, u143704, u143753, u143760, u143772, u143825, u143832, u143844, u143880, u143938, u143959, u143996, u144008, u144010, u144041, u144086, u144108, u144138, u144155, u144164, u144166, u144174, u144281, u144451, u144454, u144457, u144463, u144469, u144575, u144580, u144599, u144612, u144635, u144671, u144701, u144734, u144754, u144784, u144792, u144819, u144840, u145007, u145040, u145041, u145076, u145110, u145115, u145147, u145148, u145197, u145305, u145438, u145449, u145476, u145524, u145572, u145614, u145633, u145653, u145685, u145690, u145696, u145723, u145724, u145760, u145768, u145801, u145893, u145918, u145937, u145961, u145972, u145991, u146027, u146042, u146051, u146058, u146108, u146151, u146177, u146186, u146209, u146255, u146284, u146314, u146366, u146415, u146449, u146453, u146475, u146490, u146508, u146546, u146562, u146568, u146571, u146627, u146643, u146655, u146659, u146665, u146671, u146677, u146681, u146704, u146732, u146778, u146783, u146801, u146850, u146860, u146902, u146907, u146940, u146954, u147000, u147091, u147210, u147250, u147286, u147352, u147355, u147377, u147421, u147424, u147437, u147459, u147479, u147530, u147547, u147571, u147603, u147649, u147652, u147674, u147692, u147728, u147756, u147769, u147787, u147790, u147916, u147944, u147962, u147985, u147988, u148026, u148077, u148107, u148123, u148131, u148132, u148150, u148178, u148230, u148270, u148346, u148388, u148399, u148429, u148464, u148477, u148487, u148509, u148531, u148558, u148592, u148594, u148616, u148621, u148641, u148666, u148713, u148721, u148724, u148736, u148744, u148759, u148868, u148913, u148992, u149007, u149034, u149067, u149082, u149147, u149250, u149286, u149289, u149295, u149353, u149401, u149415, u149512, u149580, u149593, u149598, u149599, u149607, u149642, u149840, u149875, u149881, u149919, u149921, u149939, u149941, u149990, u150001, u150011, u150012, u150018, u150039, u150078, u150104, u150116, u150118, u150158, u150169, u150176, u150180, u150182, u150183, u150185, u150186, u150190, u150234, u150240, u150273, u150285, u150317, u150337, u150345, u150351, u150364, u150389, u150400, u150413, u150434, u150444, u150445, u150452, u150503, u150519, u150526, u150574, u150586, u150670, u150699, u150740, u150748, u150758, u150780, u150796, u150817, u150879, u150915, u150968, u150973, u150974, u150976, u150980, u150988, u151004, u151015, u151025, u151075, u151080, u151091, u151094, u151105, u151114, u151137, u151155, u151159, u151175, u151179, u151180, u151185, u151222, u151388, u151397, u151400, u151481, u151515, u151532, u151563, u151565, u151611, u151615, u151658, u151687, u151721, u151741, u151748, u151752, u151816, u151899, u151953, u152036, u152046, u152061, u152062, u152092, u152109, u152222, u152264, u152266, u152345, u152346, u152374, u152470, u152558, u152566, u152567, u152580, u152606, u152620, u152621, u152644, u152669, u152823, u152858, u152876, u152885, u152963, u152973, u153050, u153111, u153131, u153171, u153202, u153288, u153304, u153316, u153341, u153371, u153396, u153405, u153409, u153410, u153411, u153427, u153458, u153476, u153551, u153558, u153578, u153628, u153632, u153636, u153695, u153703, u153820, u153869, u153902, u153903, u153919, u154033, u154052, u154083, u154147, u154195, u154224, u154227, u154232, u154245, u154255, u154281, u154333, u154433, u154439, u154541, u154560, u154593, u154662, u154703, u154720, u154723, u154732, u154734, u154743, u154751, u154753, u154759, u154764, u154767, u154768, u154777, u154778, u154784, u154793, u154818, u154991, u155061, u155065, u155068, u155086, u155150, u155189, u155201, u155227, u155306, u155314, u155396, u155440, u155446, u155505, u155526, u155626, u155630, u155668, u155672, u155686, u155695, u155700, u155731, u155735, u155739, u155749, u155764, u155769, u155811, u155844, u155852, u155868, u155877, u155973, u156001, u156014, u156015, u156018, u156067, u156145, u156162, u156170, u156216, u156246, u156291, u156371, u156383, u156407, u156412, u156419, u156444, u156483, u156524, u156528, u156537, u156538, u156546, u156551, u156599, u156681, u156706, u156719, u156754, u156772, u156864, u156867, u156927, u156951, u156953, u156974, u156975, u157063, u157120, u157149, u157154, u157160, u157172, u157183, u157193, u157200, u157263, u157291, u157338, u157386, u157420, u157455, u157458, u157466, u157484, u157497, u157526, u157670, u157705, u157724, u157769, u157778, u157782, u157844, u157847, u157853, u157857, u157874, u157875, u157879, u157882, u157888, u157893, u157924, u157936, u158160, u158199, u158214, u158281, u158282, u158293, u158299, u158303, u158305, u158335, u158390, u158393, u158460, u158471, u158555, u158606, u158645, u158654, u158659, u158701, u158728, u158759, u158777, u158829, u158849, u158912, u158928, u158933, u159003, u159069, u159200, u159250, u159349, u159382, u159391, u159415, u159458, u159460, u159489, u159525, u159526, u159550, u159553, u159557, u159558, u159561, u159579, u159584, u159585, u159606, u159609, u159611, u159616, u159671, u159678, u159729, u159792, u159815, u159850, u159870, u159877, u159902, u159937, u159938, u159939, u159948, u159956, u160044, u160062, u160076, u160116, u160118, u160123, u160198, u160201, u160262, u160263, u160339, u160350, u160357, u160364, u160367, u160369, u160371, u160391, u160393, u160403, u160422, u160467, u160497, u160502, u160516, u160541, u160561, u160592, u160594, u160607, u160614, u160637, u160652, u160680, u160694, u160712, u160782, u160805, u160836, u160932, u160936, u160951, u160972, u161031, u161047, u161054, u161078, u161084, u161088, u161102, u161135, u161142, u161166, u161168, u161186, u161216, u161228, u161229, u161267, u161295, u161313, u161320, u161321, u161361, u161451, u161499, u161524, u161531, u161544, u161565, u161569, u161598, u161617, u161684, u161700, u161717, u161792, u161805, u161934, u161936, u161950, u161999, u162055, u162070, u162075, u162080, u162081, u162127, u162163, u162169, u162196, u162284, u162294, u162298, u162304, u162320, u162325, u162358, u162378, u162394, u162463, u162544, u162561, u162602, u162606, u162613, u162624, u162630, u162654, u162655, u162658, u162671, u162681, u162684, u162699, u162739, u162742, u162743, u162754, u162839, u162880, u162891, u162933, u162953, u162963, u162975, u163022, u163045, u163049, u163070, u163168, u163181, u163189, u163224, u163301, u163341, u163595, u163604, u163636, u163668, u163686, u163702, u163710, u163759, u163777, u163792, u163878, u163879, u164010, u164099, u164113, u164138, u164189, u164221, u164360, u165239, u165261, u165326, u165411, u165449, u165460, u165462, u165476, u165492, u165655, u165721, u165820, u165873, u165895, u165919, u165924, u165931, u165937, u165951, u166004, u166032, u166057, u166074, u166120, u166121, u166164, u166167, u166168, u166236, u166251, u166257, u166273, u166377, u166440, u166481, u166486, u166520, u166577, u166651, u166685, u166715, u166790, u166819, u166917, u166932, u166941, u166950, u166959, u166963, u167009, u167019, u167061, u167076, u167206, u167240, u167283, u167293, u167358, u167388, u167436, u167627, u167676, u167684, u167714, u167715, u167739, u167901, u167941, u168089, u168157, u168216, u168421, u168436, u168592, u168616, u168623, u168741, u168825, u168859, u168875, u168909, u168919, u169209, u169265, u169360, u169379, u169515, u169517, u169535, u169555, u169565, u169578, u169595, u169626, u169643, u169646, u169664, u169721, u169756, u169794, u169823, u169931, u170800, u170832, u171050, u171458, u171469, u171497, u171536, u172027, u172059, u172151, u172265, u172278, u172307, u172324, u172332, u172451, u172641, u172669, u172679, u172785, u172795, u172844, u172969, u173006, u173064, u173071, u173110, u173134, u173170, u173203, u173206, u173232, u173325, u173332, u173334, u173339, u173416, u173449, u173461, u173472, u173509, u173533, u173577, u173664, u173672, u173678, u173776, u173778, u173883, u173991, u174297, u174298, u174700, u175128, u175173, u175203, u175334, u175420, u175451, u176079, u176889, u176926, u176947, u176954, u176999, u177019, u177027, u177053, u177056, u177066, u177074, u177124, u177152, u177200, u177245, u177265, u177278, u177279, u177280, u177293, u177297, u177304, u177305, u177331, u177338, u177345, u177367, u177402, u177430, u177441, u177444, u177482, u177624, u177654, u177776, u177780, u177893, u177894, u177919, u177995, u177997, u178008, u178013, u178042, u178045, u178057, u178093, u178167, u178178, u178197, u178310, u178322, u178334, u178401, u178425, u178437, u178439, u178442, u178461, u178465, u178474, u178485, u178493, u178499, u178504, u178520, u178523, u178577, u178611, u178623, u178626, u178654, u178686, u178693, u178773, u179004, u179033, u179259, u179298, u179322, u179348, u179391, u179435, u179474, u179487, u179500, u179532, u179580, u179590, u179682, u179704, u179764, u179767, u179772, u179789, u179911, u179912, u179921, u179932, u179998, u180096, u180113, u180135, u180182, u180402, u180406, u180416, u180441, u180529, u180614, u180636, u180639, u180642, u180669, u180674, u180681, u180693, u180739, u180766, u180906, u180912, u180920, u180927, u181001, u181006, u181032, u181036, u181038, u181052, u181093, u181113, u181156, u181171, u181181, u181183, u181184, u181188, u181195, u181208, u181234, u181240, u181261, u181343, u181435, u181467, u181519, u181538, u181572, u181585, u181594, u181634, u181675, u181677, u181683, u181720, u181725, u181821, u181838, u181848, u181892, u181900, u181907, u181908, u181926, u181938, u181940, u181961, u182056, u182110, u182175, u182233, u182237, u182238, u182248, u182264, u182266, u182280, u182283, u182289, u182309, u182326, u182328, u182342, u182343, u182388, u182410, u182413, u182417, u182451, u182476, u182502, u182544, u182551, u182552, u182553, u182555, u182559, u182560, u182567, u182584, u182596, u182609, u182611, u182639, u182642, u182663, u182702, u182716, u182728, u182736, u182764, u182798, u182934, u182957, u182987, u182994, u182999, u183024, u183032, u183045, u183067, u183081, u183082, u183106, u183108, u183119, u183132, u183149, u183157, u183167, u183198, u183237, u183361, u183395, u183477, u183481, u183502, u183509, u183515, u183525, u183545, u183562, u183581, u183663, u183668, u183733, u183766, u183785, u183853, u183864, u183869, u183874, u183911, u183974, u184058, u184136, u184161, u184180, u184214, u184231, u184315, u184387, u184393, u184398, u184438, u184439, u184491, u184513, u184517, u184519, u184538, u184547, u184548, u184558, u184566, u184574, u184603, u184618, u184632, u184654, u184665, u184666, u184674, u184684, u184703, u184730, u184735, u184742, u184877, u185053, u185082, u185117, u185256, u185286, u185340, u185359, u185378, u185480, u185581, u185582, u185604, u185623, u185626, u185630, u185640, u185641, u185667, u185669, u185779, u185800, u185807, u185871, u185898, u185944, u18595, u185964, u185997, u186013, u186035, u186041, u186045, u186046, u186074, u186088, u186111, u186141, u186144, u186147, u186165, u186168, u186280, u186299, u186305, u186316, u186323, u186326, u186350, u186356, u186359, u186362, u186404, u186410, u186427, u186448, u186451, u186466, u186468, u186471, u186529, u186595, u186600, u186617, u186621, u186626, u186643, u186646, u186696, u186801, u186877, u186905, u186927, u186930, u186960, u186987, u187089, u187152, u187155, u187163, u187177, u187190, u187274, u187313, u187608, u187640, u187746, u187880, u187908, u187976, u188159, u188256, u188262, u188268, u188286, u188354, u188358, u188360, u188374, u188390, u188391, u188405, u188461, u188517, u188551, u188567, u188597, u188666, u188679, u188705, u188707, u188762, u188797, u188808, u188877, u188898, u188937, u188939, u188950, u188962, u188968, u189042, u189070, u189092, u189098, u189131, u189142, u189180, u189196, u189203, u189238, u189240, u189247, u189270, u189279, u189285, u189299, u189308, u189356, u189402, u189413, u189415, u189416, u189436, u189504, u189540, u189550, u189583, u189589, u189677, u189691, u189725, u189744, u189810, u189817, u189867, u189871, u189904, u189907, u189910, u189922, u189960, u190008, u190108, u190133, u190186, u190202, u190237, u190327, u190335, u190336, u190338, u190343, u190345, u190380, u190403, u190486, u190500, u190556, u190627, u190678, u190717, u190745, u190790, u190849, u190857, u190895, u190912, u190957, u190965, u191073, u191076, u191101, u191248, u191358, u191371, u191435, u191454, u191496, u191627, u19168, u191770, u191778, u191779, u191876, u191882, u191890, u191891, u191918, u191930, u191938, u192017, u192078, u192130, u192145, u192222, u192250, u192271, u192276, u192287, u192362, u192419, u192535, u192545, u192562, u192782, u192784, u192844, u192845, u192857, u192875, u192882, u192888, u19290, u192971, u193277, u193350, u193384, u193401, u193405, u193462, u193490, u193513, u193544, u193748, u193774, u193775, u193787, u193801, u193824, u193828, u193843, u193848, u193861, u193913, u193972, u193987, u194052, u194086, u194094, u194107, u194127, u194153, u194180, u194229, u194232, u194233, u194234, u194255, u194267, u194271, u194283, u194291, u194302, u194308, u194341, u194420, u194427, u194454, u194455, u194456, u194459, u194461, u194507, u194549, u194550, u194562, u194599, u194634, u194638, u194684, u194692, u194695, u194704, u194717, u194782, u194816, u194836, u194892, u194949, u194986, u194993, u195146, u195159, u195196, u195227, u195236, u195259, u195330, u195375, u195405, u195562, u195570, u195577, u195630, u195651, u195655, u195671, u195679, u195707, u195823, u195848, u195854, u195932, u195947, u195999, u196038, u196118, u196120, u196175, u196199, u196279, u196289, u196371, u196471, u196499, u196525, u196590, u196597, u196598, u196619, u196688, u196718, u196963, u197059, u197061, u197088, u197093, u197192, u197228, u197229, u197251, u197270, u197294, u197310, u197361, u197447, u197454, u197675, u197748, u197783, u197907, u197989, u198082, u198105, u198526, u198558, u198695, u198701, u198730, u198796, u198869, u198875, u199121, u199153, u199185, u199257, u199261, u199270, u199278, u199295, u199348, u199393, u199558, u199582, u199652, u199663, u199678, u199685, u199726, u199976, u200015, u200022, u200030, u200049, u200230, u200238, u200256, u200267, u200296, u200333, u200336, u200341, u200375, u200378, u200564, u200578, u200584, u200589, u200593, u200650, u200678, u200703, u200705, u200723, u200770, u200773, u200786, u200789, u200824, u200825, u200863, u200870, u200878, u200899, u200962, u200967, u200972, u200998, u201045, u201048, u201080, u201083, u201096, u201118, u201213, u201220, u201225, u201238, u201244, u201250, u201266, u201279, u201283, u201327, u201329, u201335, u201343, u201357, u201367, u201372, u201373, u201444, u201491, u201503, u201506, u201507, u201531, u201542, u201589, u201667, u201781, u201824, u201867, u201880, u201892, u201958, u201992, u202049, u202054, u202064, u202109, u202123, u202125, u202138, u202169, u202176, u202178, u202183, u202185, u202201, u202217, u202224, u202226, u202239, u202247, u202248, u202253, u202322, u202354, u202359, u202374, u202381, u202383, u202391, u202401, u202405, u202448, u202456, u202531, u202604, u202620, u202623, u202681, u202721, u202732, u202768, u202772, u202778, u202783, u202978, u203000, u203005, u203019, u203033, u203034, u203035, u203090, u203118, u203119, u203126, u203137, u203175, u203176, u203203, u203242, u203257, u203320, u203331, u203337, u203356, u203368, u203404, u203452, u203580, u203632, u203653, u203701, u203713, u203756, u203795, u203797, u203879, u203890, u203953, u203956, u204004, u204027, u204046, u204050, u204062, u204082, u204092, u204093, u204117, u204118, u204180, u204193, u204208, u204229, u204272, u204296, u204303, u204322, u204334, u204335, u204348, u204355, u204357, u204379, u204382, u204390, u204394, u204404, u204435, u204469, u204470, u204472, u204480, u204524, u204565, u204589, u204691, u204708, u20479, u204793, u204817, u204884, u204898, u204905, u204909, u204913, u204914, u205020, u205049, u205053, u205055, u205068, u205117, u205122, u205127, u205152, u205162, u205219, u205228, u205267, u205269, u205302, u205305, u205309, u205310, u205333, u205350, u205383, u205406, u205441, u205475, u205489, u205494, u205499, u205525, u205547, u205550, u205553, u205574, u205584, u205595, u205619, u205623, u205753, u205784, u205793, u205812, u205813, u205821, u205861, u205894, u205902, u205927, u205998, u206006, u206010, u206039, u206083, u206109, u206114, u206166, u206202, u206203, u206205, u206245, u206279, u206293, u206305, u206314, u206320, u206382, u206408, u206416, u206433, u206446, u206448, u206457, u206483, u206490, u206503, u206558, u206647, u206654, u206659, u206676, u206696, u206698, u206752, u206762, u206789, u206794, u206802, u206821, u206828, u206829, u206830, u206831, u206885, u206949, u206967, u206982, u206987, u207007, u207012, u207082, u207130, u207135, u207162, u207203, u207215, u207238, u207278, u207290, u207291, u207349, u207383, u207389, u207418, u207423, u207452, u207495, u207502, u207516, u207517, u207552, u207562, u207605, u207606, u207609, u207613, u207616, u207642, u207643, u207651, u207654, u207672, u207674, u207676, u207677, u207684, u207693, u207700, u207706, u207749, u207750, u207756, u207781, u207800, u207808, u207811, u207830, u207831, u207839, u20784, u207908, u207923, u207955, u207973, u207991, u207993, u208007, u208017, u208019, u208020, u208022, u208036, u208079, u208093, u208103, u208118, u208157, u208193, u208237, u208239, u208262, u208275, u208297, u208298, u208303, u208329, u208368, u208378, u208418, u208424, u208446, u208464, u208466, u208495, u208502, u208541, u208546, u208574, u208577, u208589, u208610, u208611, u208615, u208619, u208623, u208648, u208717, u208728, u208741, u208761, u208762, u208764, u208768, u208779, u208780, u208790, u208798, u208801, u208808, u208843, u208871, u208875, u208898, u208900, u208973, u209025, u209027, u209049, u209067, u209078, u209108, u209122, u209155, u209171, u209214, u209223, u209244, u209257, u209306, u209307, u209324, u209343, u209345, u209349, u209368, u209371, u209401, u209408, u209412, u209448, u209473, u209522, u209523, u209552, u209564, u209577, u209587, u209595, u209607, u209649, u209659, u209710, u209724, u209736, u209750, u209751, u209758, u209759, u209767, u209806, u209826, u209830, u209839, u209864, u209884, u209922, u209940, u209983, u210001, u210015, u210074, u210095, u210108, u210118, u210128, u210142, u210152, u210170, u210204, u210206, u210222, u210233, u210244, u210352, u210398, u210472, u210510, u210545, u210571, u210574, u210580, u210587, u210598, u210601, u210606, u210610, u210618, u210620, u210622, u210644, u210655, u210752, u210759, u210768, u210789, u210811, u210824, u210827, u210835, u210860, u210938, u210940, u210947, u211002, u211032, u211038, u211058, u211063, u211066, u211068, u211084, u211085, u211088, u211120, u211138, u211156, u211202, u211233, u211238, u211270, u211317, u211356, u211357, u211365, u211373, u211461, u211473, u211491, u211493, u211509, u211518, u211524, u211538, u211600, u211630, u211631, u211766, u211818, u211847, u211848, u211910, u211939, u211986, u212052, u212056, u212078, u212108, u212130, u212282, u212304, u212316, u212319, u212331, u212335, u212362, u212418, u212522, u212592, u212607, u212628, u212634, u212655, u212662, u212690, u212770, u212803, u212811, u212821, u212851, u212867, u212879, u212913, u212917, u212950, u212954, u212973, u212979, u212996, u213025, u213030, u213041, u213056, u213080, u213120, u213126, u213152, u213188, u213226, u213234, u213330, u213350, u213376, u213397, u213408, u213429, u213433, u213442, u213453, u213507, u213511, u213515, u213526, u213570, u213600, u213611, u213646, u213655, u213661, u213701, u213741, u213745, u213757, u213767, u213773, u213780, u213784, u213806, u213809, u213828, u213854, u213862, u213868, u213889, u213903, u213921, u213937, u213939, u213951, u213956, u213963, u214029, u214053, u214059, u214081, u214093, u214101, u214110, u214124, u214136, u214146, u214156, u214157, u214159, u214165, u214175, u214179, u214186, u214327, u214340, u214397, u214402, u214404, u214428, u214434, u214441, u214458, u214463, u214473, u214489, u214555, u214565, u214578, u214602, u214621, u214705, u214721, u214752, u214816, u214861, u214887, u214924, u214949, u214961, u214978, u214981, u214983, u214996, u215007, u215010, u215036, u215099, u215153, u215154, u215157, u215180, u215192, u215202, u215210, u215211, u215285, u215318, u215403, u215443, u215452, u215485, u215546, u215552, u215579, u215580, u215604, u215611, u215817, u215957, u215965, u216118, u216123, u216172, u216244, u216262, u216297, u216319, u216321, u216325, u216333, u216375, u216384, u216445, u216447, u216596, u216621, u216729, u216733, u216755, u216776, u216810, u216829, u216968, u216990, u217011, u217192, u217207, u217530, u217545, u217809, u217905, u217912, u217925, u217973, u218030, u218178, u218236, u218284, u218331, u218336, u218339, u218438, u218525, u218532, u218749, u218758, u218934, u219033, u219076, u219137, u219188, u219194, u219209, u219210, u219249, u219331, u219332, u219353, u219361, u219726, u220101, u220120, u220136, u220165, u220188, u220301, u220339, u220348, u220354, u220373, u220423, u220434, u220438, u220455, u220457, u220486, u220506, u220618, u220621, u220649, u220654, u220659, u220664, u220680, u220696, u220715, u220720, u220742, u220747, u220767, u220786, u220795, u220801, u220807, u220840, u220897, u220902, u220924, u220929, u220976, u220985, u221003, u221051, u221080, u221119, u221152, u221170, u221179, u221211, u221272, u221285, u221308, u221333, u221374, u221394, u221457, u221520, u221536, u221541, u221567, u221631, u221652, u221729, u221742, u221766, u221780, u221793, u221871, u221892, u222007, u222010, u222030, u222060, u222062, u222065, u222077, u222116, u222120, u222212, u222273, u222310, u222311, u222341, u222348, u222371, u222381, u222385, u222388, u222401, u222419, u222428, u222478, u222490, u222527, u222556, u222648, u222693, u222787, u222807, u222822, u222831, u222894, u222964, u222978, u222990, u223026, u223069, u223075, u223122, u223138, u223166, u223180, u223233, u223337, u223351, u223357, u223371, u223375, u223424, u223428, u223463, u223492, u223498, u223547, u223598, u223599, u223603, u223610, u223621, u223622, u223628, u223638, u223651, u223653, u223658, u223685, u223705, u223716, u223729, u223737, u223766, u223845, u223866, u223879, u223885, u223980, u223995, u224131, u224167, u224246, u224258, u224309, u224335, u224451, u224453, u224482, u224597, u224637, u224687, u224693, u224705, u224708, u224734, u224742, u224761, u224766, u224789, u224793, u224811, u224852, u224865, u224868, u224876, u224879, u224880, u224881, u224884, u224897, u224902, u224926, u224929, u224949, u225028, u225037, u225045, u225054, u225094, u225098, u225099, u225105, u225114, u225115, u225116, u225117, u225150, u225204, u225259, u225282, u225292, u225304, u225317, u225336, u225342, u225430, u225436, u225438, u225442, u225443, u225527, u225542, u225549, u225582, u225593, u225617, u225618, u225624, u225643, u225652, u225704, u225706, u225772, u225819, u225829, u225848, u225970, u225978, u226027, u226043, u226046, u226049, u226090, u226129, u226133, u226173, u226196, u226203, u226221, u226222, u226231, u226237, u226245, u226296, u226319, u226431, u226435, u226486, u226509, u226557, u226620, u226640, u226645, u226651, u226698, u226702, u22675, u226751, u226778, u226868, u226926, u226933, u226970, u226980, u226993, u227018, u227060, u227068, u227069, u227075, u227088, u227101, u227149, u227152, u227153, u227159, u227163, u227186, u227197, u227201, u227211, u227234, u227246, u227258, u227301, u227374, u227382, u227391, u227433, u227450, u227468, u227508, u227523, u227525, u227546, u227553, u227575, u227605, u227627, u227659, u227664, u227686, u227721, u227811, u227822, u227830, u227838, u227860, u227864, u227883, u227922, u227944, u227961, u227967, u227984, u228096, u228158, u228233, u228265, u228377, u228426, u228487, u228528, u228543, u228546, u228582, u228588, u228589, u228594, u228671, u228733, u228749, u228752, u228757, u228758, u228761, u228767, u228797, u228821, u228830, u228835, u228916, u228917, u228933, u228939, u228948, u228955, u228958, u228960, u228961, u228965, u228972, u228993, u229023, u229034, u229037, u229095, u229106, u229123, u229133, u229136, u229158, u229172, u229204, u229269, u229301, u229317, u229332, u229354, u229366, u229369, u229409, u229419, u229421, u229433, u229441, u229513, u229539, u229556, u229566, u229577, u229624, u229629, u229711, u229720, u229754, u229763, u229764, u229780, u229806, u229825, u229829, u229876, u229910, u229921, u229965, u229975, u229978, u229985, u230021, u230114, u230190, u230200, u230201, u230229, u230254, u230257, u230323, u23033, u230473, u230486, u230544, u230603, u230645, u230649, u230666, u230749, u230759, u230772, u230842, u230919, u230932, u230946, u230976, u230980, u231020, u231062, u231104, u231111, u231127, u231148, u231158, u231170, u231181, u231200, u231203, u231286, u231323, u231362, u231429, u231443, u231499, u231502, u231520, u231589, u231632, u231641, u231649, u231713, u231722, u231759, u231760, u231778, u231824, u231828, u231837, u231841, u231866, u231891, u231899, u231964, u231966, u231979, u231986, u231992, u232000, u232002, u232024, u232108, u232216, u232250, u232319, u232321, u232323, u232353, u232358, u232545, u232559, u232590, u232650, u232658, u232661, u232662, u232664, u232676, u232691, u232745, u232767, u232780, u23279, u232801, u232812, u232882, u232883, u232884, u232935, u232987, u233000, u233002, u233010, u233059, u233108, u233155, u233160, u233161, u233163, u233198, u233199, u233203, u233214, u233224, u233240, u233251, u233283, u233291, u233300, u233328, u233342, u233361, u233396, u233437, u233458, u233468, u233553, u233568, u233592, u233600, u233612, u233624, u233753, u233754, u233787, u233803, u233813, u233821, u233875, u233883, u233906, u233951, u233959, u234025, u234045, u234137, u234191, u234229, u234317, u234335, u234365, u234367, u234402, u234409, u234463, u234486, u234499, u234561, u234563, u234601, u234667, u234675, u234676, u234721, u234770, u234784, u234787, u234846, u234883, u235070, u235340, u235525, u235530, u235623, u235626, u235658, u235659, u235684, u235970, u236077, u236416, u236433, u236770, u237354, u237403, u237482, u237510, u237557, u237604, u237625, u237635, u237641, u237645, u237672, u237716, u237722, u237725, u237786, u237836, u237842, u237851, u237876, u237901, u237912, u237938, u237989, u237998, u238031, u238040, u238045, u238050, u238106, u238116, u238147, u238158, u238228, u238231, u238273, u238322, u238347, u238353, u238372, u238376, u238486, u238496, u238525, u238542, u238558, u238571, u238572, u238579, u238581, u238602, u238649, u238661, u238683, u238703, u238738, u238748, u238761, u238798, u238811, u238838, u238860, u238912, u238986, u238992, u239067, u239068, u239178, u239361, u239412, u239429, u239449, u239510, u239511, u239516, u239552, u239582, u239655, u239794, u239827, u239831, u239981, u239986, u240024, u240057, u240117, u240119, u240121, u240149, u240160, u240163, u240183, u240207, u240234, u240260, u240272, u240344, u240369, u240397, u240401, u240443, u240452, u240493, u240516, u240519, u240521, u240524, u240527, u240531, u240543, u240553, u240573, u24059, u240611, u240642, u240697, u240756, u240798, u240827, u240871, u240874, u240875, u240877, u240880, u240888, u240893, u240900, u240910, u240929, u240946, u240980, u240984, u240993, u241011, u241058, u241062, u241083, u241084, u241107, u241139, u241152, u241222, u241246, u241289, u241342, u241345, u241353, u241356, u241365, u241371, u241390, u241398, u241453, u241582, u241645, u241686, u241731, u241733, u241744, u241752, u241779, u241787, u241800, u241801, u241802, u241818, u241870, u241909, u241975, u241976, u242029, u242036, u242054, u242072, u242079, u242161, u242208, u242227, u242230, u242237, u242276, u242278, u242352, u242473, u242524, u242584, u242631, u242736, u242842, u242859, u242890, u242950, u243000, u243017, u243046, u243095, u243138, u243139, u243281, u243298, u243325, u243413, u243422, u243463, u243489, u243520, u243532, u243538, u243575, u243584, u243631, u243639, u243653, u243661, u243668, u243688, u243754, u243786, u243789, u243790, u243800, u243804, u243807, u243811, u243832, u243899, u243928, u243950, u243956, u243962, u243976, u243987, u243990, u243998, u243999, u244005, u244010, u244038, u24406, u244086, u244212, u24423, u244252, u244312, u244330, u244370, u244389, u244396, u244401, u244445, u244449, u244485, u244529, u244635, u244661, u244672, u244697, u244713, u244749, u244805, u244812, u244814, u244834, u244847, u244873, u244915, u244950, u244961, u244965, u245054, u245117, u245154, u245186, u245201, u245209, u245319, u245366, u245386, u245406, u245433, u245451, u245458, u245529, u245637, u245646, u245647, u245681, u245682, u245716, u245755, u245836, u245859, u245868, u245876, u245898, u245947, u245962, u245991, u246001, u246009, u246216, u246273, u246274, u246281, u246328, u246337, u246344, u246472, u246482, u246516, u246517, u246544, u246571, u246605, u246679, u246706, u246707, u246708, u246734, u246735, u246748, u246826, u246859, u246864, u246875, u246898, u246899, u246903, u246908, u246975, u246992, u247090, u247093, u247132, u247157, u247162, u247170, u247191, u247206, u247237, u247307, u247384, u247448, u247452, u247493, u247676, u247688, u247692, u247726, u247728, u247805, u247812, u247835, u247842, u247876, u247902, u247930, u247952, u247953, u247964, u248005, u248036, u248044, u248069, u248079, u248096, u248175, u248206, u248221, u248226, u248229, u248235, u248248, u248282, u248317, u248366, u248411, u248449, u248463, u248471, u248552, u248589, u248629, u248690, u248715, u248769, u248805, u248835, u248850, u248889, u248899, u248924, u248946, u249015, u249051, u249098, u249105, u249110, u249125, u249128, u249138, u249142, u249167, u249184, u249237, u249366, u249397, u249401, u249414, u249416, u249469, u249472, u249497, u249502, u249527, u249535, u249547, u249550, u249566, u249602, u249606, u249617, u249630, u249640, u249648, u249661, u249673, u249722, u249741, u249813, u249907, u249946, u250029, u250075, u250251, u250346, u250369, u250427, u250431, u250441, u250473, u250561, u250587, u250614, u250626, u250630, u250751, u250850, u250903, u250916, u250929, u251021, u251029, u251074, u251147, u251175, u251206, u251243, u251360, u251440, u251492, u251508, u251552, u251574, u251813, u251928, u252009, u252018, u252020, u252041, u252087, u252123, u252127, u252134, u252164, u252183, u252184, u252200, u252313, u252385, u252407, u252423, u252480, u252495, u252515, u252541, u252559, u252595, u252623, u252655, u252661, u252666, u252682, u252756, u252816, u252948, u252949, u252964, u253052, u253111, u253304, u253340, u253385, u253395, u253413, u253500, u253506, u253520, u253595, u253721, u253825, u253829, u253864, u254492, u254573, u254633, u254653, u254666, u254687, u254691, u254705, u254707, u254712, u254728, u254734, u254752, u254802, u254803, u254894, u254928, u254945, u254964, u255006, u255029, u255092, u255155, u255163, u255179, u255185, u255240, u255249, u255252, u255281, u255286, u255309, u255310, u255313, u255322, u255324, u255353, u255364, u255367, u255368, u255380, u255448, u255514, u255529, u255553, u255558, u255565, u255580, u255586, u255606, u255611, u255618, u255645, u255671, u255720, u255751, u255790, u255826, u255850, u255858, u255874, u255889, u255891, u255913, u255933, u255945, u255953, u255979, u255986, u256030, u256049, u256050, u256122, u256130, u256145, u256209, u256224, u256238, u256250, u256282, u256315, u256354, u256375, u256397, u256398, u256421, u256443, u256451, u256464, u256468, u256517, u256540, u256558, u256639, u256670, u256671, u256681, u256693, u256720, u256756, u256810, u256903, u256943, u257044, u257094, u257116, u257145, u257172, u257175, u257201, u257252, u257359, u257364, u257427, u257474, u257518, u257597, u257620, u257628, u257636, u257638, u257676, u257684, u257760, u257789, u257928, u257946, u257973, u257977, u257978, u257995, u258003, u258020, u258021, u258025, u258104, u258153, u258201, u258225, u258231, u258234, u258243, u258278, u258299, u258332, u258336, u258340, u258394, u258424, u258445, u258452, u258512, u258591, u258731, u258732, u258744, u258754, u258796, u258892, u259025, u259027, u259044, u259064, u259074, u259084, u259085, u259095, u259133, u259137, u259155, u259156, u259162, u259179, u259180, u259237, u259244, u259250, u259329, u259358, u259359, u259396, u259412, u259425, u259440, u259452, u259453, u259460, u259473, u259476, u259478, u259500, u259514, u259570, u259715, u259737, u259744, u259759, u259761, u259767, u259774, u259775, u259783, u259789, u259805, u259825, u259832, u259901, u259914, u259937, u259944, u259990, u259993, u260030, u260031, u260108, u260145, u260150, u260181, u260216, u260283, u260332, u260339, u260340, u260365, u260417, u260421, u260427, u260443, u260508, u260514, u260525, u260533, u260588, u260592, u260594, u260597, u260603, u260618, u260640, u260645, u260734, u260740, u260774, u260775, u260777, u260787, u260838, u260869, u260871, u260875, u260895, u260907, u260930, u260948, u260984, u261010, u261041, u261043, u261075, u261082, u261092, u261111, u261114, u261119, u261131, u261137, u261149, u261173, u261193, u261216, u261218, u261243, u261263, u261272, u261281, u261301, u261320, u261336, u261342, u261358, u261361, u261363, u261399, u261473, u261494, u261495, u261503, u261506, u261508, u261517, u261529, u261559, u261588, u261590, u261607, u261610, u261618, u261691, u261716, u261896, u261902, u261903, u261910, u261927, u261950, u261952, u261970, u261973, u261994, u262004, u262006, u262022, u262034, u262035, u262042, u262049, u262055, u262082, u262084, u262090, u262112, u262114, u262116, u262128, u262149, u262221, u262241, u262266, u262272, u262277, u262312, u262341, u262343, u262349, u262354, u262373, u262413, u262418, u262424, u262427, u262437, u262441, u262462, u262467, u262540, u262565, u262601, u262627, u262646, u262664, u262676, u262695, u262698, u262733, u262760, u262772, u262781, u262786, u262787, u262790, u262839, u262875, u262924, u262985, u262989, u262993, u263020, u263021, u263044, u263070, u263075, u263081, u263103, u263113, u263120, u263126, u263140, u263196, u263233, u263236, u263277, u263312, u263315, u263319, u263326, u263394, u263396, u263404, u263431, u263445, u263451, u263514, u263551, u263625, u263658, u263668, u263670, u263672, u263674, u263675, u263756, u263757, u263760, u263773, u263781, u263793, u263802, u263823, u263848, u263849, u263867, u263886, u263903, u263905, u263916, u263927, u263936, u263939, u263981, u264072, u264089, u264096, u264133, u264171, u264213, u264243, u264275, u264277, u264286, u264292, u264330, u264377, u264388, u264415, u264475, u264496, u264514, u264557, u264588, u264666, u264708, u264717, u264724, u264834, u265003, u265026, u265052, u265090, u265201, u265250, u265526, u265540, u265708, u265782, u26579, u265812, u265837, u265893, u265922, u265932, u265967, u266043, u266095, u266245, u266269, u266304, u266467, u266549, u26656, u266681, u266698, u266738, u266856, u266897, u266908, u266970, u266999, u267063, u267064, u267150, u267204, u267211, u267216, u267221, u267227, u267254, u267256, u267279, u267287, u267290, u267316, u267334, u267349, u267435, u267438, u267440, u267451, u267465, u267466, u267483, u267484, u267485, u267507, u267549, u267599, u267602, u267603, u267626, u267628, u267657, u267662, u267675, u267688, u267702, u267707, u267743, u267752, u267790, u267794, u267835, u267849, u267967, u267975, u268000, u268001, u268023, u268044, u268057, u268064, u268091, u268107, u268228, u268229, u268271, u268273, u268343, u268372, u268377, u268390, u268443, u268472, u268548, u268575, u268576, u268655, u268686, u268689, u268728, u268741, u268820, u268844, u268901, u268904, u268910, u268942, u268972, u268986, u268999, u269060, u269065, u269075, u269076, u269083, u269096, u269102, u269125, u269178, u269201, u269249, u269257, u269263, u269272, u269277, u269329, u269962, u270329, u270396, u270626, u270630, u27070, u271054, u27110, u271178, u271253, u271643, u271665, u271692, u271792, u271812, u271813, u271856, u271951, u272101, u272142, u272146, u272155, u272414, u272422, u272424, u272436, u272439, u272508, u272524, u272603, u272628, u272642, u272649, u272669, u272702, u272709, u272740, u272744, u272810, u272815, u272832, u272834, u273113, u273401, u273456, u273465, u273470, u27349, u273615, u273699, u273747, u273809, u273908, u273917, u273972, u27399, u27403, u274219, u274220, u274230, u274238, u274244, u274302, u274370, u274384, u274389, u274397, u274403, u274418, u274425, u274456, u274469, u274512, u274524, u274527, u274530, u274682, u27470, u274718, u274867, u274890, u274947, u274954, u274955, u274956, u274971, u275026, u275039, u275059, u275166, u275181, u275208, u275242, u275302, u275357, u275392, u275437, u275451, u275466, u275491, u275519, u275620, u275629, u275777, u275808, u275837, u275841, u275849, u275897, u275946, u275960, u275980, u276001, u276035, u276046, u276071, u276072, u276078, u276093, u276122, u276127, u276132, u276136, u276253, u276262, u276273, u276277, u276278, u276295, u276314, u276336, u276455, u276464, u276473, u276480, u276483, u276490, u276512, u276524, u276655, u276737, u276774, u277371, u277453, u277508, u277528, u277584, u277618, u277709, u277724, u277786, u277825, u277886, u277897, u277939, u277941, u278015, u278144, u278145, u278178, u278195, u278208, u278233, u278261, u278355, u278400, u278441, u278465, u278466, u278567, u278596, u278692, u278706, u278717, u278834, u27889, u278915, u278979, u278981, u278984, u279031, u279179, u279186, u279191, u27926, u279307, u279310, u279437, u279468, u279565, u279642, u279646, u279671, u279679, u279680, u279728, u279731, u279747, u279776, u279798, u279802, u279831, u279837, u279838, u279859, u279874, u279889, u279905, u279932, u279953, u279959, u279963, u279984, u280006, u280030, u280068, u280079, u280132, u280135, u280154, u280173, u280200, u280234, u280241, u280286, u280296, u280309, u280327, u280335, u280353, u280360, u280374, u280391, u280403, u280475, u280478, u280559, u280571, u280598, u280628, u280644, u280665, u280687, u280700, u280729, u280750, u280766, u280782, u280814, u280843, u280863, u280877, u280891, u280959, u280966, u280972, u280973, u281013, u281016, u281024, u281031, u281139, u281230, u281250, u281262, u281279, u281288, u281381, u281399, u281408, u281419, u281422, u281475, u281528, u281531, u281544, u281559, u281585, u281587, u281615, u281802, u281970, u281975, u282059, u282068, u282089, u282119, u28215, u282176, u282251, u28235, u282359, u282360, u282382, u282409, u282513, u282517, u282576, u282684, u282701, u282711, u282726, u282854, u282862, u282901, u282974, u283053, u283064, u283127, u283209, u283226, u283236, u283246, u283248, u283258, u283271, u283282, u283298, u283316, u283327, u283328, u283406, u283412, u283429, u283453, u283459, u283628, u283640, u283669, u283678, u283683, u283758, u283772, u283773, u283778, u283820, u283868, u283954, u283978, u283992, u283994, u284059, u284071, u28412, u284138, u284233, u284257, u284272, u284503, u284505, u284511, u284669, u28469, u284710, u284818, u28483, u284866, u284872, u284890, u284927, u285012, u285047, u285048, u285067, u285158, u285178, u285194, u285211, u285239, u285278, u285325, u285486, u285525, u285531, u285533, u285552, u285555, u285561, u285578, u285583, u285674, u285680, u285739, u285786, u285787, u285866, u285867, u285884, u285885, u285912, u285927, u286029, u286053, u286062, u286110, u286111, u286289, u286306, u286320, u286322, u286323, u286332, u286337, u286461, u286476, u286572, u286625, u286648, u286676, u286677, u286759, u286817, u286832, u286863, u286884, u286909, u286913, u286917, u286921, u286923, u286957, u287016, u287053, u287059, u287091, u287110, u287148, u287170, u287183, u287189, u287196, u287202, u287210, u287214, u287215, u287218, u287226, u287244, u287262, u287269, u287319, u287358, u287398, u287461, u287475, u287480, u287509, u287708, u287730, u287784, u287833, u287843, u287849, u287856, u287867, u287888, u287902, u287917, u288064, u288112, u288184, u288335, u288337, u288339, u288356, u288359, u288409, u288417, u288481, u288539, u288577, u288632, u288642, u288647, u288671, u288691, u288884, u289002, u289044, u289050, u289053, u289057, u289075, u289078, u289084, u289166, u28920, u289241, u289271, u289305, u289313, u289317, u289375, u289389, u289396, u289427, u289449, u289545, u289622, u289643, u289689, u289700, u289716, u289782, u289816, u289869, u289882, u290038, u290087, u290111, u290167, u290191, u290219, u290245, u290261, u290293, u290300, u290335, u290361, u290463, u290499, u290523, u290559, u290618, u290664, u290672, u290687, u290891, u290934, u290972, u290977, u290983, u291071, u291104, u291185, u291193, u291219, u291223, u291224, u291308, u291323, u291335, u291363, u291402, u291406, u291422, u291534, u291586, u291598, u291659, u291674, u291682, u291735, u291777, u291810, u291837, u291839, u291885, u291887, u291989, u29199, u292001, u292073, u292074, u292075, u292083, u292153, u292156, u292178, u292181, u292188, u292190, u29222, u29224, u292251, u292304, u292320, u292331, u292382, u292422, u292430, u292431, u292473, u292476, u292573, u292600, u292622, u292642, u292660, u292699, u292778, u292839, u292848, u292853, u292873, u292887, u292893, u292911, u292924, u292979, u293001, u293052, u293065, u293152, u293277, u293295, u293310, u29341, u293458, u293462, u293463, u293469, u293479, u293508, u293521, u293540, u293641, u293694, u293711, u293773, u293849, u293890, u293926, u293940, u293956, u293973, u294016, u294127, u294143, u294161, u294177, u294215, u294365, u294391, u294431, u294434, u294436, u294446, u29445, u29448, u294495, u294508, u294544, u294545, u294546, u294599, u294601, u294635, u294642, u294670, u294689, u294875, u294926, u294934, u294965, u295039, u295048, u295057, u295074, u295084, u295088, u295095, u295151, u295223, u29523, u295239, u295299, u295319, u295378, u295380, u295382, u295395, u295415, u29545, u295453, u295480, u295550, u295635, u295673, u295697, u295738, u295752, u295782, u295784, u295799, u295806, u295837, u29584, u295856, u295879, u295880, u295893, u295897, u295948, u295981, u295989, u295995, u296016, u296072, u296117, u296170, u296241, u296258, u296289, u296296, u296300, u296343, u296385, u296614, u29662, u296622, u29923, u29967, u29992, u29995, u30048, u30228, u30246, u30369, u30451, u30579, u30603, u30633, u30660, u30694, u31060, u31570, u31644, u32112, u32146, u32277, u32650, u32815, u32827, u32863, u32894, u32912, u32990, u32998, u33010, u33040, u33054, u33061, u33657, u33674, u33692, u33694, u33712, u33741, u33754, u34116, u34125, u34223, u34243, u34260, u34509, u34524, u34568, u34572, u34776, u34802, u34817, u34855, u34886, u34981, u35511, u35515, u35591, u35622, u35633, u35665, u36474, u36707, u37219, u37267, u37310, u37376, u37466, u37546, u37763, u37809, u37923, u37924, u37971, u37982, u38009, u38028, u38224, u38231, u38479, u38609, u38610, u38639, u38851, u38949, u38968, u39027, u39031, u39032, u39033, u39035, u39090, u39115, u39132, u39215, u39358, u39410, u39537, u39553, u39587, u39597, u39722, u39734, u40093, u40224, u41481, u41486, u41555, u41616, u41627, u41719, u41728, u41743, u41954, u41966, u42093, u42530, u42784, u42785, u42789, u42818, u42868, u42942, u42982, u43063, u43080, u43081, u43505, u43507, u43629, u43848, u44003, u44056, u44065, u44112, u44955, u45159, u45282, u45348, u45416, u45421, u45439, u45459, u45685, u45768, u45802, u45806, u45857, u46191, u46293, u46374, u46383, u46423, u46572, u46575, u46599, u46674, u46687, u46784, u46799, u46849, u46887, u46909, u46953, u46955, u46980, u47101, u47148, u47154, u47230, u47293, u47319, u47368, u47388, u47554, u47640, u47822, u47843, u47941, u47986, u48039, u48241, u48266, u48287, u48309, u48320, u48496, u48510, u48517, u48634, u48647, u48692, u48695, u48696, u48704, u48732, u48734, u48756, u48757, u48773, u48830, u48937, u49117, u49153, u49193, u49222, u49224, u49504, u49536, u49566, u49575, u49751, u49821, u49894, u49972, u50028, u50052, u50066, u50115, u50140, u50169, u50220, u50294, u50406, u50586, u50712, u50774, u50799, u50806, u50960, u50966, u51002, u51047, u51073, u51084, u51105, u51110, u51119, u51137, u51156, u51161, u51179, u51225, u51227, u51308, u51339, u51374, u51378, u51431, u51444, u51483, u51540, u51544, u51681, u51690, u51767, u51806, u51826, u51829, u51862, u51933, u51936, u51941, u51983, u52010, u52018, u52029, u52150, u52156, u52197, u52330, u52356, u52394, u52424, u52440, u52503, u52541, u52544, u52549, u52608, u52621, u52623, u52683, u52707, u52713, u52761, u52892, u52903, u52943, u52968, u53076, u53093, u53118, u53131, u53142, u53165, u53225, u53234, u53298, u53304, u53347, u53348, u53403, u53413, u53440, u53466, u53532, u53624, u53721, u53843, u53844, u53892, u53993, u54000, u54018, u54037, u54102, u54115, u54157, u54203, u54220, u54250, u54293, u54303, u54327, u54391, u54425, u54524, u54566, u54577, u54602, u54670, u54746, u54794, u54840, u54865, u54920, u54985, u54999, u55170, u55206, u55247, u55360, u55385, u55412, u55479, u55495, u55536, u55600, u55642, u55733, u55807, u55820, u55840, u55903, u56007, u56023, u56051, u56065, u56098, u56180, u56229, u56263, u56278, u56317, u56332, u56336, u56446, u56452, u56475, u56545, u56554, u56560, u56607, u56681, u56722, u56724, u56747, u56748, u56766, u56775, u56784, u56793, u56809, u56868, u56900, u57012, u57035, u57057, u57062, u57098, u57135, u57156, u57184, u57196, u57258, u57268, u57286, u57330, u57420, u57494, u57595, u57609, u57645, u57677, u57685, u57764, u57828, u57837, u57917, u57979, u58024, u58048, u58069, u58081, u58124, u58151, u58182, u58198, u58209, u58224, u58278, u58294, u58327, u58341, u58387, u58418, u58438, u58450, u58469, u58485, u58515, u58553, u58566, u58679, u58691, u58709, u58736, u58753, u58797, u58915, u58990, u59029, u59030, u59059, u59111, u59178, u59182, u59213, u59236, u59307, u59321, u59380, u59389, u59392, u59576, u59598, u59675, u59679, u59706, u59736, u59858, u59916, u59960, u59974, u59980, u60034, u60087, u60112, u60287, u60401, u60429, u60505, u60559, u60583, u60638, u60639, u60654, u60665, u60689, u60717, u60816, u60994, u61021, u61085, u61115, u61171, u61185, u61186, u61204, u61337, u61388, u61435, u61488, u61491, u61542, u61560, u61572, u61647, u61758, u61762, u61803, u61820, u61855, u61864, u61876, u61898, u61923, u62016, u62033, u62052, u62057, u62092, u62103, u62114, u62150, u62191, u62296, u62485, u62490, u62568, u62636, u62648, u62683, u62735, u62782, u62819, u62823, u62958, u63103, u63233, u63311, u63312, u63394, u63411, u63506, u63598, u63601, u63636, u63640, u63676, u63679, u63708, u63734, u63737, u63959, u63987, u64053, u64083, u64095, u64097, u64255, u64272, u64408, u64446, u64502, u64549, u64550, u64632, u64652, u64753, u64754, u64770, u64828, u64859, u64892, u64901, u64906, u64937, u64954, u64977, u65106, u65124, u65237, u65258, u65278, u65412, u65598, u65703, u65722, u65727, u65739, u65756, u65875, u65925, u65955, u66091, u66121, u66183, u66192, u66303, u66327, u66359, u66437, u66503, u66522, u66617, u66620, u66626, u66637, u66661, u66697, u66730, u66752, u66754, u66861, u66867, u66901, u66931, u66965, u67028, u67039, u67098, u67171, u67186, u67337, u67374, u67423, u67464, u67505, u67523, u67535, u67543, u67572, u67587, u67661, u67679, u67720, u67721, u67777, u67786, u67803, u67855, u67857, u67860, u67913, u67922, u67962, u68071, u68145, u68160, u68181, u68214, u68388, u68402, u68411, u68419, u68433, u68485, u68501, u68505, u68523, u68556, u68564, u68632, u68706, u68728, u68771, u68862, u68910, u68984, u69016, u69065, u69133, u69147, u69185, u69252, u69261, u69276, u69279, u69338, u69361, u69388, u69465, u69497, u69593, u69627, u69671, u69680, u69783, u69830, u69842, u69889, u69939, u70021, u70047, u70050, u70101, u70162, u70191, u70244, u70272, u70348, u70525, u70542, u70563, u70568, u70585, u70603, u70606, u70629, u70654, u70753, u70787, u70809, u70826, u70891, u70936, u70967, u70973, u71033, u71104, u71141, u71153, u71235, u71286, u71318, u71501, u71635, u71682, u71699, u71771, u71790, u71815, u71893, u71958, u71964, u71968, u71997, u72200, u72222, u72257, u72289, u72331, u72421, u72427, u72458, u72507, u72547, u72564, u72592, u72636, u72665, u72705, u72761, u72804, u72866, u73024, u73035, u73044, u7310, u73207, u73215, u73280, u73413, u73481, u73518, u73533, u73554, u73604, u73708, u73730, u73752, u73771, u73806, u73854, u73899, u74003, u74069, u74080, u74274, u74302, u74332, u74356, u74423, u74436, u74488, u74524, u74593, u74678, u74693, u74696, u74710, u74719, u74764, u74939, u74942, u75019, u75026, u75130, u75213, u75289, u75326, u75350, u75372, u75403, u75420, u75473, u75486, u75494, u75496, u75579, u75596, u75597, u75602, u75635, u75683, u75695, u75728, u75729, u75781, u75818, u75864, u75876, u75904, u75929, u75998, u76072, u76118, u76141, u76291, u76413, u76493, u76519, u76544, u76554, u76590, u76699, u76750, u76812, u76884, u76895, u76926, u76930, u76955, u76972, u77167, u77170, u77172, u77187, u77198, u77204, u77284, u77300, u77315, u77316, u77317, u77340, u77386, u77395, u77474, u77485, u77690, u77797, u77807, u77832, u78005, u78006, u78012, u78057, u78070, u78175, u78192, u78246, u78268, u78355, u78373, u78380, u78427, u78496, u78500, u78544, u78576, u78579, u78698, u78730, u78733, u78749, u78862, u78901, u78924, u79015, u79041, u79074, u79157, u79255, u79257, u79270, u79292, u79314, u79342, u79384, u79389, u79396, u79415, u79473, u79504, u79529, u79614, u79619, u79620, u79644, u79697, u79785, u79791, u79851, u79928, u80031, u80257, u80328, u80339, u80375, u80392, u80396, u80454, u80511, u80667, u80673, u80680, u80686, u80775, u80817, u80838, u80880, u80889, u80895, u80907, u80920, u80926, u81010, u81029, u81120, u81124, u81132, u81169, u81203, u81214, u81229, u81256, u81267, u81276, u81284, u81317, u81405, u81475, u81498, u81549, u81584, u81613, u81678, u81768, u81791, u81862, u81866, u81874, u81937, u82164, u82189, u82228, u82301, u82316, u82338, u82362, u82516, u82528, u82546, u82584, u82594, u82613, u82614, u82645, u82647, u82843, u82859, u82927, u82961, u83003, u83018, u83025, u83032, u83120, u83145, u83148, u83183, u83205, u83288, u83312, u83331, u83339, u83444, u83463, u83491, u83492, u83573, u83610, u83635, u83661, u83754, u83816, u83863, u83878, u83991, u84029, u84071, u84091, u84186, u84201, u84363, u84364, u84387, u84461, u84485, u84671, u84697, u84701, u84710, u84753, u84789, u84813, u84857, u84898, u84952, u84963, u84977, u84979, u84986, u84990, u85000, u85007, u85014, u85037, u85065, u85125, u85149, u85301, u85340, u85373, u85408, u85576, u85608, u85629, u85762, u85767, u85796, u85889, u85909, u85923, u85937, u85940, u85963, u85969, u85986, u85987, u86079, u86099, u86113, u86129, u86215, u86429, u86489, u86552, u86631, u86639, u86681, u86688, u86740, u86856, u86861, u86935, u86955, u86957, u86959, u86966, u86980, u86992, u86993, u87037, u87139, u87205, u87234, u87243, u87269, u87286, u87321, u87362, u87399, u87410, u87431, u87462, u87511, u87544, u87573, u87577, u87668, u87676, u87678, u87690, u87754, u87795, u87799, u87800, u87810, u87841, u87843, u87868, u87905, u87928, u87935, u88012, u88055, u88073, u88175, u88210, u88334, u88362, u88382, u88392, u88414, u88423, u88451, u88475, u88495, u88507, u88509, u88513, u88532, u88558, u88592, u88595, u88600, u88662, u88738, u88742, u88790, u88829, u88848, u88875, u88941, u89047, u89089, u89094, u8911, u89111, u89118, u89166, u89278, u89324, u89390, u89414, u89429, u89459, u89463, u89469, u89530, u89568, u89586, u89601, u89641, u89755, u89782, u89813, u89842, u89869, u89968, u90043, u90066, u90072, u90091, u90133, u90177, u90189, u90289, u90310, u90363, u90364, u90399, u90433, u90474, u90476, u90521, u90526, u90613, u90618, u90651, u90670, u90698, u90750, u90769, u90788, u90800, u90801, u90809, u90818, u90837, u90850, u90863, u90929, u90991, u91024, u91031, u91032, u91048, u91083, u91085, u91132, u91161, u91210, u91218, u91234, u91245, u91254, u91287, u91292, u91320, u91360, u91371, u91414, u91458, u91574, u91584, u91592, u91604, u91642, u91644, u91658, u91675, u91744, u91862, u91906, u91950, u92053, u92105, u92128, u92136, u92153, u92155, u92163, u92185, u92194, u92211, u92255, u92279, u92284, u92301, u92305, u92330, u92340, u92351, u92396, u92416, u92568, u92709, u92763, u92765, u92775, u92786, u92801, u92827, u92829, u92865, u92915, u92989, u92998, u93006, u93035, u93134, u93170, u93208, u93277, u93280, u93291, u93313, u93358, u93395, u93456, u93483, u93519, u93572, u93598, u93711, u93783, u93821, u93841, u93863, u93871, u93875, u93877, u93889, u93912, u93962, u93999, u94049, u94095, u94121, u94134, u94163, u94176, u94207, u94221, u94322, u94359, u94362, u94474, u94524, u94532, u94552, u94714, u94760, u94815, u94816, u95011, u95041, u95084, u95119, u95127, u95151, u95375, u95384, u95406, u95432, u95498, u95506, u95522, u95523, u95535, u95536, u95558, u95596, u95634, u95661, u95692, u95749, u95754, u95763, u95794, u95844, u95913, u95926, u95942, u95976, u95987, u95992, u95993, u96001, u96009, u96118, u96149, u96171, u96184, u96239, u96259, u96330, u96341, u96420, u96424, u96443, u96457, u96465, u96511, u96613, u96665, u96704, u96705, u96794, u96878, u96881, u96891, u96900, u96956, u97055, u97062, u97088, u97329, u97396, u97421, u97425, u97532, u97572, u97641, u97687, u97693, u97775, u97782, u97810, u97840, u97844, u97894, u98002, u98013, u98140, u98182, u98192, u98338, u98365, u98479, u98560, u98611, u98632, u98651, u98715, u98744, u98754, u98762, u98794, u98851, u98856, u98885, u98915, u99026, u99224, u99228, u99256, u99261, u99277, u99278, u99283, u99323, u99364, u99371, u99390, u99442, u99468, u99501, u99520, u99564, u99619, u99622, u99654, u99665, u99700, u99703, u99727, u99812, u99830, u99835, u99838, u99840, u99845, u99846, u99876, u99898, u99920, u99936, u99941, u99949
[2025-01-11 17:12:15,334]    INFO [extract] Capturing records from BUS file extracted_seqs/tmp/output.bus to extracted_seqs/tmp/output_extracted_u100074.bus with capture list extracted_seqs/tmp/pull_out_reads_transcript_ids_temp.txt
[2025-01-11 17:12:15,334]   DEBUG [extract] bustools capture -o extracted_seqs/tmp/output_extracted_u100074.bus -c extracted_seqs/tmp/pull_out_reads_transcript_ids_temp.txt -e extracted_seqs/tmp/matrix_no_mm.ec -t extracted_seqs/tmp/transcripts.txt --transcripts extracted_seqs/tmp/output.bus
[2025-01-11 17:12:15,739]   DEBUG [extract] Parsing transcripts .. done
[2025-01-11 17:12:15,739]   DEBUG [extract] Parsing ECs .. done
[2025-01-11 17:12:15,739]   DEBUG [extract] Parsing capture list .. done
[2025-01-11 17:12:16,840]   DEBUG [extract] Read in 3086 BUS records, wrote 323 BUS records
[2025-01-11 17:12:17,943]   DEBUG [extract] extracted_seqs/tmp/output_extracted_u100074.bus passed validation
[2025-01-11 17:12:17,943]    INFO [extract] Sorting BUS file extracted_seqs/tmp/output_extracted_u100074.bus to extracted_seqs/tmp/output_extracted_u100074_sorted.bus
[2025-01-11 17:12:17,943]   DEBUG [extract] bustools sort -o extracted_seqs/tmp/output_extracted_u100074_sorted.bus -T tmp -t 8 -m 2G --flags extracted_seqs/tmp/output_extracted_u100074.bus
[2025-01-11 17:12:18,044]   DEBUG [extract] Warning: Number of threads cannot be greater than or equal to 2. Setting number of threads to 2
[2025-01-11 17:12:20,252]   DEBUG [extract] all fits in buffer
[2025-01-11 17:12:21,553]   DEBUG [extract] Read in 323 BUS records
[2025-01-11 17:12:21,554]   DEBUG [extract] reading time 1.1e-05s
[2025-01-11 17:12:21,554]   DEBUG [extract] sorting time 2.1e-05s
[2025-01-11 17:12:21,554]   DEBUG [extract] writing time 0.000122s
[2025-01-11 17:12:21,554]    INFO [extract] Extracting BUS file extracted_seqs/tmp/output_extracted_u100074_sorted.bus to extracted_seqs/u100074
[2025-01-11 17:12:21,554]   DEBUG [extract] bustools extract -o extracted_seqs/u100074 -f /content/SRR11777734_short.fastq -N 1 extracted_seqs/tmp/output_extracted_u100074_sorted.bus
[2025-01-11 17:12:34,783]   DEBUG [extract] Read in 323 BUS records
[2025-01-11 17:12:34,807]   DEBUG [main] Removing `extracted_seqs/tmp` directory
CPU times: user 10.1 s, sys: 1.37 s, total: 11.5 s
Wall time: 21min 17s

# Grab the reads extracted by kb extract for each virus ID
def get_all_sequences(fastq_file):
    """Extract all sequences from a gzipped FASTQ file."""
    sequences = []
    with gzip.open(fastq_file, "rt") as f:
        for record in SeqIO.parse(f, "fastq"):
            sequences.append(str(record.seq))
    return sequences

extracted_seqs_dict = {}
for vir_id in vir_ids_to_extract:
    # Grab the reads from the fastq file returned by kb extract
    seqs = get_all_sequences(f"{outfolder_extract}/{vir_id}/1.fastq.gz")
    # Save the extracted reads for each virus ID in a dictionary
    extracted_seqs_dict[vir_id] = seqs

BLAST extracted reads

BLAST the first read extracted for the first virus ID (here, u100074) using `gget blast <https://pachterlab.github.io/gget/en/blast.html>`__ to see which known reference sequences it aligns to:

seq_idx = 0
raw_read = extracted_seqs_dict[vir_ids_to_extract[0]][seq_idx]

blast_df = gget.blast(raw_read)
blast_df

INFO:gget.utils:Sequence recognized as nucleotide sequence.
INFO:gget.utils:BLAST will use program 'blastn' with database 'nt'.
INFO:gget.utils:BLAST initiated with search ID S59M77GX016. Estimated time to completion: 30 seconds.
INFO:gget.utils:Retrieving results...
/usr/local/lib/python3.10/dist-packages/gget/gget_blast.py:327: FutureWarning: Passing literal html to 'read_html' is deprecated and will be removed in a future version. To read from a literal string, wrap it in a 'StringIO' object.
  results_df = pd.read_html(str(dsc_table))[0]

	Description	Scientific Name	Common Name	Taxid	Max Score	Total Score	Query Cover	E value	Per. Ident	Acc. Len	Accession
0	Severe acute respiratory syndrome coronavirus ...	Severe acute respiratory syndrome coronavirus 2	NaN	2697049	124	124	100%	9.000000e-25	96.05%	29890	OV107737.2
1	Severe acute respiratory syndrome coronavirus ...	Severe acute respiratory syndrome coronavirus 2	NaN	2697049	124	124	100%	9.000000e-25	96.05%	29655	OQ122586.1
2	Severe acute respiratory syndrome coronavirus ...	Severe acute respiratory syndrome coronavirus 2	NaN	2697049	124	124	100%	9.000000e-25	96.05%	29715	OQ579903.1
3	Severe acute respiratory syndrome coronavirus ...	Severe acute respiratory syndrome coronavirus 2	NaN	2697049	124	124	100%	9.000000e-25	96.05%	29718	OP415689.1
4	Severe acute respiratory syndrome coronavirus ...	Severe acute respiratory syndrome coronavirus 2	NaN	2697049	124	124	100%	9.000000e-25	96.05%	29841	OX424285.1
5	Severe acute respiratory syndrome coronavirus ...	Severe acute respiratory syndrome coronavirus 2	NaN	2697049	124	124	100%	9.000000e-25	96.05%	29850	OX858193.1
6	Severe acute respiratory syndrome coronavirus ...	Severe acute respiratory syndrome coronavirus 2	NaN	2697049	124	124	100%	9.000000e-25	96.05%	29890	OV024669.1
7	Severe acute respiratory syndrome coronavirus ...	Severe acute respiratory syndrome coronavirus 2	NaN	2697049	124	124	100%	9.000000e-25	96.05%	29747	OQ697877.1
8	Severe acute respiratory syndrome coronavirus ...	Severe acute respiratory syndrome coronavirus 2	NaN	2697049	124	124	100%	9.000000e-25	96.05%	29868	OZ080834.1
9	Severe acute respiratory syndrome coronavirus ...	Severe acute respiratory syndrome coronavirus 2	NaN	2697049	124	124	100%	9.000000e-25	96.05%	29817	OW805761.1
10	Severe acute respiratory syndrome coronavirus ...	Severe acute respiratory syndrome coronavirus 2	NaN	2697049	124	124	100%	9.000000e-25	96.05%	29747	OQ698290.1
11	Severe acute respiratory syndrome coronavirus ...	Severe acute respiratory syndrome coronavirus 2	NaN	2697049	124	124	100%	9.000000e-25	96.05%	29823	OK581539.1
12	Severe acute respiratory syndrome coronavirus ...	Severe acute respiratory syndrome coronavirus 2	NaN	2697049	124	124	100%	9.000000e-25	96.05%	29873	OY803461.1
13	Severe acute respiratory syndrome coronavirus ...	Severe acute respiratory syndrome coronavirus 2	NaN	2697049	124	124	100%	9.000000e-25	96.05%	29652	OR650313.1
14	Severe acute respiratory syndrome coronavirus ...	Severe acute respiratory syndrome coronavirus 2	NaN	2697049	124	124	100%	9.000000e-25	96.05%	29890	OY500377.1
15	Severe acute respiratory syndrome coronavirus ...	Severe acute respiratory syndrome coronavirus 2	NaN	2697049	124	124	100%	9.000000e-25	96.05%	29874	OR938031.1
16	Severe acute respiratory syndrome coronavirus ...	Severe acute respiratory syndrome coronavirus 2	NaN	2697049	124	124	100%	9.000000e-25	96.05%	29897	OX775923.1
17	Severe acute respiratory syndrome coronavirus ...	Severe acute respiratory syndrome coronavirus 2	NaN	2697049	124	124	100%	9.000000e-25	96.05%	29769	OL963202.1
18	Severe acute respiratory syndrome coronavirus ...	Severe acute respiratory syndrome coronavirus 2	NaN	2697049	124	124	100%	9.000000e-25	96.05%	29876	OX776611.1
19	Severe acute respiratory syndrome coronavirus ...	Severe acute respiratory syndrome coronavirus 2	NaN	2697049	124	124	100%	9.000000e-25	96.05%	29769	OL901140.1
20	Severe acute respiratory syndrome coronavirus ...	Severe acute respiratory syndrome coronavirus 2	NaN	2697049	124	124	100%	9.000000e-25	96.05%	29867	OZ059218.1
21	Severe acute respiratory syndrome coronavirus ...	Severe acute respiratory syndrome coronavirus 2	NaN	2697049	124	124	100%	9.000000e-25	96.05%	29715	OQ467913.1
22	Severe acute respiratory syndrome coronavirus ...	Severe acute respiratory syndrome coronavirus 2	NaN	2697049	124	124	100%	9.000000e-25	96.05%	29824	OM914127.1
23	Severe acute respiratory syndrome coronavirus ...	Severe acute respiratory syndrome coronavirus 2	NaN	2697049	124	124	100%	9.000000e-25	96.05%	29890	OU878122.1
24	Severe acute respiratory syndrome coronavirus ...	Severe acute respiratory syndrome coronavirus 2	NaN	2697049	124	124	100%	9.000000e-25	96.05%	29873	OW071185.1
25	Severe acute respiratory syndrome coronavirus ...	Severe acute respiratory syndrome coronavirus 2	NaN	2697049	124	124	100%	9.000000e-25	96.05%	29778	OY667194.1
26	Severe acute respiratory syndrome coronavirus ...	Severe acute respiratory syndrome coronavirus 2	NaN	2697049	124	124	100%	9.000000e-25	96.05%	29885	OX860078.1
27	Severe acute respiratory syndrome coronavirus ...	Severe acute respiratory syndrome coronavirus 2	NaN	2697049	124	124	100%	9.000000e-25	96.05%	29897	OY296656.1
28	Severe acute respiratory syndrome coronavirus ...	Severe acute respiratory syndrome coronavirus 2	NaN	2697049	124	124	100%	9.000000e-25	96.05%	29890	OV027327.1
29	Severe acute respiratory syndrome coronavirus ...	Severe acute respiratory syndrome coronavirus 2	NaN	2697049	124	124	100%	9.000000e-25	96.05%	29890	OV089785.1
30	Severe acute respiratory syndrome coronavirus ...	Severe acute respiratory syndrome coronavirus 2	NaN	2697049	124	124	100%	9.000000e-25	96.05%	29890	OV269224.1
31	Severe acute respiratory syndrome coronavirus ...	Severe acute respiratory syndrome coronavirus 2	NaN	2697049	124	124	100%	9.000000e-25	96.05%	29890	OV319883.1
32	Severe acute respiratory syndrome coronavirus ...	Severe acute respiratory syndrome coronavirus 2	NaN	2697049	124	124	100%	9.000000e-25	96.05%	29867	OY802686.1
33	Severe acute respiratory syndrome coronavirus ...	Severe acute respiratory syndrome coronavirus 2	NaN	2697049	124	124	100%	9.000000e-25	96.05%	29844	OX431773.1
34	Severe acute respiratory syndrome coronavirus ...	Severe acute respiratory syndrome coronavirus 2	NaN	2697049	124	124	100%	9.000000e-25	96.05%	29652	OQ681834.1
35	Severe acute respiratory syndrome coronavirus ...	Severe acute respiratory syndrome coronavirus 2	NaN	2697049	124	124	100%	9.000000e-25	96.05%	29715	OQ654309.1
36	Severe acute respiratory syndrome coronavirus ...	Severe acute respiratory syndrome coronavirus 2	NaN	2697049	124	124	100%	9.000000e-25	96.05%	29812	OV225377.1
37	Severe acute respiratory syndrome coronavirus ...	Severe acute respiratory syndrome coronavirus 2	NaN	2697049	124	124	100%	9.000000e-25	96.05%	29720	OR603704.1
38	Severe acute respiratory syndrome coronavirus ...	Severe acute respiratory syndrome coronavirus 2	NaN	2697049	124	124	100%	9.000000e-25	96.05%	29721	OQ727759.1
39	Severe acute respiratory syndrome coronavirus ...	Severe acute respiratory syndrome coronavirus 2	NaN	2697049	124	124	100%	9.000000e-25	96.05%	29844	OX677911.1
40	Severe acute respiratory syndrome coronavirus ...	Severe acute respiratory syndrome coronavirus 2	NaN	2697049	124	124	100%	9.000000e-25	96.05%	29859	OK383762.1
41	Severe acute respiratory syndrome coronavirus ...	Severe acute respiratory syndrome coronavirus 2	NaN	2697049	124	124	100%	9.000000e-25	96.05%	29867	OY797923.1
42	Severe acute respiratory syndrome coronavirus ...	Severe acute respiratory syndrome coronavirus 2	NaN	2697049	124	124	100%	9.000000e-25	96.05%	29867	OY800051.1
43	Severe acute respiratory syndrome coronavirus ...	Severe acute respiratory syndrome coronavirus 2	NaN	2697049	124	124	100%	9.000000e-25	96.05%	29787	OL884912.1
44	Severe acute respiratory syndrome coronavirus ...	Severe acute respiratory syndrome coronavirus 2	NaN	2697049	124	124	100%	9.000000e-25	96.05%	29823	OV239519.1
45	Severe acute respiratory syndrome coronavirus ...	Severe acute respiratory syndrome coronavirus 2	NaN	2697049	124	124	100%	9.000000e-25	96.05%	29694	ON588468.1
46	Severe acute respiratory syndrome coronavirus ...	Severe acute respiratory syndrome coronavirus 2	NaN	2697049	124	124	100%	9.000000e-25	96.05%	29890	OV182074.1
47	Severe acute respiratory syndrome coronavirus ...	Severe acute respiratory syndrome coronavirus 2	NaN	2697049	121	121	100%	1.000000e-23	94.74%	29902	OY941257.1
48	Severe acute respiratory syndrome coronavirus ...	Severe acute respiratory syndrome coronavirus 2	NaN	2697049	121	121	100%	1.000000e-23	94.74%	29171	OP402614.1
49	Severe acute respiratory syndrome coronavirus ...	Severe acute respiratory syndrome coronavirus 2	NaN	2697049	121	121	100%	1.000000e-23	94.74%	29338	OL489374.1

It looks like reads mapped to virus ID u100074 might have also belonged to SARS-CoV2 (these reads mapped to u100074 instead of u1—for which the predicted taxonomy is SARS-CoV), potentially suggesting the presence of multiple SARS-CoV2 variants or viral evolution in the cardiomyocytes.

For inspiration regarding further analysis, each figure from our manuscript can be reproduced using these notebooks.