This is a test version of Biostars. For the public version, visit https://www.biostars.org.
How to calculate percent coverage of each read for nanopore

Hello, I am assigned to investigate the read the cover the reference more than 80%. I know that samtools can be used to calculate coverage and depth for the interested regions but what about each read? For example, I have read 1-10 that cover region A, I want to know their percent coverage for all of each read. Expect output: read 1 70% cover, leftmost position 1 (like sam format column) Is there any tools or script that can be used for this? Thanks in advance.

alignment nanopore sequencing
• 5,122 views • link
written by sonsunjirachote • sonsunjirachote

It seems like several concepts may be confused. "Reads" don't have a coverage, they're just one sequenced molecule. Instead, mapped reads can be characterized by their "length" and "mismatch rate". So are you trying to get the coverage of the sequence, get reads where >80% match the reference, get reads that cover 80% of the reference (read length), or something else entirely?

• 0 views • link
written by Jautis Jautis

It seems like several concepts may be confused. "Reads" don't have a coverage, they're just one sequenced molecule. Instead, mapped reads can be characterized by their "length" and "mismatch rate". So are you trying to get the coverage of the sequence, get reads where >80% match the reference, get reads that cover 80% of the reference (read length), or something else entirely?

• 0 views • link
written by Jautis Jautis

Yes, I want the reads that cover 80% of my interested region.

• 0 views • link
written by sonsunjirachote • sonsunjirachote

1 answer

Have a look at something like sambamba depth to obtain coverage stats (for a bam mapped to a reference genome)

From a bash script I use frequently:

    input=$i
    sec_input=${input%%.bam}

    window=100000
    overlap=50000
    covMax=999999999
    threads=1

sambamba depth window -t $threads --max-coverage=$covMax --window-size=$window --overlap $overlap -c 0.00001 ${sec_input}.bam > ${sec_input}_cov_window.txt &

• 0 views • link
written by colindaven colindaven

Reads aligned to your reference of interest in BAM format.

• 0 views • link
written by GenoMax GenoMax

I have tried with my sorted bam file but the message show that sambamba-depth: All files must be coordinate-sorted. Processing reference #3317 (tufA) Processing reference #3318 (fusA) Then it generated empty file, only the header column show in file.

• 0 views • link
written by sonsunjirachote • sonsunjirachote

It looks like you aligned to genes individually ? The standard approach is to align the reads to the complete genome.

• 0 views • link
written by colindaven colindaven

Oh, I see. the region that I interested is the gene. Anyways thanks.

• 0 views • link
written by sonsunjirachote • sonsunjirachote

You could use samtools idxstats your.bam to get number of reads that are aligned to each gene in your reference.

• 0 views • link
written by GenoMax GenoMax

I think you need to adjust your question. What did you do ? Which organism ? What are your research questions. Normally, align to the genome, then view the coordinate sorted bam in IGV.

For gene positions, get a GFF3 or GTF file and view this together with your reads. Good luck

• 0 views • link
written by colindaven colindaven

Thanks for advice! I have question about mapping to reference genome (FASTA) with gene annotation file (GTF or GFF) what is difference of output that they generate in SAM file when compare to mapping with mapping to reference genome (FASTA) only. Thanks in advance!

• 0 views • link
written by sonsunjirachote • sonsunjirachote

Log in to answer this question.

More posts like this