This is a test version of Biostars. For the public version, visit https://www.biostars.org.
How do I consider that my reads might map to the 3' UTR (3'-tag RNA seq)?

I used a 3'-tag library kit and sequenced RNA samples. I have read that because of the library preparation method, the reads will be biased towards the 3' end, and may even fully/mostly map to the 3'-UTR. I have a T2T genome (although I don't know how complete the annotation is), but I don't know how to utilize the 3' UTR during the assignment step. Currently ~30% of my reads are unassigned - unmapped (featurecounts).

The gff version of the annotation has 3'-UTR:

cut -f3 annotation.gff3 | sort | uniq -c
output:  319956 CDS
 390123 exon
  63954 five_prime_UTR
  44375 gene
 323672 intron
   2939 lnc_RNA
  63512 mRNA
  64065 three_prime_UTR

But when converting to a gtf I only have exon and cds: cut -f3 annotation.gtf | sort | uniq -c

output: 319956 CDS 390123 exon

I am using STAR and featurecounts to map my reads. I also read that I could use bedtools slop to extend features downstream, which may also address this issue.

Any thoughts or advice?

rna-seq

1 answer

The records in your gff3 are overlapping - CDS, three_prime_UTR and five_primer_UTR are all components of the exon. If you were to visualise it, it would look like this:

           gene |>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>|
     transcript |>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>|
           exon |>>>>>>>|    |>>>>|     |>>>>>>>>>>>>>>>>|
            CDS      |>>|    |>>>>|     |>>|
 five_prime_utr |>>>>|
three_prime_utr                            |>>>>>>>>>>>>>|

Your conversion program has carried over only the exon and CDS entries, but actaully, that is all that is neccessary, as you can get all the information you need from the these:

           exon |>>>>>>>|    |>>>>|     |>>>>>>>>>>>>>>>>|
            CDS      |>>|    |>>>>|     |>>|

Almost all quantification programs will only use exon entries in counting, which is fine, it means that as long as the annotation is correct you will capture all the reads from the 3' UTR, as the 3' UTR is part of the exons.

I would guess your unassigned reads come from three places:

  1. Missmapping: Reads have aligned in intergenic space that should have been mapped to a gene.
  2. Genomic background: carry over in the library prep process.
  3. oligo-dT priming at intronic homopolyer tracks (3'-end sequencing uses an oligo-dT primer to capture polyA tails, but it turns out that polyA tails arn't the only places you find runs of As in the genome.
  4. Missannotaiton: Even in the human genome, the ends of genes are not perfectly annotated. Nor is every poly-adenylated lncRNA in the annotation. You could be either getting reads from unannotated 3' ends of annotated genes, or even reads from entirely unannotated genes. If you've got them, we've previously used the read2s from 3'-end sequencing to reannotated the UTRs of species with less than complete 3'-end annotation.

Log in to answer this question.