Compare BAM files
1
1
Entering edit mode
7.6 years ago
Cindy ▴ 70

I want to compare two bam files to see whether they are same, how many reads they overlapped, and how many reads are unique. bam1 file was got several years ago and then I convert it to bam1_r1_fq and bam1_r2_fq and then remap to the new reference genome and got new sam file. Then I sorted, added header etc and got the bam2 file. When I use

samtools view -f 64 -F 2304 bam1 | cut -f 1,3,4 | LC_ALL=C sort -t ' ' -k 1,1  > bam1.txt   
samtools view -f 64 -F 2304 bam2 | cut -f 1,3,4 | LC_ALL=C sort -t ' ' -k 1,1  > bam2.txt       
join -t ' ' -1 1 -2 1 bam1.txt bam2.txt > bam.comparaison

I got error like this :

join: bam1.txt:4077: is not sorted: HWI-D00222:175:H7C06ADXX:2:1101:1250:20054  MT  10679
join: bam2.txt:4077  is not sorted: HWI-D00222:175:H7C06ADXX:2:1101:1250:20054  MT  10679

Then I want to try "cmpbams" that @Pierre Lindenbaum wrote, and it showed

mkdir -p lib/com/github/samtools/htsjdk/2.9.1/ && curl -Lk  -o "lib/com/github/samtools/htsjdk/2.9.1/htsjdk-2.9.1.jar" "http://central.maven.org/maven2/com/github/samtools/htsjdk/2.9.1/htsjdk-2.9.1.jar"
/bin/bash: curl: command not found
maven.mk:102: recipe for target 'lib/com/github/samtools/htsjdk/2.9.1/htsjdk-2.9.1.jar' failed
make: *** [lib/com/github/samtools/htsjdk/2.9.1/htsjdk-2.9.1.jar] Error 127

How can I fix one of this problem or both? Thanks ahead.

sequencing • 2.6k views
ADD COMMENT
3
Entering edit mode
7.6 years ago
GenoMax 148k

Since many aligners work in non-deterministic mode the files will not be 100% identical (more so now due to different reference genomes, differences that may be there due to actual alignment logic, if you used different aligners/different versions of the same aligner etc).

It also appears that you are missing the program curl for your OS. Install it using appropriate means.

ADD COMMENT
0
Entering edit mode

Yes, I need to install other programs. It works now. Thanks

ADD REPLY

Login before adding your answer.

Traffic: 2210 users visited in the last hour
Help About
FAQ
Access RSS
API
Stats

Use of this site constitutes acceptance of our User Agreement and Privacy Policy.

Powered by the version 2.3.6