[Aman Mahajan]

SAM Format

I have converted my alignment file to SAM format using samtools. I am not sure if it is correct. It looks like this. Any idea about it?..Thanks

HWI-ST141_0380:5:1101:1213:2108#TGACCA 65 199172 35 30 50M * 0 0 CGTCTATAAACCTCTCACATGCTGGTGGATACTTGGTCATGGAAAAATCC hhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhh NM:i:2 MD:Z:40C1C7
HWI-ST141_0380:5:1101:1176:2150#TGACCA 81 57103 30 30 50M * 0 0 TCTCCGTCGTCGCAGACGAGTTCATGGCCAGGGACATCACATTCGAGAAC

[son_nexg]

Hi I am experiencing problems using .sam output by Tophat with Cufflink.
My first question is : What options should I use to generate a .bam file with Tophat ?

If thats not possible, could someone please help me fix the error I am getting when I try to run Cufflink on the .sam file:

You are using Cufflinks v1.0.3, which is the most recent release.
[bam_header_read] EOF marker is absent.
[bam_header_read] invalid BAM binary header (this is not a BAM file).
File Sample_1001645_9/accepted_hits.sam doesn't appear to be a valid BAM file, trying SAM...
[15:21:47] Loading reference annotation.
[15:21:51] Inspecting reads and determining fragment length distribution.
> Processing Locus chrY:59330251-59343488 [************************ ] 97%
Error: this SAM file doesn't appear to be correctly sorted!
current hit is at chrX:145287, last one was at chrM:16469
Cufflinks requires that if your file has SQ records in
the SAM header that they appear in the same order as the chromosomes names
in the alignments.
If there are no SQ records in the header, or if the header is missing,
the alignments must be sorted lexicographically by chromsome
name and by position.




I also tried to sort the sam output using samtools but that also gives and error that '@SQ header is missing' hence can not sort it using samtools!!

Thank you.

[ben.weisburd]

I was seeing this error also ("Error: Segment join failed with err = 1").
Turned out the problem was that I was running more than one tophat process simultaneously in the same directory (I wanted to process multiple fastq files in parallel). When I launched the processes in separate directories, everything worked fine.

-Ben

I just downloaded the latest Tophat-1.0.10, and found there is an segment join failed error at the last.
And Tophat cannot generate the output files.
But the previous version did not have this error.
Is this a new incompatibility or anything else?

[Fri Jul 31 13:58:04 2009] Beginning TopHat run (v1.0.10)
-----------------------------------------------
[Fri Jul 31 13:58:04 2009] Preparing output location ./tophat_out/
[Fri Jul 31 13:58:04 2009] Checking for Bowtie index files
[Fri Jul 31 13:58:04 2009] Checking for reference FASTA file
[Fri Jul 31 13:58:04 2009] Checking for Bowtie
Bowtie version: 0.10.0.0
[Fri Jul 31 13:58:04 2009] Checking reads
seed length: 50bp
format: fastq
quality scale: --phred33-quals
[Fri Jul 31 13:59:06 2009] Reading known junctions from GFF file
Splitting reads into 2 segments
[Fri Jul 31 14:00:21 2009] Mapping reads against h_sapiens_asm with Bowtie
[Fri Jul 31 14:05:34 2009] Mapping reads against h_sapiens_asm with Bowtie
[Fri Jul 31 14:10:42 2009] Searching for junctions via segment mapping
[Fri Jul 31 14:22:10 2009] Retrieving sequences for splices
[Fri Jul 31 14:24:18 2009] Indexing splices
Warning: Encountered reference sequence with only gaps
Warning: Encountered reference sequence with only gaps
.....

[Fri Jul 31 14:24:41 2009] Mapping reads against segment_juncs with Bowtie
[Fri Jul 31 14:26:35 2009] Mapping reads against segment_juncs with Bowtie
[Fri Jul 31 14:28:30 2009] Joining segment hits
[FAILED]
Error: Segment join failed with err = 1

[b1nu]

My bad. Hadn't specified the '-S' option to samtools to indicate that input is in SAM format. Works perfectly now.

I am using tophat 1.0.13 & bowtie 0.12.2, but still get the error
"[bam_pileup] fail to read the header: non-exisiting file or wrong format."

The tophat output does mark the quality scale as 'phred33'.

[Tue Feb 23 13:18:51 2010] Checking reads
	seed length:	 42bp
	format:		 fastq
	quality scale:	 phred33 (default)

[b1nu]

I am using tophat 1.0.13 & bowtie 0.12.2, but still get the error
"[bam_pileup] fail to read the header: non-exisiting file or wrong format."

The tophat output does mark the quality scale as 'phred33'.

