/
binBamReads
executable file
·46 lines (39 loc) · 1.24 KB
/
binBamReads
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
#!/usr/bin/env ruby
if (ARGV.size != 1)
STDERR.printf("Usage: %s bam-file\n", File.basename($0))
exit(1)
end
bam, rest = ARGV
goodPaired = 0
pairMap1 = Hash.new
pairMap2 = Hash.new
`samtools view #{bam}`.split("\n").each do |line|
name, flag, contig = line.chomp.split("\t")
if (name[1] == "1")
pairMap1[name] = contig
pairMap2[name] = "*" if (!pairMap2[name])
else
name[1] = "1"
pairMap2[name] = contig
pairMap1[name] = "*" if (!pairMap1[name])
end
end
fname = File.basename(bam, ".bam")
mismatch = File.new(fname + ".mismatch", "w")
unpaired = File.new(fname + ".unpaired", "w")
unmapped = File.new(fname + ".unmapped", "w")
pairMap1.keys.sort.each do |key|
if (pairMap1[key] != pairMap2[key] && pairMap1[key] != "*" && pairMap2[key] != "*")
key2 = key.dup
key2[1] = "2"
mismatch.printf("%s\t%s\t%s\t%s\n", key, key2, pairMap1[key],pairMap2[key])
elsif(pairMap1[key] == "*" && pairMap2[key] == "*")
unmapped.printf("%s\t%s\t%s\t%s\n", key, key2, pairMap1[key],pairMap2[key])
elsif(pairMap1[key] == "*" || pairMap2[key] == "*")
unpaired.printf("%s\t%s\t%s\t%s\n", key, key2, pairMap1[key],pairMap2[key])
end
end
mismatch.close
unpaired.close
unmapped.close
STDERR.printf("There are %d pairs", pairMap1.keys.size)