Add rules for generating gene expression & count tables

Description

For downstream analysis (most importantly for differential gene expression analysis) it would be desirable to generate summary tables for expression (TPM) and projected raw counts for each feature (gene, transcript, transcript region and/or exon).

Possible implementations

A single table could be created individually for each sample and feature type. Alternatively or in addition, a single common table could be created for all samples in a given run, each for expression and projected raw counts and for each of the feature types.

Additional context

In this regard it might be worthwhile to discuss the implementation of a results database that allows, e.g., to check the expression of genes of interest across all samples that have been run in the past.

Related issues

It is possible that #23 (closed) already addresses this issue. #30 (closed) includes a rule for collecting Salmon transcript expression estimates for all samples in a run and compiling a common table.

changed the description

Should be partially fixed by #30 (closed)

changed the description

Thanks, I referenced the wrong issue. Fixed it now. The thing is that #30 (closed) is still open and doesn't include any additional information apart from the rule. I think there's still some room for discussion here on what we would really like/need as an output for downstream analysis (e.g., what is our level of analysis: genes, transcripts, exons, CDS, UTRs, all of it?), what MultiQC and other tools possibly provides etc. That's why I opened this issue. But perhaps it's better to copy the rule in here as a comment and close #30 (closed) to bundle everything here.

Rule to collect all TPM values from Salmon, courtesy of @gypas:

################################
### Filter salmon (TPM)
################################################################################

rule filter_salmon_TPM:
    input:
        salmon_genes_out = os.path.join(config["output_dir"], "{sample}", "salmon", "quant_reads", "quant.genes.sf"),
    output:
        salmon_filtered = os.path.join(config["output_dir"], "{sample}", "salmon", "quant_reads", "quant.genes.TPM.tsv")
    log:
        os.path.join(config["local_log"], "filter_salmon_TPM_{sample}.log")
    run:
        df = pd.read_csv(input.salmon_genes_out, header=0, sep="\t")
        df = df[["Name", "TPM"]].copy()
        df.set_index("Name", inplace=True)
        df.to_csv(output.salmon_filtered, header=True, sep="\t", index=True)

################################################################################
### Collect all salmon reads TPM
################################################################################

rule collect_salmon_reads_TPM:
    input:
        salmon_filtered = expand(os.path.join(config["output_dir"], "{sample}", "salmon", "quant_reads", "quant.genes.TPM.tsv"), sample=get_samples())
    output:
        salmon_filtered_all = os.path.join(config["output_dir"], "salmon", "quant_reads", "quant.genes.TPM.tsv")
    params:
        sample_name = expand("{sample}", sample=get_samples())
    log:
        os.path.join(config["local_log"], "collect_salmon_reads_TPM.log")
    run:
        df_merge = pd.DataFrame()
        for i in range(len(input.salmon_filtered)):
            df = pd.read_csv(input.salmon_filtered[i], header=0, sep="\t")
            df.columns = ["Name", params.sample_name[i]]
            if df_merge.empty:
                df_merge = df.copy()
            else:
                df_merge = pd.merge(df_merge, df, on="Name")
        df_merge.to_csv(output.salmon_filtered_all, header=True, sep="\t", index=False)

Done.

Thanks. We can also use the inbuilt option from salmon:

salmon quantmerge --help

And here is an example rule:

rule salmon_quantmerge:
    input:
        salmon_in = expand(os.path.join(config["output_dir"], "{sample}", "salmon", "quant_reads", "quant.sf"), sample=get_samples())
    output:
        salmon_out = os.path.join(config["output_dir"], "summary", "salmon", "quantmerge", "tpm.tsv")
    params:
        salmon_dir = expand(os.path.join(config["output_dir"], "{sample}", "salmon", "quant_reads"), sample=get_samples()),
        sample_name_list = expand("{sample}", sample=get_samples())
    log:
        os.path.join(config["local_log"], "salmon_quantmerge.log")
    threads:    1
    singularity:
        "docker://zavolab/salmon:0.11.0"
    shell:
        "(salmon quantmerge \
        --quants {params.salmon_dir} \
        --names {params.sample_name_list} \
        --column tpm \
        --output {output.salmon_out}) &> {log}"

assigned to @devagy74

added 1 deleted label

changed milestone to %First usable version

removed milestone

@gypas Give us a chance to finish that first milestone We can open another one for missing rules

added To Do label

changed milestone to %v0.1.0 release

changed title from Generate gene expression & count tables to Add rules for generating gene expression & count tables

mentioned in merge request !30 (merged)

Partly resolved by !30 (merged) (TPM but not counts)

mentioned in merge request !33 (merged)

closed via merge request !33 (merged)

mentioned in commit f31261b5

mentioned in issue #104 (closed)