# Converting BUS format into sparse matrix

Lambda Moses

#### 2022-04-26

# Contents

## 0.1 Introduction

The Barcode, UMI, Set (BUS) format is a new way to represent pseudoalignments of reads from RNA-seq. Files of this format can be efficiently generated by the command line tool `kallisto bus`. With `kallisto bus` and this package, we go from the fastq files to the sparse matrix used for downstream analysis such as with Seurat within half an hour, while Cell Ranger would take hours.

In this vignette, we convert an 10x 1:1 mouse and human cell mixture dataset from the BUS format to a sparse matrix. To see how the BUS format can be generated from fastq file, as well as more in depth vignettes, see the website of this package.

Note that this vignette is deprecated and is kept for historical reasons as it was implemented when `kallisto | bustools` was experimental. The functionality of `make_sparse_matrix` has been implemented more efficiently in the command line tool `bustools`. Please use the updated version of `bustools` and if you wish, the wrapper `kb` instead.

## 0.2 Download the dataset

    # The dataset package
    library(TENxBUSData)
    library(BUSpaRse)
    library(Matrix)
    library(zeallot)
    library(ggplot2)

    TENxBUSData(".", dataset = "hgmm100")
    #> snapshotDate(): 2022-04-19
    #> see ?TENxBUSData and browseVignettes('TENxBUSData') for documentation
    #> loading from cache
    #> The downloaded files are in F:\biocbuild\bbs-3.15-bioc\tmpdir\RtmpeESxMO\Rbuild1f60400f61f5\BUSpaRse\vignettes/out_hgmm100
    #> [1] "F:\\biocbuild\\bbs-3.15-bioc\\tmpdir\\RtmpeESxMO\\Rbuild1f60400f61f5\\BUSpaRse\\vignettes/out_hgmm100"

## 0.3 Convert to sparse matrix

First, we map transcripts, as in the `kallisto` index, to the corresponding genes.

    tr2g <- transcript2gene(species = c("Homo sapiens", "Mus musculus"), type = "vertebrate",
                            kallisto_out_path = "./out_hgmm100", ensembl_version = 99,
                            write_tr2g = FALSE)
    #> Querying biomart for transcript and gene IDs of Homo sapiens
    #> Querying biomart for transcript and gene IDs of Mus musculus

    head(tr2g)
    #>          transcript              gene gene_name chromosome_name
    #> 1 ENST00000434970.2 ENSG00000237235.2     TRDD2              14
    #> 2 ENST00000415118.1 ENSG00000223997.1     TRDD1              14
    #> 3 ENST00000448914.1 ENSG00000228985.1     TRDD3              14
    #> 4 ENST00000632684.1 ENSG00000282431.1     TRBD1               7
    #> 5 ENST00000390583.1 ENSG00000211923.1  IGHD3-10              14
    #> 6 ENST00000431440.2 ENSG00000232543.2  IGHD4-11              14

Here we make both the gene count matrix and the TCC matrix.

    c(gene_count, tcc) %<-% make_sparse_matrix("./out_hgmm100/output.sorted.txt",
                                   tr2g = tr2g, est_ncells = 1e5,
                                   est_ngenes = nrow(tr2g))
    #> Reading matrix.ec
    #> Processing ECs
    #> Matching genes to ECs
    #> Reading data
    #> Constructing gene count matrix
    #> Constructing TCC matrix

## 0.4 Remove empty droplets

Here we have a sparse matrix with genes in rows and cells in columns.

    dim(gene_count)
    #> [1]  26937 151449

This dataset should only have about 100 cells, but here we get over 100,000. In fact, most of the barcodes correspond to empty droplets; they can be removed by filtering out barcodes with too few UMI.

    tot_counts <- Matrix::colSums(gene_count)
    summary(tot_counts)
    #>     Min.  1st Qu.   Median     Mean  3rd Qu.     Max.
    #>     1.00     1.00     2.00    24.77     5.00 64041.00

    df1 <- get_knee_df(gene_count)
    infl1 <- get_inflection(df1)

    knee_plot(df1, infl1)

    gene_count <- gene_count[, tot_counts > infl1]
    dim(gene_count)
    #> [1] 26937   122

Then this sparse matrix can be used in Seurat for downstream analysis.

Likewise, we can remove empty droplets from the TCC matrix.

    dim(tcc)
    #> [1] 129371 157659

