Finding exons and introns sequences

Question

[elisheva]

Hello everybody!!
I have to analyze exons' and intron's sequences of many organisms.
My question is: what is the most efficient way to retrieve all those sequences in fasta format?
Or in other words: which Database holds accessible information about exons and introns sequences?
I thought to download gff files of all the organisms, filter exons, and introns and then get their sequences by using bedtools get fasta.
But this process requires to download many genomes and seems not to be effective.
Any suggestion for this purpose?

[vkkodali_ncbi]

How many genomes are you working with? Will you be downloading all intron and exon sequences for each genome? If so, where are the getting the coordinates from? For RefSeq data, you can download gff3 files, parse them for intron and exon coordinates and use edirect to download sequences in fasta format. But edirect would be an inefficient way to do this if you want to download sequences for all of the introns and exons. Downloading the entire genome sequence in fasta format to disk first would be much more efficient.

[elisheva]

I am working on 80 organisms from Ensembl db. So this is not practical to download their full genome. Besides, on their gff3/gtf files, there are no introns at all ):

[Pierre Lindenbaum]

Besides, on their gff3/gtf files, there are no introns at all

prokaryotes ?

[elisheva]

of course not. I am looking only on mammals.

Answer 1

[vkkodali_ncbi]

You can extract the intron features from a GFF3 file using the script shown here:

	#!/usr/bin/env python
	import csv
	import os
	import sys
	import gzip
	from collections import defaultdict
	import argparse
	# See http://stackoverflow.com/questions/14207708/ioerror-errno-32-broken-pipe-python
	from signal import signal, SIGPIPE, SIG_DFL
	signal(SIGPIPE, SIG_DFL)
	parser = argparse.ArgumentParser(description ="""This script parses input
	gff3 file and generates a table of introns""")
	parser.add_argument('-i', '--infile', help="input file; default STDIN")
	parser.add_argument('-o', '--outfile', help="output file; default STDOUT")
	parser.add_argument('-z', '--gzip', help="input is gzip compressed", \
	action = 'store_true')
	args = parser.parse_args()
	if args.infile:
	if args.gzip:
	gff3_file = gzip.open(args.infile, 'rt')
	else:
	gff3_file = open(args.infile, 'r')
	else:
	gff3_file = sys.stdin
	if args.outfile:
	introns_file = open(args.outfile, 'w')
	else:
	introns_file = sys.stdout
	def get_exons(gff3_file):
	tbl = csv.reader(gff3_file, delimiter = '\t')
	exons_dict = defaultdict(list)
	for line in tbl:
	if not line[0].startswith('#'):
	[
	chrom,
	feat_source,
	feat_type,
	start,
	stop,
	score,
	strand,
	phase,
	attribs
	] = line
	if (feat_type == "exon"
	and 'transcript_id=' in attribs):
	new_attribs = process_attribs(attribs)
	if 'GeneID' in new_attribs:
	gene_id = new_attribs['GeneID']
	elif 'gene' in new_attribs:
	gene_id = new_attribs['gene']
	else:
	gene_id = 'UNKNOWN'
	tx = new_attribs['transcript_id']
	start, stop = int(start), int(stop)
	exons_dict[(chrom, strand, gene_id, tx)].append((start, stop))
	elif(feat_type == "exon"
	and 'Parent=transcript' in attribs):
	attribs = attribs.split(';')
	new_attribs = {}
	for attrib in attribs:
	attrib = attrib.split('=')
	new_attribs[attrib[0]] = attrib[1]
	gene_id = new_attribs['Parent'].strip('transcript:')
	tx = new_attribs['Parent'].strip('transcript:')
	start, stop = int(start), int(stop)
	exons_dict[(chrom, strand, gene_id, tx)].append((start, stop))
	gff3_file.close()
	return exons_dict
	def process_attribs(attribs):
	new_attribs = {}
	attribs = list(filter(None, attribs.split(';'))) ## removes empty strings, needed because some gff3 lines have ";;"
	for attrib in attribs:
	k, v = attrib.split('=')
	if k == 'Dbxref':
	xrefs = v.split(',')
	for xref in xrefs:
	terms = xref.split(':')
	new_attribs[terms[-2]] = terms[-1]
	else:
	new_attribs[k] = v
	return new_attribs
	def get_introns(exons_dict):
	introns_dict = defaultdict(set)
	for tx_info, exons in exons_dict.items():
	if len(exons) > 1:
	[
	chrom,
	strand,
	gene_id,
	tx
	] = tx_info
	exons = sorted(exons)
	i = 0
	while (i + 1) < len(exons):
	introns = (exons[i][1] + 1, exons[i+1][0] - 1)
	i = i + 1
	introns_dict[(chrom, strand, gene_id, tx)].add(introns)
	return introns_dict
	def tabulate_introns(introns_dict, introns_file):
	tbl = csv.writer(introns_file,
	delimiter = '\t',
	lineterminator = os.linesep)
	tbl.writerow(['#chrom',
	'intron_start',
	'intron_end',
	'strand',
	'gene_id',
	'tx_acc',
	'intron_num',
	'intron_ct'])
	for tx_info, introns in introns_dict.items():
	[
	chrom,
	strand,
	gene_id,
	tx_acc
	] = tx_info
	if strand == '+':
	introns = sorted(introns)
	elif strand == '-':
	introns = sorted(introns, reverse = True)
	num_introns = len(introns)
	for intron in introns:
	tbl.writerow([chrom,
	intron[0],
	intron[1],
	strand,
	gene_id,
	tx_acc,
	introns.index(intron) + 1,
	num_introns])
	introns_file.close()
	exons_dict = get_exons(gff3_file)
	introns_dict = get_introns(exons_dict)
	tabulate_introns(introns_dict, introns_file)

