Question: Error in samtools uniquely mapped reads file
0
gravatar for Whoknows
4.4 years ago by
Whoknows740
Tehran,Iran
Whoknows740 wrote:

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.

rna-seq samtools tophat • 1.9k views
ADD COMMENTlink modified 4.4 years ago by Manvendra Singh2.1k • written 4.4 years ago by Whoknows740
1
gravatar for Devon Ryan
4.4 years ago by
Devon Ryan91k
Freiburg, Germany
Devon Ryan91k wrote:
  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.
ADD COMMENTlink written 4.4 years ago by Devon Ryan91k

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?

ADD REPLYlink written 4.4 years ago by Whoknows740

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

ADD REPLYlink written 4.4 years ago by Devon Ryan91k

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

ADD REPLYlink written 4.4 years ago by Whoknows740
0
gravatar for Manvendra Singh
4.4 years ago by
Manvendra Singh2.1k
Berlin, Germany
Manvendra Singh2.1k wrote:

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

 

ADD COMMENTlink written 4.4 years ago by Manvendra Singh2.1k
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: 1474 users visited in the last hour