Hi friends
I want to calculate average of read count / coverage for all genes, I mean creating a list which shows number of aligned reads for each location.
For better understanding, making a proportion for all genes location and an average list of read count for the locations like :
Location = Upsrtream 0-200 200-400 400-600 ... Downstream
Average read counts= 20reads 15reads 32reads ....
How can I make this list from bam or wig file?
3 answers
I recently found myself looking for the exact info. I too decided to go with bedtools. This is my final one liner assuming a pre-sorted bam file.
bedtools coverage -sorted -d -g human_g1k_v37.genome -a genes.bed -b my_sorted.bam \
| sort -k4,4 \
| bedtools groupby -g 4 -c 8 -o mean \
| sort -k1,1
Genome file:
$ cat human_g1k_v37.genome
1 249250621
2 243199373
3 198022430
4 191154276
5 180915260
6 171115067
7 159138663
8 146364022
9 141213431
10 135534747
11 135006516
12 133851895
13 115169878
14 107349540
15 102531392
16 90354753
17 81195210
18 78077248
19 59128983
20 63025520
21 48129895
22 51304566
X 155270560
Y 59373566
MT 16569
Interval file:
$ head genes.bed
#chrom start end name score strand
1 2337184 2337293 PEX10 . -
1 2337902 2338078 PEX10 . -
1 2338138 2338414 PEX10 . -
1 2339870 2340317 PEX10 . -
1 2341789 2341910 PEX10 . -
1 2342177 2342327 PEX10 . -
Ouput:
ABCA12 193.5
ABCA3 195.35
ABCB1 12.67
ABCC8 176.47
ACADL 226.97
ACOX1 178.44
ADA 149.33
ADAMTS13 187.34
ADGRV1 196.86
AFF2 163.39
Hi again,\ I use different way to find gene base coverages
bedtools bamtobed -i input_bam | bedtools coverage -mean -b - -a input_bed > output.file\ it gave me mean depth for each region that included in bed file.
Then I calculated length for each region in the bedtools output. After that I grouped file according to gene and obtained each genes total length and I divided depth by related genes length then I summed each value according to related genes that I obtain. While it sounds complicated, it's actually easy.
I hope this help guys who want to calculate gene-based coverage.
Thanks.
How about using bedtools?
First, sort the bamfile:
samtools sort bamfile Output_sorted
Then run bedtools:
bedtools genomecov -ibam Output_sorted.bam
This will print something like (taken from the manual):
chr1 554304 554309 5
chr1 554309 554313 6
chr1 554313 554314 1
chr1 554315 554316 6
chr1 554316 554317 5
chr1 554317 554318 1
chr1 554318 554319 2
chr1 554319 554321 6
chr1 554321 554323 1
chr1 554323 554334 7
So from 554304 to 554309 there are 5 reads aligning etc.
If you have a second bed file of gene coordinates you can then report the intersection between both files with bedtools intersect
I wrote https://github.com/lindenb/jvarkit/wiki/BamStats05 ; It needs a BED and a BAM; The output looks like
#chrom start end gene length mincov maxcov mean nocoveragepb percentcovered
Try NGS: RNA Analysis > htseq-count. The SAMtools and BEDtool items mentioned in this thread area also wrapped for Galaxy.
Log in to answer this question.
-
Count Matrix normalisation for downstream analysis and for creating heatmap of targeted genes
written by Vikram •
-
Multiple replicates with one SoapDenovo-Trans run?
written by nina.maryn •
-
lower exon percentage
written by anubratadas •
-
What causes strange oscillations in the count distribution of Salmon and STAR aligned gene tables?
written by Gabriel
-
Segmentation fault (core dumped) in SOAP
written by choi.yisoo.hi •
-
Deletions in a consensus sequence.
written by kspata
-
how to find the locations of protein products of a certain list of genes?
written by Yingzi Zhang
-
Tophat/cuffdiff versus CLC for mapping RNA-seq reads
written by muppetleague •
-
Which Of The Genes Are Enriched With Repeat Elements
written by roll
-
Rna-Seq Raw Counts
written by Gregor Rot
you can get the read count for each coordinate from bam using samtools view option.
samtools view bamfile chr10:3100260-3110000 | wc -l
this u can loop in command line for all the genes u have
or just
samtools view -c bamfile chr10:3100260-3110000but your solution just count the reads, it doesn't consider the deletions, the clipping, etc...