Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Jump to bottom

Added new option: samtools depth -m <int> maximum coverage depth [2147483647] #375

Closed
wants to merge 1 commit into from
Closed

Added new option: samtools depth -m <int> maximum coverage depth [2147483647] #375

Changes from all commits
Commits
Show all changes
1 commit
Select commit
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
10 changes: 8 additions & 2 deletions bam2depth.c
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -69,7 +69,7 @@ int read_file_list(const char *file_list,int *n,char **argv[]);

int main_depth(int argc, char *argv[])
{
int i, n, tid, beg, end, pos, *n_plp, baseQ = 0, mapQ = 0, min_len = 0, status = EXIT_SUCCESS, nfiles;
int i, n, tid, beg, end, pos, *n_plp, baseQ = 0, mapQ = 0, min_len = 0, status = EXIT_SUCCESS, nfiles, max_depth = -1;
const bam_pileup1_t **plp;
char *reg = 0; // specified region
void *bed = 0; // BED data structure
Expand All @@ -79,7 +79,7 @@ int main_depth(int argc, char *argv[])
bam_mplp_t mplp;

// parse the command line
while ((n = getopt(argc, argv, "r:b:q:Q:l:f:")) >= 0) {
while ((n = getopt(argc, argv, "r:b:q:Q:l:f:m:")) >= 0) {
switch (n) {
case 'l': min_len = atoi(optarg); break; // minimum query length
case 'r': reg = strdup(optarg); break; // parsing a region requires a BAM header
Expand All @@ -90,6 +90,7 @@ int main_depth(int argc, char *argv[])
case 'q': baseQ = atoi(optarg); break; // base quality threshold
case 'Q': mapQ = atoi(optarg); break; // mapping quality threshold
case 'f': file_list = optarg; break;
case 'm': max_depth = atoi(optarg); break; // maximum coverage depth
}
}
if (optind == argc && !file_list) {
Expand All @@ -102,6 +103,7 @@ int main_depth(int argc, char *argv[])
fprintf(stderr, " -q <int> base quality threshold\n");
fprintf(stderr, " -Q <int> mapping quality threshold\n");
fprintf(stderr, " -r <chr:from-to> region\n");
fprintf(stderr, " -m <int> maximum coverage depth [2147483647]\n");
fprintf(stderr, "\n");
return 1;
}
Expand Down Expand Up @@ -164,6 +166,10 @@ int main_depth(int argc, char *argv[])

// the core multi-pileup loop
mplp = bam_mplp_init(n, read_bam, (void**)data); // initialization
if (0 < max_depth)
bam_mplp_set_maxcnt(mplp,max_depth); // set maximum coverage depth
else
bam_mplp_set_maxcnt(mplp,2147483647); // set default maximum coverage depth to MAX INT32 = (2^31)-1
n_plp = calloc(n, sizeof(int)); // n_plp[i] is the number of covering reads from the i-th BAM
plp = calloc(n, sizeof(bam_pileup1_t*)); // plp[i] points to the array of covering reads (internal in mplp)
while (bam_mplp_auto(mplp, &tid, &pos, n_plp, plp) > 0) { // come to the next covered position
Expand Down