are there some software to compute depth？

Question

[Tao]

Are there some software to compute depth if a site was covered by read1 read2 or both note as one?

Answer 1

[seidel]

There are many tools that can do this. A few that come to mind are (1) bedtools coverage, if you're comfortable using linux, this is a straightforward way to proceed, and (2) R and the GenomicRanges library. If you know R, you can import bam files, calculate coverage with a single function call, and have a direct read out of base depth at any position, or easily quantify reads across features. (3) samtools depth. See your question on samtools. There are several other ways of doing this that perhaps others can mention. You didn't specify any particular preferred skill set or environment.

Answer 2

[Alex Reynolds]

If you want more statistics than just counts of reads that overlap gene annotations (or other sites, as determined by you), perhaps the following answer will help.

Starting with a generic BAM file, it can be converted to BED via bam2bed.

The 13th column (QUAL) can be converted from its printable characters back to an average Phred quality score, and then swapped with the fifth (score) column, to following the UCSC specification for BED files.

This numerical value can subsequently be used for mapping, applying score statistical functions with bedmap (shown further below):

$ bam2bed < reads.bam | awk -v FS="\t" -v OFS="\t" 'BEGIN{ for (n=0; n<256; n++) ord[sprintf("%c",n)]=n }{ old=$5; l=split($13,q,""); s=0; for(n=1; n<=l; n++) { s+=ord[q[n]]-33; } $5=s/l; $13=old; print $0; }' > reads.bed

You can start with genes as described in another answer in this thread, which are in a file called genes.bed.

Now you can map your BED file containing genes-of-interest to the reads from the BAM file:

$ bedmap --skip-unmapped --delim '\t' --echo --count --bases-uniq --echo-ref-size --bases-uniq-f --mean genes.bed reads.bed > answer.bed

The file answer.bed will contain summary statistics for reads mapped to genes and look something like this, unique to the makeup of your genes and reads:

$ head answer.bed
chr1    199064  200374  ... 404 1310    1310    1.000000    25
chr1    6853080 6870040 ... 204 9044    16960   0.533255    20
chr1    8004008 8027135 ... 507 22308   23127   0.964587    20
chr1    8061658 8098833 ... 366 37175   37175   1.000000    25
chr1    8174215 8211962 ... 369 34253   37747   0.907436    35
chr1    8873384 8890703 ... 434 17319   17319   1.000000    50
chr1    8874708 8890922 ... 420 16214   16214   1.000000    45
chr1    15833039    15850691    ... 271 17652   17652   1.000000    50
chr1    16643023    16645543    ... 266 2520    2520    1.000000    45
chr1    16665212    16667822    ... 281 2610    2610    1.000000    30

The first three columns are the intervals defined by genes. These and metadata columns come from using --echo.

Columns in between will be metadata about the gene (symbol, strand, etc.) and depend on the source of gene annotations.

The last five columns come from the remaining options (--count, --bases-uniq, etc.) and are as follows:

• The number of reads in the BAM file that overlap an interval in the BED file by one or more bases
• The number of bases in the interval in the BED file that have coverage by reads from the BAM file
• The length of the interval in the BED file
• The fraction of bases in the interval in the BED file that have coverage by reads from the BAM file
• The (arithmetic) mean of Phred quality scores of reads over the region

The order of the specified bedmap operands --count, --bases-uniq, --echo-ref-size, --bases-uniq-f, --mean write the aforementioned five column values in that same order.

Other operations are available and can be added to the bedmap command, if useful. Additionally, if one base of read coverage is too lenient, additional operands are available to further constrain the overlap criterion.

Run bedmap --help and take a look at the "Overlap Options" and other sections, or take a look at the online documentation for a more complete walkthrough.