Error in samtools uniquely mapped reads file

Question

[Whoknows]

Hi

I used samtools to extract uniquely mapped reads from a Tophat accepted_hits.bam file by command :

samtools view q 10 accepted_hits.bam > UniqueHits.bam

but during Cuffdiff procedure it makes this error

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.

1. Is my samtools extraction command makes Sam file in output?? And if yes

2. Can I use below command for creating bam output?? (and then Sorted bam from that)

   samtools view -b -q 10 accepted_hits.bam > Unique.bam
   samtools sort Unique.bam  Unique.sorted
   

Thanks

Answer 1

[Devon Ryan]

1. No, you meant samtools view -bq 10 accepted_hits.bam > UniqueHits.bam. However, don't do this! Cuffdiff is designed explicitly to use ambiguous alignments.
2. Unique.bam is already sorted.

[Whoknows]

My accepted_hits.bam size is 2.2GB but Unique.bam is 8.6 GB so I'm not sure it could be correct,

If Unique.bam is sorted, why cuffdiff makes error?

[Devon Ryan]

Without the -b flag, that's a SAM file. Tophat defaults to sorting alignments, so unless you told it not to do that then accepted_hits.bam is already sorted (and filtering a sorted file with samtools will always produce a sorted file).

[Whoknows]

I'm not sure maybe it because of my genome sequence headers, my genome is Wheat and they might have long header for each sequence

Answer 2

[Manvendra Singh]

Try this,

samtools view -h accepted_hits.bam | grep -E '^@|NH:i:1' | samtools view -Sb - > uniquely_mapped.bam

I think that it should help