<p>Dear all,</p>

<p>I have bedg file created from bedtools -</p>

<pre>
chr1 1  10 num_of_reads

chr1 10 20 num_of_reads

chr1 20 30 num_of_reads

.

.

.</pre>

<p>3th column is information about overlapping reads at this region. I would like to add to 4th column information about GC content in reads from my BAM file. Do you have any idea how to do that? Thank you for any help.</p>

<p>So result should be:</p>

<pre>
ch start stop num_of_reads cg_content</pre>

Here's a python solution. You need to have the [pysam][1] library installed (`pip install pysam`):

```
import pysam
from sys import argv

bam_file = argv[1]
bed_file = argv[2]

bam = pysam.AlignmentFile(bam_file, 'rb')
with open(bed_file) as bf:
    for line in bf:
        line_parts = line.strip().split()
        chr = line_parts[0]
        start = int(line_parts[1])
        end = int(line_parts[2])
        read_data = bam.fetch(chr, start, end)
        total_bases = 0
        gc_bases = 0
        for read in read_data:
            seq = read.query_sequence
            total_bases += len(seq)
            gc_bases += len([x for x in seq if x == 'C' or x == 'G'])
        if total_bases == 0:
            gc_percent = 'No Reads'
        else:
            gc_percent = '{0:.2f}%'.format(float(gc_bases)/total_bases * 100)
        print '{0}\t{1}'.format(line.strip(), gc_percent)
```

I didn't understand your comment to my question. If you want to print the total number of GC bases in the reads of a coordinate set, just change the last line to:

    print '{0}\t{1}'.format(line.strip(), gc_bases)

You run the program like this:

    python [script_name] [bam_file] [bed_file]

It prints the modified BED data to STDOUT.

 [1]: https://github.com/pysam-developers/pysam