This is a test version of Biostars. For the public version, visit https://www.biostars.org.
Error in samtools uniquely mapped reads file

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

tophat samtools rna-seq

2 answers

  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.

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?

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

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

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

Log in to answer this question.