0
$\begingroup$

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?

Improve this question
$\endgroup$
6
  • $\begingroup$ 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? $\endgroup$
    – swbarnes2
    Commented Feb 4, 2021 at 23:09
  • $\begingroup$ @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. $\endgroup$
    – annabelperry
    Commented Feb 5, 2021 at 3:09
  • 1
    $\begingroup$ @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 $\endgroup$
    – annabelperry
    Commented Feb 5, 2021 at 3:10
  • 1
    $\begingroup$ 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 $\endgroup$
    – swbarnes2
    Commented Feb 5, 2021 at 20:04
  • $\begingroup$ 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). $\endgroup$
    – Maximilian Press
    Commented Feb 6, 2021 at 14:57

1 Answer 1

Reset to default
1
$\begingroup$

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

Improve this answer
$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.