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.

bowtie2-build  

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. 




9 comments:

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.
divinespellhome@gmail.com
WhatsApp him now +1(267)527-9481

Anonymous said...

POWERFUL LOTTERY SPELL CASTER DR GBOJIE 2018/2019
i am very grateful sharing this great testimonies with you, The best thing that has ever happened in my life is how i win the lottery. I am a woman who believe that one day i will win the lottery.finally my dreams came through when i email Dr gbojie . and tell him i need the lottery numbers. i have come a long way spending money on ticket just to make sure i win. But i never know that winning was so easy until the day i meant the spell caster online which so many people has talked about that he is very great in casting lottery spell, so i decide to give it a try.I contacted this man and he did a spell and he gave me the winning lottery numbers. But believe me when the draws were out i was among winners. i win 1.900.000 million Dollar. Dr. gbojie truly you are the best, with these man you can will millions of money through lottery. i am so very happy to meet these man, i will forever be grateful to you. Email him for your own winning lottery numbers gbojiespiritualspelltemple@gmail.com. OR call him +2349066410185.or check out his website :http://gbojiespiritualtemple.website2.me

Lucy Morgan said...

God is Good!
I promised God that I would share my testimony on this website. I had all the signs of  STD Virus but I was not too sure as to which one. I did a lot of online research and scared myself straight for a whole week before going to see the nurse. She took one look at my genital part and first said that it could just be the anatomy of my body, then she said it looked like genital warts and that I MAY have herpes. I was devastated. She gave me some medicine for the herpes and some cream for the warts. I was also tested for every single STD including herpes. I went home and cried searching the web for all sorts of cures for herpes and awaiting my results. I saw a post whereby Dr. Kham cured Herpes and other diseases, I copied his contacts out and added him on whats app via +2348159922297. The next day my test result was ready and I confirmed  Herpes positive. I told Dr. Kham about my health problems and he assured me of a cure. He prepared his herbal medicine and sent it to me. I took it for 14 days (2 weeks). Before the completion of the 14 days in which I completed the dose, the Blisters and Warts that were on my body were cleared. I went back for check-up and I was told I'm free from the virus. Dr. Kham cures all types of diseases and viruses with the help of his herbal medicine. You can reach Dr. Kham via his email address on (dr.khamcaregiver@gmail.com) or call him on +2348159922297 and whats app him on same number, but You can Also Visit his website to know more about him at > https://drkhamcaregiver.wixsite.com/drkhamcaregiverherba  for more info.  

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.

Tamara Barrow said...

Greetings to the general public, I want to tell about how I was cured of
HIV disease by a Doctor called Dr. voodoo . I was browsing the Internet
searching for remedy on HIV and I saw a comment of people talking about how
Doctor voodoo cured them. I was scared because I never believed in the
Internet but i was convince to give him a try because i having no hope of
been cured of HIV so I decided to contact him with his email that was
listed on the comment voodoospelltemple66@gmail.com when I contacted
him he gave me hope and send a Herbal medicine to me that I took and
it seriously worked for me, am a free person now without problem,
my HIV result came out negative. I pray for you Dr. voodoo God will give
you everlasting life, you shall not die before your time for being sincere
and great men. Am so happy, you can also contact him if you have any
problem Email: voodoospelltemple66@gmail.com
or Add him up on WhatsAp +2348140120719

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.