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?
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:
- Missmapping: Reads have aligned in intergenic space that should have been mapped to a gene.
- Genomic background: carry over in the library prep process.
- 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.
- 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.
What program did you use for the conversion? I suggest trying AGAT (https://github.com/NBISweden/AGAT ) to do the conversion.