Getting consensus/assembly from RAD-seq data locus, non-model organism

Question

[Andres]

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,

Answer 1

[Philipp Bayer]

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