Skip to content

Commit 81cfd58

Browse files
yangao07yangao07
committed
fix corrupted VCF (#10)
1 parent 687fd44 commit 81cfd58

5 files changed

Lines changed: 44 additions & 62 deletions

File tree

‎src/bam_utils.c‎

Lines changed: 6 additions & 32 deletions
Original file line numberDiff line numberDiff line change
@@ -685,13 +685,6 @@ int collect_digar_from_eqx_cigar(bam_chunk_t *chunk, int read_i, const struct ca
685685
digar->qual = (uint8_t*)malloc(qlen * sizeof(uint8_t));
686686
for (int i = 0; i < qlen; ++i) {
687687
digar->qual[i] = bam_get_qual(read)[i];
688-
if (digar->qual[i] > 99) {
689-
fprintf(stderr, "Warning: read %s has base quality %d > 99, set to 99\n", bam_get_qname(read), digar->qual[i]);
690-
digar->qual[i] = 99;
691-
} else if (digar->qual[i] < 0) {
692-
fprintf(stderr, "Warning: read %s has invalid base quality %d\n", bam_get_qname(read), digar->qual[i]);
693-
digar->qual[i] = 0;
694-
}
695688
chunk->qual_counts[digar->qual[i]]++;
696689
}
697690
int _n_digar = 0, _m_digar = 2 * n_cigar; digar1_t *_digars = (digar1_t*)malloc(_m_digar * sizeof(digar1_t));
@@ -843,13 +836,6 @@ int collect_digar_from_cs_tag(bam_chunk_t *chunk, int read_i, const struct call_
843836
digar->qual = (uint8_t*)malloc(qlen * sizeof(uint8_t));
844837
for (int i = 0; i < qlen; ++i) {
845838
digar->qual[i] = bam_get_qual(read)[i];
846-
if (digar->qual[i] > 99) {
847-
fprintf(stderr, "Warning: read %s has base quality %d > 99, set to 99\n", bam_get_qname(read), digar->qual[i]);
848-
digar->qual[i] = 99;
849-
} else if (digar->qual[i] < 0) {
850-
fprintf(stderr, "Warning: read %s has invalid base quality %d\n", bam_get_qname(read), digar->qual[i]);
851-
digar->qual[i] = 0;
852-
}
853839
chunk->qual_counts[digar->qual[i]]++;
854840
}
855841
int _n_digar = 0, _m_digar = 2 * n_cigar; digar1_t *_digars = (digar1_t*)malloc(_m_digar * sizeof(digar1_t));
@@ -1015,13 +1001,6 @@ int collect_digar_from_MD_tag(bam_chunk_t *chunk, int read_i, const struct call_
10151001
digar->qual = (uint8_t*)malloc(qlen * sizeof(uint8_t));
10161002
for (int i = 0; i < qlen; ++i) {
10171003
digar->qual[i] = bam_get_qual(read)[i];
1018-
if (digar->qual[i] > 99) {
1019-
fprintf(stderr, "Warning: read %s has base quality %d > 99, set to 99\n", bam_get_qname(read), digar->qual[i]);
1020-
digar->qual[i] = 99;
1021-
} else if (digar->qual[i] < 0) {
1022-
fprintf(stderr, "Warning: read %s has invalid base quality %d\n", bam_get_qname(read), digar->qual[i]);
1023-
digar->qual[i] = 0;
1024-
}
10251004
chunk->qual_counts[digar->qual[i]]++;
10261005
}
10271006
int _n_digar = 0, _m_digar = 2 * n_cigar; digar1_t *_digars = (digar1_t*)malloc(_m_digar * sizeof(digar1_t));
@@ -1199,18 +1178,8 @@ int collect_digar_from_ref_seq(bam_chunk_t *chunk, int read_i, const struct call
11991178
digar->qual = (uint8_t*)malloc(qlen);
12001179
for (int i = 0; i < qlen; ++i) {
12011180
digar->qual[i] = bam_get_qual(read)[i];
1202-
if (digar->qual[i] > 99) {
1203-
fprintf(stderr, "Warning: read %s has base quality %d > 99, set to 99\n", bam_get_qname(read), digar->qual[i]);
1204-
digar->qual[i] = 99;
1205-
} else if (digar->qual[i] < 0) {
1206-
fprintf(stderr, "Warning: read %s has invalid base quality %d\n", bam_get_qname(read), digar->qual[i]);
1207-
digar->qual[i] = 0;
1208-
}
12091181
chunk->qual_counts[digar->qual[i]]++;
12101182
}
1211-
if (strcmp("SRR25029837.8497337", bam_get_qname(read)) == 0) {
1212-
fprintf(stderr, "DBG: ref_seq: %.*s\n", (int)(ref_end-ref_beg+1), ref_seq);
1213-
}
12141183
int _n_digar = 0, _m_digar = 2 * n_cigar; digar1_t *_digars = (digar1_t*)malloc(_m_digar * sizeof(digar1_t));
12151184
int rlen = bam_cigar2rlen(n_cigar, cigar); int tlen = chunk->whole_ref_len;
12161185
xid_queue_t *q = init_xid_queue(rlen, max_s, win);
@@ -1356,7 +1325,7 @@ int bam_chunk_init0(bam_chunk_t *chunk, int n_reads) {
13561325
chunk->ref_seq = NULL;
13571326
chunk->low_comp_cr = NULL;
13581327
// intermediate
1359-
chunk->qual_counts = (int*)calloc(100, sizeof(int));
1328+
chunk->qual_counts = (int*)calloc(256, sizeof(int));
13601329
chunk->is_skipped = (uint8_t*)calloc(n_reads, sizeof(uint8_t));
13611330
chunk->n_clean_agree_snps = (int*)malloc(n_reads * sizeof(int));
13621331
chunk->n_clean_conflict_snps = (int*)malloc(n_reads * sizeof(int));
@@ -1520,6 +1489,11 @@ void get_bam_chunk_reg_ref_seq0(faidx_t *fai, bam_chunk_t *chunk, hts_pos_t beg,
15201489

15211490
int len;
15221491
chunk->ref_seq = faidx_fetch_seq(fai, chunk->tname, ref_beg, ref_end, &len); // ref_beg & ref_end: 0-based
1492+
if (chunk->ref_seq == NULL || len <= 0) {
1493+
_err_error("Failed to fetch reference sequence for %s:%" PRId64 "-%" PRId64 " from the fasta file.\n", chunk->tname, beg, end);
1494+
_err_error_exit("Please make sure the reference genome and the alignment file match.\n");
1495+
1496+
}
15231497
chunk->ref_beg = ref_beg+1; // 1-based
15241498
chunk->ref_end = ref_beg+len; // 1-based
15251499

‎src/call_var_main.c‎

Lines changed: 3 additions & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -497,6 +497,7 @@ static void collect_regions(call_var_pl_t *pl, call_var_opt_t *opt, int n_region
497497
}
498498
// check if region(s) is provided
499499
if (n_regions > 0) {
500+
opt->only_autosome = 0, opt->only_autosome_XY=0;
500501
if (opt->reg_bed_fn != NULL) {
501502
_err_error("Both region(s) and region bed file are provided. Only region(s) will be used.\n");
502503
}
@@ -511,6 +512,7 @@ static void collect_regions(call_var_pl_t *pl, call_var_opt_t *opt, int n_region
511512
return collect_regions_from_region_list(opt, pl, iter, n_regions, regions);
512513
}
513514
} else if (opt->reg_bed_fn != NULL) { // check if region bed file is provided
515+
opt->only_autosome = 0, opt->only_autosome_XY=0;
514516
collect_regions_from_bed_file(opt, pl);
515517
if (pl->n_reg_chunks == 0) {
516518
_err_error("Failed to parse provided region bed file: %s\n", opt->reg_bed_fn);
@@ -677,7 +679,7 @@ static void call_var_usage(void) {//main usage
677679
fprintf(stderr, "\n");
678680

679681
fprintf(stderr, "Options:\n");
680-
fprintf(stderr, " Intput:\n");
682+
fprintf(stderr, " Input:\n");
681683
fprintf(stderr, " --hifi HiFi reads [Default]\n");
682684
fprintf(stderr, " --ont ONT reads [False]\n");
683685
fprintf(stderr, " --region-file FILE region bed file [NULL]\n");

‎src/call_var_main.h‎

Lines changed: 1 addition & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -68,7 +68,7 @@
6868
#define LONGCALLD_MIN_SOMATIC_TE_ALT_DP 1 // >= 1 read supporting somatic TE variant
6969
// fisher is NOT used for PacBio-HiFi
7070
#define LONGCALLD_MIN_SOMATIC_FISHER_PVAL 0.05
71-
#define LONGCALLD_STRAND_BIAS_PVAL_ONT 0.05
71+
#define LONGCALLD_STRAND_BIAS_PVAL_ONT 0.01
7272
// BLT50: 1/10/3 works better
7373
// #define LONGCALLD_SOMATIC_BETA_ALPHA 2 // beta prior for somatic variant calling
7474
// #define LONGCALLD_SOMATIC_BETA_BETA 10 // beta prior for somatic variant calling

‎src/collect_var.c‎

Lines changed: 12 additions & 5 deletions
Original file line numberDiff line numberDiff line change
@@ -265,12 +265,18 @@ int collect_cand_vars(const call_var_opt_t *opt, bam_chunk_t *chunk, int n_var_s
265265
int var_is_strand_bias(cand_var_t *var, const call_var_opt_t *opt) {
266266
int for_alt_cov = var->strand_to_alle_covs[0][1]; // forward strand alt coverage
267267
int rev_alt_cov = var->strand_to_alle_covs[1][1]; // reverse strand alt coverage
268+
int for_ref_cov = var->strand_to_alle_covs[0][0]; // forward strand ref coverage
269+
int rev_ref_cov = var->strand_to_alle_covs[1][0]; // reverse strand ref coverage
268270
int expected_alt_cov = (for_alt_cov + rev_alt_cov) / 2; // expected alt coverage
269271
if (expected_alt_cov == 0) return 0; // no alt coverage, no strand bias
270272
// check if strand bias is significant
271273
float fisher_p = fisher_exact_test(for_alt_cov, rev_alt_cov, expected_alt_cov, expected_alt_cov, opt);
272-
if (fisher_p < opt->strand_bias_pval) return 1; // significant strand bias
273-
else return 0; // no significant strand bias
274+
// float fisher_p = fisher_exact_test(for_alt_cov, rev_alt_cov, for_ref_cov, rev_ref_cov, opt);
275+
if (fisher_p < opt->strand_bias_pval) {
276+
// fprintf(stderr, "ONT-StrandBias: %" PRId64 " for_alt_cov=%d, rev_alt_cov=%d, for_ref_cov=%d, rev_ref_cov=%d, fisher p-value=%.5f\n",
277+
// var->pos, for_alt_cov, rev_alt_cov, for_ref_cov, rev_ref_cov, fisher_p);
278+
return 1; // significant strand bias
279+
} else return 0; // no significant strand bias
274280
// strand bias: 1) >=3 vs 0, 2) >= 3 folds
275281
// int min_fold = 3;
276282
// if (var->alle_covs[1] < min_fold) return 0;
@@ -1005,10 +1011,11 @@ void collect_digars_from_bam(bam_chunk_t *chunk, const struct call_var_pl_t *pl)
10051011
if (ret < 0) chunk->is_skipped[i] = BAM_RECORD_WRONG_MAP;
10061012
}
10071013
// print chunk->qual_counts
1008-
int valid_quals[100], n_valid_quals = 0; // XXX MAX_QUAL = 99
1014+
int n_all_quals = 256;
1015+
int valid_quals[256], n_valid_quals = 0; // XXX MAX_QUAL = 255
10091016
int64_t n_total_counts = 0;
1010-
for (int i = 0; i < 100; ++i) n_total_counts += chunk->qual_counts[i];
1011-
for (int i = 0; i < 100; ++i) {
1017+
for (int i = 0; i < n_all_quals; ++i) n_total_counts += chunk->qual_counts[i];
1018+
for (int i = 0; i < n_all_quals; ++i) {
10121019
if (chunk->qual_counts[i] <= 0) continue; // skip 0 counts
10131020
if (chunk->qual_counts[i] >= 0.001 * n_total_counts) valid_quals[n_valid_quals++] = i;
10141021
// if (LONGCALLD_VERBOSE >= 0) {

‎src/vcf_utils.c‎

Lines changed: 22 additions & 23 deletions
Original file line numberDiff line numberDiff line change
@@ -20,7 +20,7 @@ void write_vcf_header(bam_hdr_t *hdr, struct call_var_opt_t *opt) {
2020
bcf_hdr_t *vcf_hdr = bcf_hdr_init("w");
2121
if (!vcf_hdr) _err_error_exit("Could not allocate VCF header.\n");
2222
// File format
23-
bcf_hdr_append(vcf_hdr, "##fileformat=VCFv4.3");
23+
// bcf_hdr_append(vcf_hdr, "##fileformat=VCFv4.3");
2424

2525
// Get current date
2626
time_t t = time(NULL);
@@ -144,14 +144,14 @@ int write_var_to_vcf(var_t *vars, const struct call_var_opt_t *opt, char *chrom)
144144
buffer = (char*)realloc(buffer, buf_m * sizeof(char));
145145
}
146146
for (int j = 0; j < var.ref_len; j++)
147-
len += snprintf(buffer + len, sizeof(buffer) - len, "%c", "ACGTN"[var.ref_bases[j]]);
147+
len += snprintf(buffer + len, buf_m - len, "%c", "ACGTN"[var.ref_bases[j]]);
148148

149149
// Write ALT
150-
len += snprintf(buffer + len, sizeof(buffer) - len, "\t");
150+
len += snprintf(buffer + len, buf_m - len, "\t");
151151
for (int j = 0; j < var.n_alt_allele; j++) {
152152
for (int k = 0; k < var.alt_len[j]; k++)
153-
len += snprintf(buffer + len, sizeof(buffer) - len, "%c", "ACGTN"[var.alt_bases[j][k]]);
154-
if (j < var.n_alt_allele - 1) len += snprintf(buffer + len, sizeof(buffer) - len, ",");
153+
len += snprintf(buffer + len, buf_m - len, "%c", "ACGTN"[var.alt_bases[j][k]]);
154+
if (j < var.n_alt_allele - 1) len += snprintf(buffer + len, buf_m - len, ",");
155155
}
156156

157157
// Structural Variant (SV) annotation
@@ -174,21 +174,21 @@ int write_var_to_vcf(var_t *vars, const struct call_var_opt_t *opt, char *chrom)
174174
}
175175

176176
// Write QUAL, FILTER, INFO
177-
len += snprintf(buffer + len, sizeof(buffer) - len, "\t%d\tPASS\t", var.QUAL);
178-
if (var.is_somatic) len += snprintf(buffer + len, sizeof(buffer) - len, "SOMATIC;");
179-
if (var.te_seq_i >= 0) len += snprintf(buffer + len, sizeof(buffer) - len, "MEI;");
180-
len += snprintf(buffer + len, sizeof(buffer) - len, "END=%" PRId64 "", var.pos + var.ref_len - 1);
177+
len += snprintf(buffer + len, buf_m - len, "\t%d\tPASS\t", var.QUAL);
178+
if (var.is_somatic) len += snprintf(buffer + len, buf_m - len, "SOMATIC;");
179+
if (var.te_seq_i >= 0) len += snprintf(buffer + len, buf_m - len, "MEI;");
180+
len += snprintf(buffer + len, buf_m - len, "END=%" PRId64 "", var.pos + var.ref_len - 1);
181181
if (is_sv) {
182-
len += snprintf(buffer + len, sizeof(buffer) - len, ";%s;%s", SVTYPE, SVLEN);
182+
len += snprintf(buffer + len, buf_m - len, ";%s;%s", SVTYPE, SVLEN);
183183
if (var.tsd_len > 0) {
184-
len += snprintf(buffer + len, sizeof(buffer) - len, ";TSD=");
185-
for (int i = 0; i < var.tsd_len; ++i) len += snprintf(buffer + len, sizeof(buffer) - len, "%c", "ACGTN"[var.tsd_seq[i]]);
186-
len += snprintf(buffer + len, sizeof(buffer) - len, ";TSDLEN=%d;POLYALEN=%d;TSDPOS1=%" PRId64 "", var.tsd_len, var.polya_len, var.tsd_pos1);
187-
if (var.tsd_pos2 > 0) len += snprintf(buffer + len, sizeof(buffer) - len, ";TSDPOS2=%" PRId64 "", var.tsd_pos2);
184+
len += snprintf(buffer + len, buf_m - len, ";TSD=");
185+
for (int i = 0; i < var.tsd_len; ++i) len += snprintf(buffer + len, buf_m - len, "%c", "ACGTN"[var.tsd_seq[i]]);
186+
len += snprintf(buffer + len, buf_m - len, ";TSDLEN=%d;POLYALEN=%d;TSDPOS1=%" PRId64 "", var.tsd_len, var.polya_len, var.tsd_pos1);
187+
if (var.tsd_pos2 > 0) len += snprintf(buffer + len, buf_m - len, ";TSDPOS2=%" PRId64 "", var.tsd_pos2);
188188
}
189-
if (var.te_seq_i >= 0) len += snprintf(buffer + len, sizeof(buffer) - len, ";REPNAME=%c%s", "+-"[var.te_is_rev], opt->te_seq_names[var.te_seq_i]);
189+
if (var.te_seq_i >= 0) len += snprintf(buffer + len, buf_m - len, ";REPNAME=%c%s", "+-"[var.te_is_rev], opt->te_seq_names[var.te_seq_i]);
190190
}
191-
len += snprintf(buffer + len, sizeof(buffer) - len, "\t");
191+
len += snprintf(buffer + len, buf_m - len, "\t");
192192

193193
// Write FORMAT and Genotype Data
194194
int gt1 = var.GT[0], gt2 = var.GT[1];
@@ -198,19 +198,19 @@ int write_var_to_vcf(var_t *vars, const struct call_var_opt_t *opt, char *chrom)
198198
if (gt1 > gt2) { int tmp = gt1; gt1 = gt2; gt2 = tmp; }
199199
}
200200
if (is_hom || var.PS == 0)
201-
len += snprintf(buffer + len, sizeof(buffer) - len, "GT:DP:AD:GQ\t%d%c%d:%d:", gt1, gt_seperator, gt2, var.DP);
201+
len += snprintf(buffer + len, buf_m - len, "GT:DP:AD:GQ\t%d%c%d:%d:", gt1, gt_seperator, gt2, var.DP);
202202
else
203-
len += snprintf(buffer + len, sizeof(buffer) - len, "GT:DP:AD:GQ:PS\t%d%c%d:%d:", gt1, gt_seperator, gt2, var.DP);
203+
len += snprintf(buffer + len, buf_m - len, "GT:DP:AD:GQ:PS\t%d%c%d:%d:", gt1, gt_seperator, gt2, var.DP);
204204

205205
for (int j = 0; j < 1 + var.n_alt_allele; j++) {
206-
if (j > 0) len += snprintf(buffer + len, sizeof(buffer) - len, ",");
207-
len += snprintf(buffer + len, sizeof(buffer) - len, "%d", var.AD[j]);
206+
if (j > 0) len += snprintf(buffer + len, buf_m - len, ",");
207+
len += snprintf(buffer + len, buf_m - len, "%d", var.AD[j]);
208208
}
209209

210210
if (is_hom || var.PS == 0)
211-
len += snprintf(buffer + len, sizeof(buffer) - len, ":%d\n", var.GQ);
211+
len += snprintf(buffer + len, buf_m - len, ":%d\n", var.GQ);
212212
else
213-
len += snprintf(buffer + len, sizeof(buffer) - len, ":%d:%" PRId64 "\n", var.GQ, var.PS);
213+
len += snprintf(buffer + len, buf_m - len, ":%d:%" PRId64 "\n", var.GQ, var.PS);
214214

215215
// Write to htsFile
216216
if (out_vcf->format.compression!=no_compression) {
@@ -222,7 +222,6 @@ int write_var_to_vcf(var_t *vars, const struct call_var_opt_t *opt, char *chrom)
222222
_err_error_exit("Could not write to VCF file.\n");
223223
}
224224
}
225-
// fprintf(stdout, "%s", buffer);
226225
n_output_vars++;
227226
}
228227
free(buffer);

0 commit comments

Comments
 (0)