Use of RNA-seq strand information in EGAPX pipeline

[BitaoQiu]

Dear EGAPX developers,

I would like to ask whether the EGAPX genome annotation pipeline takes the strand information of RNA-seq data into account.

From my understanding, STAR alignment does not require specifying strand information. However, for annotation tasks (e.g., assigning the correct orientation of ncRNAs or defining 3′ UTRs when they overlap with other genes), I would expect strand information to improve accuracy.

Could you please clarify whether strand information is incorporated into the EGAPX pipeline when stranded RNA-seq data are provided, or if I am missing something here?

Best regards,
Bitao

[murphyte]

Could you please clarify whether strand information is incorporated into the EGAPX pipeline when stranded RNA-seq data are provided, or if I am missing something here?

EGAPx includes a step rnaseq_divide_by_strandedness which analyzes a sampling of reads from each run for whether the spliced reads behave as if they were generated using a stranded sequencing method or not. That is, strandedness is computed from the alignments themselves so there's no need to provide this information in the configuration.

For runs that are stranded, all reads are treated as having that expected orientation in the downstream processing steps. Otherwise, only reads that are spliced or have polyA tails are considered to have a defined orientation, and the remaining unspliced reads are treated as ambiguously oriented. Ambiguously oriented alignments can contribute to models on either strand, whereas alignments with a defined orientation can only be used for model generation on that strand.

Here's a visual of two overlapping 3' UTRs where stranded RNA-seq helps to resolve them:
https://www.ncbi.nlm.nih.gov/gdv/browser/genome/?cfg=NCID_1_30521447_130.14.18.128_9146_1757349955_2071810285
model generation is using the RNA-seq alignments themselves, whereas this visual is transforming the alignments into coverage graphs, but you can see how in the top graph (the 3rd track) it's not really clear where the alignments for the left gene end, whereas that's a bit more clear with the bottom stranded track. Model generation is also trying to account for polyadenylation evidence seen in the RNA-seq evidence, which can result in a longer or shorter 3' UTR than what might be apparent in the coverage graphs.

[BitaoQiu]

@murphyte Thanks a lot for the quick and detailed explanation, much appreciated!

[Alteroldis]

Hi, @murphyte, I wanted to check how egapx recognized the strand-specificity in my libraries and looked at where the rnaseq_divide_by_strandedness process was writing in the working directory. I ended up finding only one file, run.strandedness, but it only contained a header like this:
#Run mate total_oount sense_count antisense_count strandedness
It is normal? All my rna-seq runs have RF type of strand-specificity.

[BitaoQiu]

