8 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 ...
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, ...
user avatar
  • 11.9k
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 ...
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 ...
user avatar
  • 19.2k
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 ...
user avatar
  • 19.2k
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.
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 ...
user avatar
  • 19.2k
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.
user avatar
  • 11.9k
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 ...
user avatar
  • 5,814
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): ...
user avatar
  • 11.9k
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 ...
user avatar
  • 1,638
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....
user avatar
  • 19.2k
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 ...
user avatar
  • 146
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 ...
user avatar
  • 1,244
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 ...
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: ...
user avatar
  • 1,332
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 ...
user avatar
  • 19.2k
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 ...
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 ...
user avatar
  • 1,214
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 ...
user avatar
  • 873
2 votes

Bedtools wrongly indicates zero coverage

If this is your output : ...
user avatar
1 vote

Understanding exercise on file coverage with question on summary statistics

Based on the description in your question, I don't think this assumption is correct: I believe what they refer in this exercise is that the target file should be something like the whole chr sizes of ...
user avatar
  • 8,161
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.
user avatar
  • 5,814
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/...
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 ...
user avatar
  • 4,622
1 vote

Coverage required

According to the book you mentionned Introduction to Genomics SecondEdition by Arthur MLesk p.113 ...
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 ...
user avatar
  • 19.2k
1 vote

Extract mapping coverage from GTF files

So consider this, you know how many reads map to each of your genes (1). Does that tell you anything about (2) the number of reads that mapped against non-gene regions (3) the number of reads that did ...
user avatar
1 vote

Double-counting coverage of overlapped read pairs

Bases from overlapping reads of the same fragment provide extra information, as their base qualities come from different sequencing cycles on the sequencer. This information can be used by downstream ...
user avatar
  • 2,096
1 vote

How many reads do I need for hybrid assembly

I have never used dbg2olc. However other OLC assemblers managed to do a solid assembly from 20x of corrected long reads. But obviously, the more the better. In you case of allotetraploid, you need to ...
user avatar
  • 5,327

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