SEQanswers

Go Back   SEQanswers > Bioinformatics > Bioinformatics



Similar Threads
Thread Thread Starter Forum Replies Last Post
1000 Genomes VCF Files ashkot Bioinformatics 8 02-25-2013 03:12 PM
confused by all the 1000 Genomes fastq files efoss Bioinformatics 3 10-11-2011 07:43 AM
bam files of exomes from 1000 genomes jorjial Bioinformatics 3 06-06-2011 04:16 AM
using BWA to align SOLiD fastq files from 1000 Genomes tgenahmet Bioinformatics 1 10-15-2009 06:19 PM
.bam alignment format from 1000 Genomes Project CellsDividing Bioinformatics 2 02-06-2009 09:21 AM

Reply
 
Thread Tools
Old 03-28-2012, 02:17 PM   #1
mez
Junior Member
 
Location: georgia

Join Date: Nov 2011
Posts: 5
Default 1000 genomes site BAM Files

Please,

I would like anyone who has had experience working with the .bam alignment files from the 1000 genomes site to help me out with this situation. I have been trying for long to get a complete FASTA sequence of the sample DNA from the .bam files for a given genomic region. I had used the samtools view, sort, index as well as pile up commands to get a pile up file and eventually extracted the consensus DNA sequence from the fourth column of the pile up file, creating the FASTA file.

The problem I am having with this approach, is that the raw DNA sequence seem to be obviously shorter in each case than the DNA sequence from the reference genome, and for different samples I get different DNA sequence length. I do not know if the issue is because of the low coverage sequence approach of the 1000 genomes site, or if I am doing something wrong.

Thanks and will appreciate any guidance.
mez is offline   Reply With Quote
Old 03-28-2012, 03:07 PM   #2
swbarnes2
Senior Member
 
Location: San Diego

Join Date: May 2008
Posts: 912
Default

You'd think that there woud be an out-of-the-box answer to making fastas from .bams, but there isn't.

The first thing I see wrong with your approach is that pileup will simply not show any line if it is totally uncovered. Are you padding your fasta with N's when there is no corresponding line in the pileup?

Older versions of samtools came with a samtools.pl that had a utility to make consensus fastq of pileups. You could try using that. perl scripts are flat text files, you can read them like an ordinary text file.

Search around seqanswers a bit, I know I've written previously on how to convert a perl script that writes a fastq into one that writes a fasta. It's just a matter of changing a few characters in the script.
swbarnes2 is offline   Reply With Quote
Old 03-28-2012, 03:11 PM   #3
Alex Renwick
Member
 
Location: Houston, Texas

Join Date: Jul 2011
Posts: 44
Default

Quote:
Originally Posted by mez View Post
... extracted the consensus DNA sequence from the fourth column of the pile up file, creating the FASTA file.
...
By "fourth column" I suppose you mean "third column", since the fourth column of the pileup format is the count of reads covering the site. See pileup format.

The sequence you get from the third column is the reference sequence, not a consensus sequence. It will be shorter because sites with no read coverage won't appear.

See http://samtools.sourceforge.net/mpileup.shtml for instructions on using samtools to generate a consensus sequence:

Code:
samtools mpileup -uf ref.fa aln.bam | bcftools view -cg - | vcfutils.pl vcf2fq > cns.fq
Alex Renwick is offline   Reply With Quote
Old 03-28-2012, 07:05 PM   #4
nexgengirl
Member
 
Location: Maryland

Join Date: Apr 2010
Posts: 31
Default

Check out this tool:

http://picard.sourceforge.net/comman...tml#SamToFastq

It will convert a sam or a bam to a fastq. That would get you at least half way to a fasta file.

I'm pretty sure there have been other posts on how to convert from fastq to fasta+qual so I would check on seqanswers for more about that.

Good luck!
nexgengirl is offline   Reply With Quote
Old 03-29-2012, 03:41 AM   #5
laura
Senior Member
 
Location: Cambridge UK

Join Date: Sep 2008
Posts: 151
Default

If you are after all the sequence for a particular individual from the 1000 genomes project you are much better using our fastq files rather than our bam files

ftp://ftp.1000genomes.ebi.ac.uk/vol1....sequence_data
ftp://ftp.1000genomes.ebi.ac.uk/vol1/ftp/sequence.index

If you want just the sequence which has been aligned to a particular region then a solution like picards SamtoFastq which has already been suggested would be best

If you want the reference genome annotated for all the variants in a particular individual then you should look at solutions like

http://www.broadinstitute.org/gsa/wi...e_reference%29

or

vcfutils.pl vcf2fq from samtools

amd our vcf files which you can get from

ftp://ftp.1000genomes.ebi.ac.uk/vol1...ease/20110521/
laura 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 09:38 PM.


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