SEQanswers

Go Back   SEQanswers > Bioinformatics > Bioinformatics



Similar Threads
Thread Thread Starter Forum Replies Last Post
Samtools error: the alignment was not sorted Desiree Wilson Bioinformatics 8 03-28-2012 02:02 PM
Fastest way to extract differing positions from each alignment in a BAM file CHRYSES Bioinformatics 5 12-14-2011 11:28 AM
BWA produces odd alignment results dandyrilla Bioinformatics 2 11-27-2011 11:28 PM
BWA: specifying SAM/BAM file header fields before read alignment? nora Bioinformatics 3 12-04-2010 09:11 PM

Reply
 
Thread Tools
Old 04-13-2015, 10:19 AM   #1
alch
Junior Member
 
Location: CA

Join Date: Apr 2015
Posts: 6
Default Alignment produces BAM file with sorted reads, but cannot see alignment

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.
alch is offline   Reply With Quote
Old 04-13-2015, 10:29 AM   #2
GenoMax
Senior Member
 
Location: East Coast USA

Join Date: Feb 2008
Posts: 7,076
Default

Did you check to see if chromosome names of your reference and the names that are in your BAM files are still matching (after all that processing)? If you look at the headers of the BAM file that should show you the names there. e.g. chr1 needs to stay chr1, it can't become anything else. Hopefully you did not have spaces in the names.

Note: This may sound stupid but have you "zoomed in" far enough in IGV. You won't see the alignments until a certain point.

Last edited by GenoMax; 04-13-2015 at 10:31 AM.
GenoMax is offline   Reply With Quote
Old 04-13-2015, 11:12 AM   #3
Richard Finney
Senior Member
 
Location: bethesda

Join Date: Feb 2009
Posts: 700
Default

Can you dump the reads using "samtools view" ?

What is the command you used in samtools?

What does the bam header look like?
Richard Finney is offline   Reply With Quote
Old 04-13-2015, 04:28 PM   #4
alch
Junior Member
 
Location: CA

Join Date: Apr 2015
Posts: 6
Default

The commands I used in samtools are:
Code:
samtools view -u -S alignment.sam > alignment.bam
samtools sort alignment.bam alignment.bam.sorted
samtools index alignment.bam.sorted.bam
The headers appear to match:

This is the genome fastq file:
>Chr1 CHROMOSOME dumped from ADB: Jun/20/09 14:53; last updated: 2009-02-02
[Chr1 sequence goes here]
>Chr2 CHROMOSOME dumped from ADB: Jun/20/09 14:54; last updated: 2009-02-02
[Chr2 sequence goes here]
[etc.]

This are the headers I get for both the unsorted and sorted bam files:
Code:
samtools view -H alignment.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
Code:
samtools view -H alignment.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
It does appear that I can dump the reads using samtools view.


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?
alch is offline   Reply With Quote
Old 04-13-2015, 11:36 PM   #5
dpryan
Devon Ryan
 
Location: Freiburg, Germany

Join Date: Jul 2011
Posts: 3,480
Default

BAM is a compressed binary format, it's not human readable. Perhaps your problem is that "Chr1" should be "chr1" (and so on).
dpryan is offline   Reply With Quote
Old 04-14-2015, 10:13 AM   #6
alch
Junior Member
 
Location: CA

Join Date: Apr 2015
Posts: 6
Default

dpryan,
I got the headers from the BAM file through the "view" command in samtools. The chromosome names seem to match the reference exactly (e.g. Chr1). Will making them lowercase make a difference? Wouldn't I have to redo the alignment if I want to change to lowercase?
alch is offline   Reply With Quote
Old 04-14-2015, 10:18 AM   #7
GenoMax
Senior Member
 
Location: East Coast USA

Join Date: Feb 2008
Posts: 7,076
Default

I think "CHROMOSOME dumped from ADB: Jun/20/09 14:53; last updated: 2009-02-02" part in the fasta header is throwing IGV off.

Can you try with a genome file that has this part removed from chromosome names? Names should be just Chr1, Chr2 etc.
GenoMax is offline   Reply With Quote
Old 04-14-2015, 10:44 AM   #8
alch
Junior Member
 
Location: CA

Join Date: Apr 2015
Posts: 6
Default

I just deleted the ends from the fastq genome file to get >Chr1, >Chr2, etc. Reloading this in IGV didn't make a difference. Do I need to reindex/realign?
alch is offline   Reply With Quote
Old 04-14-2015, 10:46 AM   #9
dpryan
Devon Ryan
 
Location: Freiburg, Germany

Join Date: Jul 2011
Posts: 3,480
Default

You might have to reindex the fasta file (actually, you might have to just remake the genome in IGV/tablet/etc., but this is quick). You do not need to realign (you'd only need to do that if you changed the actual sequence).
dpryan is offline   Reply With Quote
Old 04-14-2015, 01:47 PM   #10
alch
Junior Member
 
Location: CA

Join Date: Apr 2015
Posts: 6
Default

Hi guys,
As it turns out, I realized that I can view the reads in the sorted BAM file. However, that's only if 1) I zoom in enough and 2) I know exactly where to scroll to within the chromosome. Is there a way to instantly locate reads? IGV and other viewers I've tried do not show any preview of where reads are and I'm just mindlessly scrolling until I find reads appear in the window.

Thanks.
alch is offline   Reply With Quote
Old 04-14-2015, 02:09 PM   #11
GenoMax
Senior Member
 
Location: East Coast USA

Join Date: Feb 2008
Posts: 7,076
Default

How many reads did you have and what was the alignment %? Do you know roughly what fold coverage to expect?

If this was a re-sequencing experiment you should have uniform coverage across the genome so reads should be easy to locate. If this was an RNA- or ChIP-seq experiment then you may have to find areas where there are genes/motifs that are enriched.
GenoMax is offline   Reply With Quote
Old 04-14-2015, 05:59 PM   #12
alch
Junior Member
 
Location: CA

Join Date: Apr 2015
Posts: 6
Default

About 42,000 reads were aligned. I am sequencing to find transposon tags. They should only appear about 1-2 times in the genome. Is there a better downstream analysis program to do what I need to do? I just need to know a summary of the alignments: so, where a particular sequence(s) aligned (e.g. chromosome), the # of reads of that particular sequence(s), and the sequence itself.
alch is offline   Reply With Quote
Reply

Thread Tools

Posting Rules
You may not post new threads
You may not post replies
You may not post attachments
You may not edit your posts

BB code is On
Smilies are On
[IMG] code is On
HTML code is Off




All times are GMT -8. The time now is 08:09 AM.


Powered by vBulletin® Version 3.8.9
Copyright ©2000 - 2020, vBulletin Solutions, Inc.
Single Sign On provided by vBSSO