Question: RNA-Seq Analysis: To sort or not to Sort?
0
gravatar for kanika.151
3.5 years ago by
kanika.15160
United States
kanika.15160 wrote:

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

 

sort rna-seq sam myposts htseq • 2.9k views
ADD COMMENTlink modified 3.5 years ago by h.mon25k • written 3.5 years ago by kanika.15160
2
gravatar for Sam
3.5 years ago by
Sam2.3k
New York
Sam2.3k wrote:

Why don't you use samtools sort directly?

Something like samtools sort sam > sort.sam

ADD COMMENTlink written 3.5 years ago by Sam2.3k

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

ADD REPLYlink written 3.5 years ago by kanika.15160

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

ADD REPLYlink written 3.5 years ago by mkulecka300

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?

ADD REPLYlink written 3.5 years ago by kanika.15160
2
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
ADD REPLYlink written 3.5 years ago by mkulecka300

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 ...

ADD REPLYlink written 3.5 years ago by kanika.15160
1

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.

ADD REPLYlink modified 3.5 years ago • written 3.5 years ago by Carlo Yague4.4k

Yes it did! Thanks :)

ADD REPLYlink written 3.5 years ago by kanika.15160
1
gravatar for h.mon
3.5 years ago by
h.mon25k
Brazil
h.mon25k wrote:

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.

ADD COMMENTlink modified 3.5 years ago • written 3.5 years ago by h.mon25k
Please log in to add an answer.

Help
Access

Use of this site constitutes acceptance of our User Agreement and Privacy Policy.
Powered by Biostar version 2.3.0
Traffic: 1799 users visited in the last hour