view raw gff3_to_introns.py hosted with ❤ by GitHub

Run it as follows:

./gff3_to_introns.py -i Homo_sapiens.GRCh38.94.chr.gff3 -o introns_file.tsv

The output file will have introns in the following format:

#chrom  intron_start  intron_end  strand  gene_id          tx_acc           intron_num  intron_ct
1       12228         12612       +       ENST00000456328  ENST00000456328  1           2
1       12722         13220       +       ENST00000456328  ENST00000456328  2           2
1       12058         12178       +       ENST00000450305  ENST00000450305  1           5
1       12228         12612       +       ENST00000450305  ENST00000450305  2           5

A few notes:

• This script works fine for RefSeq GFF3 and Ensembl GFF3 files. I did not test with others though.
• For Ensembl GFF3 files, the gene_id column does not have the actual gene_id; it has the transcript accession instead.
• The coordinates are 1-based just like they are for GFF3 files. BED files are 0-based.

Once you have the coordinates for introns, you should be able to use bedtools getfasta to fetch the exact sequence.

Answer 2

[Vitis]

I think Ensembl API provides an entry to tap into the sequence database and fetch the fasta sequences in a programmatic way. It requires a bit PERL programming, but isn't too bad. For the organisms that are not in Ensembl, you'll have to download the fasta and gff files and fetch the sequences locally, probably by BioPerl.

See this:

https://uswest.ensembl.org/info/docs/api/index.html

and this:

https://bioperl.org/howtos/SeqIO_HOWTO.html

[elisheva]

I tried it. But I have to say it's quite unclear for someone who is totally new to perl. Besides, I tried only to get transcript Ids of one chromosome and it takes about 3 minutes, so I guess for downloading full sequences it will take too much time.

[Vitis]

Ensembl could be slow, depending on your connection speed. But working with all exons and introns of 80 species is also a very very big endeavor. Downloading 80 genomes and extracting sequences based on GFF3 could still be an option if I'm doing this, because at least the extraction part could be very fast and efficient, see this:

https://bioperl.org/howtos/Local_Databases_HOWTO.html

[elisheva]

Thank you for your response. But as I mentioned above, gff3 doesn't include introns coordinates for some reason. Therefore, I can't see how it will be helpfull. And one more thing, assuming I do have the coordinates and the complete genomes, I guess bedtools getfasta will more efficient for this case than bioperl

[Vitis]

I think intron coordinates could be inferred from GFF3 with some carefully designed calculations. Introns are just sequences between the exons. You only need to be careful with the tricky ones involving UTRs.