6
$\begingroup$

People talk about deep sequencing. Is there any way to calculate how deep the sequencing is ? What should be the optimum depth to get reliable data ?

I am doing whole genome sequencing of a virus genome which is 10 kb long. I got 80000 reads from Illumina sequencing. Is there any way to tell about my sequencing depth ?

$\endgroup$
8
$\begingroup$

I found a post useful for this topic.

It explains the difference of coverage and depth. It also has a useful explanation on how to calculate coverage and depth.

enter image description here

Here is a copy of what the link says just encase the post is removed:

Depth of coverage

How strong is a genome "covered" by sequenced fragments (short reads)?

Per-base coverage is the average number of times a base of a genome is sequenced. The coverage depth of a genome is calculated as the number of bases of all short reads that match a genome divided by the length of this genome. It is often expressed as 1X, 2X, 3X,... (1, 2, or, 3 times coverage).

Breadth of coverage

How much of a genome is "covered" by short reads? Are there regions that are not covered, even not by a single read?

Breadth of coverage is the percentage of bases of a reference genome that are covered with a certain depth. For example: 90% of a genome is covered at 1X depth; and still 70% is covered at 5X depth.

SAMtools: get breadth of coverage

Get breadth of coverage of a reference genome.

How many bases of my reference genome are covered by my short reads sample at a level of 5X or higher?

 # set shell variables to your file-names

 REF_GENOME_FILE=my_ref_genome.fna
 WGS_SAMPLE=my_short_read_sample.fastq 
 MIN_COVERAGE_DEPTH=5
  1. Map short reads against reference genome (bowtie2)

    # create bowtie2 index database (database name: refgenome)
    
    bowtie2-build ${REF_GENOME_FILE} refgenome
    ls
    refgenome.1.bt2
    refgenome.2.bt2
    refgenome.3.bt2
    refgenome.4.bt2
    refgenome.rev.1.bt2
    refgenome.rev.2.bt2
    
    
    # map reads and sort bam file
    bowtie2 -x refgenome --no-unal -U ${WGS_SAMPLE} -S - -p 12 | \
    samtools view -bS - | \
    samtools sort -m 5G - mapping_result_sorted.bam
    
    -U ${WGS_SAMPLE} consider all reads as unpaired reads
    -S -             bowtie output in SAM format, write to standard out
    -x refgenome     use reference genome bowtie2-database
    --no-unal        suppress SAM records for reads that failed to align
    -p 12            use 12 processors
    
  2. Get breadth of coverage (SAMtools)

    # create samtools index
    samtools index mapping_result_sorted.bam
    
    # get total number of bases covered at MIN_COVERAGE_DEPTH or higher
    samtools mpileup mapping_result_sorted.bam | awk -v 
    X="${MIN_COVERAGE_DEPTH}" '$4>=X' | wc -l
    32876
    
    # get length of reference genome
    bowtie2-inspect -s refgenome | awk '{ FS = "\t" } ; BEGIN{L=0}; {L=L+$3}; 
    END{print L}'
    45678
    

Result: breadth of reference genome coverage total number of covered bases: 32876 (with >= 5X coverage depth) Depth of coverage (average per-base coverage): 0.719 X (32876 ÷ 45678) (total number of covered bases divided by reference genome length) percent: 71.9% (0.719 × 100) (coverage breadth in percent)

$\endgroup$
  • 1
    $\begingroup$ Nice picture! Could you include the relevant parts of the linked page in the answer to make the answer self sustained (and avoid having to go one link further). $\endgroup$ – llrs Oct 31 '17 at 7:39
  • $\begingroup$ Thanks for the answer Samantha. This exactly answers my question. $\endgroup$ – L R Joshi Nov 1 '17 at 3:30
7
$\begingroup$

Sequencing depth is typically calculated as the number of total bases sequenced divided by the number of bases in the target genome. An Illumina sequencing run with 2x125 bp reads and 500 million read pairs sequenced would be a sequencing depth of about 40X (assuming my calculations are correct for a 3 billion base genome).

The sequencing depth depends on how long each read is, as well as whether it is a paired-end sequencing run. 80,000 read pairs with each having 125 bp would be a coverage of 2000X for a 10 kbp viral genome (assuming the sample were entirely composed of that virus).

The optimum depth changes depending on the application. A genotyping-by-sequencing analysis, combined with imputation, could work well with 96 samples at a 1X coverage (or possibly less). On the other end of the scale, a cell-free DNA analysis to identify structural variation in metastasing cancer may need a few thousand times coverage of the host DNA.

$\endgroup$
7
$\begingroup$

There are several questions in your post I'll try to answer each one:

Is there any way to calculate how deep the sequencing is ?

See gringer's answer. TLDR: The depth of the sequencing is how many times each position has been sequenced.

What should be the optimum depth to get reliable data ?

The optimum depth depends on what you want to do with that data. Usually, 30X or 50X is done, but in targeted sequencing you can sequence a lot more reaching 1000X (which would be really deep).

Is there any way to tell about my sequencing depth ?

If you know the length of the region sequenced and how many bases you have in the reads, you can calculate the general sequencing depth.

$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.