Skip to content

riboWaltz with yeast data #92

Description

@sopenaml

Hi,

I've been trying to run riboWaltz with yeast data with limited success. The issue as it has been discussed in this one and this one is the lack of UTR information in the GTF. I'm working with a custom GTF and Fasta files, generated by the lab I'm collaborating with. I have done alignment using STAR and producing "_Aligned.toTranscriptome.sorted" files. I have then created an annotation file as follows:

gtf_df    <- as.data.frame(import(gtf_file))
exon_gtf  <- subset(gtf_df, type == "exon" & !is.na(transcript_id))
tx_lengths <- aggregate(width ~ transcript_id, data = exon_gtf, FUN = sum)
colnames(tx_lengths) <- c("transcript", "l_tr")

### record the gene_id for each transcript so we can join with the GFF.
tx_gene <- unique(exon_gtf[, c("transcript_id", "gene_id")])
colnames(tx_gene) <- c("transcript", "gene_id")
tx_lengths <- merge(tx_lengths, tx_gene, by = "transcript", all.x = TRUE)

### get gene lengths from GFF
gff_df  <- as.data.frame(import(gff_file))
cds_df  <- subset(gff_df, type == "CDS")
if (nrow(cds_df) == 0) {
  stop("No CDS features found in GFF. Cannot build riboWaltz annotation.")
}
cds_df$gene_id <- as.character(cds_df$Parent)
cds_lengths <- aggregate(width ~ gene_id, data = cds_df, FUN = sum)
colnames(cds_lengths) <- c("gene_id", "l_cds")

### merge on gene_id

annot_df <- merge(tx_lengths, cds_lengths, by = "gene_id", all.x = TRUE)

###  If a transcript has no CDS match, default CDS length to full transcript length
annot_df$l_cds[is.na(annot_df$l_cds)] <- annot_df$l_tr[is.na(annot_df$l_cds)]

### (If UTRs aren't explicit, setting them to 0 forces riboWaltz to read from position 1)
annot_df$l_utr5 <- 0
annot_df$l_utr3 <- annot_df$l_tr - annot_df$l_cds - annot_df$l_utr5

final_columns <- c("transcript", "l_tr", "l_utr5", "l_cds", "l_utr3")
annotation_dt <- as.data.table(annot_df[, final_columns])

The problem with this annotation is that all UTRs are 0 and psite doesn't like that, but if I artificially I add:

annotation_dt[, l_utr5 := 30]
annotation_dt[, l_tr := l_utr5 + l_cds + l_utr3]

psite function runs ok. However, the lack of 3UTR causes an issue with functions such region_psite, as it selects transcripts with particular UTR length. I've tested altering the function to make l_utr3=0 and the function runs. My biggest issue is that when running metaprofile_psite, I get a metaprofile plot that is off centered and the signal in the 3' UTR is virtually gone (see file attached). I imagine is due to the artificially added UTR length, although I can't make full sense of it, as the peak is around nt, despite my artificial 5UTR value is 30 and "best offset: 12nts from the 5' end"

My questions are, is my approach valid with regards to the 5'UTR? I'm only after the PO per read length to use elsewhere, or do I have to create a new gtf with UTRs added to it, realigned and rerun the analysis, as it was suggested in previous issues. Since my alignment is to the genome, despite generating a transcriptome outputs I would have thought the realignment is not required, but just double checking here. Thank you for your help!

Image

Activity

  1. fabiolauria commented on Jul 7, 2026

    @fabiolauria
    Collaborator

    Hi there,
    thanks for using RiboWaltz, and thanks for taking the time to go through the previous issues before opening a new one.

    I think the safest and most robust approach would still be to modify the GTF by adding the artificial UTRs there, then realign the reads. This way, everything downstream remains fully consistent, and you avoid having to patch annotations or intermediate objects afterwards, with the risk of introducing inconsistencies.

    That said, if your only goal is to identify the P-site offsets, I think your approach of modifying annotation_dt could potentially be acceptable. The main issue I can think of, and which might explain the behavior you're observing, is related to the reads imported with bamtolist(). If I understood correctly, you only modified the annotation, so the 5' UTR and the transcript lengths are different from the original annotation (and I would also recommend extending the 3' UTR so that the total transcript length increases by 60 nt). However, the read coordinates (end5 and end3) obtained from the BAM are still referenced to the transcript coordinates without the artificial UTRs.

    To keep everything consistent, you would also need to shift of 30 nt end5 and end3 accordingly after importing the BAM files. Otherwise, the annotation and the read coordinates no longer refer to the same transcript coordinate system, which could explain why functions such as metaprofile_psite() produce shifted profiles.

    I hope this helps, and please let me know if I've misunderstood any part of your workflow.

    Best
    Fabio

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment

Metadata

Metadata

Assignees

No one assigned

    Labels

    No labels
    No labels

    Type

    No type

    Projects

    No projects

      Milestone

      No milestone

      Relationships

      None yet

      Development

      No branches or pull requests

      Issue actions