Non-standard nucleotides in my fast converted from a bam

Question

[biogirl]

Hi,

I have converted my sorted aligned BAM file (paired-end short reads aligned to reference) to a consensus fasta file using the smalls mpileup, vcfutils and smalls seqtk pipeline, but I have some non standard nucleotides (i.e. not A,C,G,T) in my fasta file. I assume that these are artifacts of the BAM file. For instance, I have normal nucleotide sequence, and then I'll have a 'W', or a 'K'. Are these following the IUB/IUPAC nomenclature, or are they something different?

Thanks

[Devon Ryan]

I actually wonder if those were inserted downstream of the BAM file, likely by vcfutils or seqtk. In order to have IUPAC codes in a BAM file either (1) the original read had to contain that code of (2) the aligner had to produce it. Both of those seem unlikely.

[biogirl]

Very good point, thank you. I'll investigate the conversion process.

[biogirl]

It seems it happens in the conversion of the BAM to the fastq. I guess this could be due to SNPs/ploidy and it can't decide what the consensus nucleotide is.

[Devon Ryan]

It looks like vcfutils is adding the IUPAC codes then: Questions Regarding Consensus Sequence Calling With Samtools / Bcftools / Vcfutils.Pl

[biogirl]

That is incredibly handy, thanks for posting that link!

[Devon Ryan]

No problem, you eliminated the other likely possibility!

Answer 1

[Matt Shirley]

You mention that you are making a consensus FASTA file. In addition to Devon's comment, there is another source of IUPAC ambiguous base codes: your reference FASTA file. I'll wager a bet that you're using hg38:

[Matt Shirley]

Note that based on the comment linked below, the UCSC version of hg38 replaces all ambiguous bases with 'N'.

GRCh38 BAM with hg38 VCFs

[biogirl]

Hi, I have checked my reference already, and there aren't any ambiguous bases in there (it's not hg38 either!). I'm currently checking the pipeline to see where it spits out these bases.