How to calculate GC content of reads that mapped to a specific gene?

Question

[bioinfo]

Hello,

I aligned my data with STAR and I have the BAM files. I would like to calculate the GC content of reads that mapped to a specific gene. I have found this thread Gc Content From Bam but this does not seem to specify to do this just for a specific gene.

Thank you

Answer 1

[bk11]

First create bam file for your gene-

samtools view -b input.bam "ChrX:10000-40000" > your_gene.bam

And then

samtools view  your_gene.bam | cut -f 10 | fold -w 1 | awk '($1=="G" || $1=="C") {N++;} END {print (N/(1.0*NR)*100);}'

OR you can do in using BBMap Suite

reformat.sh -Xmx2g in=your_gene.bam out=your_gene.fa | stats.sh -Xmx2g in=your_geneIN.fa

[bioinfo]

Thank you so much. I used samtools and it worked perfectly.