Script to report GC content of the island returned by SICER

Question

[Soumitra Pal]

Hi,

I need to find the GC content of the islands returned by the SICER tool. Does anybody know about such a script? I want to reuse if it exists rather than reinvent.

Any code will be fine, however python will be preferred.

Another way to ask my question. I need to solve the following problem.

Input: 1) Genome, 2) Intervals on genome Output: Percent of GC in those intervals.

Thank you. Soumitra

[Soumitra Pal]

Thank you @EagleEye.

Answer 1

[WouterDeCoster]

So you have a bed file containing intervals of interest and a fasta of the genome? Sounds rather trivial.

You can use bedtools getfasta to slice out the fragments of the genome you need based on the bed file. Then you have fragments for which you need to calculate the GC content.

I mashed a bit of code together from a related script to produce the following:

Takes in a fasta file, returns table of identifier, GC content and length (bonus).

Optionally you could write a script which takes the entire genome, slices out the parts you want and calculates on that.

[Soumitra Pal]

Thank you WouterDeCoster for your quick reply. I will use bedtools and your code, as you suggested. I hope the bedtools will handle 'N' in the genome file and also multiple chromosomes. Let me explore that. Thanks again.

[WouterDeCoster]

I hope the bedtools will handle 'N' in the genome file and also multiple chromosomes

That should work without issues.

[Soumitra Pal]

Yes, it worked without any issues. Thank you.

I also found out that 'nuc' instead of 'getfasta' from bedtools can also be used. It directly gives GC contents and other information. However, I will still use your code as it will help me do some more analysis.

Answer 2

[EagleEye]

Check out this post

A: GC percent detection for several number of genes

[Soumitra Pal]

Thank you too EagleEye.