Teichlab / bracer

BraCeR - reconstruction of B cell receptor sequences from single-cell RNAseq data
Other
41 stars 22 forks source link

Add option to skip transcriptome quantification #17

Open dcroote opened 6 years ago

dcroote commented 6 years ago

For those interested in only the assembled sequences, this option can substantially reduce runtime by skipping the Kallisto build transcriptome index step. Note Kallisto is still used so as to not break anything downstream, but the index is built from only the assembled sequences rather than entire transcriptome and assembled sequences.

For bracer test adding --no_transcriptome_quant decreased the runtime from 311s to 60s (80% decrease) on my machine.

mstubb commented 6 years ago

Thanks for this. Can you show that this gives equivalent quantification for the assembled sequences when compared with making an index from the whole transcriptome? It was always my concern that having so few sequences in the index would violate assumptions in the Kallisto/Salmon quantification models and so would give spurious quants.

For an alternative way to solve this same problem, have a look at the -small_index option for TraCeR (https://github.com/Teichlab/tracer#options). Perhaps we should try to bring that into BraCeR and unify some functionality. Interested to hear your thoughts.

PS - apologies for not accepting your CI pull request. I've been travelling and haven't had chance to look at it. Hoping to do so this week.

dcroote commented 6 years ago

@mstubb - the quantification will be different. My thought is that this option would be for those who are performing mapping / quantification by another method e.g. STAR + htseq and are not interested in the pseudoalignment transcriptome counts. We could explicitly write TPM: N/A or similar via the get_summary core function to make the output explicit in stating that the TPM will not be reliable if --no_transcriptome_quant is used

mstubb commented 6 years ago

That makes sense although worth noting that the quantifications are used to filter the reconstructed sequences in cases where more than two recombinants are detected for a particular locus (eg 3 Igh sequences). The two most highly expressed are assumed to be the 'correct' ones. If we implemented this we'd need to show that the ranking by expression of the recombinants was the same even if the TPMs were wildly different. Or, only allow this mode to be used without attempting any filtering.

dcroote commented 6 years ago

I agree that filtering is useful and agree that we would want to demonstrate that ignoring the transcriptome doesn't affect recombinant selection.

To verify this I've created a test where I've added the following to the fastq reads of cell1 within test_data:

The TPM results from filtered_BCR_seqs/filtered_BCRs.txt:

assembly IGHV3-30 IGKV2-28 IGLV1-40 IGLV2-14
with transcriptome 126425.0 33325.5 134683.0 77786.6
without transcriptome 290882.0 116700.0 387686.0 204732.0

The raw TPM numbers are much higher without the transcriptome as expected, but when normalized by their TPM sum they are quite comparable:

assembly IGHV3-30 IGKV2-28 IGLV1-40 IGLV2-14
with transcriptome 0.34 0.09 0.36 0.21
without transcriptome 0.29 0.12 0.39 0.2