annotate deseq-hts_1.0/mex/get_bam_properties.cpp @ 11:cec4b4fb30be draft default tip

DEXSeq version 1.6 added
author vipints <vipin@cbio.mskcc.org>
date Tue, 08 Oct 2013 08:22:45 -0400
parents 94a108763d9e
children
Ignore whitespace changes - Everywhere: Within whitespace: At end of lines:
rev   line source
0
94a108763d9e deseq-hts version 1.0 wraps the DESeq 1.6.0
vipints
parents:
diff changeset
1 /*
94a108763d9e deseq-hts version 1.0 wraps the DESeq 1.6.0
vipints
parents:
diff changeset
2 * This program is free software; you can redistribute it and/or modify
94a108763d9e deseq-hts version 1.0 wraps the DESeq 1.6.0
vipints
parents:
diff changeset
3 * it under the terms of the GNU General Public License as published by
94a108763d9e deseq-hts version 1.0 wraps the DESeq 1.6.0
vipints
parents:
diff changeset
4 * the Free Software Foundation; either version 3 of the License, or
94a108763d9e deseq-hts version 1.0 wraps the DESeq 1.6.0
vipints
parents:
diff changeset
5 * (at your option) any later version.
94a108763d9e deseq-hts version 1.0 wraps the DESeq 1.6.0
vipints
parents:
diff changeset
6 *
94a108763d9e deseq-hts version 1.0 wraps the DESeq 1.6.0
vipints
parents:
diff changeset
7 * Written (W) 2009-2011 Regina Bohnert
94a108763d9e deseq-hts version 1.0 wraps the DESeq 1.6.0
vipints
parents:
diff changeset
8 * Copyright (C) 2009-2011 Max Planck Society
94a108763d9e deseq-hts version 1.0 wraps the DESeq 1.6.0
vipints
parents:
diff changeset
9 */
94a108763d9e deseq-hts version 1.0 wraps the DESeq 1.6.0
vipints
parents:
diff changeset
10
94a108763d9e deseq-hts version 1.0 wraps the DESeq 1.6.0
vipints
parents:
diff changeset
11
94a108763d9e deseq-hts version 1.0 wraps the DESeq 1.6.0
vipints
parents:
diff changeset
12 #include <stdio.h>
94a108763d9e deseq-hts version 1.0 wraps the DESeq 1.6.0
vipints
parents:
diff changeset
13 #include <stdlib.h>
94a108763d9e deseq-hts version 1.0 wraps the DESeq 1.6.0
vipints
parents:
diff changeset
14 #include <signal.h>
94a108763d9e deseq-hts version 1.0 wraps the DESeq 1.6.0
vipints
parents:
diff changeset
15 #include <ctype.h>
94a108763d9e deseq-hts version 1.0 wraps the DESeq 1.6.0
vipints
parents:
diff changeset
16 #include <assert.h>
94a108763d9e deseq-hts version 1.0 wraps the DESeq 1.6.0
vipints
parents:
diff changeset
17
94a108763d9e deseq-hts version 1.0 wraps the DESeq 1.6.0
vipints
parents:
diff changeset
18 #include <vector>
94a108763d9e deseq-hts version 1.0 wraps the DESeq 1.6.0
vipints
parents:
diff changeset
19 using std::vector;
94a108763d9e deseq-hts version 1.0 wraps the DESeq 1.6.0
vipints
parents:
diff changeset
20 #include <string>
94a108763d9e deseq-hts version 1.0 wraps the DESeq 1.6.0
vipints
parents:
diff changeset
21 using std::string;
94a108763d9e deseq-hts version 1.0 wraps the DESeq 1.6.0
vipints
parents:
diff changeset
22 #include <algorithm>
94a108763d9e deseq-hts version 1.0 wraps the DESeq 1.6.0
vipints
parents:
diff changeset
23 using std::find;
94a108763d9e deseq-hts version 1.0 wraps the DESeq 1.6.0
vipints
parents:
diff changeset
24 using std::min;
94a108763d9e deseq-hts version 1.0 wraps the DESeq 1.6.0
vipints
parents:
diff changeset
25
94a108763d9e deseq-hts version 1.0 wraps the DESeq 1.6.0
vipints
parents:
diff changeset
26 #include <mex.h>
94a108763d9e deseq-hts version 1.0 wraps the DESeq 1.6.0
vipints
parents:
diff changeset
27
94a108763d9e deseq-hts version 1.0 wraps the DESeq 1.6.0
vipints
parents:
diff changeset
28
94a108763d9e deseq-hts version 1.0 wraps the DESeq 1.6.0
vipints
parents:
diff changeset
29 char *get_string(const mxArray *prhs);
94a108763d9e deseq-hts version 1.0 wraps the DESeq 1.6.0
vipints
parents:
diff changeset
30
94a108763d9e deseq-hts version 1.0 wraps the DESeq 1.6.0
vipints
parents:
diff changeset
31 typedef unsigned int uint32_t;
94a108763d9e deseq-hts version 1.0 wraps the DESeq 1.6.0
vipints
parents:
diff changeset
32 typedef unsigned char uint8_t;
94a108763d9e deseq-hts version 1.0 wraps the DESeq 1.6.0
vipints
parents:
diff changeset
33
94a108763d9e deseq-hts version 1.0 wraps the DESeq 1.6.0
vipints
parents:
diff changeset
34 /*
94a108763d9e deseq-hts version 1.0 wraps the DESeq 1.6.0
vipints
parents:
diff changeset
35 * [read_len num_reads] = get_bam_properties(fname, path_samtools, contig_name)
94a108763d9e deseq-hts version 1.0 wraps the DESeq 1.6.0
vipints
parents:
diff changeset
36 *
94a108763d9e deseq-hts version 1.0 wraps the DESeq 1.6.0
vipints
parents:
diff changeset
37 * -- input --
94a108763d9e deseq-hts version 1.0 wraps the DESeq 1.6.0
vipints
parents:
diff changeset
38 * prhs[0] file name of paired reads in BAM format (sorted by read id)
94a108763d9e deseq-hts version 1.0 wraps the DESeq 1.6.0
vipints
parents:
diff changeset
39 * prhs[1] path to samtools
94a108763d9e deseq-hts version 1.0 wraps the DESeq 1.6.0
vipints
parents:
diff changeset
40 * prhs[2] contig name
94a108763d9e deseq-hts version 1.0 wraps the DESeq 1.6.0
vipints
parents:
diff changeset
41 *
94a108763d9e deseq-hts version 1.0 wraps the DESeq 1.6.0
vipints
parents:
diff changeset
42 * -- output --
94a108763d9e deseq-hts version 1.0 wraps the DESeq 1.6.0
vipints
parents:
diff changeset
43 * plhs[0] length of read
94a108763d9e deseq-hts version 1.0 wraps the DESeq 1.6.0
vipints
parents:
diff changeset
44 * plhs[1] number of unique reads
94a108763d9e deseq-hts version 1.0 wraps the DESeq 1.6.0
vipints
parents:
diff changeset
45 */
94a108763d9e deseq-hts version 1.0 wraps the DESeq 1.6.0
vipints
parents:
diff changeset
46 void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[]) {
94a108763d9e deseq-hts version 1.0 wraps the DESeq 1.6.0
vipints
parents:
diff changeset
47 // checks for the right number of arguments
94a108763d9e deseq-hts version 1.0 wraps the DESeq 1.6.0
vipints
parents:
diff changeset
48 if (nrhs !=3 || nlhs > 2) {
94a108763d9e deseq-hts version 1.0 wraps the DESeq 1.6.0
vipints
parents:
diff changeset
49 mexErrMsgTxt("number of input and output args should be 3 and 2\nUSAGE:\n [read_len, num_reads] = get_bam_properties(fname, path_samtools, contig_name)\n");
94a108763d9e deseq-hts version 1.0 wraps the DESeq 1.6.0
vipints
parents:
diff changeset
50 return;
94a108763d9e deseq-hts version 1.0 wraps the DESeq 1.6.0
vipints
parents:
diff changeset
51 }
94a108763d9e deseq-hts version 1.0 wraps the DESeq 1.6.0
vipints
parents:
diff changeset
52
94a108763d9e deseq-hts version 1.0 wraps the DESeq 1.6.0
vipints
parents:
diff changeset
53 signal(SIGCHLD, SIG_IGN); // avoid zombies
94a108763d9e deseq-hts version 1.0 wraps the DESeq 1.6.0
vipints
parents:
diff changeset
54
94a108763d9e deseq-hts version 1.0 wraps the DESeq 1.6.0
vipints
parents:
diff changeset
55 // read input arguments
94a108763d9e deseq-hts version 1.0 wraps the DESeq 1.6.0
vipints
parents:
diff changeset
56 char *fname = get_string(prhs[0]);
94a108763d9e deseq-hts version 1.0 wraps the DESeq 1.6.0
vipints
parents:
diff changeset
57 char *path_samtools = get_string(prhs[1]);
94a108763d9e deseq-hts version 1.0 wraps the DESeq 1.6.0
vipints
parents:
diff changeset
58 char *contig_name = get_string(prhs[2]);
94a108763d9e deseq-hts version 1.0 wraps the DESeq 1.6.0
vipints
parents:
diff changeset
59 char command[10000];
94a108763d9e deseq-hts version 1.0 wraps the DESeq 1.6.0
vipints
parents:
diff changeset
60
94a108763d9e deseq-hts version 1.0 wraps the DESeq 1.6.0
vipints
parents:
diff changeset
61 sprintf(command, "%s./samtools view %s %s 2>/dev/null", path_samtools, fname, contig_name);
94a108763d9e deseq-hts version 1.0 wraps the DESeq 1.6.0
vipints
parents:
diff changeset
62 //printf("%s\n", command);
94a108763d9e deseq-hts version 1.0 wraps the DESeq 1.6.0
vipints
parents:
diff changeset
63
94a108763d9e deseq-hts version 1.0 wraps the DESeq 1.6.0
vipints
parents:
diff changeset
64 // get number of unique reads
94a108763d9e deseq-hts version 1.0 wraps the DESeq 1.6.0
vipints
parents:
diff changeset
65 int status;
94a108763d9e deseq-hts version 1.0 wraps the DESeq 1.6.0
vipints
parents:
diff changeset
66 uint32_t num_unique_reads = 0;
94a108763d9e deseq-hts version 1.0 wraps the DESeq 1.6.0
vipints
parents:
diff changeset
67 char command2[10000];
94a108763d9e deseq-hts version 1.0 wraps the DESeq 1.6.0
vipints
parents:
diff changeset
68 sprintf(command2, "%s | cut -f 1 | sort -u | wc -l", command);
94a108763d9e deseq-hts version 1.0 wraps the DESeq 1.6.0
vipints
parents:
diff changeset
69 FILE* fp = popen(command2, "r");
94a108763d9e deseq-hts version 1.0 wraps the DESeq 1.6.0
vipints
parents:
diff changeset
70 if (fp == NULL) {
94a108763d9e deseq-hts version 1.0 wraps the DESeq 1.6.0
vipints
parents:
diff changeset
71 mexErrMsgTxt("Error using popen\n");
94a108763d9e deseq-hts version 1.0 wraps the DESeq 1.6.0
vipints
parents:
diff changeset
72 }
94a108763d9e deseq-hts version 1.0 wraps the DESeq 1.6.0
vipints
parents:
diff changeset
73 int num_scans = 1;
94a108763d9e deseq-hts version 1.0 wraps the DESeq 1.6.0
vipints
parents:
diff changeset
74 num_scans = fscanf(fp, "%d", &num_unique_reads);
94a108763d9e deseq-hts version 1.0 wraps the DESeq 1.6.0
vipints
parents:
diff changeset
75 if (num_scans != 1) {
94a108763d9e deseq-hts version 1.0 wraps the DESeq 1.6.0
vipints
parents:
diff changeset
76 rewind(fp);
94a108763d9e deseq-hts version 1.0 wraps the DESeq 1.6.0
vipints
parents:
diff changeset
77 char ret[1000];
94a108763d9e deseq-hts version 1.0 wraps the DESeq 1.6.0
vipints
parents:
diff changeset
78 fgets(ret, 1000, fp);
94a108763d9e deseq-hts version 1.0 wraps the DESeq 1.6.0
vipints
parents:
diff changeset
79 fprintf(stdout, "%s", ret);
94a108763d9e deseq-hts version 1.0 wraps the DESeq 1.6.0
vipints
parents:
diff changeset
80 mexErrMsgTxt("Could not determine number of reads\n");
94a108763d9e deseq-hts version 1.0 wraps the DESeq 1.6.0
vipints
parents:
diff changeset
81 }
94a108763d9e deseq-hts version 1.0 wraps the DESeq 1.6.0
vipints
parents:
diff changeset
82 status = pclose(fp);
94a108763d9e deseq-hts version 1.0 wraps the DESeq 1.6.0
vipints
parents:
diff changeset
83 //printf("%i", num_unique_reads);
94a108763d9e deseq-hts version 1.0 wraps the DESeq 1.6.0
vipints
parents:
diff changeset
84
94a108763d9e deseq-hts version 1.0 wraps the DESeq 1.6.0
vipints
parents:
diff changeset
85 // select reads for given positions and strand
94a108763d9e deseq-hts version 1.0 wraps the DESeq 1.6.0
vipints
parents:
diff changeset
86 int num_rows_selected = min((int) num_unique_reads, 100);
94a108763d9e deseq-hts version 1.0 wraps the DESeq 1.6.0
vipints
parents:
diff changeset
87 sprintf(command, "%s | head -n %i | cut -f 1-11", command, num_rows_selected);
94a108763d9e deseq-hts version 1.0 wraps the DESeq 1.6.0
vipints
parents:
diff changeset
88 fp = popen(command, "r");
94a108763d9e deseq-hts version 1.0 wraps the DESeq 1.6.0
vipints
parents:
diff changeset
89 if (fp == NULL) {
94a108763d9e deseq-hts version 1.0 wraps the DESeq 1.6.0
vipints
parents:
diff changeset
90 mexErrMsgTxt("Error using popen\n");
94a108763d9e deseq-hts version 1.0 wraps the DESeq 1.6.0
vipints
parents:
diff changeset
91 }
94a108763d9e deseq-hts version 1.0 wraps the DESeq 1.6.0
vipints
parents:
diff changeset
92 /* SAM format
94a108763d9e deseq-hts version 1.0 wraps the DESeq 1.6.0
vipints
parents:
diff changeset
93 1: read id, 2: flag, 3: reference name, 4: start (1-based, incl.), 5: mapping quality,
94a108763d9e deseq-hts version 1.0 wraps the DESeq 1.6.0
vipints
parents:
diff changeset
94 6: CIGAR, 7: mate reference name, 8: mate start (1-based, incl.), 9: insert size, 10: read, 11: quality
94a108763d9e deseq-hts version 1.0 wraps the DESeq 1.6.0
vipints
parents:
diff changeset
95 12+: additional tags
94a108763d9e deseq-hts version 1.0 wraps the DESeq 1.6.0
vipints
parents:
diff changeset
96 */
94a108763d9e deseq-hts version 1.0 wraps the DESeq 1.6.0
vipints
parents:
diff changeset
97 uint32_t read_idx = 0, row_idx = 0, num_col = 0;
94a108763d9e deseq-hts version 1.0 wraps the DESeq 1.6.0
vipints
parents:
diff changeset
98 uint32_t flag = 0, start_pos = 0, map_score = 0, mate_end_pos = 0, num_matches = 0, num_del = 0, num_ins = 0, ins_size = 0;
94a108763d9e deseq-hts version 1.0 wraps the DESeq 1.6.0
vipints
parents:
diff changeset
99 char ri [1000], read_contig_name [1000], cg [1000], mate_read_id [1000], read [1000], read_qual [1000];
94a108763d9e deseq-hts version 1.0 wraps the DESeq 1.6.0
vipints
parents:
diff changeset
100 string last_read_id;
94a108763d9e deseq-hts version 1.0 wraps the DESeq 1.6.0
vipints
parents:
diff changeset
101 vector<uint32_t> block_lengths, block_starts;
94a108763d9e deseq-hts version 1.0 wraps the DESeq 1.6.0
vipints
parents:
diff changeset
102 vector<string> read_ids;
94a108763d9e deseq-hts version 1.0 wraps the DESeq 1.6.0
vipints
parents:
diff changeset
103 vector<string>::iterator it;
94a108763d9e deseq-hts version 1.0 wraps the DESeq 1.6.0
vipints
parents:
diff changeset
104
94a108763d9e deseq-hts version 1.0 wraps the DESeq 1.6.0
vipints
parents:
diff changeset
105 uint32_t read_len = 0;
94a108763d9e deseq-hts version 1.0 wraps the DESeq 1.6.0
vipints
parents:
diff changeset
106 bool empty_line = true;
94a108763d9e deseq-hts version 1.0 wraps the DESeq 1.6.0
vipints
parents:
diff changeset
107 int num_rows = 0;
94a108763d9e deseq-hts version 1.0 wraps the DESeq 1.6.0
vipints
parents:
diff changeset
108 while(empty_line && num_rows < num_rows_selected) {
94a108763d9e deseq-hts version 1.0 wraps the DESeq 1.6.0
vipints
parents:
diff changeset
109 num_col = fscanf(fp, "%s\t%i\t%s\t%i\t%i\t%s\t%s\t%i\t%i\t%s\t%s", &ri, &flag, &read_contig_name, &start_pos, &map_score, &cg, &mate_read_id, &mate_end_pos, &ins_size, &read, &read_qual);
94a108763d9e deseq-hts version 1.0 wraps the DESeq 1.6.0
vipints
parents:
diff changeset
110 if (num_col != 11) {
94a108763d9e deseq-hts version 1.0 wraps the DESeq 1.6.0
vipints
parents:
diff changeset
111 mexErrMsgTxt("error reading SAM line\n");
94a108763d9e deseq-hts version 1.0 wraps the DESeq 1.6.0
vipints
parents:
diff changeset
112 }
94a108763d9e deseq-hts version 1.0 wraps the DESeq 1.6.0
vipints
parents:
diff changeset
113
94a108763d9e deseq-hts version 1.0 wraps the DESeq 1.6.0
vipints
parents:
diff changeset
114 string cigar = (string) cg;
94a108763d9e deseq-hts version 1.0 wraps the DESeq 1.6.0
vipints
parents:
diff changeset
115 // ignore lines with reads w/o mapping information
94a108763d9e deseq-hts version 1.0 wraps the DESeq 1.6.0
vipints
parents:
diff changeset
116 if (start_pos == 0 || cigar.compare("*")==0) {
94a108763d9e deseq-hts version 1.0 wraps the DESeq 1.6.0
vipints
parents:
diff changeset
117 continue;
94a108763d9e deseq-hts version 1.0 wraps the DESeq 1.6.0
vipints
parents:
diff changeset
118 }
94a108763d9e deseq-hts version 1.0 wraps the DESeq 1.6.0
vipints
parents:
diff changeset
119 // parse CIGAR
94a108763d9e deseq-hts version 1.0 wraps the DESeq 1.6.0
vipints
parents:
diff changeset
120 uint last_c = 0;
94a108763d9e deseq-hts version 1.0 wraps the DESeq 1.6.0
vipints
parents:
diff changeset
121 string last_str;
94a108763d9e deseq-hts version 1.0 wraps the DESeq 1.6.0
vipints
parents:
diff changeset
122 num_matches = 0;
94a108763d9e deseq-hts version 1.0 wraps the DESeq 1.6.0
vipints
parents:
diff changeset
123 char *end = NULL;
94a108763d9e deseq-hts version 1.0 wraps the DESeq 1.6.0
vipints
parents:
diff changeset
124 uint32_t tmp_nm = 0, tmp_nd = 0, tmp_ni = 0;
94a108763d9e deseq-hts version 1.0 wraps the DESeq 1.6.0
vipints
parents:
diff changeset
125 uint32_t last_block_start = 0, last_block_length = 0, last_intron_len = 0;
94a108763d9e deseq-hts version 1.0 wraps the DESeq 1.6.0
vipints
parents:
diff changeset
126 block_lengths.clear(); block_starts.clear();
94a108763d9e deseq-hts version 1.0 wraps the DESeq 1.6.0
vipints
parents:
diff changeset
127
94a108763d9e deseq-hts version 1.0 wraps the DESeq 1.6.0
vipints
parents:
diff changeset
128 for (uint c = 0; c < cigar.size(); c++) {
94a108763d9e deseq-hts version 1.0 wraps the DESeq 1.6.0
vipints
parents:
diff changeset
129 switch (cigar[c]) {
94a108763d9e deseq-hts version 1.0 wraps the DESeq 1.6.0
vipints
parents:
diff changeset
130 case 'M':
94a108763d9e deseq-hts version 1.0 wraps the DESeq 1.6.0
vipints
parents:
diff changeset
131 last_str = cigar.substr(last_c, c-last_c);
94a108763d9e deseq-hts version 1.0 wraps the DESeq 1.6.0
vipints
parents:
diff changeset
132 tmp_nm = strtoul(last_str.c_str(), &end, 10);
94a108763d9e deseq-hts version 1.0 wraps the DESeq 1.6.0
vipints
parents:
diff changeset
133 if (*end != '\0')
94a108763d9e deseq-hts version 1.0 wraps the DESeq 1.6.0
vipints
parents:
diff changeset
134 mexErrMsgTxt("error: number of mismatches\n");
94a108763d9e deseq-hts version 1.0 wraps the DESeq 1.6.0
vipints
parents:
diff changeset
135 end = NULL;
94a108763d9e deseq-hts version 1.0 wraps the DESeq 1.6.0
vipints
parents:
diff changeset
136 last_block_length += tmp_nm;
94a108763d9e deseq-hts version 1.0 wraps the DESeq 1.6.0
vipints
parents:
diff changeset
137 num_matches += tmp_nm;
94a108763d9e deseq-hts version 1.0 wraps the DESeq 1.6.0
vipints
parents:
diff changeset
138 last_c = c + 1;
94a108763d9e deseq-hts version 1.0 wraps the DESeq 1.6.0
vipints
parents:
diff changeset
139 break;
94a108763d9e deseq-hts version 1.0 wraps the DESeq 1.6.0
vipints
parents:
diff changeset
140 case 'I':
94a108763d9e deseq-hts version 1.0 wraps the DESeq 1.6.0
vipints
parents:
diff changeset
141 last_str = cigar.substr(last_c, c-last_c);
94a108763d9e deseq-hts version 1.0 wraps the DESeq 1.6.0
vipints
parents:
diff changeset
142 tmp_ni = strtoul(last_str.c_str(), &end, 10);
94a108763d9e deseq-hts version 1.0 wraps the DESeq 1.6.0
vipints
parents:
diff changeset
143 if (*end != '\0')
94a108763d9e deseq-hts version 1.0 wraps the DESeq 1.6.0
vipints
parents:
diff changeset
144 mexErrMsgTxt("error: number of insertions\n");
94a108763d9e deseq-hts version 1.0 wraps the DESeq 1.6.0
vipints
parents:
diff changeset
145 end = NULL;
94a108763d9e deseq-hts version 1.0 wraps the DESeq 1.6.0
vipints
parents:
diff changeset
146 num_ins += tmp_ni;
94a108763d9e deseq-hts version 1.0 wraps the DESeq 1.6.0
vipints
parents:
diff changeset
147 last_c = c + 1;
94a108763d9e deseq-hts version 1.0 wraps the DESeq 1.6.0
vipints
parents:
diff changeset
148 break;
94a108763d9e deseq-hts version 1.0 wraps the DESeq 1.6.0
vipints
parents:
diff changeset
149 case 'D':
94a108763d9e deseq-hts version 1.0 wraps the DESeq 1.6.0
vipints
parents:
diff changeset
150 last_str = cigar.substr(last_c, c-last_c);
94a108763d9e deseq-hts version 1.0 wraps the DESeq 1.6.0
vipints
parents:
diff changeset
151 tmp_nd = strtoul(last_str.c_str(), &end, 10);
94a108763d9e deseq-hts version 1.0 wraps the DESeq 1.6.0
vipints
parents:
diff changeset
152 if (*end != '\0')
94a108763d9e deseq-hts version 1.0 wraps the DESeq 1.6.0
vipints
parents:
diff changeset
153 mexErrMsgTxt("error: number of deletions\n");
94a108763d9e deseq-hts version 1.0 wraps the DESeq 1.6.0
vipints
parents:
diff changeset
154 end = NULL;
94a108763d9e deseq-hts version 1.0 wraps the DESeq 1.6.0
vipints
parents:
diff changeset
155 num_del += tmp_nd;
94a108763d9e deseq-hts version 1.0 wraps the DESeq 1.6.0
vipints
parents:
diff changeset
156 last_block_length += tmp_nd;
94a108763d9e deseq-hts version 1.0 wraps the DESeq 1.6.0
vipints
parents:
diff changeset
157 last_c = c + 1;
94a108763d9e deseq-hts version 1.0 wraps the DESeq 1.6.0
vipints
parents:
diff changeset
158 break;
94a108763d9e deseq-hts version 1.0 wraps the DESeq 1.6.0
vipints
parents:
diff changeset
159 case 'N':
94a108763d9e deseq-hts version 1.0 wraps the DESeq 1.6.0
vipints
parents:
diff changeset
160 last_str = cigar.substr(last_c, c-last_c);
94a108763d9e deseq-hts version 1.0 wraps the DESeq 1.6.0
vipints
parents:
diff changeset
161 last_intron_len = strtoul(last_str.c_str(), &end, 10);
94a108763d9e deseq-hts version 1.0 wraps the DESeq 1.6.0
vipints
parents:
diff changeset
162 end = NULL;
94a108763d9e deseq-hts version 1.0 wraps the DESeq 1.6.0
vipints
parents:
diff changeset
163 last_c = c + 1;
94a108763d9e deseq-hts version 1.0 wraps the DESeq 1.6.0
vipints
parents:
diff changeset
164 break;
94a108763d9e deseq-hts version 1.0 wraps the DESeq 1.6.0
vipints
parents:
diff changeset
165 case 'S':
94a108763d9e deseq-hts version 1.0 wraps the DESeq 1.6.0
vipints
parents:
diff changeset
166 break;
94a108763d9e deseq-hts version 1.0 wraps the DESeq 1.6.0
vipints
parents:
diff changeset
167 case 'H':
94a108763d9e deseq-hts version 1.0 wraps the DESeq 1.6.0
vipints
parents:
diff changeset
168 break;
94a108763d9e deseq-hts version 1.0 wraps the DESeq 1.6.0
vipints
parents:
diff changeset
169 case 'P':
94a108763d9e deseq-hts version 1.0 wraps the DESeq 1.6.0
vipints
parents:
diff changeset
170 break;
94a108763d9e deseq-hts version 1.0 wraps the DESeq 1.6.0
vipints
parents:
diff changeset
171 default:
94a108763d9e deseq-hts version 1.0 wraps the DESeq 1.6.0
vipints
parents:
diff changeset
172 break;
94a108763d9e deseq-hts version 1.0 wraps the DESeq 1.6.0
vipints
parents:
diff changeset
173 }
94a108763d9e deseq-hts version 1.0 wraps the DESeq 1.6.0
vipints
parents:
diff changeset
174 if (cigar[c] == 'N' || c==cigar.size()-1) {
94a108763d9e deseq-hts version 1.0 wraps the DESeq 1.6.0
vipints
parents:
diff changeset
175 block_starts.push_back(last_block_start);
94a108763d9e deseq-hts version 1.0 wraps the DESeq 1.6.0
vipints
parents:
diff changeset
176 last_block_start = last_block_start + last_block_length + last_intron_len;
94a108763d9e deseq-hts version 1.0 wraps the DESeq 1.6.0
vipints
parents:
diff changeset
177 last_intron_len = 0;
94a108763d9e deseq-hts version 1.0 wraps the DESeq 1.6.0
vipints
parents:
diff changeset
178 block_lengths.push_back(last_block_length);
94a108763d9e deseq-hts version 1.0 wraps the DESeq 1.6.0
vipints
parents:
diff changeset
179 last_block_length = 0;
94a108763d9e deseq-hts version 1.0 wraps the DESeq 1.6.0
vipints
parents:
diff changeset
180 }
94a108763d9e deseq-hts version 1.0 wraps the DESeq 1.6.0
vipints
parents:
diff changeset
181 }
94a108763d9e deseq-hts version 1.0 wraps the DESeq 1.6.0
vipints
parents:
diff changeset
182 read_len = 0;
94a108763d9e deseq-hts version 1.0 wraps the DESeq 1.6.0
vipints
parents:
diff changeset
183 for (uint n = 0; n < block_lengths.size(); n++) {
94a108763d9e deseq-hts version 1.0 wraps the DESeq 1.6.0
vipints
parents:
diff changeset
184 read_len += block_lengths[n];
94a108763d9e deseq-hts version 1.0 wraps the DESeq 1.6.0
vipints
parents:
diff changeset
185 }
94a108763d9e deseq-hts version 1.0 wraps the DESeq 1.6.0
vipints
parents:
diff changeset
186 empty_line = false;
94a108763d9e deseq-hts version 1.0 wraps the DESeq 1.6.0
vipints
parents:
diff changeset
187 } // end of stream parsing
94a108763d9e deseq-hts version 1.0 wraps the DESeq 1.6.0
vipints
parents:
diff changeset
188
94a108763d9e deseq-hts version 1.0 wraps the DESeq 1.6.0
vipints
parents:
diff changeset
189 status = pclose(fp);
94a108763d9e deseq-hts version 1.0 wraps the DESeq 1.6.0
vipints
parents:
diff changeset
190
94a108763d9e deseq-hts version 1.0 wraps the DESeq 1.6.0
vipints
parents:
diff changeset
191 if (empty_line)
94a108763d9e deseq-hts version 1.0 wraps the DESeq 1.6.0
vipints
parents:
diff changeset
192 mexErrMsgTxt("Could not determine read length\n");
94a108763d9e deseq-hts version 1.0 wraps the DESeq 1.6.0
vipints
parents:
diff changeset
193
94a108763d9e deseq-hts version 1.0 wraps the DESeq 1.6.0
vipints
parents:
diff changeset
194 plhs[0] = mxCreateDoubleScalar((double) read_len);
94a108763d9e deseq-hts version 1.0 wraps the DESeq 1.6.0
vipints
parents:
diff changeset
195 plhs[1] = mxCreateDoubleScalar((double) num_unique_reads);
94a108763d9e deseq-hts version 1.0 wraps the DESeq 1.6.0
vipints
parents:
diff changeset
196
94a108763d9e deseq-hts version 1.0 wraps the DESeq 1.6.0
vipints
parents:
diff changeset
197 return;
94a108763d9e deseq-hts version 1.0 wraps the DESeq 1.6.0
vipints
parents:
diff changeset
198 }
94a108763d9e deseq-hts version 1.0 wraps the DESeq 1.6.0
vipints
parents:
diff changeset
199
94a108763d9e deseq-hts version 1.0 wraps the DESeq 1.6.0
vipints
parents:
diff changeset
200
94a108763d9e deseq-hts version 1.0 wraps the DESeq 1.6.0
vipints
parents:
diff changeset
201 char *get_string(const mxArray *prhs) {
94a108763d9e deseq-hts version 1.0 wraps the DESeq 1.6.0
vipints
parents:
diff changeset
202 char *buf;
94a108763d9e deseq-hts version 1.0 wraps the DESeq 1.6.0
vipints
parents:
diff changeset
203 int buflen;
94a108763d9e deseq-hts version 1.0 wraps the DESeq 1.6.0
vipints
parents:
diff changeset
204 if (!prhs)
94a108763d9e deseq-hts version 1.0 wraps the DESeq 1.6.0
vipints
parents:
diff changeset
205 mexErrMsgTxt("get_string called with NULL pointer arg");
94a108763d9e deseq-hts version 1.0 wraps the DESeq 1.6.0
vipints
parents:
diff changeset
206 if (!mxIsChar(prhs))
94a108763d9e deseq-hts version 1.0 wraps the DESeq 1.6.0
vipints
parents:
diff changeset
207 mexErrMsgTxt("input is not a string");
94a108763d9e deseq-hts version 1.0 wraps the DESeq 1.6.0
vipints
parents:
diff changeset
208 if (mxGetM(prhs) != 1)
94a108763d9e deseq-hts version 1.0 wraps the DESeq 1.6.0
vipints
parents:
diff changeset
209 mexErrMsgTxt("input is not a row vector");
94a108763d9e deseq-hts version 1.0 wraps the DESeq 1.6.0
vipints
parents:
diff changeset
210 buflen = mxGetN(prhs) + 1;
94a108763d9e deseq-hts version 1.0 wraps the DESeq 1.6.0
vipints
parents:
diff changeset
211 buf = (char*) malloc(buflen);
94a108763d9e deseq-hts version 1.0 wraps the DESeq 1.6.0
vipints
parents:
diff changeset
212 /* copy the string from prhs into buf and add terminating NULL char */
94a108763d9e deseq-hts version 1.0 wraps the DESeq 1.6.0
vipints
parents:
diff changeset
213 if (mxGetString(prhs, buf, buflen))
94a108763d9e deseq-hts version 1.0 wraps the DESeq 1.6.0
vipints
parents:
diff changeset
214 mexErrMsgTxt("not enough space");
94a108763d9e deseq-hts version 1.0 wraps the DESeq 1.6.0
vipints
parents:
diff changeset
215 return buf;
94a108763d9e deseq-hts version 1.0 wraps the DESeq 1.6.0
vipints
parents:
diff changeset
216 }