I am trying to call consensus by mapping PE data onto assembled genome by following the below mentioned steps in that order:
- mapping reads onto genome (bwa mem)
- sam to bam to sorted bam (samtools)
- sorted bam to vcf (samtools mpileup)
- vcf to fastq (samtools-bcftools vcfutils)
- fastq to fasta (seqret)
Is it expected to get the a significantly less number of fasta records in the consensus fasta file as compared to the number of scaffolds in the assembly fasta file? The interesting part is that the total number of bases in the assembly fasta file and the consensus fasta file is almost equivalent. That means that the data was not lost. This is again supported by a significantly good alignment of the PE reads onto the assembled genome.
So, the question is, what exactly happened while calling the consensus?
probably already the first step would have many reads discarded as unmapped.... look into the stats of each step and/or into number of alignments in BAM (eg: wc -l my.fastq (divided by 4), wc -l my.sam, etc)
Dear Michal
98% of my reads mapped to the genome which suggests not many reads were discarded.
- Badges
- Science topic
- Similar topics
- Biological Science
- Physiology
- Physiological Processes
- Stress