determining read depth for every position in a sequence
I am trying to create a table for a specific sequence, about 1900bp long and have for each nucleotide the read depth achieved by sequencing. I saw you can check read depth by:
samtools view file.bam chr 1:100-500 | wc
but that would require 1900 commands (not good...) I am sure that there is a smarter way to do this - can anyone offer any help? Thanks!
2 answers
you should use genomeCoverageBed from bedtools
bedtools genomecov computes histograms (default), per-base reports (-d) and BEDGRAPH (-bg) summaries of feature coverage (e.g., aligned sequences) for a given genome.
that is I think exactly what you need
From samtools:
samtools depth input.bam >input.coverage
Refer to this post:
Log in to answer this question.
More posts like this
-
Read depth per base DNA seq
written by Sabeen •
-
Count number of unique alignments to each chromosome from sam file
written by sg197 •
-
How split regions in half to get at most 1000 bp length?
written by star
-
Combine multiple samtools view
written by rah •
-
determining read-depth of given size of genomic window via samtools
written by ella •
-
Using Samtools depth to get number of forward and reverse reads per base from a .bam file
written by jspainhour3 •
-
Samtools mpileup, coverage and view provide different read count
written by george •
-
How to display a VCF/BCF file or stream as a paginated table in a python web framework (e.g. Django…
written by William
-
Multithreading Python3 + Samtools
written by glihm
-
Discrepancy In Samtools Mpileup/Depth And Bedtools Genomecoveragebed Counts
written by Surya Saha
