SEQanswers

Go Back   SEQanswers > General



Similar Threads
Thread Thread Starter Forum Replies Last Post
Extract base from bam file empyrean Bioinformatics 4 07-03-2012 03:47 AM
Per base sequence coverage from sam/bam file? ewilbanks Bioinformatics 7 06-06-2012 01:03 PM
NEw to Chip-seq and have .bam/.sam/.bam.bai files... then what? NGS newbie Bioinformatics 11 05-25-2011 07:48 AM
Base coverage from Bam ElMichael Bioinformatics 4 12-01-2010 09:18 AM
Extracting base coverage and quality values from BAM files unagaswamy Bioinformatics 9 11-11-2010 06:55 AM

Reply
 
Thread Tools
Old 10-29-2013, 07:05 AM   #1
Coru S
Junior Member
 
Location: Germany

Join Date: Oct 2013
Posts: 4
Default How to extract base coverage from bam/bai-files

Hi Everybody,

I'm still quite new to the world of next-generation sequencing. I started to work with the Ion Torrent PGM several months ago.
I've got a question related to the analysis of the coverage depth of different nucleotides/dels/insertions on the base positions of my reads.

I normally use the bam and corresponding bai-files generated by the Torrent Suite Software and load them into the IGV-Genome Browser.
There, one can get information about single base position coverage by mouseover (see screenshot).

One can see the total base count, count of A/T/G/Cs (+ and - reads), and count of dels and insertions.
I would like to have this data in a table, to use it in an Excel spreadsheet.

Can anybody tell me, how I can extract this data from my bam/bai-files for a defined region of interest? I'd appreciate any hints.
The most preferred output would be a data table like this:
e.g.
chrom position total count count of:A+ A- T+ T- C+ C- G+ G- del+ del- ins+ ins-
chr13 33211005 50 20 21 2 3 0 0 0 0 2 0 1 1

I've searched the forums for related threads, but I only found some bioinformatic information that I couldn't understand. I'm just a biochemist, and I have no big skills in bioinformatics. E.g., I don't know about Linux OS and how to use software or commands there.
Any help is much appreciated.

Thank you very much in advance.
Greetings
Coru
Coru S is offline   Reply With Quote
Old 10-29-2013, 08:07 AM   #2
r_j_p
Junior Member
 
Location: New York

Join Date: Aug 2013
Posts: 4
Default

How big is your region of interest, and do you want one row for every position? If you do it could be a huge spreadsheet that Excel might struggle with.

A VCF file, (for example as generated by samtools), has similar information for sites with SNPs & indels, which you can then filter to help distinguish true variants from sequencing errors.
r_j_p is offline   Reply With Quote
Old 10-29-2013, 12:48 PM   #3
gringer
David Eccles (gringer)
 
Location: Wellington, New Zealand

Join Date: May 2011
Posts: 824
Default

I would use bedtools (which can use BAM files as input). Use genomecov to get the coverage all across the genome:

http://bedtools.readthedocs.org/en/l...genomecov.html

Or just coverage if you have a specific gff/bed file with the region that you want to determine the coverage for:

http://bedtools.readthedocs.org/en/l.../coverage.html

edit: sorry, that's just for the base-pair coverage (independent of the base). I suspect that per-nucleotide counts in a particular region is straying into "needs a custom script to work it out from the samtools mpileup output" territory.

Last edited by gringer; 10-29-2013 at 12:53 PM.
gringer is offline   Reply With Quote
Old 10-31-2013, 08:28 AM   #4
geneart
Member
 
Location: DC area

Join Date: Sep 2011
Posts: 42
Default

Hello all,
I have a similar related question to extract chromosome positions and depth per region. So I used bedtools to get bam to bed conversion:
I got :
Chr1 11419 11425 HWUSI-# 1 + 1
Chr1 13877 13892 HWUSI-# 0 - 2
Chr1 14714 14721 HWUSI-# 1 + 3
Chr1 16155 16163 HWUSI-# 1 - 4
Chr1 18239 18249 HWUSI-# 1 + 5
Now I looked up the manual for bedtools and it indicates that col5 as scores 1 to 1000 . I tried looking more into it, but cannot understand what those 1 and 0s stand for specifically.I know they are scores and that's it. I also have downstream in this column 32, 42 etc and not jut 1 or 0.I have cut and pasted just the top 5 outputs.
Also the last column looks like it is assigned a number for each cluster. Is this correct?
Please can some one clarify this for me?
Also if I need to get just the first 4 columns as is and also get the seq itself in the next column and counts of tags in the next columns, can I use bedtools to do it? If so , how? or do I need to write a script so that I have:
Chr# start end name seq # of tags pos/neg strand cluster#. from my bam files?
Please any help will do
Thanks in advance
geneart.
geneart is offline   Reply With Quote
Old 12-13-2013, 03:33 AM   #5
Coru S
Junior Member
 
