Jan 31, 2018

Genome Coverage from BAM file

There are many excellent tools for analysis of Next Gen Sequencing data in the standard BAM alignment format so I was surprised how difficult it was for me to get a nice graph of genome coverage.  This will be trivial for a lot of hard core bioinformatics coders, so just move along if you are bored/annoyed.

I needed to check the evenness of coverage across intervals of a bacterial genome that we were re-sequencing for various experimental reasons. I aligned my FASTQ to a reference genome from GenBank using Bowtie2. There are several nice tools in the SAMTools and BEDTools kits that produce either a base by base coverage count or a histogram of coverage showing how many bases are covered how deeply.  I wanted a map at 1 Kb resolution. It took a while to figure out that I first need to make a BED file of intervals from my genome - with correct names for the 'chromosomes' that match the SAM header, and then use 'samtools bedcov' to get my intervals. Then a simple graph from Excel or R shows the coverage per interval along the genome.

Here are the steps (as much for me to remember as for usefulness to anyone else)

1) Create a Bowtie2 index of the reference genome (can be from GenBank, or can be a de novo assembly of contigs created locally from this FASTQ data).   Reference_in must be FASTA format.
bt2_base is the name you will call the index.


2) Align the FASTQ file(s) to the Reference. bt2-idx is the name of the index created in the previous step. There are a ton of options for stringency of alignment, format of input data, etc. I set this to use 16 CPU threads.  I usually like to leave out the unaligned reads, which can reduce the output file size somewhat. [Of course, the unalinged reads are the goal when you use Bowtie as a filter to remove human from microbiome data or any other sort of contaminant screen.] 

bowtie2 -p 16 
--no-unal -x -1 -2 -S

3) Is is annoying that Bowtie produces output in SAM format, and 99% of the time, the very first thing you have to do is convert to sorted BAM. Note that samtools sort puts its own .bam on the end, so if your are not careful you will get files named output.bam.bam

samtools view -bS output.sam | samtools sort - file_sorted

4) Create a 'genome' file for your reference genome. This is just a tab delimited file that names the chromosomes (or individudal sequences that may be in a multi-FASTA file, such as contig names).
It is super easy to mess this up. The best way is to view the top of your output.sam file created by bowtie2. The lines that start with @ are your chromosome headers, and they very helpfully already show the length of each one.  This is a bit of a pain if you  have a genome with lots of contigs, but a little 'cut' and 'paste' in bash or Excel will get you there.

Here is what mine looked like:

@HD     VN:1.0  SO:unsorted
@SQ     SN:MKZW02000001.1       LN:5520555
@SQ     SN:MKZW02000002.1       LN:248293
@PG     ID:bowtie2      PN:bowtie2      VN:2.2.7        CL:"/local/apps/bowtie2/2.2.7/bowtie2-align-s --wrapper basic-0 -p 8 -x Kluy

and here is the genome file I made:

MKZW02000001.1  5520555
MKZW02000002.1  248293

5) Make a set of intervals with bedtools makewindows. I wanted 1 Kb intervals, so I use -w 1000.
The result is a simple BED file with one line for each 1Kb window of the genome. 

bedtools makewindows -g genome.txt -w 1000 > genome_1k.bed

$ more genome_1k.bed

MKZW02000001.1  0       1000
MKZW02000001.1  1000    2000
MKZW02000001.1  2000    3000
MKZW02000001.1  3000    4000

6) Use samtools bedcov to count the total number of bases in the BAM file that are located in each of the intervals ('sum of per base read depths per BED region'). This works much faster than any other coverage tool that I have tested. 

samtools bedcov genome_1k.bed kv_sorted.bam > kv_1k.cov

7) Divide the sum of coverage by the window size (/1000 in my case), and plot the average coverage per window as a scatter plot, using the end of each interval as the X axis and the coverage as the Y. 
Histogram will only work nicely if you have very few intervals. In my case, high and low coverage outlier intervals are easily visible. 


Tim said...

Seen deeptools bamCoverage? https://deeptools.readthedocs.io/en/latest/content/tools/bamCoverage.html

Unknown said...

With a grateful hearth , I want to give my sincere appreciation to Dr Sambo He is the best , he's so wonderful and helping. Few weeks ago , I contacted him for a lotto lucky number, without delay , he gave me the real lucky winning numbers I played and won $100,000,000. Contact him for lucky winning numbers and your story will change for good . Thank you Dr Sambo I will forever be grateful for your kind gesture. Contact him today through his email he can help you.
WhatsApp him now +1(267)527-9481

Anacyte Laboratories said...

Very nice posting. Your article is quite informative. Thanks for the same. Our service also helps you to market your products with various marketing.

Click Here:- single cell rna sequencing

Donna said...

I am Doctor Paul I got affected with HIV in the process of attending to my HIV patient I tried all I can to get cured but all to no avail, until I saw a post in a health forum about a herbalist man who prepare herbal medication to cure all kind of diseases including HIV virus, at first I doubted if it was real but decided to give it a try, when I contact this herbalist via his email Blessedlovetemple@gmail.com and he prepared a HIV herbal cure and sent it to me via fed-ex delivery company service, when I received this herbal cure, he gave me step by directions on how to apply it, when I applied it as instructed, I was totally cured of this deadly disease within 5 days of usage, I am now free from the deadly disease called HIV, all thanks to Dr Mark. Contact this great herbal spell caster. Kindly contact him. Blessedlovetemple@gmail.com
He cures all kinds of sickness or diseases such as: 1. HERPES VIRUS 2. LASSA FEVER 3. GONORRHEA 4. HIV/AID 5. EX BACK.
Thanks Dr Mark for saving my life.

Ric Clayton said...

I really want to thank Dr Emu for saving my marriage. My wife really treated me badly and left home for almost 3 month this got me sick and confused. Then I told my friend about how my wife has changed towards me. Then she told me to contact Dr Emu that he will help me bring back my wife and change her back to a good woman. I never believed in all this but I gave it a try. Dr Emu casted a spell of return of love on her, and my wife came back home for forgiveness and today we are happy again. If you are going through any relationship stress or you want back your Ex or Divorce husband you can contact his whats app +2347012841542 or email emutemple@gmail.com website: Https://emutemple.wordpress.com/

Unknown said...

I want to share this great testimony to the world on how Dr VOODOO cure me from HSV1&2 with his herbs, I was nervous when i first contact him about the cure for HSV but i decided to give him a try because i was desperate to get cured and be free from HSV. Dr VOODOO prepared the remedy and sent it to me through UPS,which I use just the way Dr VOODOO instructed me and thank God today I am a beneficiary to this cure and I went back to the hospital after 7 days of taking the herbs and I tested HSV1&2 Negative. So I will tell you all who are looking for a cure to his/her HSV1&2 that Dr VOODOO took research before he could finally get the solution to it and a lot of people are benefiting from him right now. He also cured my friend from HPV. Dr VOODOO heals with natural herbs. Please I urge you to contact him now through his email address: voodoospelltemple66@gmail.com or call and WhatsApp him on +2348140120719. He is capable of curing HIV/AIDS, HERPES, HPV, HSV1&2, CANCER of all kinds, DIABETES and so many other infections.