0
$\begingroup$

I have two resources of .bam files. One is generated by our lab (1 sample = 1 bam). One is downloaded online (again 1 sample = 1 bam). For the downloaded samples the chromosomes are labelled: chr1, chr2, chr3 etc For our lab samples, the chromosomes are labelled: 1, 2, 3 etc. I want to generate a single VCF file of variants across all samples.

I'm using bcftools:

bcftools mpileup -Ov -f ref.fasta -b samples.txt | bcftools call -mv -o bamMge.vcf

However, I get no calls and the repeated error: [E::faidx_adjust_position] The sequence "1" was not found

I have two questions:

  1. Is my strategy correct (i.e. bcftools mpileup)? Do I need to incorporate extra steps (NB I also have matching gVCF files)

  2. How can I either (i) alter the chromosome labeling of one of the subsets of bam files, or, (ii) use a mapping file to match chromosome labels during the mpileup run?

Thanks

$\endgroup$
2
$\begingroup$

The main issue here is that your bam files have different chromosome labels. (i.e. 1,2,3 vs. chr1,chr2,chr3) as you mentioned.

This suggests that the data was mapped to a slightly different version of the reference genome. As a first layer, you might want to double check that these reference chromosomes are cross compatible (e.g. both correspond to hg19 for instance). To find more information on the reference genome used, you can check the header information using "samtools view -H "

Assuming that the reference genomes are cross compatible, you can try converting it to sam, adding the missing "chr" label in a pipe, and converting it back to a bam file. However, this will also require modification of the header file which can be challenging in itself.

The easiest (and more rigorous) approach is probably just to convert the batch of bam files with the wrong chromosome labels back into fastq files. You will then want to remap them using the exact same reference genome as the other batch of bam files. After which, you should be able to perform variant calling on the combined batch of bam files.

| improve this answer | |
$\endgroup$
1
$\begingroup$

I am not familiar with getting variants from bam files, so I will just answer the second part. From the error message, it seems that it is necessary to have consistent chromosome labels across datasets. It is also a good idea to do so because this keeps your datasets consistent and makes things a lot easier. To change the chromosomes labels, you can write your own scripts (I use Python, for example). You can also find several solutions on BioStars (https://www.biostars.org/p/98582/).

| improve this answer | |
$\endgroup$
  • $\begingroup$ Editing BAMs or VCFs to change chromosome names is really dangerous if you're not 100% sure you know what you're doing. I've seen scripts that screwed up the file contents in subtle enough ways to mess up downstream analyses. Even if you avoid that, you can end up combining data aligned to different genome reference builds, which also ends up affecting analysis results. It's much safer to realign from scratch like KT8 describes. $\endgroup$ – Geraldine_VdAuwera Apr 23 at 16:00

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.