This is a test version of Biostars. For the public version, visit https://www.biostars.org.
Alignment produces BAM file with sorted reads, but cannot see alignment in viewer programs.

Hi all,

This is my first time performing NGS alignment. I am using command line tools. What I did so far, was to first use cat to combine all chromosomes of my reference genome. I then indexed the genome with bwa (-a bwtsw). I then de-multiplexed and sorted into separate FASTQ files by barcode using FASTX Toolkit's Barcode Splitter. Next, I used FASTX trimmer to remove 75bp of barcode, primer, and transposon sequence (I am over trimming, but I started with 300 bases). I then aligned with bwa, trying both the mem and samse algorithms. Then, with samtools, I converted to BAM format with the view command, sorted, and indexed. In viewer programs such as Tablet, IGV, and BamView, I either saw no reads/blank or, in the case of Tablet, I saw that it had categorized reads by chromosome, but clicking on a chromosome showed nothing. This is true for the output of both algorithms. Not sure if it is worth noting, but in Tablet, the mismatches column for each chromosome all show a "?". What could go wrong along the way?

Thanks for your help.

alignment next-gen sequencing
• 3,849 views • link
written by ac • ac

I think it may be worth noting that I just found out that I can view the reads if I open the SAM file, but not the BAM file. Maybe something went wrong during the conversion?

• 0 views • link
written by ac • ac

It might be helpful to see the actual commands you executed to convert SAM to BAM.

• 0 views • link
written by Kamil Kamil

This is the code I used:

samtools view -u -S alignment.sam > alignment.bam
samtools sort alignment.bam alignment.bam.sorted
samtools index alignment.bam.sorted.bam
• 0 views • link
written by ac • ac

1 answer

I might guess that your genome browser expects chromosome names like "chr1" and your BAM file has chromosome names like "1", or vice versa. You'll want to make sure that the genome browser is using exactly the same reference sequence as your read mapper.

Check the chromosome names in the FASTA file that is used by your genome browser:

head -n1 genome.fa

Also, check the chromosome names in your BAM file header:

samtools view -H file.bam

And column 3 in your BAM file:

samtools view file.bam | head
• 1 views • link
written by Kamil Kamil

Here are the results:

head -n1 genome.fas
>Chr1 CHROMOSOME dumped from ADB: Jun/20/09 14:53; last updated: 2009-02-02

samtools view -H alignmentfile.bam
@SQ    SN:Chr1    LN:30427671
@SQ    SN:Chr2    LN:19698289
@SQ    SN:Chr3    LN:23459830
@SQ    SN:Chr4    LN:18585056
@SQ    SN:Chr5    LN:26975502
@SQ    SN:chloroplast    LN:154478
@SQ    SN:mitochondria    LN:366924
samtools view -H alignmentfile.bam.sorted.bam
@HD    VN:1.3    SO:coordinate
@SQ    SN:Chr1    LN:30427671
@SQ    SN:Chr2    LN:19698289
@SQ    SN:Chr3    LN:23459830
@SQ    SN:Chr4    LN:18585056
@SQ    SN:Chr5    LN:26975502
@SQ    SN:chloroplast    LN:154478
@SQ    SN:mitochondria    LN:366924
samtools view file.bam | head

Appears to give what I think are 10 or so alignments? I truncated the rest of the output.

M01440:92:000000000-AEGP7:1:1101:9444:1247    16    Chr5    5962715    32    50S176M    *    0    0    [...]
• 0 views • link
written by ac • ac

Log in to answer this question.

More posts like this