I'm working on a project where I am analyzing the performance of an alignment workflow. My goal is to find regions in the resulting BAM file where there are outstanding discrepancies or anything that indicates my assembly/alignment has "mistakes". I'm working with human genomes from the 1000 genomes project, they are high coverage Illumina paired end reads, each read is 100bp.

My workflow so far:

  • input paired end FASTQ files (human) into SPAdes to create contigs.fasta file
  • align paired end FASTQ files to contigs.fasta file with minimap2 to get SAM file
  • convert SAM to BAM with samtools
  • sort and index BAM file
  • load BAM file into IGV and load contigs.fasta as the reference

In the IGV GUI, I have colored alignments by insert size and pair orientation. This has revealed a bunch of insertions, tandem duplications, and inversions.

I am wondering:

  • Do these colored reads indicate that the assembler made "mistakes" or could they be there just by chance?
  • Do SNPs indicate anything about how SPAdes/minimap2/samtools performed?
  • Is there a better way to sort/color the alignments based on what I'm trying to do?

Thank you so much to anyone who has any advice!

EDIT: here is a screenshot of the alignment. For clarification, there are 8 genomes of interest, hence 8 tracks


Background Integrative genome viewer (IGV) is a well known Broad Institute tool and available here.

The most recent papers are: James T. Robinson, Helga Thorvaldsdóttir, Aaron M. Wenger, Ahmet Zehir, Jill P. Mesirov; Variant Review with the Integrative Genomics Viewer. Cancer Res 1 November 2017; 77 (21): e31–e34. https://doi.org/10.1158/0008-5472.CAN-17-0337

igv.js: an embeddable JavaScript implementation of the Integrative Genomics Viewer (IGV) James T. Robinson, Helga Thorvaldsdóttir, Douglass Turner, Jill P. Mesirov bioRxiv 2020.05.03.075499; doi: https://doi.org/10.1101/2020.05.03.075499

  • $\begingroup$ seeing some browser shots would help. In general, a few weird reads is no big deal. A lot of reads showing systematic issues might be more of a problem. It would also help to know what kind of reads. We trust Illumina more for base-level accuracy, and long reads more for contiguity. What species? $\endgroup$ Nov 18, 2021 at 14:37
  • $\begingroup$ It's a human genome from the 1000 genomes project! And thank you, I just added a photo for clarification $\endgroup$ Nov 18, 2021 at 16:02
  • $\begingroup$ thanks- can you share some of the track labels? I assume some of that stuff is gene tracks (e.g. green bars), but there are some coverage plots that don't show the read pileups, which makes me worry that I don't understand what's going on. $\endgroup$ Nov 18, 2021 at 17:32
  • 2
    $\begingroup$ perhaps more to the point, there are 3 pertinent pieces of information you can get from this approach: (1) aligning reads against GRCh38, (2) reads against your assemblies, (3) your assemblies against GRCh38 (off the top of my head). I would suggest doing all 3 of these approaches, and comparing them. 1 and 3 can be visualized at the same time, which might be particularly helpful in finding e.g. regions of ambiguity due to repetitive regions. You appear to be doing 2, which might be a good way to identify regions to check by 1 and 3. $\endgroup$ Nov 18, 2021 at 17:51
  • $\begingroup$ What do you mean by track labels? this site makes it seem like it's just the name of the file I'm using. In this case all my files are from de novo assemblies, I only got the FASTQ files off the internet $\endgroup$ Nov 18, 2021 at 19:45


Your Answer

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

Browse other questions tagged or ask your own question.