# Diploid Consensus Sequences

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.

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:


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?

• What you are asking for doesn't make sense. How on earth is any software supposed to determine if a discrepancy is a polymorphism or an error if you remove all the supporting reads? Feb 4, 2021 at 23:09
• @swbarnes2 It's not that I want to simply remove the supporting reads - I should clarify my intent. My end goal is to use the consensus genome as input in a homebrew microsat-calling software. I don't want this software to call sequencing errors as microsats, so I want to merge all reads of the same sequence into a single consensus read based on base quality scores, generating two sets of sequences which are the pipeline's best guess at the individual's haplomes. I will edit my original post to make this clear. Feb 5, 2021 at 3:09
• @swbarnes2 I want to do what is described in this article, but I don't have PacBio long-reads to compliment my Illumina short-reads: doi.org/10.1093/bioinformatics/bty279 Feb 5, 2021 at 3:10
• FWIW, I do not think that removing lines from a bam file is what you want to do at all. I think you need to be calling variants with phasing data, and using that vcf to make a modified fasta. A quick google turned this up biostars.org/p/225262 Feb 5, 2021 at 20:04
• This is a little hard to answer as is, could you give an example of the kind of data you want to get out? Do you want to phase the haplotypes, like in the linked paper? I think swbarnes2 has the right of it, that calling variants is the way to go. Also- duplicates are usually defined as reads with the same coordinates; unless you have very high coverage this is unlikely to remove more than a very small fraction of informative reads (they are more likely PCR or optical duplicates, i.e. artifacts). Feb 6, 2021 at 14:57