RNA-Seq Analysis: To sort or not to Sort?

Question

[kanika.151]

Hello,

I am trying my hand at RNA-Seq Analysis by going TopHat-> HTSeq-> edgeR

After TopHat, conversion of bam file to sam is recommended then to sort is recommended.

This is how I converted:

#convert bam file to sam file
samtools view -h -o out.sam in.bam

and this is how I sorted:

sort -s -k my_file.sam > my_file_sorted.sam

Although, my sorted.sam file kept on giving me errors while running through HTSeq such as error when reading sam/bam file raised in count.py:84 for another file it would say 'seq' and 'qualstr' do not have the same length.

My HTSeq versions are uptodate. Also, I have PySam installed.

But, when I ran unsorted bam files or sam files through HTSeq then they were getting processed.

I want to know the significance of sorting a file and also why is HTSeq processing unsorted files?

my HTSeq command is as follows:

nohup time -p python -m HTSeq.scripts.count -f bam -s yes --idattr=ID hits.bam anno.gff3 &>mylog &

The only time I got an error for unsorted file was when they had separate pair end and single end reads in one file which separated by using

samtools view -bf 1 foo.bam > pair.bam
samtools view -bF 1 foo.bam > single.bam

Answer 1

[Sam]

Why don't you use samtools sort directly?

Something like samtools sort sam > sort.sam

[kanika.151]

Okay. I will try that too. Can it be used to sort bam files too?

[mkulecka]

Yes, it can. By the way, did you check if maybe your original files are already sorted? Some programmes sort files for you.

[kanika.151]

All ready sorted? as in? we ran through NGS QC then through TopHat...I don't remember sorting them...

How do I check for it?

[mkulecka]

cat your_file.sam | grep SO

If you get something like "SO:coordinate" then your files are already sorted.

With bam files:

samtools view -H your_bam_file.bam | grep SO

[kanika.151]

How does samtools work?

the command I used was:

samtools sort -n in.bam sort.sam

but it is giving me separate files in output:

sort.sam.001.bam
sort.sam.002.bam ...

[Carlo Yague]

Hi !

You just need to wait. samtools sort first creates several sorted files that are merged in one final file at the end of the process. It can take some time.

[kanika.151]

Yes it did! Thanks :)

[GB]

Hi Carlo Yague, would you know how long samtools sort usually take? I feel like it would depend on how big my bam file is (437 GB), but is >3 days normal? Should I not be worried as long as the number of temp files is increasing/ there are temp output files? Thank you!

[Carlo Yague]

Hi, to be honest, I never tried to sort a file that large, but yeah, the larger the file, the longer it takes... In such a case, it becomes especially important to optimize the -@ and -m parameters of samtools (number of threads and available memory per thread). Hope this helps. Carlo.

Answer 2

[h.mon]

To properly account for paired reads, htseq-count expects them to be sorted, either by name (-r name parameter, the default for htseq-count) or by position (-r pos). There is a bug with -r pos which causes htseq-count to crash, I think it is still unfixed.

I think htseq-count does not check if the reads are sorted, if you run it on an unsorted file (or sorted by position), it will not account for read pairing and you will get approximately double counts per feature - a bit less, because sometimes just one read of the pair maps. When I did this mistake, I got between 1.8x and 1.9x more counts per gene, on average.

P.S. edit: for a RNAseq differential expression analysis, you should consider Subread+featureCounts, or STAR with --quantMode GeneCounts, both replacing TopHat+HTSeq, then proceeding to edgeR as before.