My goal is to generate a diploid consensus genome. By "diploid consensus", I mean that I want to merge all supporting reads for a genomic region into a single haplome such that the resulting sam file contains TWO sets of reads, one from each haplome.
I have Illumina short-read files from the NCBI's SRA. These short-reads were obtained from a diploid species. After alignment, I want to remove all duplicate reads and sequencing error reads so the sam files have only two homologs for each genomic region.
I am worried that, by simply using Picard's MarkDuplicates with the "REMOVE_DUPLICATES" option set to "TRUE", I will retain only the homolog with the highest quality and output a haploid version of my genome. Here is my current workflow:
## Map Illumina reads to the reference genome. Include read IDs
bwa mem -M -R '@RG\tID:BioSample.LibraryName\tSM:BioSample\tLB:LibraryName\tPU:SingleFlowcell.SingleLane\tPL:ILLUMINA\t' reference.fasta SRR[number]_1.fastq SRR[number]_2.fastq > SRR[number].sam
## Sort Illumina reads by coordinate
java -jar $EBROOTPICARD/picard.jar SortSam I=SRR[number].sam O=sorted_SRR[number].sam CREATE_INDEX=true SORT_ORDER=coordinate TMP_DIR=$TMPDIR
## Remove duplicate reads from sam file
## *Note: How do I remove sampling error reads?
java -jar $EBROOTPICARD/picard.jar MarkDuplicates REMOVE_DUPLICATES=true REMOVE_SEQUENCING_DUPLICATES=true I= 'sorted_SRR[number].sam' O= nodup_SRR[number].sam METRICS_FILE= metrics_SRR[number].txt TMP_DIR=$TMPDIR
## Convert duplicate-free sam file to bam file
samtools view -S -b nodup_SRR[number].sam > nodup_SRR[number].bam
## Re-align reads around indels to reduce slippage.
gatk LeftAlignIndels -R reference.fasta -I nodup_SRR[number].bam -O aligned_SRR[number].bam
Does anyone know of a program which will whittle my sam files down to ONLY two homologs per region?