Dear all,
I have bedg file created from bedtools -
chr1 1 10 num_of_reads chr1 10 20 num_of_reads chr1 20 30 num_of_reads . . .
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.
So result should be:
ch start stop num_of_reads cg_content
4 answers
Here's a python solution. You need to have the pysam 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.
Thank you for help - it is installed. But when I try to run I have following error:
Traceback (most recent call last):
File "cg_content.py", line 10, in <module>
bam = pysam.AlignmentFile(bam_file.bam, 'rb')
AttributeError: 'module' object has no attribute 'AlignmentFile'
Do you know where could be a problem?
Dear Dan, thank you so much - I did not installed pyrex and Sphinx - right now it looks it work perfectly fine. May I ask - in some intervals there is no reads = 0 - I have at the end error bug:
gc_percent = float(gc_bases)/total_bases * 100
ZeroDivisionError: float division by zero
Is there any solution how to avoid in python to divide by zero? Very very thank you for your time and help!!!
Just wanted to add a link to this website, which can help debugging AlignmentFile() issues:
You could combine the bedtools map tool with awk to do this. Specifically, the example below uses map to count the number of overlapping reads for each interval (-c 10 -o count) and concatenate all of the nucleotide sequences in each overlapping BAM alignment (-c 10 -o collapse). Thus the fourth column of the output is the count of overlapping reads and the fifth column is the concatenated sequence from each read. Awk is then used to compute the GC fraction from this concatenated sequence.
➜ cat test.bed
20 0 1000000
20 1000000 2000000
20 2000000 3000000
➜ bedtools map -a test.bed \
-b <(curl -s ftp://ftp-trace.ncbi.nih.gov/1000genomes/ftp/data/HG00096/alignment/HG00096.chrom20.ILLUMINA.bwa.GBR.low_coverage.20120522.bam) \
-c 10,10 \
-o count,concat \
| awk -v OFS="\t" '{n=length($5); gc=gsub("[gcGC]", "", $5); print $1,$2,$3,$4,gc/n}'
20 0 1000000 47199 0.45724
20 1000000 2000000 47929 0.445155
20 2000000 3000000 49902 0.445785
•
written
by
Aaronquinlan
In order for the map tool to work correctly, all files must use the same chromosome sorting order. Based on what you describe, it sounds like the order of chromosomes in your BED file differes from the BAM file. Our next release will be able to detect these changes automagically.
I am really amazed by the simplicity and awesomeness of this solution.
Just one change I would like to make it to work would be
| awk -v OFS="\t" '{n=length($6); gc=gsub("[gcGC]", "", $6); print $1,$2,$3,$4,gc/n}'
As $5 is the count of the number of reads that overlapped with the bed region
This doesn't seem to be the case.
Below is the output of the first portion of the one-liner before the AWK command.
chrX 301493 301678 2 ACGGCATACGCGATAACGCTATGTGTCTGGAGTTAAGAG
chrX 305864 306261 2 TGACTCGGCCTTAGCCTCCCGAGTAGCTGGGACTACAGA
Grabbing field #5 would be the correct move here. There isn't a field 6 at all.
Building off my previous answer (How To Calculate A Region‘ Gc Content In Reference Genome？？), you could first convert the BAM reads into FASTA with samtools and awk, then pipe that result to a second awk script that converts FASTA to a BED-like result that includes columns with per-FASTA element GC content:
$ samtools view reads.bam \
| awk ' \
BEGIN { \
OFS="\t"; \
} \
{ \
print ">"$3"-"($4-1)"-"($4+$9)"\n"$10 \
}' - \
| awk ' \
BEGIN { \
FS=""; \
cg=0; \
t=0; \
} \
{ \
if ($1 != ">") { \
for (i = 1; I <= NF; i++) { \
t++; \
if ($i ~ /[CGcg]/) { \
cg++;
} \
} \
} \
else { \
if (t > 0) { \
print e"\t"cg"\t"t"\t"(cg/t); \
cg = 0; \
t = 0; \
} \
h = substr($0, 2); \
split(h, a, "-"); \
e = a[1]"\t"a[2]"\t"a[3]; \
} \
} \
END { \
print e"\t"cg"\t"t"\t"(cg/t); \
}' - \
> answer.bed
You'll need to modify this (untested) pipeline to get the exact columns you want, but I think this should get you about 95% of the way there.
•
written
by
Alex Reynolds
It looks like you modified your question since I answered it. However, you can take the result from this pipeline and do a bedmap --mean operation against a BED file of pre-made windows (say, 10 bases apart), to quickly and efficiently calculate the average GC content over each window.
Log in to answer this question.
-
Concatenate information columns from multiple identical bed files.
written by nicole.zeltser •
-
Merge overlapping and adjacent features in the BED file having the same label in the name (4-th) co…
written by Denis
-
Matching genomic intervals in two bed files
written by arsala521 •
-
Bedtools intersect with index
written by landscape95
-
Counting summary info (count of reads and average GC) for segment with specific length
written by Korsocius
-
Get gVCF file from genome coordinates
written by Paul
-
Normalization reads between samples for aneuploidy detection
written by Paul
-
Identify Overlapping And Non Overlapping Regions For Paired-End Data
written by ancient_learner
-
How To Get Annotation For Bed File From Another Bed File
written by k.nirmalraman
-
Cluster Intervals By Gene Id In A Bed File
written by Rad
Should the cg_content column be expressed as a percentage of bases which are G or C?
better would be just a number - cg content - I now how to calculate separately from bam -
but I need to calculate according BED file.
It is important consider that such measurements will be biased by coverage