Location: Germany

Join Date: Oct 2013
Posts: 4
Default

Hi All,

Thank you very much for your comments and hints.
I found a computer scientist who wrote a small program that can do the job.

Greetings
Rayk
Coru S is offline   Reply With Quote
Old 01-14-2014, 02:52 AM   #6
gray_cambs
Junior Member
 
Location: Cambridge UK

Join Date: Jan 2014
Posts: 1
Default

Corus S
Would you be able to post the code as I also would be interested in doing this search??

Cheers

Gray
gray_cambs is offline   Reply With Quote
Old 01-15-2014, 07:31 AM   #7
Coru S
Junior Member
 
Location: Germany

Join Date: Oct 2013
Posts: 4
Default

Well, I don't have the source code, just the final files of the program.

But I can ask the guy who wrote it, if he is willing to share the code.

I'll let you know...
Coru S is offline   Reply With Quote
Old 01-17-2014, 07:11 AM   #8
Coru S
Junior Member
 
Location: Germany

Join Date: Oct 2013
Posts: 4
Default

Hi Gray,

I contacted him, but unfortunately he doesn't want to reveal the code or the final files in the public at the moment, since the program is still beta status.

When he has extended the program he will make it publicly available in the final version.

Greetings
Coru S is offline   Reply With Quote
Old 01-17-2014, 01:52 PM   #9
gringer
David Eccles (gringer)
 
Location: Wellington, New Zealand

Join Date: May 2011
Posts: 824
Default

Quote:
I contacted him, but unfortunately he doesn't want to reveal the code or the final files in the public at the moment, since the program is still beta status.
Okay, well I've realised you can do this just using samtools and common Linux command-line programs. Here are base counts:
Code:
#forward
samtools view -F 0x10 sampled_lib19_extended.bam YAL005C_mRNA-51+51bp:500-3000 | awk -F '\t' '{print $10}' | fold -w 1 | sort | uniq -c
#reverse
samtools view -f 0x10 sampled_lib19_extended.bam YAL005C_mRNA-51+51bp:500-3000 | awk -F '\t' '{print $10}' | fold -w 1 | sort | uniq -c
#generic
samtools view [filter] input.bam sequence:bpRange | awk -F '\t' '{print $10}' | fold -w 1 | sort | uniq -c
[These one-liners can be easily modified to put sequence names at the start of the lines if necessary]

For a quick INDEL count, you can use the CIGAR column (#6):

Code:
#generic (use -f 0x10 or -F 0x10 for forward/reverse filter as necessary)
samtools view [filter] input.bam sequence:bpRange | awk -F '\t' '{print $6}' | fold -w 1 | grep '[ID]' | sort | uniq -c
gringer is offline   Reply With Quote
Old 01-27-2014, 12:07 PM   #10
yl01
Member
 
Location: Germany

Join Date: Aug 2012
Posts: 25
Default

Igvtools (which is a command line tools associated with IGV) count can do this task easily. It produces a text file in wig format. However it may not easy to open it in Excel as desired.
yl01 is offline   Reply With Quote
Old 01-03-2019, 03:50 PM   #11
Katherine_B
Junior Member
 
Location: Toronto

Join Date: Sep 2016
Posts: 4
Default

Quote:
Originally Posted by yl01 View Post
Igvtools (which is a command line tools associated with IGV) count can do this task easily. It produces a text file in wig format. However it may not easy to open it in Excel as desired.
Could you possibly share the igvtools command line to do this? I have been trying to figure out how to get the base coverage values at each position from my BAM files. I'd want to know which base is at a position and how many reads produce that base so I can create some custom scripts using that information.
I tried:

igv count --bases input_sorted.bam output_basecounts.wig

But my .wig files are empty. Should I be saving the output differently??

Last edited by Katherine_B; 01-03-2019 at 03:53 PM. Reason: Forgot to include the command I used
Katherine_B is offline   Reply With Quote
Old 01-03-2019, 04:53 PM   #12
HESmith
Senior Member
 
Location: Bethesda MD

Join Date: Oct 2009
Posts: 502
Default

Use BEDtools 'genomecov' command with '-d' flag to obtain per-base depth of coverage.
HESmith is offline   Reply With Quote
Reply

Thread Tools

Posting Rules
You may not post new threads
You may not post replies
You may not post attachments
You may not edit your posts

BB code is On
Smilies are On
[IMG] code is On
HTML code is Off




All times are GMT -8. The time now is 10:12 PM.


Powered by vBulletin® Version 3.8.9
Copyright ©2000 - 2019, vBulletin Solutions, Inc.
Single Sign On provided by vBSSO