scrnaseq icon indicating copy to clipboard operation
scrnaseq copied to clipboard

Published 20 hours ago • verified publisher icon nf-core

Output 10x counts

Open kafkasl opened this issue 2 years ago • 13 comments

This PR adds support to generate 10x count files as output (features.tsv, barcodes.tsv, and matrix.mtx) as part of the pipeline.

Issue: https://github.com/nf-core/scrnaseq/issues/66

PR checklist

  • [x] This comment contains a description of changes (with reason).
  • [ ] If you've fixed a bug or added code that should be tested, add tests!
  • [ ] If you've added a new tool - have you followed the pipeline conventions in the contribution docs- [ ] If necessary, also make a PR on the nf-core/scrnaseq branch on the nf-core/test-datasets repository.
  • [x] Make sure your code lints (nf-core lint).
  • [x] Ensure the test suite passes (nextflow run . -profile test,docker --outdir <OUTDIR>).
  • [ ] Usage Documentation in docs/usage.md is updated.
  • [ ] Output Documentation in docs/output.md is updated.
  • [x] CHANGELOG.md is updated.
  • [x] README.md is updated (including new tool citations and authors/contributors).

kafkasl avatar Oct 11 '22 16:10 kafkasl

nf-core lint overall result: Passed :white_check_mark:

Posted for pipeline commit aae4f5e

+| ✅ 158 tests passed       |+

:white_check_mark: Tests passed:

Run details

  • nf-core/tools version 2.6
  • Run at 2022-11-17 14:11:35

github-actions[bot] avatar Oct 11 '22 16:10 github-actions[bot]

Two questions:

  • isn't this essentially throwing away sample information? I see you are only storing the cell index in the barcodes.tsv.
  • I'm afraid this will be quite slow when the merged dataset get's large (i.e. > 100k cells). This is one of the reasons to prefer binary formats like h5 over text representation. Maybe not run this by default, or at least make it easy to switch off.

grst avatar Oct 12 '22 06:10 grst

Regarding 1, the barcodes file contains the actual barcodes list, it looks like this:

AAACCTGAGCTAGTGG
AAACCTGCATCGACGC
...

Concerning the second point, do you know which part is the one making it slow? I agree to not run it by default. I'll add a new parameter.

kafkasl avatar Oct 12 '22 08:10 kafkasl

Regarding 1, the barcodes file contains the actual barcodes list, it looks like this

yes, so how will you actually know which barcode comes from which biological sample?

Concerning the second point, do you know which part is the one making it slow? I agree to not run it by default. I'll add a new parameter.

Just writing a text file vs. writing a binary file. Tbh, I don't really have experience with writing these files, but reading feat/barcode/mtx is way slower than reading the h5 files produced by cellranger.

grst avatar Oct 12 '22 09:10 grst

Regarding 1, the barcodes file contains the actual barcodes list, it looks like this

yes, so how will you actually know which barcode comes from which biological sample?

Because those 3 files will be generated on a per-sample basis. This is the folder structure:

├── Sample_X_matrix
│   ├── barcodes.tsv
│   ├── features.tsv
│   └── matrix.mtx
└── Sample_Y_matrix
    ├── barcodes.tsv
    ├── features.tsv
    └── matrix.mtx

Concerning the second point, do you know which part is the one making it slow? I agree to not run it by default. I'll add a new parameter.

Just writing a text file vs. writing a binary file. Tbh, I don't really have experience with writing these files, but reading feat/barcode/mtx is way slower than reading the h5 files produced by cellranger.

Yeah, unfortunately for some integrations we need those files in text format.

I have added the parameter --output_10x and it's disabled by default.

kafkasl avatar Oct 12 '22 09:10 kafkasl

As mentioned on slack, no new functionality should be needed to produce these files - Maybe some publishDir options need to be tuned to get them into the right places in the output directory.

All the aligners produce mtx files already and we use that python script to build h5ad files from them.In the example output of the workflow:

grst avatar Oct 13 '22 05:10 grst

I got your point. However, I have checked those files and there are two problems with them.

  1. The matrix is transposed so it will failed when trying to load them in R
  2. The genes files contains the ensemble IDs but we want that file to include also the gene names (which we are doing in our code).

My goal is to format correctly those files independently of the aligner so we can add another step to upload them automatically for downstream analysis. I am wary of adding those "R formatting" options into the actual aligners' steps to avoid polluting them with formatting. How would you do that?

kafkasl avatar Oct 13 '22 10:10 kafkasl

