alevin: Almost no overlapping cell barcodes between GEX and HTO libraries

[ivanek]

alevin (single-cell mode)

Describe the bug

After mapping sample libraries with alevin I got very low number (only 22) of cell barcodes present in both gene expression data (GEX, 7954 cells reported) and cell hashing libraries (HTO, 97305 cells reported). When mapping with cellranger v4.0.0 the overlap is reasonable (GEX 6844 cells reported, HTO 5990 cells reported, overlap of 5966 cell barcodes).

For GEX, all cells identified by cellranger were also reported by salmon alevin. The only difference is on HTO libraries. Given the fact, that cellranger is able to identify correctly cell barcodes in HTO libraries, there is most likely a bug in alevin.

To Reproduce

The GEX library was mapped with following command (the full path was replaced by relative path in this bug report):

salmon alevin -l ISR \
-1 ./fastq/BSSE_QGF_150046_HTJYKDRXX_1_24_h_timepoint_SI-TT-A1__S1_L001_R1_001.fastq.gz \
    ./fastq/BSSE_QGF_150046_HTJYKDRXX_2_24_h_timepoint_SI-TT-A1__S1_L002_R1_001.fastq.gz \
-2 ./fastq/BSSE_QGF_150046_HTJYKDRXX_1_24_h_timepoint_SI-TT-A1__S1_L001_R2_001.fastq.gz \
     ./fastq/BSSE_QGF_150046_HTJYKDRXX_2_24_h_timepoint_SI-TT-A1__S1_L002_R2_001.fastq.gz \
--chromiumV3 -i ./salmon_alevin_gencode_v35/gencode.v35.annotation.expanded.sidx \
-p 8 --dumpFeatures -o ./salmon_alevin_output/BSSE_QGF_150046_GEX \
--tgMap ./salmon_alevin_gencode_v35/gencode.v35.annotation.expanded.tx2gene.tsv

LOG: salmon_alevin_output/BSSE_QGF_150046_GEX/alevin/alevin.log:

[2020-10-14 22:55:44.218] [alevinLog] [info] Found 612271 transcripts(+194 decoys, +39 short and +0 duplicate names in the index)
[2020-10-14 22:55:45.097] [alevinLog] [info] Filled with 612310 txp to gene entries
[2020-10-14 22:55:45.231] [alevinLog] [info] Found all transcripts to gene mappings
[2020-10-14 22:55:45.298] [alevinLog] [info] Processing barcodes files (if Present)

[2020-10-14 22:59:24.696] [alevinLog] [info] Done barcode density calculation.
[2020-10-14 22:59:24.696] [alevinLog] [info] # Barcodes Used: 151479984 / 151479984.
[2020-10-14 22:59:30.432] [alevinLog] [info] Knee found left boundary at  6957
[2020-10-14 22:59:32.208] [alevinLog] [info] Gauss Corrected Boundary at  6957
[2020-10-14 22:59:32.208] [alevinLog] [info] Learned InvCov: 103.418 normfactor: 2347.53
[2020-10-14 22:59:32.208] [alevinLog] [info] Total 7958(has 1000 low confidence) barcodes
[2020-10-14 22:59:34.312] [alevinLog] [info] Done True Barcode Sampling
[2020-10-14 22:59:34.896] [alevinLog] [info] Total 7.07491% reads will be thrown away because of noisy Cellular barcodes.
[2020-10-14 22:59:35.370] [alevinLog] [info] Done populating Z matrix
[2020-10-14 22:59:35.415] [alevinLog] [info] Total 69392 CB got sequence corrected
[2020-10-14 22:59:35.430] [alevinLog] [info] Done indexing Barcodes
[2020-10-14 22:59:35.430] [alevinLog] [info] Total Unique barcodes found: 3616446
[2020-10-14 22:59:35.430] [alevinLog] [info] Used Barcodes except Whitelist: 68653
[2020-10-14 22:59:36.166] [alevinLog] [info] Done with Barcode Processing; Moving to Quantify

[2020-10-14 22:59:36.166] [alevinLog] [info] parsing read library format
[2020-10-14 23:34:55.585] [alevinLog] [info] Starting optimizer

