Visualizing Coverage Of Ngs Reads Across Reference Genome
3
2
Entering edit mode
9.2 years ago
nicolazilio ▴ 20

I am interested in finding out at how well my genome is covered (depth of coverage) by the reads I got out on an Illumina run to determine if there are loci that are under-represented, or missing (i.e. large deletion) or over-represented. Is there any script that can do that? I already aligned the reads using BWA.

Thank you.

coverage sequencing • 6.1k views
ADD COMMENT
4
Entering edit mode
9.2 years ago
KCC ★ 4.0k

You probably have a file in SAM format. Convert to BAM (samtools):

samtools view -bS -F 4 data.sam > data.bam

Convert from BAM to BED (bedtools command):

bamToBed -i data.bam > data.bed

BED to Sorted BED (unix system command):

sort -k1,1 data.bed > data.sorted

Sorted BED to BedGraph (bedtools command):

genomeCoverageBed -bga -i data.sorted -g chrom.sizes > data.bedGraph

Convert to bigwig for fast viewing (UCSC utilities):

bedGraphToBigWig data.bedGraph chrom.sizes data.bw

Finally, use a program like IGV to view it. On a Mac, it's as easy as selecting the right genome, then drag and drop the bigwig file into the window. You can now explore your underrepresented features.

Note that the chrom.sizes file is a tab delimited file with the format: <chromosome name><tab><number of bases in chromosome>

ADD COMMENT
2
Entering edit mode

You can skip all the bed stuff and just load a sorted and indexed bam file directly into IGV. It shows a coverage track along with the reads, so it should give you the same result with fewer steps.

ADD REPLY
0
Entering edit mode
ADD REPLY
3
Entering edit mode
9.2 years ago

If you want to eyeball coverage, one easy way is with IGV, from the Broad institute.

ADD COMMENT
0
Entering edit mode

IGV works well, but make sure you have a coverage file. This can be created using igvtools (should be under "Run --> Run igvtools...")

ADD REPLY
0
Entering edit mode

Like I said, you can eyeball coverage. Just throw the indexed .bam in there.

ADD REPLY
0
Entering edit mode
9.2 years ago
nicolazilio ▴ 20

Thank you, that worked!

ADD COMMENT

Login before adding your answer.

Traffic: 2169 users visited in the last hour
Help About
FAQ
Access RSS
API
Stats

Use of this site constitutes acceptance of our User Agreement and Privacy Policy.

Powered by the version 2.3.6