<p>How to separate reads, one for each reference haplotype.</p>
<ol>
<li>Map all the reads against the reference genome</li>
<li>Use GATK to call SNPs</li>
<li>Use GATK tool [FastaAlternateReferenceMaker][1] to build the alternative reference</li>
<li>You now have the two &quot;haploid references&quot; (quotes needed because you have to take into account that they will not necessarily reflect the exact haplotypic structure, it could be a mix of the two haplotypes)</li>
<li>Align the reads against the two references, and attribute each read to the reference to which they map better. My colleagues developed for a different purpose (detection of Allelic Imbalance) a pipeline that includes the step for determining to which of two references a read belongs. The tool is called [BayesASE][2]</li>
</ol>
<p>In general you may study all the work that has been done in the field of allelic imbalance and allele specific expression, where this problem has been encountered and solved.
[1]: <a href="https://gatk.broadinstitute.org/hc/en-us/articles/360037594571-FastaAlternateReferenceMaker" rel="nofollow noreferrer">https://gatk.broadinstitute.org/hc/en-us/articles/360037594571-FastaAlternateReferenceMaker</a>
[2]: <a href="https://github.com/McIntyre-Lab/BayesASE" rel="nofollow noreferrer">https://github.com/McIntyre-Lab/BayesASE</a></p>

CC BY-SA 4.0
2021-02-08T08:14:11.463Z
2021-02-08T08:14:11.463Z
Fabio Marroni
<p>My goal is to generate a diploid consensus genome. By &quot;diploid consensus&quot;, 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.</p>
<p>I have Illumina short-read files from the NCBI's SRA. These short-reads were obtained from a diploid species.</p>
<p>I am worried that, by simply using Picard's MarkDuplicates with the &quot;REMOVE_DUPLICATES&quot; option set to &quot;TRUE&quot;, I will retain only the homolog with the highest quality and output a haploid version of my genome.
Here is my current workflow:</p>
<pre><code>
## 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 &gt; SRR[number].sam  

## Sort Illumina reads by coordinate
java -jar <span class="math-container">$EBROOTPICARD/picard.jar SortSam I=SRR[number].sam O=sorted_SRR[number].sam CREATE_INDEX=true SORT_ORDER=coordinate TMP_DIR=$</span>TMPDIR

## Remove duplicate reads from sam file
## *Note: How do I remove sampling error reads?
java -jar <span class="math-container">$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=$</span>TMPDIR

## Convert duplicate-free sam file to bam file
samtools view -S -b nodup_SRR[number].sam &gt; 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 
</code></pre>
<p>Does anyone know of a program which will whittle my sam files down to ONLY two homologs per region?</p>

CC BY-SA 4.0
2021-02-04T23:09:38.717Z
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?
swbarnes2
CC BY-SA 4.0
2021-02-05T03:09:45.813Z
@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.
annabelperry
CC BY-SA 4.0
2021-02-05T03:10:25.737Z
@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: https://doi.org/10.1093/bioinformatics/bty279
annabelperry
CC BY-SA 4.0
2021-02-05T20:04:49.293Z
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  https://www.biostars.org/p/225262/
swbarnes2
CC BY-SA 4.0
2021-02-06T14:57:22.663Z
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).
Maximilian Press
CC BY-SA 4.0
2021-02-12T01:24:43.053Z
@MaximilianPress Ah I have failed to clarify my end goal. I am designing a software for microsatellite detection from short-reads (no such software currently exists which meets my needs). I only want to call variants in "real" reads rather than duplicates. However, I have found an alternative to diploid consensus creation - I am instead going to run my software on short-reads (including duplicates) and use variant frequencies to filter out false positives.
annabelperry
<sam><markduplicate>
Diploid Consensus Sequences
# Diploid Consensus Sequences

(Asked by: username_0 on 2021-02-03)

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:

## 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?

**Tags:** sam, markduplicate
> **Comments:**
> **username_1**: 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?  
> **username_0**: @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.  
> **username_0**: @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: https://doi.org/10.1093/bioinformatics/bty279  
> **username_1**: 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  https://www.biostars.org/p/225262/  
> **username_2**: 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).  
> **username_0**: @MaximilianPress Ah I have failed to clarify my end goal. I am designing a software for microsatellite detection from short-reads (no such software currently exists which meets my needs). I only want to call variants in "real" reads rather than duplicates. However, I have found an alternative to diploid consensus creation - I am instead going to run my software on short-reads (including duplicates) and use variant frequencies to filter out false positives.  

---

## 1 answers

### Answer by username_3 on 2021-02-08. Score: 1

How to separate reads, one for each reference haplotype.

1. Map all the reads against the reference genome
2. Use GATK to call SNPs
3. Use GATK tool [FastaAlternateReferenceMaker][1] to build the alternative reference
4. You now have the two "haploid references" (quotes needed because you have to take into account that they will not necessarily reflect the exact haplotypic structure, it could be a mix of the two haplotypes)
5. Align the reads against the two references, and attribute each read to the reference to which they map better. My colleagues developed for a different purpose (detection of Allelic Imbalance) a pipeline that includes the step for determining to which of two references a read belongs. The tool is called [BayesASE][2]

In general you may study all the work that has been done in the field of allelic imbalance and allele specific expression, where this problem has been encountered and solved.
[1]: <https://gatk.broadinstitute.org/hc/en-us/articles/360037594571-FastaAlternateReferenceMaker>
[2]: <https://github.com/McIntyre-Lab/BayesASE>


---

