1000 genomes site BAM Files

[mez]

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.

[swbarnes2]

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.

[Alex Renwick]

... 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:

samtools mpileup -uf ref.fa aln.bam | bcftools view -cg - | vcfutils.pl vcf2fq > cns.fq

[nexgengirl]

Check out this tool:

http://picard.sourceforge.net/command-line-overview.shtml#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!

[laura]

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/wiki/index.php/Creating_a_Fasta_Alternate_Reference_%28incorporating_variants_into_the_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/