What is the reference sequence? ( can I find it in a .bai index?)

[keebs42]

Working with the new 1000 genomes (pilot2) release, there are .bai files associated with the new .bam files. I've used 'samtools pileup' to create the pileup files, but the reference column is populated with 'N'. I am trying to determine the reference base (without simply looking at the ncbi human reference and comparing positions), but can't seem to find out how. Is there a way to a) use the .bai file in the process of making the pileup so the reference base column is filled correctly, or b) extract the reference directly from the .bai file?

Sorry if I've missed something obvious, but I can't find anything in the samtools documentation that answers my question.

Thanks,
Jonathan

[keebs42]

Going to try setting the -f flag when running samtools pileup..

[JackyH]

Hi,

have you managed to solve this problem at all? I am seeing exactly the same and can't quite figure out where I've gone wrong.

Thanks,

Jacky

[keebs42]

I ended up finding downloading the NCBI human reference build 36 in fasta format (split by chromosome), and then using the -f flag when creating the pileup. The reference column seems to be correct after doing this.

samtools pileup -f ref.fasta alignment.bam > alignment.pu

Hope that helps!
Jonathan

[JackyH]

Mh, that's curious. I have tried that with RefSeq as a reference but I still don't see the reference base. Maybe it's an issue with the format of the Fasta header (contains a colon in my case).

Thanks for your fast reply anyhow!
Jacky

[RAS]

Still can't get correct reference sequence column

Hi,

I still can't get a pileup file where the reference sequence shows bases instead of "N"s. I'd like to create pileup files of sequences from the 1,000 genomes project aligned with NCBI human reference sequences. I am using the -f flag--indicating that the reference sequence is in FASTA format--and also need to use the -c flag--indicating that the pileup file should have the consensus sequence for the original .bam file. In the main samtools-0.1.7a folder of a GNU/Linux computer, I've typed many variants of the following:

./samtools pileup -cf /ifs/scratch/.../humanReferenceGenome/UCSCBuild36/chrX.fa /ifs/scratch/.../fatherAlignment/NA12891.chromX.ILLUMINA.bwa.CEU.high_coverage.20100517.bam > NA12891.ChrX.UCSC36.pileup

The resulting pileup file contains the NA12891 consensus sequence. I've tried using a number of reference sequences, including builds 36.1, 36.2, 36.3 and 37 of the NCBI reference genome and build 36 of the UCSC reference genome, in the hope that one of these reference sequences would also appear in the pileup file. I would very much appreciate any suggestions.

Thanks,
Rebecca