# Make a fasta out of mapped reads without taking into account the reference

I've aligned reads onto an haploid reference genome. I'd like to have a consensus sequence of all my aligned reads that doesn't take into account the reference genome I used (i.e. I just want the nucleotides from my reads at each base (taking into account frequency and quality) but if there is no read a one position I rather have N/? than the reference nucleotide. Is there a way to do this ? Thanks

vSNP outputs a VCF file that includes positions with no coverage, represented as "N". Then there is a tool to merge a VCF file into the FASTA reference vsnp_merge_vcf_into_fasta.py. It will merge those "N" regions into the FASTA reference. There are options for considering frequency and quality. It's limited to just SNPs and those regions along the reference with no coverage. It does not merge small indels or unmapped reads. Pilon is a better tool for applying indels, but I'm not aware of it letting regions with no coverage be represented as "N".

vSNP_step1.py -r1 *_R1*fastq.gz -r2 *_R2*fastq.gz -r ref.fasta


Use *zc.vcf output to merge into FASTA

vsnp_merge_vcf_into_fasta.py -f ref.fasta -v *_zc.vcf


This is similar to this question, which links to this solution.

The difference here is that you want to go through and replace low coverage (maybe zero-coverage?) positions with Ns rather than using the reference. You can likely accomplish this using samtools depth to measure coverage and then writing a script that will go through each position in your reference, check its coverage, and string-replace with an N. Probably a bash one-liner could do this if someone very devoted to bash one-liners came up with one (I don't care for them, but other people are more creative than I). Very quick and dirty pseudocode might look like this:

### Someone could come up with something better if they spent >3min on it.
### should be easy to get it working if you use biopython
# read_and_preprocess() stands for whatever you need to do to get data in the right shape
# reference_with_variants.fasta comes from the solution linked above

threshold = 0
fasta_string = ''
for position in reference_fasta:
if coverage_data[position] <= threshold:
fasta_string += "N"
else:
fasta_string += reference_fasta[position]

# nonexistent fn that writes a fasta.
write_as_fasta(fasta_string)


I am a little skeptical of this approach because it seems like if you do not trust the reference sequence to stand in for zero-coverage areas, then you probably shouldn't trust it to align to and call variants against. Read mapping can be a rather noisy process, especially when you are using a more distantly related reference. I think this is related to the @swbarnes2 comment.

I would suggest instead indicating the regions of low/zero coverage by putting them in lowercase acgt with your covered regions uppercase ACGT, which preserves the coverage information that you care about but still allows you access to the reference information. It's not a perfect solution, but I think it makes the most out of the marginal data you're working with, and it is trivial to change the above example code to do it.

Hope that helps-

Have you googled de novo assembly?

• I can’t do a de-novo assembly as I have really low coverage. I’m working with really short reads from ancient DNA. – LauraR Jan 31 '20 at 20:07
• You can make a consensus fasta from your bam, but that uses the reference. If your read depth is that low, can you be confident of discrepancies between your reads and the reference? – swbarnes2 Jan 31 '20 at 20:36
• No, but that’s another problem. Thanks. The thing I want to do can also be useful when you want to reassemble only some targeted genes from medium coverage data I guess. – LauraR Feb 1 '20 at 8:22