Actually, maybe this could be simplified and at the same time used to solve https://github.com/nf-core/scrnaseq/issues/159:

  • Add gene symbols already in anndata (https://github.com/nf-core/scrnaseq/issues/159)
  • Add an additional flag --export_mtx to @fmalmeida's mtx_to_h5ad.py script, to enable mtx export on demand.
  • Output mtx directly from that python script. Like that we can profit from the file already being loaded into memory and there is very little complexity added (no more processes)

At the same time we should think about how we organize the output directory. We currently have the mtx_conversions output directory for each aligner with the h5ad and rds files in them. This is fine, but I think the name is not ideal and it needs some substructure, e.g.

aligner_cellranger/matrices
|
| - per_sample
|   | 
|    - h5ad
|    - mtx
|    - seurat
  - merged
    |
     - h5ad
     - seurat

@apeltzer, @fmalmeida, @kafkasl, what do you think?

grst avatar Oct 14 '22 08:10 grst

Well, I agree that having it standard so they can used in the same way afterwards is a good idea. I gues, for example, instead of adding this export on the pyScript, we can actaully add a step before, which does this standardization on files, and "saves" it and then the other conversions come from these standardised matrices and we work on re-shaping the conversion modules and scripts to use that new (standard) formatting.

Maybe this makes more sense, no?

At the same time we should think about how we organize the output directory. We currently have the mtx_conversions output directory for each aligner with the h5ad and rds files in them. This is fine, but I think the name is not ideal and it needs some substructure, e.g.

On this second comment, I totally agree that we should reshape it and I have no comments on it. I liked the structure proposed.

fmalmeida avatar Oct 14 '22 09:10 fmalmeida

export on the pyScript, we can actaully add a step before, which does this standardization on files, and "saves" it and then the other conversions come from these standardised matrices and we work on re-shaping the conversion modules and scripts to use that new (standard) formatting.

It depends a bit what needs to be done. Reading/writing mtx files in a Python script is way slower than h5(ad) files. Not reading (i.e. using the data already in memory) is even better. So purely from a runtime perspective, it is beneficial to read whatever output files the aligners create once and then write all desired outputs.

That being said, transposing a mtx matrix could probably done on the command line using awk, which would be very fast.

grst avatar Oct 14 '22 10:10 grst

I agree that we should not read often but once and output what is necessary / standardize this a bit. The conversion modules also need to add versions for example, so would be great to do all of the above to make sure we're also following best practices 👍

If that even closes more issues even better 🥳

apeltzer avatar Oct 14 '22 10:10 apeltzer

I think what you mentioned @grst sounds good. We have an idea of how to add the gene symbols and the --export_mtx parameter. We will have to look into how to change the output structure. I don't know what @apeltzer refers to with "adding versions".

kafkasl avatar Oct 17 '22 07:10 kafkasl

The current mxt_... modules do not output the required versions.yml which should be done too 👍🏻 Thats something that was missed but is also highlighted by the nf-core lint tool, so we should adhere to the standard of generating this when the conversion is done :-)

apeltzer avatar Oct 17 '22 07:10 apeltzer

@grst @apeltzer I've modified the PR to do the following:

  1. Standardize the mtx_conversions outputs. They now look like this:
├── mtx_conversions
│   ├── Sample_X
│   │   ├── Sample_X_matrix.h5ad
│   │   ├── Sample_X_matrix.rds
│   │   ├── barcodes.tsv
│   │   ├── features.tsv
│   │   └── matrix.mtx
│   ├── Sample_Y
│   │   ├── Sample_Y_matrix.h5ad
│   │   ├── Sample_Y_matrix.rds
│   │   ├── barcodes.tsv
│   │   ├── features.tsv
│   │   └── matrix.mtx
  1. Export the 10x counts format files, barcodes.tsv, features.tsv, and matrix.mtx. This has been done inside an existing module from data already in memory so there shouldn't be any performance impact other than writing to disk.
  2. I've enriched the features.tsv files with the gene names, by default they only have the gene IDs. I extracted them from txp2gene for kallisto, and the geneInfo.tab file for star & cellranger. For alevin, we haven't managed to find where to get that translation info so far.
  3. I added the export the 10x counts param --export_mtx that you suggested but setting it to false breaks the downstream process mtx_to_seurat which depends on this matrix counts being exported so I think it should not be added.

as a side note, the prettier command and nf-core schema build have conflicting formatting for nextflow_schema.json

kafkasl avatar Nov 11 '22 08:11 kafkasl

@grst I've addressed all your comments. The most important ones:

  • mtx_to_h5ad now exports a versions file
  • cellrange_mtx_to_h5ad has been deleted and now uses also mtx_to_h5ad
  • improved the parsing of t2g-like files with your suggestions

kafkasl avatar Nov 15 '22 09:11 kafkasl