add CLEAN tag

Author: yangao07

README.md

@@ -14,7 +14,7 @@
 
 * Fix variant calling close to start/end boundaries of chromosomes
 * Add `--out-sv-rnames` and `--out-som-sv-rnames` to output SV-supporing read names (tag: `SVREADS`) in FORMAT field of VCF
-<!-- * Add --STR  -->
+* Add `CLEAN` INFO tag for variants in clean regions, i.e., SNPs or simple small indels (≤5bp) in non-repetitive regions which are generally more reliable, to help with downstream filtering and benchmarking
 
 
 ## Getting Started

src/call_var_main.h

@@ -109,7 +109,7 @@ typedef struct {
     uint8_t type; // BAM_CINS/BAM_CDEL/BAM_CDIFF
     hts_pos_t pos, PS; // phase set
     uint8_t *ref_bases; int ref_len;
-    uint8_t is_somatic, is_sv;
+    uint8_t is_somatic, is_sv, is_clean; // clean: SNP or simple indel in non-repetitive region
     uint8_t *tsd_seq; int tsd_len, polya_len; hts_pos_t tsd_pos1, tsd_pos2; // target site duplication, 2 TSDs for DEL
     int te_seq_i, te_is_rev; // char *rep_name, *rep_family, *rep_class; use te_seq_i to retrieve TE sequence info
 

src/collect_var.c

@@ -1412,6 +1412,7 @@ int make_variants(const call_var_opt_t *opt, bam_chunk_t *chunk, var_t **_var) {
         }
         var->vars[i].n_alt_allele = 0;
         var->vars[i].is_sv = 0;
+        var->vars[i].is_clean = (chunk->var_i_to_cate[cand_i] & LONGCALLD_CAND_GERMLINE_CLEAN_VAR_CATE) != 0;
         for (int hap=1; hap <= 2; ++hap) {
             int hap_alle = hap == 1 ? hap1_alle : hap2_alle;
             if (hap_alle != 0) { // alt allele

src/vcf_utils.c

@@ -58,6 +58,7 @@ void write_vcf_header(bam_hdr_t *hdr, struct call_var_opt_t *opt) {
     // INFO fields
     bcf_hdr_append(vcf_hdr, "##INFO=<ID=END,Number=1,Type=Integer,Description=\"End position of the variant described in this record\">");
     bcf_hdr_append(vcf_hdr, "##INFO=<ID=SOMATIC,Number=0,Type=Flag,Description=\"Somatic/mosaic variant\">");
+    bcf_hdr_append(vcf_hdr, "##INFO=<ID=CLEAN,Number=0,Type=Flag,Description=\"Clean-region variant (SNP or simple indel in non-repetitive region)\">");
     bcf_hdr_append(vcf_hdr, "##INFO=<ID=SVTYPE,Number=1,Type=String,Description=\"Type of structural variant\">");
     bcf_hdr_append(vcf_hdr, "##INFO=<ID=SVLEN,Number=A,Type=Integer,Description=\"Difference in length between REF and ALT alleles\">");
     // TSD info
@@ -177,6 +178,7 @@ int write_var_to_vcf(var_t *vars, const struct call_var_opt_t *opt, bam_chunk_t
 
         // Write QUAL, FILTER, INFO
         len += snprintf(buffer + len, buf_m - len, "\t%d\tPASS\t", var.QUAL);
+        if (var.is_clean) len += snprintf(buffer + len, buf_m - len, "CLEAN;");
         if (var.is_somatic) len += snprintf(buffer + len, buf_m - len, "SOMATIC;");
         if (var.te_seq_i >= 0) len += snprintf(buffer + len, buf_m - len, "MEI;");
         len += snprintf(buffer + len, buf_m - len, "END=%" PRId64 "", var.pos + var.ref_len - 1);