Hi I have RAD-seq libraries from a non-model organism. I mapped those reads (with BWA) to a gene from a related species. I need to get a "consensus" from the mapped reads. The thing is that the coverage is uneven 30-300x for some areas and 0 for others (That is the nature of RAD-seq data).
What would be the most reasonable/acceptable way to get a consensus/aseembly from tose reads. Taking into account that at any given locus i have like 300bp of a consistent high coverage region. The gene in question is about 900 bp.
I want that consensus so that i can incorporate that sequence into a phylogenetic gene tree from an old paper were they have a lot a sequences from the genera i'm working with.
I've already tried to use bcftools consensus but i feel that it relies too much in the consensus sequence. ¿Is there any good few-reads small-contig assembler?. Or what do you recommend?
Any help is highly appreciated, best regards,
1 answer
Armin Scheben has a nice protocol on RADSeq data analysis; he suggests using Stacks for this. https://github.com/ascheben/RAD_analysis_workflow
He has this line at the end of the Stacks run:
populations -P ${datadir} -M ${popmap} -t ${threads} --vcf
You can alternatively have populations also output consensus fastas via --fasta-loci instead of --vcf
Log in to answer this question.