This is a test version of Biostars. For the public version, visit https://www.biostars.org.
Gc Content From Bam

Dear all, I have one problem, I would like to count the content of GC nucleotids from bam file for each read. I think I have to count it from $10, but I do not know how? Could you help me please? Thank you

gc bam awk perl
• 8,385 views • link
updated by Ram Ram • written by filipzembol filipzembol

2 answers 

awk '{ n=length($10); print $10, gsub(/[GCCgcs]/,"",$10)/n;}' your.sam

• 179 views • link
updated by Ram Ram • written by Pierre Lindenbaum Pierre Lindenbaum

I am not sure if it is good, because every time I receive the error awk: (FILENAME=filip.sam FNR=1) fatal: division by zero attempted, than there is the problem with "n" or I do not know.

• 0 views • link
updated by Ram Ram • written by filipzembol filipzembol

I tested my awk script without the sam header.

• 0 views • link
updated by Ram Ram • written by Pierre Lindenbaum Pierre Lindenbaum

Oh ok, because I use to convert samtools view -h (with header). Without header is working, than thank you so much for your help.

• 0 views • link
updated by Ram Ram • written by filipzembol filipzembol 
samtools view [bam] | perl -lane '$N =()= $F[9] =~ /GC/gi; print $N;'

The odd looking =()= operator instantly converts the match from array to scalar context.

Perl one liner 'a' splits by tabs, 'l' prints a new line, 'n' operates on all lines and 'e' indicates expression.

• 301 views • link
updated by Ram Ram • written by hardingnj hardingnj

An alternate Perl part:

perl -p -lae '$_ =()= $F[9] =~ /GC/gi'

• 0 views • link
written by Kenosis Kenosis

I am not sure, but I try the perl code and I recieve different number from the awk code and awk code has a right numbers of GC count (I did not use /n)..

• 0 views • link
updated by Ram Ram • written by filipzembol filipzembol

Log in to answer this question.

More posts like this
  • Read File Format .Bam
    written by hudarembang • hudarembang
    <p>My research is related to the manufacture of snp database. constraints that occurs is how to read the file. <a href='http://samtools.sourceforge.net/SAM1.pdf'>bam</a>. in this case, I …
  • Separate Bam File
    written by filipzembol filipzembol
    <p>Dear all, I have one problem I need to count number of rows and GC content in <a href='http://samtools.sourceforge.net/SAM1.pdf'>bam</a> file for individual segments, but I …
  • Pcr Free Variant Calls
    written by tommy.carstensen • tommy.carstensen
    <p>If I have a PCR free high coverage <a href='http://samtools.sourceforge.net/SAM1.pdf'>bam</a> and a low coverage <a href='http://samtools.sourceforge.net/SAM1.pdf'>bam</a> and I downsample the PCR free high coverage <a …
  • Filtration Of Reads With Length Lower Than 30 From Bam
    written by filipzembol filipzembol
    <p>Dear all, I have one question how could I filtrate the reads from <a href='http://samtools.sourceforge.net/SAM1.pdf'>bam</a> file, which have length of read lower than 30 bp. …
  • Compute Binned Gc-Normalized Read Counts From Bam File
    written by Christian Christian
    <p>Binned read counts can be readily obtained from <a href='http://samtools.sourceforge.net/SAM1.pdf'>BAM</a> files using for example 'samtools depth' or <a href='https://code.google.com/p/bedtools/'>bedtools</a>. However, I would be interested in …
  • Merging Bam Files
    written by rickyflintoff • rickyflintoff
    <p>I have two <a href='http://samtools.sourceforge.net/SAM1.pdf'>BAM</a> files that have a different sample in read group (@RG) section. I would like the sample to be the same. …
  • How Do I Pull Out One Read (Given Its Id) From A Bam?
    written by zam.iqbal.genome zam.iqbal.genome
    <p>Hi there</p> <p>Does anyone know how I can pull just one read from a <a href='http://samtools.sourceforge.net/SAM1.pdf'>BAM</a> file, when I know the read-id (and sequence) of …
  • Extracting Read Count For Each Gene/Exon From Rna-Seq Bam Files
    written by Sahel Sahel
    <p>Hello all,</p> <p>I have some RNA-seq <a href='http://samtools.sourceforge.net/SAM1.pdf'>bam</a> files which I would like to submit to edgeR for differential expression analysis. My problem is how …
  • Soap Result Converted To Bam Not Working In Samtool
    written by Sarwar • Sarwar
    <p>I converted the output of soap2 which is map file to <a href='http://samtools.sourceforge.net/SAM1.pdf'>sam</a> file using map2sam.pl and then converted this <a href='http://samtools.sourceforge.net/SAM1.pdf'>sam</a> to <a href='http://samtools.sourceforge.net/SAM1.pdf'>bam</a> …
  • How To Read Bam.Ebd, Bam.Ebi Format Dataset
    written by Lerong • Lerong
    <p>I have the sequencing data divided into <a href='http://samtools.sourceforge.net/SAM1.pdf'>bam</a>.ebd and <a href='http://samtools.sourceforge.net/SAM1.pdf'>bam</a>.ebi format for each individual like:</p> <pre><code>HG00096.mapped.ILLUMINA.bwa.GBR.low_coverage.20101123.bam.ebd HG00096.mapped.ILLUMINA.bwa.GBR.low_coverage.20101123.bam.ebi HG00102.mapped.SOLID.bfast.GBR.low_coverage.20101123.bam.ebd HG00102.mapped.SOLID.bfast.GBR.low_coverage.20101123.bam.ebi </code></pre> <p>Anyone know how …