TopHat v1.2.0 sort header

[Jon_Keats]

I never used Tophat v1.1.1 which listed a fix for the sam sort header (see below) but in newest version, TopHat v1.2.0 all my sam files have sort headers of "sorted" not "coordinate". Oddly, parsing the file through Samtools sort does not fix the problem but parsing it through Picard does. Also the sam headers are not listed in numeric order:

Such as:
chr1
chr11
chr12
...
chr2
chr20

Not:
chr1
chr2
chr3
chr4


Anyone else seeing these minor issues?


----Previous release notes-----

TopHat 1.1.1 release 10/11/2010

This release of TopHat includes some fixes related to Colorspace read mapping.

* Negative quality values are now handled correctly.
* Comments at the beginning of csfasta files no longer trigger an error.
* --integer-quals no longer conflicts with -i
* The header in TopHat BAM files now correctly lists the sort order as coordinate, with group order reference

[colindaven]

Picard seems to do a better job than Samtools of putting whether the BAM has been sorted in the SAM header. I've been working with a group of computer scientists who picked up on this one, so I have changed from Samtools to Picard for SAM/BAM conversion and sorting.

[cjp]

Hi Jon,

You can use picard ReorderSam:

http://picard.sourceforge.net/command-line-overview.shtml#ReorderSam

http://www.broadinstitute.org/gsa/wiki/index.php/ReorderSam

First, TopHat gives the wrong sort order in the header, so you'll have to change that else picard will complain.

e.g.,

samtools view -H acc.bam | sed 's/sorted/unsorted/' > acc.header.sam
samtools reheader acc.header.sam acc.bam > acc_head.bam

This is the picard command that works for me (it uses the order of sequences in the reference file in the output BAM file):

java -jar /path/to/picard/jars/ReorderSam.jar I=acc_head.bam O=acc_order.bam R=/path/to/ref/human_g1k_v37.fasta

Then you may have to re-sort the BAM file. Although, if you trust TopHat's sorting, I guess you can change the sed line above to: sed 's/sorted/coordinate/'. For me, I like to add read group info anyway as a lot of software like GATK needs them to run, so I use picard AddOrReplaceReadGroups, which also allows you to sort with the SO option:

e.g.,

name="MY_SAMPLE"

java -jar /path/to/picard/jars/AddOrReplaceReadGroups.jar I=acc_order.bam O=acc_rg.bam RGID=$name RGLB=$name RGPL=ILLUMINA RGPU=$name RGSM=$name SO=coordinate

There are other options I use with picard:

TMP_DIR=/path/to/tmp VALIDATION_STRINGENCY=SILENT VERBOSITY=ERROR QUIET=true CREATE_INDEX=true

Chris

[cedance]

The post is quite old and newer versions of tophat (since 1.3.0 I guess), with collaboration from picard developers, have overcome these issues and also SAM format TLEN parameter etc...
Its better to use 1.3.1 (1.3.2 is out but still in beta) in my opinion.

[cjp]

The post is quite old and newer versions of tophat (since 1.3.0 I guess), with collaboration from picard developers, have overcome these issues and also SAM format TLEN parameter etc...
Its better to use 1.3.1 (1.3.2 is out but still in beta) in my opinion.

I still use old versions of TopHat for short reads because in the TopHat page it says this about version 1.3:

"For short reads (usually <45-bp), it is recommended that users decrease segment length (--segment-length) to about half the read length and segment mismatches (--segment-mismatches) to 0 or 1"

When I ran it on 36bp data, it was necessary to play with these settings and I got different results than I did with TopHat 1.2 - reads didn't align across splice sites and aligned in different places or across different splice sites. In the help pages, I couldn't find an explanation of why the new version needed these new parameter changes but they weren't needed in older versions of TopHat.

On long read data, I think TopHat 1.3 seems to work well.

[cedance]

Din't know that. Thanks for letting me know.Yes, that makes total sense. Fortunately, I work on 80bp paired end reads. The problem I faced with Tophat 1.2.0 is that the column 9 of SAM format = TLEN was 0 always. I would like to know the entire fragment length that's mapped.

Best,
Arun.

[cjp]

That's true - you can use picard FixMateInformation, but I don't know how well it works with reads that align overs introns:

http://picard.sourceforge.net/command-line-overview.shtml#FixMateInformation

Chris

[cjp]

Are you sure it is fixed, I just found one of my TopHat 1.3 files (the @PG line says so anyway) and it looks like this in the header (chromosomes are still in the order 1,10,11):

@HD VN:1.0 SO:coordinate
@SQ SN:1 LN:249250621
@SQ SN:10 LN:135534747
@SQ SN:11 LN:135006516
...
@SQ SN:GL000247.1 LN:36422
@SQ SN:GL000248.1 LN:39786
@SQ SN:GL000249.1 LN:38502
@SQ SN:MT LN:16569
@SQ SN:X LN:155270560
@SQ SN:Y LN:59373566
@PG ID:TopHat VN:1.3.1 CL:/home/cjp64/src/tophat-1.3.1/src/tophat -p 12 --segment-length 15 --segment-mismatches 0 -o A37_2_west -G /home/easih/gtf/hg19_ccds_08022011.gtf /home/easih/refs/human_1kg/bowtie/human_g1k_v37 /scratch/svvd2/A37/A3700002.1.f

[cedance]

cjp, its fixed in Tophat 1.3.1

TLEN field in SAM format is correctly output

I was talking about this (bold area):

5_Solexa_0503:5:75:14816:7572#0 99 SL2.40ch01 17202 255 75M = 17281 159 CGGCCGCACAGTTATTCGTGATGTCGCCATCGGATGTGGCCATAGTAATCACGGTATGTTTATTGGGGCTGCCGG CCCCCCCCCCCCCCCCCCCCCCCDCCCCCCCCCACCBCCCC@C@C@DCDCABBC?BA=@CC<@C=BBBDB@:@8@ NH:i:1 NM:i:2
5_Solexa_0503:5:75:14816:7572#0 147 SL2.40ch01 17281 255 80M = 17202 -159 TTGGGTCTTGGAGGAGGCTCTATGTCACTTGTTGGACAACTCGGTGGACAAACAGGTGGAGCCTTTAGTTACTGTTTGGA >DCCB@@?@DCDDCDBCDCACDCDAC?A=CCCCCCACCDCCACCCCCDCCCDCCCCCCCCCCCCCCCCCCCCCCCCCCCC NH:i:1 NM:i:2

[cjp]

The original post says this:

I never used Tophat v1.1.1 which listed a fix for the sam sort header (see below) but in newest version, TopHat v1.2.0 all my sam files have sort headers of "sorted" not "coordinate". Oddly, parsing the file through Samtools sort does not fix the problem but parsing it through Picard does. Also the sam headers are not listed in numeric order:

Such as:
chr1
chr11
chr12
...
chr2
chr20

Not:
chr1
chr2
chr3
chr4


Anyone else seeing these minor issues?

...

So this was why I mentioned ReorderSam. The TLEN thing is a separate issue, but I thought you said TopHat 1.3 fixed the issue of the header having the wrong chromosome order, but it doesn't seem to in my file.

[cjp]

Here is an example of short reads being aligned differently by TopHat 1.1, 1.2 and 1.3 (even though I set the segment values the same as mentioned in the TopHat home page):



Chris