This dataset should only have about 100 cells, but here we get over 100,000. In fact, most of the barcodes correspond to empty droplets; they can be removed by filtering out barcodes with too few UMI.

    tot_counts <- Matrix::colSums(tcc)
    summary(tot_counts)
    #>     Min.  1st Qu.   Median     Mean  3rd Qu.     Max.
    #>     1.00     1.00     2.00    25.84     5.00 69235.00

    df2 <- get_knee_df(tcc)
    infl2 <- get_inflection(df2)

    knee_plot(df2, infl2)

    tcc <- tcc[, tot_counts > infl2]
    dim(tcc)
    #> [1] 129371    121

    sessionInfo()
    #> R version 4.2.0 RC (2022-04-19 r82224 ucrt)
    #> Platform: x86_64-w64-mingw32/x64 (64-bit)
    #> Running under: Windows Server x64 (build 20348)
    #>
    #> Matrix products: default
    #>
    #> locale:
    #> [1] LC_COLLATE=C
    #> [2] LC_CTYPE=English_United States.utf8
    #> [3] LC_MONETARY=English_United States.utf8
    #> [4] LC_NUMERIC=C
    #> [5] LC_TIME=English_United States.utf8
    #>
    #> attached base packages:
    #> [1] stats     graphics  grDevices utils     datasets  methods   base
    #>
    #> other attached packages:
    #> [1] ggplot2_3.3.5     zeallot_0.1.0     Matrix_1.4-1      BUSpaRse_1.10.0
    #> [5] TENxBUSData_1.9.0 BiocStyle_2.24.0
    #>
    #> loaded via a namespace (and not attached):
    #>   [1] colorspace_2.0-3              rjson_0.2.21
    #>   [3] ellipsis_0.3.2                XVector_0.36.0
    #>   [5] GenomicRanges_1.48.0          farver_2.1.0
    #>   [7] bit64_4.0.5                   interactiveDisplayBase_1.34.0
    #>   [9] AnnotationDbi_1.58.0          fansi_1.0.3
    #>  [11] xml2_1.3.3                    cachem_1.0.6
    #>  [13] knitr_1.38                    jsonlite_1.8.0
    #>  [15] Rsamtools_2.12.0              dbplyr_2.1.1
    #>  [17] png_0.1-7                     shiny_1.7.1
    #>  [19] BiocManager_1.30.17           compiler_4.2.0
    #>  [21] httr_1.4.2                    assertthat_0.2.1
    #>  [23] fastmap_1.1.0                 lazyeval_0.2.2
    #>  [25] cli_3.3.0                     later_1.3.0
    #>  [27] htmltools_0.5.2               prettyunits_1.1.1
    #>  [29] tools_4.2.0                   gtable_0.3.0
    #>  [31] glue_1.6.2                    GenomeInfoDbData_1.2.8
    #>  [33] dplyr_1.0.8                   rappdirs_0.3.3
    #>  [35] Rcpp_1.0.8.3                  Biobase_2.56.0
    #>  [37] jquerylib_0.1.4               vctrs_0.4.1
    #>  [39] Biostrings_2.64.0             ExperimentHub_2.4.0
    #>  [41] rtracklayer_1.56.0            xfun_0.30
    #>  [43] stringr_1.4.0                 plyranges_1.16.0
    #>  [45] mime_0.12                     lifecycle_1.0.1
    #>  [47] restfulr_0.0.13               ensembldb_2.20.0
    #>  [49] XML_3.99-0.9                  AnnotationHub_3.4.0
    #>  [51] zlibbioc_1.42.0               scales_1.2.0
    #>  [53] BSgenome_1.64.0               hms_1.1.1
    #>  [55] promises_1.2.0.1              MatrixGenerics_1.8.0
    #>  [57] ProtGenerics_1.28.0           parallel_4.2.0
    #>  [59] SummarizedExperiment_1.26.0   AnnotationFilter_1.20.0
    #>  [61] yaml_2.3.5                    curl_4.3.2
    #>  [63] memoise_2.0.1                 sass_0.4.1
    #>  [65] biomaRt_2.52.0                stringi_1.7.6
    #>  [67] RSQLite_2.2.12                BiocVersion_3.15.2
    #>  [69] highr_0.9                     S4Vectors_0.34.0
    #>  [71] BiocIO_1.6.0                  GenomicFeatures_1.48.0
    #>  [73] BiocGenerics_0.42.0           filelock_1.0.2
    #>  [75] BiocParallel_1.30.0           GenomeInfoDb_1.32.0
    #>  [77] rlang_1.0.2                   pkgconfig_2.0.3
    #>  [79] matrixStats_0.62.0            bitops_1.0-7
    #>  [81] evaluate_0.15                 lattice_0.20-45
    #>  [83] purrr_0.3.4                   GenomicAlignments_1.32.0
    #>  [85] bit_4.0.4                     tidyselect_1.1.2
    #>  [87] magrittr_2.0.3                bookdown_0.26
    #>  [89] R6_2.5.1                      IRanges_2.30.0
    #>  [91] magick_2.7.3                  generics_0.1.2
    #>  [93] DelayedArray_0.22.0           DBI_1.1.2
    #>  [95] pillar_1.7.0                  withr_2.5.0
    #>  [97] KEGGREST_1.36.0               RCurl_1.98-1.6
    #>  [99] tibble_3.1.6                  crayon_1.5.1
    #> [101] utf8_1.2.2                    BiocFileCache_2.4.0
    #> [103] rmarkdown_2.14                progress_1.2.2
    #> [105] grid_4.2.0                    blob_1.2.3
    #> [107] digest_0.6.29                 xtable_1.8-4
    #> [109] tidyr_1.2.0                   httpuv_1.6.5
    #> [111] stats4_4.2.0                  munsell_0.5.0
    #> [113] bslib_0.3.1