Question: Nucleotide coverage if >5 reads
0
gravatar for Vca80553
4 weeks ago by
Vca805530
Stockholm
Vca805530 wrote:

Hi Everyone,

I have generated BAM files after mapping my Illumina reads to a reference genome. Now I want to know how much of the reference genome is covered (aligned/mapped) by reads (e.g. 98% of the reference genome is covered by reads). For that I used Bbmap/pileup.sh. It works fast and it is accurate.

However, the 98% I get is the % of nucleotides covered by at least 1 read. It it possible to modify this, and get % of nucleotides that were covered at least with 5 reads?

Thanx,

Sara

pileup sequence alignment • 257 views
ADD COMMENTlink modified 4 weeks ago by Pierre Lindenbaum100k • written 4 weeks ago by Vca805530
1

Take a look at the pileup.sh command options I have compiled in this thread (search for name pileup) to see if one of those options helps.

ADD REPLYlink written 4 weeks ago by genomax37k
1

Hello Vca80553!

It appears that your post has been cross-posted to another site: http://seqanswers.com/forums/showthread.php?t=78704

This is typically not recommended as it runs the risk of annoying people in both communities.

ADD REPLYlink written 4 weeks ago by Pierre Lindenbaum100k

Yes, sorry. I didn't know this. My fault. Will not happen again. Thanx.

ADD REPLYlink written 29 days ago by Vca805530

For pileup.sh

Output format (tab-delimited):
ID, Avg_fold, Length, Ref_GC, Covered_percent, Covered_bases, Plus_reads, Minus_reads, Read_GC, Median_fold, Std_Dev

Plus_reads:        Number of reads mapped to plus strand
Minus_reads:       Number of reads mapped to minus strand

Perhaps you can add those two numbers and filter on that to get the info you need.

ADD REPLYlink modified 4 weeks ago • written 4 weeks ago by genomax37k

Tagging: Brian Bushnell

ADD REPLYlink written 4 weeks ago by genomax37k

Hi,

Thanx for answering.

I am not sure I understand your answer. (I am pretty new with this, sorry)

The output format gives me something like this for each sample:

ID, Avg_fold, Length, Ref_GC, Covered_percent, Covered_bases, Plus_reads, Minus_reads, Read_GC, Median_fold, Std_Dev
S1588, 24909.1657, 7904, 100.0000, 98.7854, 7808, 1105713, 1105713, 6112, 20615.50

The Plus reads and Minus reads give me the total number of mapped reads to the whole reference genome. How can I filter from there?

ADD REPLYlink modified 4 weeks ago by genomax37k • written 4 weeks ago by Vca805530
1

Sorry. You are right. For some reason I was thinking that the output of pileup.sh is per base.

ADD REPLYlink written 4 weeks ago by genomax37k
1
gravatar for harold.smith.tarheel
4 weeks ago by
United States
harold.smith.tarheel3.9k wrote:

You can use BEDtools genomecov command:

bedtools genomecov -ibam FILE.BAM -d | awk '$4 > 4' | wc -l

Divide the output by the size of the genome to get the percentage.

ADD COMMENTlink written 4 weeks ago by harold.smith.tarheel3.9k
Please log in to add an answer.

Help
Access

Use of this site constitutes acceptance of our User Agreement and Privacy Policy.
Powered by Biostar version 2.3.0
Traffic: 1474 users visited in the last hour