Skip to main content
9 votes
Accepted

Calculating average coverage for .bam files (sequence data)

For a quick estimate you’re making it more complicated than necessary. The theoretical average coverage is $\frac{n \cdot \hat l}{N}$ where $n$ is the number of reads, $\hat l$ is the average read ...
Konrad Rudolph's user avatar
6 votes
Accepted

Coverage calculation: long reads (RNA-seq)

The read length is irrelevant when calculating the mean coverage statistic. It's simply the total number of bases sequenced divided by the target Xome length. In the example provided in the question, ...
gringer's user avatar
  • 14.2k
5 votes

Contamination on genome assembly

This could be organism-specific. We don't have a lot of info so far, so I would check a few more things: Run something like FRC_align. Check if there's a clear signal between regions flagged as ...
Bastian Schiffthaler's user avatar
5 votes
Accepted

How to best detect the "peaks" in RNA-seq data that are not assigned to any gene?

Basically what you're discovering is that there are unannotated expressed features, so your task isn't really finding peaks, but rather finding novel expressed transcripts. For that, you can use ...
Devon Ryan's user avatar
  • 19.6k
5 votes
Accepted

Double-counting coverage of overlapped read pairs

samtools depth is a simplified version of samtools mpileup, which handles overlapping regions by default. Since overlapping ...
Devon Ryan's user avatar
  • 19.6k
4 votes
Accepted

How to calculate overall reference coverage with MUMmer?

I believe all you need to do is to run dnadiff from MUMmer. That will run a comparison and output a number of useful metrics.
5heikki's user avatar
  • 156
4 votes
Accepted

Counts obtained by featureCounts seem much less than observed coverage

Given the high level of multimapping in this region, you'll need to use the -M --primary options if you want to keep many of the alignments. I would be very ...
Devon Ryan's user avatar
  • 19.6k
4 votes
Accepted

How to calculate AUC in coverage graph

Yes, just add up the coverage at each location; that's integration in a nutshell. Or, alternatively / equivalently, count the total number of reads within the region.
gringer's user avatar
  • 14.2k
4 votes

How can I calculate coverage at single bases using a bam file?

git clone https://github.com/lh3/htsbox cd htsbox && make ./htsbox pileup -cCf ref.fa aln.bam | less -S This output a VCF containing positions covered by ...
user172818's user avatar
  • 6,535
4 votes

How can I calculate coverage at single bases using a bam file?

Note: not yet tested, so there may be some additional fiddling with command line options needed The per-base depth can be obtained from samtools depth (-a includes zero-coverage positions): ...
gringer's user avatar
  • 14.2k
3 votes

If a gene is expressed at a level of 1/1200 compared to the average gene, how is probability 50:50 that we have a read mapped to it?

I actually disagree with that.. I guess I write this down as a discussion. The expected value for an average gene is 1200. For this gene at 1/1200 expression, you expect 1 read. However, because of ...
StupidWolf's user avatar
  • 1,688
3 votes
Accepted

Plotting coverage of annotation over collection of region

If one assumes that repeat elements of a given type (e.g., LINEs) don't overlap each other, then the following will work: Split your BED file by repeat element, such that you have a LINE.bed, SINE....
Devon Ryan's user avatar
  • 19.6k
3 votes
Accepted

How many reads do I need to cover the entire genome?

Let's see what your slides are claiming $\epsilon$ to be, using: $$ C = \frac{NL}{G} \approx \ln \left( \frac{G}{L\epsilon} \right) $$ We can rearrange to get $$ \epsilon \approx \frac{G}{L} \times ...
aesthete's user avatar
  • 146
3 votes
Accepted

Is it possible to, all in memory without writing a file, given a list of AlignedSegments and a bam header, build a bam, index and get coverage?

Assuming you are on an operating system that supports them, like Linux, you can try using named pipes in a ramdisk or tmpfs. Not the most elegant approach, but I don't know if what you are asking for ...
terdon's user avatar
  • 10.2k