[2020-10-14 23:34:55.742] [alevinLog] [warning] mrna file not provided; using is 1 less feature for whitelisting
[2020-10-14 23:34:55.742] [alevinLog] [warning] rrna file not provided; using is 1 less feature for whitelisting
[2020-10-14 23:42:45.591] [alevinLog] [info] Total 103442266.00 UMI after deduplicating.
[2020-10-14 23:42:45.591] [alevinLog] [info] Total 4315139 BiDirected Edges.
[2020-10-14 23:42:45.591] [alevinLog] [info] Total 33691 UniDirected Edges.
[2020-10-14 23:42:45.591] [alevinLog] [warning] Skipped 3 barcodes due to No mapped read
[2020-10-14 23:42:45.638] [alevinLog] [info] Clearing EqMap; Might take some time.
[2020-10-14 23:43:20.265] [alevinLog] [info] Starting white listing of 7954 cells
[2020-10-14 23:43:20.265] [alevinLog] [info] Starting to make feature Matrix
[2020-10-14 23:43:20.280] [alevinLog] [info] Done making feature Matrix
[2020-10-14 23:43:20.327] [alevinLog] [info] Finished white listing
[2020-10-14 23:43:20.375] [alevinLog] [info] Finished optimizer

The HTO library was mapped with following command (the full path was replaced by relative path in this bug report):

salmon alevin -l ISR \
-1 ./cellranger/fastq/lane1/BSSE_QGF_150050_S1_L001_R1_001.fastq.gz \
    ./cellranger/fastq/lane2/BSSE_QGF_150050_S1_L002_R1_001.fastq.gz \
-2 ./cellranger/fastq/lane1/BSSE_QGF_150050_S1_L001_R2_001.fastq.gz \
    ./cellranger/fastq/lane2/BSSE_QGF_150050_S1_L002_R2_001.fastq.gz \
-i ./salmon_alevin_gencode_v35/hto_index 
-o ./salmon_alevin_output/BSSE_QGF_150050_HTO \
-p 8 --citeseq \
--end 5 --umiLength 10 --barcodeLength 16 \
--featureStart 10 --featureLength 15 \
--naiveEqclass

LOG: salmon_alevin_output/BSSE_QGF_150050_HTO/alevin/alevin.log

[2020-10-20 10:05:32.442] [alevinLog] [info] set CITE-seq minScoreFraction parameter to : 0.65
[2020-10-20 10:05:32.447] [alevinLog] [info] Found 2 transcripts(+0 decoys, +0 short and +0 duplicate names in the index)
[2020-10-20 10:05:32.447] [alevinLog] [info] Filled with 2 txp to gene entries
[2020-10-20 10:05:32.447] [alevinLog] [info] Found all transcripts to gene mappings
[2020-10-20 10:05:32.453] [alevinLog] [info] Processing barcodes files (if Present)

[2020-10-20 10:07:04.690] [alevinLog] [info] Done barcode density calculation.
[2020-10-20 10:07:04.690] [alevinLog] [info] # Barcodes Used: 63789818 / 63789818.
[2020-10-20 10:07:04.910] [alevinLog] [info] Forcing to use 100000 cells
[2020-10-20 10:07:05.715] [alevinLog] [info] Throwing 2695 barcodes with < 10 reads
[2020-10-20 10:07:06.417] [alevinLog] [info] Total 97306(has 201 low confidence) barcodes
[2020-10-20 10:07:07.621] [alevinLog] [info] Done True Barcode Sampling
[2020-10-20 10:07:07.847] [alevinLog] [info] Total 2.16112% reads will be thrown away because of noisy Cellular barcodes.
[2020-10-20 10:07:12.867] [alevinLog] [info] Done populating Z matrix
[2020-10-20 10:07:12.867] [alevinLog] [info] Total 0 CB got sequence corrected
[2020-10-20 10:07:12.867] [alevinLog] [info] Done indexing Barcodes
[2020-10-20 10:07:12.867] [alevinLog] [info] Total Unique barcodes found: 893751
[2020-10-20 10:07:12.867] [alevinLog] [info] Used Barcodes except Whitelist: 0
[2020-10-20 10:07:13.571] [alevinLog] [info] Done with Barcode Processing; Moving to Quantify

[2020-10-20 10:07:13.571] [alevinLog] [info] parsing read library format
[2020-10-20 10:09:19.837] [alevinLog] [info] Starting optimizer

[2020-10-20 10:09:20.548] [alevinLog] [warning] mrna file not provided; using is 1 less feature for whitelisting
[2020-10-20 10:09:20.548] [alevinLog] [warning] rrna file not provided; using is 1 less feature for whitelisting
[2020-10-20 10:09:24.017] [alevinLog] [info] Total 43095134.00 UMI after deduplicating.
[2020-10-20 10:09:24.017] [alevinLog] [info] Total 0 BiDirected Edges.
[2020-10-20 10:09:24.017] [alevinLog] [info] Total 0 UniDirected Edges.
[2020-10-20 10:09:24.020] [alevinLog] [info] Finished optimizer

Specifically, please provide at least the following information:

