From 165f1d8cb88f04a48b7da3f4e5e5fd2a4f074972 Mon Sep 17 00:00:00 2001 From: Heng Li Date: Wed, 10 Nov 2010 17:10:21 +0000 Subject: [PATCH] * samtools-0.1.9-8 (r807) * collect indel candidates only from specified platforms (@RG-PL) * merge varFilter and filter4vcf in vcfutils.pl --- Makefile | 1 + bam2bcf.h | 4 ++- bam2bcf_indel.c | 74 ++++++++++++++++++++++++++++++++++++++++++-- bam_plcmd.c | 16 +++++++--- bamtk.c | 2 +- bcftools/vcfutils.pl | 67 +++++++++++---------------------------- 6 files changed, 106 insertions(+), 58 deletions(-) diff --git a/Makefile b/Makefile index eeae213..965491a 100644 --- a/Makefile +++ b/Makefile @@ -64,6 +64,7 @@ glf.o:glf.h sam_header.o:sam_header.h khash.h bcf.o:bcftools/bcf.h bam2bcf.o:bam2bcf.h errmod.h bcftools/bcf.h +bam2bcf_indel.o:bam2bcf.h errmod.o:errmod.h faidx.o:faidx.h razf.h khash.h diff --git a/bam2bcf.h b/bam2bcf.h index 570990e..e0673e5 100644 --- a/bam2bcf.h +++ b/bam2bcf.h @@ -17,6 +17,7 @@ typedef struct __bcf_callaux_t { char *inscns; uint16_t *bases; errmod_t *e; + void *rghash; } bcf_callaux_t; typedef struct { @@ -42,7 +43,8 @@ extern "C" { int bcf_call_combine(int n, const bcf_callret1_t *calls, int ref_base /*4-bit*/, bcf_call_t *call); int bcf_call2bcf(int tid, int pos, bcf_call_t *bc, bcf1_t *b, bcf_callret1_t *bcr, int is_SP, const bcf_callaux_t *bca, const char *ref); - int bcf_call_gap_prep(int n, int *n_plp, bam_pileup1_t **plp, int pos, bcf_callaux_t *bca, const char *ref); + int bcf_call_gap_prep(int n, int *n_plp, bam_pileup1_t **plp, int pos, bcf_callaux_t *bca, const char *ref, + const void *rghash); #ifdef __cplusplus } diff --git a/bam2bcf_indel.c b/bam2bcf_indel.c index de9df3e..2d0b5c6 100644 --- a/bam2bcf_indel.c +++ b/bam2bcf_indel.c @@ -1,13 +1,61 @@ #include #include +#include #include "bam.h" #include "bam2bcf.h" #include "ksort.h" #include "kaln.h" +#include "khash.h" +KHASH_SET_INIT_STR(rg) #define MINUS_CONST 0x10000000 #define INDEL_WINDOW_SIZE 50 +void *bcf_call_add_rg(void *_hash, const char *hdtext, const char *list) +{ + const char *s, *p, *q, *r, *t; + khash_t(rg) *hash; + if (list == 0 || hdtext == 0) return _hash; + if (_hash == 0) _hash = kh_init(rg); + hash = (khash_t(rg)*)_hash; + if ((s = strstr(hdtext, "@RG\t")) == 0) return hash; + do { + t = strstr(s + 4, "@RG\t"); // the next @RG + if ((p = strstr(s, "\tID:")) != 0) p += 4; + if ((q = strstr(s, "\tPL:")) != 0) q += 4; + if (p && q && (t == 0 || (p < t && q < t))) { // ID and PL are both present + int lp, lq; + char *x; + for (r = p; *r && *r != '\t' && *r != '\n'; ++r); lp = r - p; + for (r = q; *r && *r != '\t' && *r != '\n'; ++r); lq = r - q; + x = calloc((lp > lq? lp : lq) + 1, 1); + for (r = q; *r && *r != '\t' && *r != '\n'; ++r) x[r-q] = *r; + if (strstr(list, x)) { // insert ID to the hash table + khint_t k; + int ret; + for (r = p; *r && *r != '\t' && *r != '\n'; ++r) x[r-p] = *r; + x[r-p] = 0; + k = kh_get(rg, hash, x); + if (k == kh_end(hash)) k = kh_put(rg, hash, x, &ret); + else free(x); + } else free(x); + } + s = t; + } while (s); + return hash; +} + +void bcf_call_del_rghash(void *_hash) +{ + khint_t k; + khash_t(rg) *hash = (khash_t(rg)*)_hash; + if (hash == 0) return; + for (k = kh_begin(hash); k < kh_end(hash); ++k) + if (kh_exist(hash, k)) + free((char*)kh_key(hash, k)); + kh_destroy(rg, hash); +} + static int tpos2qpos(const bam1_core_t *c, const uint32_t *cigar, int32_t tpos, int is_left, int32_t *_tpos) { int k, x = c->pos, y = 0, last_y = 0; @@ -58,12 +106,30 @@ static inline int est_indelreg(int pos, const char *ref, int l, char *ins4) return max_i - pos; } -int bcf_call_gap_prep(int n, int *n_plp, bam_pileup1_t **plp, int pos, bcf_callaux_t *bca, const char *ref) +int bcf_call_gap_prep(int n, int *n_plp, bam_pileup1_t **plp, int pos, bcf_callaux_t *bca, const char *ref, + const void *rghash) { extern void ks_introsort_uint32_t(int, uint32_t*); int i, s, j, k, t, n_types, *types, max_rd_len, left, right, max_ins, *score, N, K, l_run, ref_type; char *inscns = 0, *ref2, *query; + khash_t(rg) *hash = (khash_t(rg)*)rghash; if (ref == 0 || bca == 0) return -1; + // mark filtered reads + if (rghash) { + N = 0; + for (s = N = 0; s < n; ++s) { + for (i = 0; i < n_plp[s]; ++i) { + bam_pileup1_t *p = plp[s] + i; + const uint8_t *rg = bam_aux_get(p->b, "RG"); + p->aux = 1; // filtered by default + if (rg) { + khint_t k = kh_get(rg, hash, (const char*)(rg + 1)); + if (k != kh_end(hash)) p->aux = 0, ++N; // not filtered + } + } + } + if (N == 0) return -1; // no reads left + } // determine if there is a gap for (s = N = 0; s < n; ++s) { for (i = 0; i < n_plp[s]; ++i) @@ -81,7 +147,7 @@ int bcf_call_gap_prep(int n, int *n_plp, bam_pileup1_t **plp, int pos, bcf_calla for (s = 0; s < n; ++s) { for (i = 0; i < n_plp[s]; ++i) { const bam_pileup1_t *p = plp[s] + i; - if (p->indel != 0) + if (p->indel != 0 && (rghash == 0 || p->aux == 0)) aux[m++] = MINUS_CONST + p->indel; j = bam_cigar2qlen(&p->b->core, bam1_cigar(p->b)); if (j > max_rd_len) max_rd_len = j; @@ -91,7 +157,9 @@ int bcf_call_gap_prep(int n, int *n_plp, bam_pileup1_t **plp, int pos, bcf_calla // squeeze out identical types for (i = 1, n_types = 1; i < m; ++i) if (aux[i] != aux[i-1]) ++n_types; - assert(n_types > 1); // there must at least one type of non-reference indel + if (n_types == 1) { // no indels + free(aux); return -1; + } types = (int*)calloc(n_types, sizeof(int)); t = 0; types[t++] = aux[0] - MINUS_CONST; diff --git a/bam_plcmd.c b/bam_plcmd.c index 54b8f91..0fbc930 100644 --- a/bam_plcmd.c +++ b/bam_plcmd.c @@ -535,7 +535,7 @@ int bam_pileup(int argc, char *argv[]) typedef struct { int max_mq, min_mq, flag, min_baseQ, capQ_thres, max_depth; - char *reg, *fn_pos; + char *reg, *fn_pos, *pl_list; faidx_t *fai; kh_64_t *hash; } mplp_conf_t; @@ -602,6 +602,8 @@ static void group_smpl(mplp_pileup_t *m, bam_sample_t *sm, kstring_t *buf, static int mpileup(mplp_conf_t *conf, int n, char **fn) { + extern void *bcf_call_add_rg(void *rghash, const char *hdtext, const char *list); + extern void bcf_call_del_rghash(void *rghash); mplp_aux_t **data; int i, tid, pos, *n_plp, beg0 = 0, end0 = 1u<<29, ref_len, ref_tid, max_depth; const bam_pileup1_t **plp; @@ -609,6 +611,7 @@ static int mpileup(mplp_conf_t *conf, int n, char **fn) bam_header_t *h = 0; char *ref; khash_t(64) *hash = 0; + void *rghash = 0; bcf_callaux_t *bca = 0; bcf_callret1_t *bcr = 0; @@ -638,6 +641,7 @@ static int mpileup(mplp_conf_t *conf, int n, char **fn) data[i]->fp = strcmp(fn[i], "-") == 0? bam_dopen(fileno(stdin), "r") : bam_open(fn[i], "r"); h_tmp = bam_header_read(data[i]->fp); bam_smpl_add(sm, fn[i], h_tmp->text); + rghash = bcf_call_add_rg(rghash, h_tmp->text, conf->pl_list); if (conf->reg) { int beg, end; bam_index_t *idx; @@ -693,6 +697,7 @@ static int mpileup(mplp_conf_t *conf, int n, char **fn) bcf_hdr_write(bp, bh); bca = bcf_call_init(-1., conf->min_baseQ); bcr = calloc(sm->n, sizeof(bcf_callret1_t)); + bca->rghash = rghash; } ref_tid = -1; ref = 0; iter = bam_mplp_init(n, mplp_func, (void**)data); @@ -731,7 +736,7 @@ static int mpileup(mplp_conf_t *conf, int n, char **fn) bcf_write(bp, bh, b); bcf_destroy(b); // call indels - if (bcf_call_gap_prep(gplp.n, gplp.n_plp, gplp.plp, pos, bca, ref) >= 0) { + if (bcf_call_gap_prep(gplp.n, gplp.n_plp, gplp.plp, pos, bca, ref, rghash) >= 0) { for (i = 0; i < gplp.n; ++i) bcf_call_glfgen(gplp.n_plp[i], gplp.plp[i], -1, bca, bcr + i); bcf_call_combine(gplp.n, bcr, -1, &bc); @@ -767,6 +772,7 @@ static int mpileup(mplp_conf_t *conf, int n, char **fn) bam_smpl_destroy(sm); free(buf.s); for (i = 0; i < gplp.n; ++i) free(gplp.plp[i]); free(gplp.plp); free(gplp.n_plp); free(gplp.m_plp); + bcf_call_del_rghash(rghash); if (hash) { // free the hash table khint_t k; for (k = kh_begin(hash); k < kh_end(hash); ++k) @@ -847,7 +853,7 @@ int bam_mpileup(int argc, char *argv[]) mplp.capQ_thres = 0; mplp.max_depth = 250; mplp.flag = MPLP_NO_ORPHAN | MPLP_REALN; - while ((c = getopt(argc, argv, "gf:r:l:M:q:Q:uaORC:BDSd:b:")) >= 0) { + while ((c = getopt(argc, argv, "gf:r:l:M:q:Q:uaORC:BDSd:b:P:")) >= 0) { switch (c) { case 'f': mplp.fai = fai_load(optarg); @@ -856,6 +862,7 @@ int bam_mpileup(int argc, char *argv[]) case 'd': mplp.max_depth = atoi(optarg); break; case 'r': mplp.reg = strdup(optarg); break; case 'l': mplp.fn_pos = strdup(optarg); break; + case 'P': mplp.pl_list = strdup(optarg); break; case 'g': mplp.flag |= MPLP_GLF; break; case 'u': mplp.flag |= MPLP_NO_COMP | MPLP_GLF; break; case 'a': mplp.flag |= MPLP_NO_ORPHAN | MPLP_REALN; break; @@ -882,6 +889,7 @@ int bam_mpileup(int argc, char *argv[]) fprintf(stderr, " -Q INT min base quality [%d]\n", mplp.min_baseQ); fprintf(stderr, " -q INT filter out alignment with MQ smaller than INT [%d]\n", mplp.min_mq); fprintf(stderr, " -d INT max per-sample depth [%d]\n", mplp.max_depth); + fprintf(stderr, " -P STR comma separated list of platforms for indels [all]\n"); fprintf(stderr, " -g generate BCF output\n"); fprintf(stderr, " -u do not compress BCF output\n"); fprintf(stderr, " -B disable BAQ computation\n"); @@ -900,7 +908,7 @@ int bam_mpileup(int argc, char *argv[]) } else mpileup(&mplp, argc - optind, argv + optind); - free(mplp.reg); + free(mplp.reg); free(mplp.pl_list); if (mplp.fai) fai_destroy(mplp.fai); return 0; } diff --git a/bamtk.c b/bamtk.c index e76a9d4..d2817b9 100644 --- a/bamtk.c +++ b/bamtk.c @@ -9,7 +9,7 @@ #endif #ifndef PACKAGE_VERSION -#define PACKAGE_VERSION "0.1.9-7 (r804)" +#define PACKAGE_VERSION "0.1.9-8 (r807)" #endif int bam_taf2baf(int argc, char *argv[]); diff --git a/bcftools/vcfutils.pl b/bcftools/vcfutils.pl index 0cf73a3..c1177ea 100755 --- a/bcftools/vcfutils.pl +++ b/bcftools/vcfutils.pl @@ -10,11 +10,10 @@ use Getopt::Std; exit; sub main { - my $version = '0.1.0'; &usage if (@ARGV < 1); my $command = shift(@ARGV); my %func = (subsam=>\&subsam, listsam=>\&listsam, fillac=>\&fillac, qstats=>\&qstats, varFilter=>\&varFilter, - hapmap2vcf=>\&hapmap2vcf, ucscsnp2vcf=>\&ucscsnp2vcf, filter4vcf=>\&filter4vcf, ldstats=>\&ldstats); + hapmap2vcf=>\&hapmap2vcf, ucscsnp2vcf=>\&ucscsnp2vcf, filter4vcf=>\&varFilter, ldstats=>\&ldstats); die("Unknown command \"$command\".\n") if (!defined($func{$command})); &{$func{$command}}; } @@ -199,8 +198,8 @@ Note: This command discards indels. Output: QUAL #non-indel #SNPs #transitions # } sub varFilter { - my %opts = (d=>1, D=>10000, a=>2, l=>10, Q=>10, s=>100, w=>10, p=>undef, F=>.001); - getopts('pd:D:l:Q:w:G:F:a:', \%opts); + my %opts = (d=>2, D=>10000, a=>2, W=>10, Q=>10, w=>10, p=>undef, 1=>1e-4, 2=>1e-100, 3=>0, 4=>1e-4); + getopts('pd:D:W:Q:w:a:1:2:3:4:', \%opts); die(qq/ Usage: vcfutils.pl varFilter [options] @@ -209,19 +208,28 @@ Options: -Q INT minimum RMS mapping quality for SNPs [$opts{Q}] -D INT maximum read depth [$opts{D}] -a INT minimum number of alternate bases [$opts{a}] -w INT SNP within INT bp around a gap to be filtered [$opts{w}] - -l INT window size for filtering adjacent gaps [$opts{l}] + -W INT window size for filtering adjacent gaps [$opts{W}] + -1 FLOAT min P-value for strand bias (given PV4) [$opts{1}] + -2 FLOAT min P-value for baseQ bias [$opts{2}] + -3 FLOAT min P-value for mapQ bias [$opts{3}] + -4 FLOAT min P-value for end distance bias [$opts{4}] -p print filtered variants + +Note: Some of the filters rely on annotations generated by SAMtools\/BCFtools. \n/) if (@ARGV == 0 && -t STDIN); # calculate the window size - my ($ol, $ow) = ($opts{l}, $opts{w}); + my ($ol, $ow) = ($opts{W}, $opts{w}); my $max_dist = $ol > $ow? $ol : $ow; # the core loop my @staging; # (indel_filtering_score, flt_tag, indel_span; chr, pos, ...) while (<>) { my @t = split; - next if (/^#/); + if (/^#/) { + print; next; + } next if ($t[4] eq '.'); # skip non-var sites + # check if the site is a SNP my $is_snp = 1; if (length($t[3]) > 1) { $is_snp = 0; @@ -238,7 +246,6 @@ Options: -Q INT minimum RMS mapping quality for SNPs [$opts{Q}] varFilter_aux(shift(@staging), $opts{p}); # calling a function is a bit slower, not much } my $flt = 0; - # parse annotations my ($dp, $mq, $dp_alt) = (-1, -1, -1); if ($t[7] =~ /DP=(\d+)/i) { @@ -258,6 +265,8 @@ Options: -Q INT minimum RMS mapping quality for SNPs [$opts{Q}] } $flt = 4 if ($dp_alt >= 0 && $dp_alt < $opts{a}); $flt = 1 if ($flt == 0 && $mq >= 0 && $mq < $opts{Q}); + $flt = 7 if ($flt == 0 && /PV4=([^,]+),([^,]+),([^,]+),([^,;\t]+)/ + && ($1<$opts{1} || $2<$opts{2} || $3<$opts{3} || $4<$opts{4})); # site dependent filters my ($rlen, $indel_score) = (0, -1); # $indel_score<0 for SNPs @@ -300,46 +309,7 @@ sub varFilter_aux { if ($first->[1] == 0) { print join("\t", @$first[3 .. @$first-1]), "\n"; } elsif ($is_print) { - print STDERR join("\t", substr("UQdDaGgsiX", $first->[1], 1), @$first[3 .. @$first-1]), "\n"; - } -} - -sub filter4vcf { - my %opts = (d=>3, D=>2000, 1=>1e-4, 2=>1e-100, 3=>0, 4=>1e-4, Q=>10, q=>3); - getopts('d:D:1:2:3:4:Q:q:', \%opts); - die(qq/ -Usage: vcfutils.pl filter4vcf [options] - -Options: -d INT min total depth (given DP or DP4) [$opts{d}] - -D INT max total depth [$opts{D}] - -q INT min SNP quality [$opts{q}] - -Q INT min RMS mapQ (given MQ) [$opts{Q}] - -1 FLOAT min P-value for strand bias (given PV4) [$opts{1}] - -2 FLOAT min P-value for baseQ bias [$opts{2}] - -3 FLOAT min P-value for mapQ bias [$opts{3}] - -4 FLOAT min P-value for end distance bias [$opts{4}]\n -/) if (@ARGV == 0 && -t STDIN); - - my %ts = (AG=>1, GA=>1, CT=>1, TC=>1); - - my @n = (0, 0); - while (<>) { - if (/^#/) { - print; - next; - } - next if (/PV4=([^,]+),([^,]+),([^,]+),([^,;\t]+)/ && ($1<$opts{1} || $2<$opts{2} || $3<$opts{3} || $4<$opts{4})); - my $depth = -1; - $depth = $1 if (/DP=(\d+)/); - $depth = $1+$2+$3+$4 if (/DP4=(\d+),(\d+),(\d+),(\d+)/); - next if ($depth > 0 && ($depth < $opts{d} || $depth > $opts{D})); - next if (/MQ=(\d+)/ && $1 < $opts{Q}); - my @t = split; - next if ($t[5] >= 0 && $t[5] < $opts{q}); - ++$n[0]; - my @s = split(',', $t[4]); - ++$n[1] if ($ts{$t[3].$s[0]}); - print; + print STDERR join("\t", substr("UQdDaGgP", $first->[1], 1), @$first[3 .. @$first-1]), "\n"; } } @@ -436,7 +406,6 @@ Command: subsam get a subset of samples fillac fill the allele count field qstats SNP stats stratified by QUAL varFilter filtering short variants - filter4vcf filtering VCFs produced by samtools+bcftools hapmap2vcf convert the hapmap format to VCF ucscsnp2vcf convert UCSC SNP SQL dump to VCF \n/); -- 2.39.2