0
$\begingroup$

I'm new to bioinformatics. I have a problem in which I have a FASTA reference genome and lots of reads in FASTQ files. Some of them could be contaminants, so I'd like to filter them out and get only the ones that align to the reference genome. I also have access to the BAM/BAM.BAI files associated with each FASTQ, although I can't know how they were aligned.

Thanks

$\endgroup$

2 Answers 2

3
$\begingroup$

Your post has multiple questions but one of them is not really a question so I'll address everything here:

1. "Some of them could be contaminants, so I'd like to filter them out and get only the ones that align to the reference genome"

Just align the reads. Contaminants should not align well and should be discarded owing to poor mapping quality score.

2. How to figure out if the FASTA file was used for the alignment

This is not really a 100% perfect solution but the BAM header should contain the exact command used for the alignment, which should contain the ref FASTA with path. You can also compare contig names and lengths between the FASTA and the BAM @SQ lines.

$\endgroup$
2
  • $\begingroup$ "Contaminants should not align well" assumption would not quite hold all the time, for example if the data is from a PDX model involving mouse as well as human reads! $\endgroup$
    – haci
    Feb 19 at 18:47
  • $\begingroup$ That's true. I was about to suggest something like Xenome or BBSplit but thought these might be simple bacterial contaminants that may not require preprocessing. $\endgroup$
    – Ram RS
    Feb 19 at 20:16
2
$\begingroup$

In general, this answer wouldn't be more correct than @Ram RS' answer but here is how I would approach:

  • I would start everything from scratch, at the end of the day you will most probably need to document the exact genome version and its associated annotation file. It looks like you do not have this information at present.

  • I don't know the experimental set-up but if you are worried about "contaminants", you can have an extra alignment step using BBSplit (link), it is capable of mapping to multiple genomes simultaneously.

  • I would give a try with nfcore/rnaseq, it abstracts away a lot of complex bioinformatics steps so that you can focus on understanding what these tools do instead of dealing with installing and running these tools. For example you can just use the --bbsplit-fasta-list argument to get rid of your contaminants.

$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

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