Entering edit mode
10.4 years ago
lkmklsmn
▴
980
Hi everyone,
I am trying to convert my paired-end RNAseq bam file into a bed file at the FRAGMENT level at a specific region of interest. To do so I believe I have to:
1. Extract reads form region of interest
samtools view -bf 0x2 file.bam chr22:50713407-50720510 > tmp.bam
2. Sort by name
samtools sort -n tmp.bam tmp_namesorted
3. Convert to bed file
bamToBed -cigar -bedpe -i tmp_namesorted.bam
But the output contains an error message:
chr22 50719392 50719581 chr22 50719529 50719604 HWI-ST1168:5:1101:2456:98694#0 255 + -
*****ERROR: -bedpe requires BAM to be sorted/grouped by query name. chr22 50714331 50715071 chr22 50715106 50715317 HWI-ST1168:5:1101:2720:13083#0 255 + -
chr22 50716598 50716673 chr22 50716883 50717034 HWI-ST1168:5:1101:2742:47756#0 255 + -
chr22 50719086 50719234 chr22 50719320 50719395 HWI-ST1168:5:1101:3091:54779#0 255 + -
chr22 50719599 50719844 chr22 50720067 50720142 HWI-ST1168:5:1101:3307:18506#0 255 + -
chr22 50719592 50719837 chr22 50719866 50720007 HWI-ST1168:5:1101:3352:29638#0 255 + -
chr22 50714319 50714393 chr22 50716357 50716432 HWI-ST1168:5:1101:3618:38910#0 255 + -
chr22 50713898 50713969 chr22 50714054 50714129 HWI-ST1168:5:1101:3740:40497#0 255 + -
chr22 50719599 50719841 chr22 50719897 50720036 HWI-ST1168:5:1101:3795:18405#0 255 + -
chr22 50713927 50714002 chr22 50714136 50714211 HWI-ST1168:5:1101:4198:53931#0 255 + -
What am I doing wrong here?
Thanks
Your bam does not look sorted.
If you look at the SAM version of the sorted file do you have pairs (2 lines) of reads with the same name?