[Tue Feb 23 13:18:51 2010] Checking reads
	seed length:	 42bp
	format:		 fastq
	quality scale:	 phred33 (default)

[Cole Trapnell]

This is a bug in 1.0.12, and is fixed in the upcoming 1.0.13. 1.0.12 was simply propogating the unconverted quality scores, where 1.0.13 converts them to Phred-33, as required by the SAM spec.

[kch]

I'm using version 1.0.10 with the "--solexa1.3-quals" flag...can you (or anyone) confirm that Tophat outputs Phred-scaled (ASCII -33 offset) base quality scores (as described in the SAM specifications)? I seem to be getting impossibly high base qualities which could be contributing to my problems in SNP filtering.

Thanks.

I am using TopHat version 1.0.12 and Bowtie version 0.11.3 with the "--solexa1.3-quals" flag as well. Similarly, I'm not sure whether the output in my accepted_hits.sam file is reporting ASCII-33 quality values. For example, in the following output:

EAS306_9079:7:69:745:1129 0 chr1 1442 1 50M * 0 0 CAGAATTGTACTGTTCTGTATCCCACCAGCAATGTCTAGGAATACCTGTT abbbbbbaabbbaababaabbbbbbbaaaabbbaaabaaaaaaaaaaaaa NM:i:1

The quality string is the same as my Illumina fastq input file and I expected the quality string to be "BCCCCCCBBCCCBBCBCBBCCCCCCCBBBBCCCBBBCBBBBBBBBBBBBB" (as converted using `maq ill2sanger`).

Any input would be appreciated, thanks!

[choy]

TopHat v1.0.12 SAM output question

I am running TopHat v1.0.12 with 40mer paired-end reads, and am having difficulty interpreting the SAM output file (accepted_hits.sam). Here is an example of a pair of output lines that correspond to the mates of a single paired-end read:

PATHBIO-SOLEXA2:2:58:734:942#0 145 chr1 554416 255 40M = 3836 0 ATAATACACACCCTCACCACTACAATCTTCCTAGGAACAA _QS_\^H^`W^Z[Ua^``M`^Jb_UNVZ\a`]baaaaV`K NM:i:1

PATHBIO-SOLEXA2:2:58:734:942#0 97 chrM 3836 255 40M = 554416 0 CCATCATGAACCTTGGCCATAATATGATTTATCTCCACAC `bbbbbbaVX`]aa`Y\X`Y]S_a\a_^XZ_bb`bbaT`` NM:i:1

My question is regarding the MNRM field (the 7th field, set to '=' here). The SAM format specification states that this field is set to the reference of the mate, and will only be set to '=' if the reference of the query sequence is the same as the reference of the mate sequence. In this case the mates map to different contigs, so why is this field set to '='?

When computing coverage (as in the coverage.wig file), do you eliminate mate pairs that have an MPOS field (the 8th field) that does not correspond to its mate?

Also, would it be possible for you to set the ISIZE parameter in the accepted_hits.sam file?

Thanks very much!

[griffon42]

Great! Looking forward to the update.

I can still email you sample data if necessary...just let me know.

Thanks again.

[Cole Trapnell]

Ah, I see where I'm going wrong. I'll have this fixed shortly.

[Cole Trapnell]

It should be, although TopHat lets Bowtie handle the actual conversion to base 33 Phred scale. Can you email me a handful of reads that wind up with bad qualities after the mapping?

[griffon42]

Thanks Cole.

Sorry for the sequential posting, but...Taking a closer look, it looks like the issue could be the base qualities returned in the 'accepted_hits.sam' output from Tophat.

I'm using version 1.0.10 with the "--solexa1.3-quals" flag...can you (or anyone) confirm that Tophat outputs Phred-scaled (ASCII -33 offset) base quality scores (as described in the SAM specifications)? I seem to be getting impossibly high base qualities which could be contributing to my problems in SNP filtering.

Thanks.

[Cole Trapnell]

I'll see what I can do for the next release

[griffon42]

Is there any way to have Tophat output an appropriate value for the "Mapping Qualities" field in the SAM format? At present, it outputs a straight "255" (which the SAM specification suggests when there is no value available), but this seems to be interfering with SNP calling in SAMtools (ie varFilter doesn't 'filter' much of anything...I end up with LOTS of false positives).

Alternatively, can anyone comment on an appropriate method of calling quality-aware, reliable cSNPs from Tophat output without a 'mapping quality' value?

Thanks!

[griffon42]

Tophat 1.0.10

Hi Cole-

The new version (1.0.10) seems to be working fine for me. Strand orientation in the output SAM file looks good, as far as I can tell.

Thanks so much for the rapid turnaround!