Selecting intergenic regions

Question

[vinayjrao]

Hi,

I have gtf files for a list of organisms. I want to extract only the intergenic sequences for each of their corresponding fasta files. Currently, I have been able to use grep -o '^..........' ecoli.fa > 10.txt to grep out the first 10 characters if the first gene starts at pos 11.

I now need to find a way to automate the process for the remaining regions, as manual inspection is impossible. How can I do the same for a huge file, for example -

Gene      Start     End
Gene1    11        157
Gene2    2878    8548
Gene3    11254  12124

The process might require file handling to subtract End of a row from the Start of the next row, and that value could be coupled with grep to give a command like grep -o '^{2877.}' ecoli.gtf | tail -f 2720 > 2720.txt

P.S. I could make the gtf file containing the intergenic lengths, but I do need help with the file handling.

Thanks in advance

[finswimmer]

Hello,

you could first create a bed file containing the coordinates for each intergenic regions with the language of your choice. Than you can use this bed file together with bedtools getfasta to get the sequences.

fin swimmer

[vinayjrao]

Thanks for the advice. It would certainly help, but to do the same, I first need to employ file handling to select the Start and End coordinates of the annotated genes. Would you happen to have a code to create the bed file?

Thanks.

[finswimmer]

How does your file with gene coordinates looks like? In your first post you've said you have a gtf file. But the example you have posted doesn't looks like one.

fin swimmer

[vinayjrao]

Sorry for the confusion on my part. I only showed the relevant columns, but it is a proper gtf file.

GeneID    Chr   Start   End   Strand    GeneInfo

[finswimmer]

@ Pierre Lindenbaum , Devon Ryan , genomax , WouterDeCoster

What happens to that thread? Who closes it and why?

fin swimmer

[GenoMax]

Only mods close threads when they are off-topic/spam etc. Users should just accept answers that are correct. That provides proper closure for the thread.

[M201570116]

Hello,

I have a GFF3 file. How can I get the intergenic regions and its positive and negative chain information?

Thank.

[GenoMax]

This is an unrelated question and should have been posted as a new thread. Also provide clear objective of what you are trying to achieve.

Answer 1

[finswimmer]

Hello again,

this doesn't look like a proper gtf file. Nevertheless here is a python script that should create a bed file for the intergenic regions. It asumes the columns you gave, that your file is sorted by coordinates and that there is header line.

Run it like this:

$ python intergenic.py genes.gtf > intergenic.bed

[vinayjrao]

Actually, I modified the gtf file I had earlier. I downloaded a gff file from NCBI, and converted it to gtf using gffread ecoli.gff -T -o ecoli.gtf. This gtf was modified to keep the columns mentioned in my previous comment.

Thanks a lot for all your help.

[vinayjrao]

I tried running the script, but line 22 sep="\t" throws an error saying invalid syntax. I haven't learnt python :(

[finswimmer]

Hello,

than you are using python2. Add this line at the beginning of the python script and run again:

from __future__ import print_function

fin swimmer

[vinayjrao]

The script is perfectly running. Thanks a lot for your time. Could you please add it as an answer so that I can accept it.

Thanks again :)

[finswimmer]

You're welcome.

I moved my post to an answer.

fin swimmer