I have a list of about 3,000 genes and I want to retrieve the genomic locus of each gene's 3' UTR. My goal is to screen all the sequence that could be part of a gene's 3' UTR, even if some of that sequence isn't always in the 3' UTR of every transcript.
I've tried the methods in a few previous answers and, for about 80% of the genes in my list, I can find the predicted start and end of the 3' UTR using Biomart or the UCSC table browser. The problem is that for each gene I get multiple results (for all the alternative transcripts of that gene), and in each result the 3' UTR starts and ends at a different place. What I would like is the site of the most upstream 3' UTR start and the most downstream 3' UTR end that have been predicted for a given gene.
Does anyone know a straightforward way to get these from UCSC or Biomart? Can I perhaps get the shortest predicted CDSend and longest predicted transcription end?
Thanks for your help!
•
written
by
max •
2 answers
What I would like is the site of the most upstream 3' UTR start and the most downstream 3' UTR end that have been predicted for a given gene.
Given the format of your posted snippet of output, you could use the following Python script to get a JSON object of gene IDs and most-upstream and most-downstream start and stop positions, respectively:
#!/usr/bin/env python
import sys
import json
def main():
c = 0
d = {}
for l in sys.stdin:
l = l.strip()
if c > 0:
try:
(gsi, tsi, gn, start, end) = l.split('\t')
if gsi not in d:
d[gsi] = {}
d[gsi]['hgnc'] = gn
if 'start' not in d[gsi] or d[gsi]['start'] < start:
d[gsi]['start'] = start
if 'end' not in d[gsi] or d[gsi]['end'] > end:
d[gsi]['end'] = end
except ValueError as ve:
pass
c += 1
sys.stdout.write("%s\n" % json.dumps(d, indent=2))
if __name__ == "__main__":
main()
On your snippet:
$ snippet2positions.py < snippet.txt
{
"ENSG00000115760": {
"start": "32842972",
"end": "32667294",
"hgnc": "BIRC6"
},
"ENSG00000102189": {
"start": "93251207",
"end": "93169786",
"hgnc": "EEA1"
},
"ENSG00000124171": {
"start": "49373281",
"end": "49370278",
"hgnc": "PARD6B"
}
}
You could write the output to a file like so:
$ snippet2positions.py < snippet.txt > positions.json
You can modify the format of this output by changing the line:
sys.stdout.write("%s\n" % json.dumps(d, indent=2))
to something else.
You could try BEDTools 'merge' to collapse overlapping UTRs into one, but I anticipate a few problems: 1) if any of the UTRs do not overlap, you'll still be left with multiples; 2) UTRs from different genes might overlap, and those would collapse as well.
Log in to answer this question.
-
identifying 3'utrs for multiple transcripts of genes from a summarized experiment.
written by RNAseqer
-
Biomart issue, why so few 3'utrs?
written by RNAseqer
-
utr sequence retrieval from UCSC
written by kim •
-
Choosing a single consensus promoter/transcript for each gene (Mouse)
written by srhic •
-
genomic features 3'utr and 5'utr
written by ginlucks •
-
retrieving refseq/ucsc transcription start site with 1kb up- and downstream flank?
written by taniasultana2004 •
-
How to retrieve "UTR" sequence from "Exon" of "Canonical Transcript" from a list of gene symbols?
written by Hyun Jin •
-
Obtaining All Cds Sequences (I.E. All Spliced Exon Variants) From Ucsc
written by Max
-
Exon Phase And Start Codon
written by Max
-
Biomart Refseq Data Problem
written by Pppp •
Can you post a snippet (
head) of your genes and UTR datasets? I have a solution in mind, but it would depend on the format of your datasets. Feel free to let us know if you'd like a hand. Cheers.Thanks for your help!
Here's a snippet of the output from Biomart.
I select dataset>'Human genes (GRCH37.p13)', Filters>'Gene Names' (here I paste a list of gene names, but I have refseq mRNA IDs too), Attributes>'Gene stable ID, Transcript stable ID, Gene name, 3' UTR start, 3' UTR end', Results>'unique results only'