1
$\begingroup$

Normally, I use bwa-mem to map my fastq reads for a bunch of species using one protein-coding gene CDS from my reference species. This was fine when I had to try it for just a group of genes...I would obtain single CDS fastq files and use them as a reference. However, we now want to obtain protein-coding sequences for ALL GENES in my exome datasets. If I have one file containing all the CDS from my reference species, can that be used in bwa-mem? The end of my custom script produces one multiple sequence alignment file when I ran it using one gene only...my concern is that bwa-mem will somehow replace certain CDS or would I just end up with over 20000+ MSAs? Which is fine because I need them to be individual MSAs for each protein-coding gene.

I use a simple python script at the end to extract base information from the mpileup files generated to create my MSAs.

Improve this question
$\endgroup$

1 Answer 1

Reset to default
1
$\begingroup$

If I have one file containing all the CDS from my reference species, can that be used in bwa-mem?

That's fine, normal usage of bwa-mem entails giving it a fasta file with all chromosomes/contigs/transcript/etc. in an organism, so your usage is no different. How your custom "make an MSA" script handles that is unknown of course.

Improve this answer
$\endgroup$
2
  • $\begingroup$ Yes unfortunately I need to figure out a new python script - I had tested it out with a small reference file containing about 10 CDS, but my final MSA was messed up it looked like it tried to merge all the CDS into one sequence. It made no sense at all. And my reference file contained a list of CDS so it was a list of 10 individual fasta sequences of diff lengths, not one fasta sequence concantenated. Maybe that was the problem? I am not sure how to make bwa-mem read each CDS line as its own reference sequence. $\endgroup$
    – CuriousDude
    Commented Sep 28, 2019 at 15:35
  • $\begingroup$ I am trying a new way out now where I assemble my reads against each chromosome from my ref species...but then I have no idea how to pull out all the individual protein-coding gene sequences. I am having trouble thinking of this differently since I have been running individual genes for the past year with no issue. $\endgroup$
    – CuriousDude
    Commented Sep 28, 2019 at 15:37

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.