2 votes
Accepted

GATK documentation for required depth to reliably call heterozygous mutation in diploid organism?

UPDATE: As an update, Sarah Walker (co-author on the poster) responded to my question on the GATK forum. She clarified with the following statement: We believe ...
Mark Ebbert's user avatar
  • 1,354
2 votes

Double-counting coverage of overlapped read pairs

One way to deal with this would be to first merge paired-end reads based upon their overlapping regions, and then map them and calculate the coverage. This way you're only counting once per unique ...
LucasBoatwright's user avatar
2 votes

Double-counting coverage of overlapped read pairs

picard CollectHsMetrics has the parameter PER_BASE_COVERAGE which take this into account. You can run it like this: ...
finswimmer's user avatar
  • 1,342
2 votes

Total read count of shotgun affected by few highly abundant sequences

That sentence in the Anders and Huber paper is a bit misleadingly worded (I don't think that was their intent). What they mean is that highly expressed and highly differentially expressed genes can ...
Devon Ryan's user avatar
  • 19.6k
2 votes

How to best detect the "peaks" in RNA-seq data that are not assigned to any gene?

Many pipelines have been set up for this such as software developed for the FANTOM Project at RIKEN. The software for this is available from our website. These were set up for analysis of CAGE ...
Tom Kelly's user avatar
  • 873
2 votes

Calculating average coverage for .bam files (sequence data)

With the current versions of samtools, you can also use: samtools coverage sample.bam The documents are here, samtools coverage – produces a histogram or table of ...
pdp10's user avatar
  • 21
2 votes

bedtools coverage - Report the depth at each position in each A feature

If you prefer to have a more concise report than a per-base one, bedtools genomecov could be a better choice. Its -bga option allows you to report consecutive bases ...
Siyuan Feng's user avatar
2 votes

Bedtools wrongly indicates zero coverage

If this is your output : ...
Paul Endymion's user avatar
2 votes
Accepted

How do you set the coverage in PacBio's Sequel II?

They just want to know the ploidy of the sample. In a diploid sample, e.g. humans, there are two haplotypes. This is mostly important for haplotype resolution (i.e. get the two full genomes in a ...
Maximilian Press's user avatar
2 votes

Low pass sequencing has been reported to detect common variants. How low can one go and get reliable data? Is 2X pass sequencing analysis possible?

This is a good idea - low pass sequencing has shown to be better than genotype arrays for things such as GWAS or QTL mapping. The primary way of processing low-pass sequencing data is to use ...
user438383's user avatar
  • 1,689
2 votes
Accepted

What should GC coverage bias plot of exome data look like?

Am I missing something? Should I, for instance, use my target regions instead of the reference genome? Yes you should do this. Why do you think it's wrong doing that? It should be the first thing to ...
Supertech's user avatar
  • 606
1 vote

How do you set the coverage in PacBio's Sequel II?

Unless you want to reconstruct all haplotypes, you can sequence to ~30X regardless of ploidy.
user172818's user avatar
  • 6,535
1 vote

how to calculate coverage of each base of mapped reads?

You are looking for bedtools genomecov, specifically -d option. See https://bedtools.readthedocs.io/en/latest/content/tools/...
Timur Shtatland's user avatar
1 vote

Metagenomic shotgun data with internal control

The problem with correlations in relative scale is know on the literature, it has been suggested to change to other metrics. So I wouldn't use it for internal control. As you can see on figure 3 of ...
llrs's user avatar
  • 4,713
1 vote

Coverage required

According to the book you mentionned Introduction to Genomics SecondEdition by Arthur MLesk p.113 ...
Pierre-Edouard Guerin's user avatar
1 vote

How to take into account alternative bwa mem mapping when computing coverage

You would want to instead use the -a option with bwa mem to output all alignments. You would be best off filtering these with a ...
Devon Ryan's user avatar
  • 19.6k

Only top scored, non community-wiki answers of a minimum length are eligible