@Alteroldis, you’re right. When I examine the code (https://github.com/ncbi/egapx/blob/main/nf/subworkflows/ncbi/rnaseq_short/bam_strandedness/main.nf), there should be more output than an empty file.

Additionally, when I checked the input for bam2asn, the run.strandedness file was also empty. This suggests that EPAGX didn’t implement the strand information in my case.

We should wait for @murphyte’s response after the government reopens....

[Alteroldis]

@BitaoQiu, and more then one file. As least three files should be in output directory of the task.
You tried to check previous version of egapx (0.4.0)?

[BitaoQiu]

@BitaoQiu, and more then one file. As least three files should be in output directory of the task. You tried to check previous version of egapx (0.4.0)?

I’m using version 0.4.0, but I’m not sure when this issue started....

[Alteroldis]

The same issue in v0.4.1. I haven't tested versions other than 0.4.0 and 0.4.1.

[Alteroldis]

The problem appears because of the sra_metadata.dat where when using a file with a list of fastqs, the run name is specified, and not the fastq file name as expected by rnaseq_divide_by_strandedness.
For example in my readlist:
rnaall_se1.1 S250015265.copy_L01_1.1.unpaired.fq.gz
Then bam file with alignments will be named genomefileprefix-S250015265.copy_L01_1.1.unpaired.fq-Aligned.out.Sorted.bam, but in sra_metadata.dat listed rnaall_se1.1.
I'll check in the next few days how this can be corrected in the egapx.py file. I think a quick solution is to create sra_metadata.dat file manually and specify the path to it in the input.yaml file with key short_reads_metadata.

[BitaoQiu]

The problem appears because of the sra_metadata.dat where when using a file with a list of fastqs, the run name is specified, and not the fastq file name as expected by rnaseq_divide_by_strandedness. For example in my readlist: rnaall_se1.1 S250015265.copy_L01_1.1.unpaired.fq.gz Then bam file with alignments will be named genomefileprefix-S250015265.copy_L01_1.1.unpaired.fq-Aligned.out.Sorted.bam, but in sra_metadata.dat listed rnaall_se1.1. I'll check in the next few days how this can be corrected in the egapx.py file. I think a quick solution is to create sra_metadata.dat file manually and specify the path to it in the input.yaml file with key short_reads_metadata.

Yes, you are right! When I provided SRR and in-house RNA data, EGAPX recognized and generated strand information for the SRR but not for my in-house data. It definitely needs to be fixed....

[BitaoQiu]

Just a quick update.

I managed to include strand information for local fastq files. I changed the local file names to SRR[0-9], and only provide *_[12].fastq files. This way, EGAPX can recognize and generate the strand information.

However, it seems that just by generating the strand information does not mean that EGAPX will use it. I checked the input for sam2asn, it seems that EGAPX will infer where the reads from by examining the *bin.bam file (merged from multiple bam files), so maybe only strand information from "true" SRR files (which contain file name in the fastq header) can be used...

[Alteroldis]

However, it seems that just by generating the strand information does not mean that EGAPX will use it.
Why did you think so?

With the readlist file contents as below and the modified egapx v0.4.1 (issue #165 ), I obtained the expected rnaseq_divide_by_strandedness task output and the workflow completed successfully:
uc.10dpd1_1.fq /path/to/fastq/uc.10dpd1_1.fq.gz
uc.10dpd1_1.fq /path/to/fastq/uc.10dpd1_2.fq.gz
The number of genes changed slightly, possibly due to trnascan task being enabled (disabled by default in the ui/assets/default_task_params.yaml, this is actually cmsearch), but the UTR length remained the same.
@BitaoQiu, did the UTR length change in your egapx-out/GNOMON/All.gnomon_evidence_stats.txt file after your tweak?

[BitaoQiu]

However, it seems that just by generating the strand information does not mean that EGAPX will use it.
Why did you think so?

With the readlist file contents as below and the modified egapx v0.4.1 (issue #165 ), I obtained the expected rnaseq_divide_by_strandedness task output and the workflow completed successfully: uc.10dpd1_1.fq /path/to/fastq/uc.10dpd1_1.fq.gz uc.10dpd1_1.fq /path/to/fastq/uc.10dpd1_2.fq.gz The number of genes changed slightly, possibly due to trnascan task being enabled (disabled by default in the ui/assets/default_task_params.yaml, this is actually cmsearch), but the UTR length remained the same. @BitaoQiu, did the UTR length change in your egapx-out/GNOMON/All.gnomon_evidence_stats.txt file after your tweak?

The reason I think EGAPX won’t utilize strand information unless derived from “real” SRR files is that after aligning and generating the strand information, EGAPX merges alignment files into bins and performs assembly using bam2asn. During this process, it uses the names of reads to infer their origins (e.g., read SRR123131.1 is from file SRR123131). Consequently, I also renamed the file names to SRR[0-9] and changed the read names with following codes in nf/subworkflows/ncbi/rnaseq_short/star_wnode/main.nf:

tmp="\${bam%.bam}.tmp.bam"   # Add sample ID prefix
samtools view -h "\$bam"  | awk -v pfx="${sampleID}." 'BEGIN{OFS="\\t"} /^@/ {print; next} { \$1 = pfx \$1; print }'  | samtools view -@ 4 -b -o "\$tmp" 
mv "\$tmp" "\$bam"

After integrating the strand information, I noticed a slight difference, but it wasn’t significant. I didn’t observe a substantial change in UTR length, with only a 2-base pair difference.