Differential expression at transcript level

[fpsanz]

Description of feature

Hello,

First of all, congratulations on your pipeline for differential abundance. I'm trying to perform a differential expression analysis at the transcript level, but I can't seem to get it to work. I've tried different parameters and input files, but I'm not succeeding. The data comes from the output of the nf-core/rnaseq pipeline. Is it possible to carry out such an analysis? Thanks,

Regards.

Fernando

[WackerO]

Hey Fernando, you should definitely be able to run the pipeline on rnaseq output, but unless you share some information about your command and the input files you used as well as the error message you received, I'm afraid I won't be able to help...

[RaverJay]

I am also interested in this. nf-core/rnaseq running star+salmon already produces transcript abundances as output, and the table has the columns when run against ensemble genome+annotation: tx gene_id <sample1> <sample2> etc

with tx having the ensembl transcript ids.

Now if you try to plug in this table + the lengths into differentialabundance it fails on validation:

-[nf-core/differentialabundance] Pipeline completed with errors-
ERROR ~ Error executing process > 'NFCORE_DIFFERENTIALABUNDANCE:DIFFERENTIALABUNDANCE:VALIDATOR (samplesheet.csv)'

Caused by:
  Process `NFCORE_DIFFERENTIALABUNDANCE:DIFFERENTIALABUNDANCE:VALIDATOR (samplesheet.csv)` terminated with an error exit status (1)

Command executed:

  validate_fom_components.R \
      --sample_metadata "samplesheet.csv" \
      --feature_metadata 'Homo_sapiens.anno.tsv' \
      --assay_files "salmon.merged.transcript_counts.tsv" \
      --contrasts_file "contrasts.csv" \
      --output_directory "study" \
      --sample_id_col 'sample' --feature_id_col 'gene_id'

<SNIP>

  [1] "Reading sample sheet at samplesheet.csv with ID col sample"
  [1] "Reading feature metadata at Homo_sapiens.anno.tsv with ID col gene_id"
  [1] "Reading assay matrix salmon.merged.transcript_counts.tsv and validating against samples and features (if supplied)"
  Error in `.rowNamesDF<-`(x, value = value) : 
    duplicate 'row.names' are not allowed
  Calls: validate_inputs ... row.names<- -> row.names<-.data.frame -> .rowNamesDF<-
  In addition: Warning message:
  non-unique values when setting 'row.names': ‘ENSG00000000003’, ‘ENSG00000000005’, ‘ENSG00000000419’, ‘ENSG00000000457’, ‘ENSG00000000460’, ‘ENSG00000000938’, ‘ENSG00000000971’, ‘ENSG00000001036’, ‘ENSG00000001084’, ‘ENSG00000001167’, ‘ENSG00000001460’, ‘ENSG00000001461’, ‘ENSG00000001497’, ‘ENSG00000001617’, ‘ENSG00000001626’, ‘ENSG00000001629’, ‘ENSG00000001630’, ‘ENSG00000001631’, ‘ENSG00000002016’, ‘ENSG00000002330’, ‘ENSG00000002549’, ‘ENSG00000002586’, ‘ENSG00000002587’, ‘ENSG00000002726’, ‘ENSG00000002745’, ‘ENSG00000002746’, ‘ENSG00000002822’, ‘ENSG00000002834’, ‘ENSG00000002919’, ‘ENSG00000002933’, ‘ENSG00000003056’, ‘ENSG00000003096’, ‘ENSG00000003137’, ‘ENSG00000003147’, ‘ENSG00000003249’, ‘ENSG00000003393’, ‘ENSG00000003400’, ‘ENSG00000003402’, ‘ENSG00000003436’, ‘ENSG00000003509’, ‘ENSG00000003756’, ‘ENSG000000 [... truncated] 
  Execution halted
  INFO:    Cleaning up image...

because the feature_id_col should not be 'gene_id', but 'tx' for this input table

Is there a way to set this parameter? Or is it set via the params in the 'Features options'?

[RaverJay]

Okay, so setting

--features_id_col             'transcript_id'
--features_type               'transcript'

made the validation succeed because features_id_col is used:

validate_fom_components.R \
    --sample_metadata "samplesheet.csv" \
    --feature_metadata 'Homo_sapiens.anno.tsv' \
    --assay_files "salmon.merged.transcript_counts.tsv" \
    --contrasts_file "contrasts.csv" \
    --output_directory "study" \
    --sample_id_col 'sample' --feature_id_col 'transcript_id'

the following filter step seems to work, yielding such a table:

