HTSeq extract part of the GTF information

[crazyhottommy]

Hi all,

I am a beginner of python and just installed the HTSeq package.

my question is:

after read in the GTF file with the HTSeq.GFF_reader function, is there an efficient way
to exact certain lines from the GTF file?

now, I have three gene lists: high expression, medium expression and low expression,
with first column the gene name, and the second column counts of RNA-seq (generated by htseq.count)

one way I can do is:
gtf_file=HTseq.GFF_reader(ifile)

genelist= open(highexpression)

for line in genelist:
line_list=line.split()
for feature in gtf_file:
if feature.name==linelist[0]: # basically check the gene_name.
do something

but it is really slow because of the nested loop.

I want to plot a TSS plot as shown here http://www-huber.embl.de/users/ander...c/tss.html#tss

but I want to group Transcription start sites according to their expression levels(high, medium and low), and then plot three curves together.
-------------------------------------------------------------
one comment on the TSS plot is that when Dr.Simon Anders (the package creator, thanks!) drew
the TSS plots, he used the start point of exon1 as TSS, it is not correct.
TSS is the transcription start site,
the start point of the exon1 is the translation start site.

to plot a TSS plot, we need the GTF file containing the 5UTR information
which can be obtained here http://genomewiki.ucsc.edu/index.php..._or_gff_format

[Simon Anders]

To start with the second one: In an Ensembl GTF file, each exon appears twice, once with type "exon", once with "CDS". The start of exon 1 is a transcription start site, the start of CDS 1 is a translation start site. (For exons without untranslated regions, the coordinates in the exon line and the CDS line are the same.)

I have never checked whether the GTF files from UCSC follow this convention, too, but given the general mess that the GFF/GTF format is, it would be no surprise if it were different.

Now, to your actual Python question. The easiest might be to first read the expression-stratified gene lists into a set:

highexpr = set()
for line in open(highexpression):
   highexpr.add( line.split()[0] )

Now, you have all the highly expressed genes in memory. Hence, when you loop through the GFF file, you do:

for feature in gtf_file:
   if feature.name in highexpr:
      do_something()

As you can see, using a container to store the information from one file in memory allows you to process your two files separately and so avoid the nested loop. This, of course is a very general and basic design pattern that you will encounter very often. (The different kinds of data containers that programming languages offer and how and when to use them is probably the next thing you need to study to further your programming skills.)

[crazyhottommy]

To start with the second one: In an Ensembl GTF file, each exon appears twice, once with type "exon", once with "CDS". The start of exon 1 is a transcription start site, the start of CDS 1 is a translation start site. (For exons without untranslated regions, the coordinates in the exon line and the CDS line are the same.)

I have never checked whether the GTF files from UCSC follow this convention, too, but given the general mess that the GFF/GTF format is, it would be no surprise if it were different.

Now, to your actual Python question. The easiest might be to first read the expression-stratified gene lists into a set:

highexpr = set()
for line in open(highexpression):
   highexpr.add( line.split()[0] )

Now, you have all the highly expressed genes in memory. Hence, when you loop through the GFF file, you do:

for feature in gtf_file:
   if feature.name in highexpr:
      do_something()

As you can see, using a container to store the information from one file in memory allows you to process your two files separately and so avoid the nested loop. This, of course is a very general and basic design pattern that you will encounter very often. (The different kinds of data containers that programming languages offer and how and when to use them is probably the next thing you need to study to further your programming skills.)

Thanks so much for your quick reply!
I am glad that your clarification about the gtf file. I usually use gtf files from UCSC.

for the python question, I actually figure it out in another way after spending some time searching online, but the basic idea is the same as yours. wow, google is the best teacher



gtf_dict={ } #make a dictionary. key is the gene_name, value is the TSS infomation
for feature in gtf_file:
gtf_dict[feature.name]= feature.iv.start_d_as_pos

tsspos=set()

for line in highexpr:
linelist=line.split()
try:
tsspos.add(gtf_dict[linelist[0]])
except:
continue # in case the key is not present

And it run much faster!

Thanks again for your nice reply.

[crazyhottommy]

Hi Simon,

I just do not want to open another thread. I want to do some Interval overlapping by HTSeq.
Let's say I have two ChIP-seq data sets, they are bed file with four columns( chr start end some_value)
and I want to get the overlapping intervals. I do not want to use nested loops...

I read http://psaffrey.wordpress.com/2011/04/

How can I do it efficiently by using the GenomicInterval class with the GenomicsInterval.overlap() method?