• Which version of salmon was used?

$salmon alevin  --version
salmon 1.3.0

• How was salmon installed (compiled, downloaded executable, through bioconda)?

• with bioconda

  • Which reference (e.g. transcriptome) was used?

• human genome hg38 (based on Gencode v35) + viral genome

  • Which read files were used?

• see code above

  • Which which program options were used?

• see code above

Expected behavior
I would expect much higher overlap of detected cell barcodes between the GEX (gene expression) and HTO (hashtag) libraries.

Desktop (please complete the following information):

• OS: CentOS 7.8, Queing System: Slurm
• Version: Linux 3.10.0-1127.8.2.el7.x86_64 #1 SMP Tue May 12 16:57:42 UTC 2020 x86_64 x86_64 x86_64 GNU/Linux

Additional context

Dear salmon alevin developers,
I am happy to share part or even complete FASTQ files for one of the samples for debugging purpose.
Thanks for your help.
Best, Robert

[k3yavi]

Hi @ivanek ,

Thanks for reaching out. Can you check if you use the following mapping file for the RNA cb, does it improves the number of common barcodes between rna and hto ?

https://github.com/10XGenomics/cellranger/raw/master/lib/python/cellranger/barcodes/translation/3M-february-2018.txt.gz

[ivanek]

Unfortunately, only very little it goes to ~400 overlapping cell barcodes.

[k3yavi]

The logs look fine to me and nothing obvious stands out from the first look. Is it possible for you to share the fastq ? I'd be happy to take a look. Also we might need the barcode tags for the hto which was used in the experiment.

[ivanek]

Hi @k3yavi ,
I have noticed, when checking the molecule_info.h5 generated by cellranger that the barcodes provided by cellranger represent the most abundant cell barcodes (using the number of UMIs) but the cell barcodes provided by salmon alevin are actually scattered along the entire dataset (with among the less abundant).
Is there a way to get the same summary data (n_umis, n_reads) from salmon alevin? I tried to use the parameter --dumpFeatures but the file featureDump.txt is empty.

[k3yavi]

I imagine you are talking about the hto data, and that's expected. Basically, with citeseq Alevin is quantitying all the barcodes we see in the experiment. The idea is that the cells which are common with the rna must me most informative and other can be ignored. I think the problem is the low number of intersecting cells. Also, it seems you have only two experiments which are hashed through hto ? Ideally it should not matter but I wonder if the very low number of hto is triggering some corner case, I have to check the data to say more. I don't think I'd need the rna data, just the lost of cb quantified by Alevin is fine but if you can share the hto fastq I can check what's going on.

[ivanek]

Why this would be expected? I thought the most abundant cell barcodes in HTO library are suppose to represent real cells and the same should be present in GEX library. The molecule_info.h5 I took also from HTO library.

[k3yavi]

May be I am not understanding the question right but when you said "the cell barcodes provided by salmon alevin are actually scattered along the entire dataset (with among the less abundant).", I presume you are talking about the 400 CB common between RNA and the ADT. My argument is we are not connecting the CB from RNA and the CB from the ADT right, there is something else in the play and I'd not trust those 400 CB for now.

[ivanek]

I meant those barcodes reported by salmon alevin (approx. 97K) are scattered among the entire dataset if you consider the reported number of UMIs by cellranger as you measure for abundance.

I am preparing the files for you and will send you a private message (email) in few minutes.

[ivanek]

I sent you email from our nextcloud instance, I hope you received it. Please keep the data confidential, especially the RNASeq data are sensitive (derived from human sample).

[k3yavi]

Thanks @ivanek for sharing the data. Unfortunately the schedule is a little busy for the afternoon but I'll get back to you later in the day. However, currently the folder for the cloud link is empty but I am guessing it's still being uploaded ?

[ivanek]

All data should be there. I am looking forward to your feedback.

[k3yavi]

All right, I think the quants look fine to me, not sure what was happening on your end but I get ~7k common cells between RNA and the ADT. I use the following script:

library(tximport)

map <- read.table("~/3M-february-2018.txt")
rnaToADT <- map$V1
names(rnaToADT) <- map$V2

hto <- tximport("~/hto_out/alevin/quants_mat.gz", "alevin")$counts
rna <- tximport("~/rna_out/alevin/quants_mat.gz", "alevin")$counts

colnames(rna) <- rnaToADT[colnames(rna)]
length(intersect(colnames(rna), colnames(hto)))
# 7482

I went further to check if the HTOs are separable and that also looks fine to me along with the RNA clustering. You can check the R script -> script.R.zip. I am closing the issue feel free to reopen if you have any further questions.