SEQanswers

Go Back   SEQanswers > Applications Forums > RNA Sequencing



Similar Threads
Thread Thread Starter Forum Replies Last Post
sorting bam files frymor Bioinformatics 23 02-10-2016 05:46 PM
Cufflinks refuses to operate on Tophat2 created bam or sam files due to sorting error amrezans Bioinformatics 1 06-24-2013 12:54 PM
sorting sam file crh Bioinformatics 2 06-16-2011 06:45 AM
Sorting SAM output from Bowtie DrD2009 Bioinformatics 9 11-10-2010 11:52 AM
Sorting large files scami Bioinformatics 3 09-20-2010 11:45 PM

Reply
 
Thread Tools
Old 06-16-2016, 09:28 AM   #1
ronaldrcutler
Member
 
Location: Virginia

Join Date: May 2016
Posts: 80
Default Sorting SAM files

Hello all. I've been fiddling around with htseq-count using SAM files that were output by hisat2. I am using paired reads here and realized as I was going through the Htseq-count Manual that I should sort my data beforehand using
Code:
samtools sort
I used this command to execute htseq-count:
Code:
htseq-count -m intersection-nonempty -i Name -s reverse -r name Sample_1_hisat2_results_sorted.sam /Volumes/cachannel/RNA_SEQ/Notch_RNASeq/9.1_Reference_Files/XENLA_UTAmayball_cdna_longest_CHRS2.gff3 >Sample_1_Hisat2_Counts.txt 2>Sample_1_Hisat2_Counts_OUTPUT_WARNINGS.txt
So now I have two htseq-count output files to compare: one where the data is sorted and one where the data is unsorted.

In the unsorted warning file I notice at the bottom: "4674733 SAM alignment pairs processed"
In the sorted warning file I notice at the bottom: "8572367 SAM alignment paris processed"

That is about 2X the amount of pairs processed for the sorted data! The count txt files for both runs look good, but I'm guessing there is much more data on counts in the sorted one. When comparing counts for specific genes between the two, the sorted appears to have around double the amount. My question is why is this and what are the consequences of sorting? Does it lead to more precise counts, or is it not even worth it?

Last edited by ronaldrcutler; 06-16-2016 at 09:36 AM.
ronaldrcutler is offline   Reply With Quote
Old 06-16-2016, 09:38 AM   #2
ronaldrcutler
Member
 
Location: Virginia

Join Date: May 2016
Posts: 80
Default

Also, another related question would be: is the SAM alignments output by hisat2 sorted?
ronaldrcutler is offline   Reply With Quote
Old 06-16-2016, 09:57 AM   #3
dpryan
Devon Ryan
 
Location: Freiburg, Germany

Join Date: Jul 2011
Posts: 3,479
Default

htseq-count wants you to "samtools sort -n", not "samtools sort". The difference is the cause of the differing results. You do not need to sort the output of hisat2 before giving it to htseq-count.

Note that since you coordinate sorted the file and then told htseq-count that it was name sorted that the results for that are...inaccurate. The file with the smaller number of processed alignments is the correct one.
dpryan is offline   Reply With Quote
Old 06-16-2016, 10:01 AM   #4
ronaldrcutler
Member
 
Location: Virginia

Join Date: May 2016
Posts: 80
Default

Thanks for the clarification, this will save a lot of time!
ronaldrcutler is offline   Reply With Quote
Old 06-28-2016, 09:51 AM   #5
ronaldrcutler
Member
 
Location: Virginia

Join Date: May 2016
Posts: 80
Default

Quote:
Originally Posted by dpryan View Post
You do not need to sort the output of hisat2 before giving it to htseq-count.
When examining the head of some SAM files I have been working with output from hisat2, I noticed that the head contains this line:
Code:
@HD	VN:1.0	SO:unsorted
I know you said hisat2 outputs sorted SAM files, so what does this mean?
ronaldrcutler is offline   Reply With Quote
Old 06-28-2016, 05:28 PM   #6
GenoMax
Senior Member
 
Location: East Coast USA

Join Date: Feb 2008
Posts: 6,780
Default

Quote:
Originally Posted by ronaldrcutler View Post
When examining the head of some SAM files I have been working with output from hisat2, I noticed that the head contains this line:
Code:
@HD	VN:1.0	SO:unsorted
I know you said hisat2 outputs sorted SAM files, so what does this mean?
You could use instead featureCounts. It is much faster and will sort the bam/sam files if needed.

Looks like HISAT2's output is unsorted.

Last edited by GenoMax; 06-28-2016 at 05:31 PM.
GenoMax is online now   Reply With Quote
Old 07-01-2016, 09:08 AM   #7
ronaldrcutler
Member
 
Location: Virginia

Join Date: May 2016
Posts: 80
Default

To follow up: sorting the sam files removed this error that I had in all of them:
Code:
Warning: Malformed SAM line: MRNM != '*' although flag bit &0x0008 set
Warning: Malformed SAM line: RNAME != '*' although flag bit &0x0004 set
Warning: Malformed SAM line: MRNM == '=' although read is not aligned.
But not this error, which was similar in all of them (however, I just ignored it):
Code:
Warning: Read ACB052:253:C76YKACXX:2:1101:2245:1957 claims to have an aligned mate which could not be found in an adjacent line.
When comparing the sorted and unsorted files using the 'diff' command, there were no differences!
ronaldrcutler is offline   Reply With Quote
Reply

Tags
hisat2, htseq count, warnings

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 03:01 AM.


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