transcript_id   EU10_R4_Ctrl_wDPO   EU11_R4_Wt_oDPO EU12_R4_Wt_wDPO EU1_R2_Ctrl_oDPO    EU2_R2_Ctrl_wDPO    EU3_R2_Wt_oDPO  EU4_R2_Wt_wDPO  EU5_R3_Ctrl_oDPO    EU6_R3_Ctrl_wDPO    EU7_R3_Wt_oDPO  EU8_R3_Wt_wDPO  EU9_R4_Ctrl_oDPO
ENST00000456328 0   0   0   0   0   0   0   0   0   0   0   1.65
ENST00000488147 112.211 78.767  74.328  63.958  41.027  119.846 73.367  53.181  60.653  47.144  54.7959.419

but it still dies on the DESeq2 process, likely because that still looks for a 'gene_id' column in the table, which is not there. 'gene_id_col' is not passed as a param to the DESeq2 process:

  ################################################
  ################################################
  ## PARSE PARAMETERS FROM NEXTFLOW             ##
  ################################################
  ################################################

  # I've defined these in a single array like this so that we could go back to an
  # optparse-driven method in future with module bin/ directories, rather than
  # the template

  # Set defaults and classes

  opt <- list(
      output_prefix = ifelse('all' == 'null', 'Ctrl_wDPO_vs_Ctrl_oDPO', 'all'),
      count_file = 'study.filtered.tsv',
      sample_file = 'samplesheet.sample_metadata.tsv',
      contrast_variable = 'condition',
      reference_level = 'Ctrl_oDPO',
      target_level = 'Ctrl_wDPO',
      blocking_variables = NULL,
      control_genes_file = '',
      transcript_lengths_file = 'salmon.merged.transcript_lengths.tsv',
      sizefactors_from_controls = FALSE,
      gene_id_col = "gene_id",
      sample_id_col = "experiment_accession",

<SNIP>

  )
  opt_types <- lapply(opt, class)

  # Apply parameter overrides

  args_opt <- parse_args('--sample_id_col "sample" --test Wald --fit_type parametric --sf_type ratio --min_replicates_for_replace 7 --use_t false --lfc_threshold 0 --alt_hypothesis greaterAbs --independent_filtering true --p_adjust_method BH --alpha 0.1 --minmu 0.5 --vs_method vst --vs_blind true --vst_nsub 1000 --shrink_lfc true --cores 1 --subset_to_contrast_samples false --blocking_variables')

Ideas?

[RaverJay]

Got it to work by renaming the column names in the counts table, lengths table AND annotation, such that the column containing the transcript_id is (wrongly) called 'gene_id'

renaming:

sed '1 s/gene_id/gene_id_orig/;s/tx/gene_id/' results_rnaseq/star_salmon/salmon.merged.transcript_counts.tsv > renamed_salmon.merged.transcript_counts.tsv
sed '1 s/gene_id/gene_id_orig/;s/tx/gene_id/' results_rnaseq/star_salmon/salmon.merged.transcript_lengths.tsv > renamed_salmon.merged.transcript_lengths.tsv
sed '1 s/gene_id/gene_id_orig/;s/transcript_id/gene_id/' results_differentialabundance/tables/annotation/Homo_sapiens.anno.tsv > renamed_Homo_sapiens.anno.tsv

running:

nextflow run nf-core/differentialabundance \
    -resume \
    -profile rnaseq,singularity \
    --input samplesheet.csv \
    --contrasts contrasts.csv \
    --matrix renamed_salmon.merged.transcript_counts.tsv \
    --transcript_length_matrix renamed_salmon.merged.transcript_lengths.tsv \
    --outdir results_transcriptsdiff_renamed \
    --features renamed_Homo_sapiens.anno.tsv

Output obviously has 'Gene ID' column name even though they are transcript ids.

Ideally there would be a way to properly parameterize which columns to use (in validation, filtering, and DESeqing) AND how to generate the annotation from the gtf

Please consider implementing it - would be a cool feature to be able to directly use transcript abundances from nf-core/rnaseq

Best

[WackerO]

Hey there!

but it still dies on the DESeq2 process, likely because that still looks for a 'gene_id' column in the table, which is not there. 'gene_id_col' is not passed as a param to the DESeq2 process:

Yep...I noticed that as well some time ago, that's a silly little bug. It's already fixed on the dev branch, so if you run that code it should work. The necessary id_col parameters should already be in the pipeline :) (There's quite a few of them unfortunately, but all of them are described on https://nf-co.re/differentialabundance/1.4.0/parameters)

[pinin4fjords]

@WackerO has indicated a fix present on the dev branch, and there has been now response by the OP since. I'm assuming this is resolved, @fpsanz please reopen if you discover that the issue persists.