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