← Back to branch summary

~ubuntu-branches/ubuntu/raring/samtools/raring

« back to all changes in this revision

Viewing changes to sam_view.c

  • Committer: Bazaar Package Importer
  • Author(s): Charles Plessy
  • Date: 2009-11-17 21:38:24 UTC
  • Revision ID: james.westby@ubuntu.com-20091117213824-dfouynpy3r7ismpj
Tags: 0.1.7a~dfsg-1
* New upstream release: new script sam2vcf.pl, and many other changes.
* Package converted to the format ‘3.0 (quilt)’ (debian/source/format).
* Wrote a manual page for razip (debian/razip.1).
* Better clean the example directory to make the source package
  buildable twice in a row (debian/rules).

Show diffs side-by-side

added added

removed removed

Lines of Context:
sam_view.c
2
2
#include <string.h>
3
3
#include <stdio.h>
4
4
#include <unistd.h>
 
5
#include <math.h>
 
6
#include "sam_header.h"
5
7
#include "sam.h"
6
8
#include "faidx.h"
7
9
 
8
10
static int g_min_mapQ = 0, g_flag_on = 0, g_flag_off = 0;
9
11
static char *g_library, *g_rg;
 
12
static int g_sol2sanger_tbl[128];
 
13
 
 
14
static void sol2sanger(bam1_t *b)
 
15
{
 
16
        int l;
 
17
        uint8_t *qual = bam1_qual(b);
 
18
        if (g_sol2sanger_tbl[30] == 0) {
 
19
                for (l = 0; l != 128; ++l) {
 
20
                        g_sol2sanger_tbl[l] = (int)(10.0 * log(1.0 + pow(10.0, (l - 64 + 33) / 10.0)) / log(10.0) + .499);
 
21
                        if (g_sol2sanger_tbl[l] >= 93) g_sol2sanger_tbl[l] = 93;
 
22
                }
 
23
        }
 
24
        for (l = 0; l < b->core.l_qseq; ++l) {
 
25
                int q = qual[l];
 
26
                if (q > 127) q = 127;
 
27
                qual[l] = g_sol2sanger_tbl[q];
 
28
        }
 
29
}
10
30
 
11
31
static inline int __g_skip_aln(const bam_header_t *h, const bam1_t *b)
12
32
{
13
33
        if (b->core.qual < g_min_mapQ || ((b->core.flag & g_flag_on) != g_flag_on) || (b->core.flag & g_flag_off))
14
34
                return 1;
15
 
        if (g_library || g_rg) {
 
35
        if (g_rg) {
16
36
                uint8_t *s = bam_aux_get(b, "RG");
17
 
                if (s) {
18
 
                        if (g_rg && strcmp(g_rg, (char*)(s + 1)) == 0) return 0;
19
 
                        if (g_library) {
20
 
                                const char *p = bam_strmap_get(h->rg2lib, (char*)(s + 1));
21
 
                                return (p && strcmp(p, g_library) == 0)? 0 : 1;
22
 
                        } return 1;
23
 
                } else return 1;
24
 
        } else return 0;
 
37
                if (s && strcmp(g_rg, (char*)(s + 1)) == 0) return 0;
 
38
        }
 
39
        if (g_library) {
 
40
                const char *p = bam_get_library((bam_header_t*)h, b);
 
41
                return (p && strcmp(p, g_library) == 0)? 0 : 1;
 
42
        }
 
43
        return 0;
25
44
}
26
45
 
27
46
// callback function for bam_fetch()
36
55
 
37
56
int main_samview(int argc, char *argv[])
38
57
{
39
 
        int c, is_header = 0, is_header_only = 0, is_bamin = 1, ret = 0, is_uncompressed = 0, is_bamout = 0;
 
58
        int c, is_header = 0, is_header_only = 0, is_bamin = 1, ret = 0, is_uncompressed = 0, is_bamout = 0, slx2sngr = 0;
40
59
        int of_type = BAM_OFDEC, is_long_help = 0;
41
60
        samfile_t *in = 0, *out = 0;
42
61
        char in_mode[5], out_mode[5], *fn_out = 0, *fn_list = 0, *fn_ref = 0;
43
62
 
44
63
        /* parse command-line options */
45
64
        strcpy(in_mode, "r"); strcpy(out_mode, "w");
46
 
        while ((c = getopt(argc, argv, "Sbt:hHo:q:f:F:ul:r:xX?T:")) >= 0) {
 
65
        while ((c = getopt(argc, argv, "Sbt:hHo:q:f:F:ul:r:xX?T:C")) >= 0) {
47
66
                switch (c) {
 
67
                case 'C': slx2sngr = 1; break;
48
68
                case 'S': is_bamin = 0; break;
49
69
                case 'b': is_bamout = 1; break;
50
70
                case 't': fn_list = strdup(optarg); is_bamin = 0; break;
96
116
        if (argc == optind + 1) { // convert/print the entire file
97
117
                bam1_t *b = bam_init1();
98
118
                int r;
99
 
                while ((r = samread(in, b)) >= 0) // read one alignment from `in'
100
 
                        if (!__g_skip_aln(in->header, b))
 
119
                while ((r = samread(in, b)) >= 0) { // read one alignment from `in'
 
120
                        if (!__g_skip_aln(in->header, b)) {
 
121
                                if (slx2sngr) sol2sanger(b);
101
122
                                samwrite(out, b); // write the alignment to `out'
 
123
                        }
 
124
                }
102
125
                if (r < -1) fprintf(stderr, "[main_samview] truncated file.\n");
103
126
                bam_destroy1(b);
104
127
        } else { // retrieve alignments in specified regions