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 ...
- 4,895
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, ...
- 15.1k
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 ...
5
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):
...
- 15.1k
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 ...
- 19.8k
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 ...
- 19.8k
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.
- 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 ...
- 19.8k
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.
- 15.1k
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 ...
- 6,605
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 ...
- 1,698
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....
- 19.8k
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 ...
- 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 ...
- 10.6k
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 ...
- 1,414
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 ...
- 316
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:
...
- 1,352
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 ...
- 19.8k
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 ...
- 873
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 ...
- 4,252
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 ...
- 1,789
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 ...
- 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 ...
- 21
2
votes
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 ...
- 636
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/...
- 841
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 ...
- 4,721
1
vote
Coverage required
According to the book you mentionned Introduction to Genomics SecondEdition by Arthur MLesk p.113
...
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 ...
- 19.8k
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 ...
Only top scored, non community-wiki answers of a minimum length are eligible