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

Question

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

Answer 1

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

Answer 2

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
coverage_data = read_and_preprocess('samtools_coverage.txt')
# reference_with_variants.fasta comes from the solution linked above
reference_fasta = read_and_preprocess('reference_with_variants.fasta')

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-

Answer 3

Have you googled de novo assembly?