add SVREADS in INFO

Author: yangao07

.gitignore

@@ -62,6 +62,7 @@ dkms.conf
 # vscode
 .lintr
 .vscode
+.claude
 
 # gprof
 gmon.out

Makefile

@@ -14,7 +14,7 @@ endif
 EXTRA_FLAGS = -Wall -Wno-misleading-indentation -Wno-unused-function #-Wno-unused-variable -Wno-alloc-size-larger-than
 
 # Define the version number
-VERSION=0.0.8
+VERSION=0.0.9
 # Get the Git commit hash
 GIT_COMMIT := $(shell git rev-parse --short HEAD 2> /dev/null)
 ifneq ($(GIT_COMMIT),)

README.md

@@ -10,25 +10,22 @@
 <!-- [![Published in Bioinformatics](https://img.shields.io/badge/Published%20in-Bioinformatics-blue.svg)](https://dx.doi.org/10.1093/bioinformatics/btaa963) -->
 <!-- [![GitHub Issues](https://img.shields.io/github/issues/yangao07/longcallD.svg?label=Issues)](https://github.com/yangao07/longcallD/issues) -->
 
-## Updates (release v0.0.8)
+## Updates (release v0.0.9)
 
-* Fix variant representation issues for overlapping variants
-* Revoke ultraLowMem mode of WFA to avoid potential variant representation issues
-* Fix output BAM errors
-* Output BAM includes all reads (even those not used for variant calling)
-* Fix a bug with --refine-aln
+* 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 INFO field of VCF
 <!-- * Add --STR  -->
 
 
 ## Getting Started
 ```sh
 # Download pre-built executables and test data (recommended)
 # Linux-x64
-wget https://github.com/yangao07/longcallD/releases/download/v0.0.8/longcallD-v0.0.8_x64-linux.tar.gz
-tar -zxvf longcallD-v0.0.8_x64-linux.tar.gz && cd longcallD-v0.0.8_x64-linux
+wget https://github.com/yangao07/longcallD/releases/download/v0.0.9/longcallD-v0.0.9_x64-linux.tar.gz
+tar -zxvf longcallD-v0.0.9_x64-linux.tar.gz && cd longcallD-v0.0.9_x64-linux
 # MacOS-arm64
-wget https://github.com/yangao07/longcallD/releases/download/v0.0.8/longcallD-v0.0.8_arm64-macos.tar.gz
-tar -zxvf longcallD-v0.0.8_arm64-macos.tar.gz && cd longcallD-v0.0.8_arm64-macos
+wget https://github.com/yangao07/longcallD/releases/download/v0.0.9/longcallD-v0.0.9_arm64-macos.tar.gz
+tar -zxvf longcallD-v0.0.9_arm64-macos.tar.gz && cd longcallD-v0.0.9_arm64-macos
 
 # PacBio HiFi reads
 ./longcallD call ./test_data/chr11_2M.fa ./test_data/HG002_chr11_hifi_test.bam --hifi > HG002_hifi_test.vcf
@@ -40,7 +37,7 @@ man ./longcallD.1
 ``` -->
 
 ## Table of Contents
-- [Updates (release v0.0.8)](#updates-release-v008)
+- [Updates (release v0.0.9)](#updates-release-v009)
 - [Getting Started](#getting-started)
 - [Table of Contents](#table-of-contents)
 - [Introduction](#introduction)
@@ -87,13 +84,13 @@ Providing the annotation sequence of common mobile elements, i.e., Alu/L1/SVA, u
 ### Pre-built executables (recommended)
 **Linux-x64**
 ```
-wget https://github.com/yangao07/longcallD/releases/download/v0.0.8/longcallD-v0.0.8_x64-linux.tar.gz
-tar -zxvf longcallD-v0.0.8_x64-linux.tar.gz
+wget https://github.com/yangao07/longcallD/releases/download/v0.0.9/longcallD-v0.0.9_x64-linux.tar.gz
+tar -zxvf longcallD-v0.0.9_x64-linux.tar.gz
 ```
 **MacOS-arm64**
 ```
-wget https://github.com/yangao07/longcallD/releases/download/v0.0.8/longcallD-v0.0.8_arm64-macos.tar.gz
-tar -zxvf longcallD-v0.0.8_arm64-macos.tar.gz
+wget https://github.com/yangao07/longcallD/releases/download/v0.0.9/longcallD-v0.0.9_arm64-macos.tar.gz
+tar -zxvf longcallD-v0.0.9_arm64-macos.tar.gz
 ```
 
 **Linux-arm64/macOS-x64**
@@ -109,9 +106,9 @@ conda install -c bioconda longcalld
 ### Build from source
 To compile longcallD from source, ensure you have **GCC/clang(9.0+)** and **zlib/libbz2/liblzma/libcurl** (for htslib) installed. 
 ```
-wget https://github.com/yangao07/longcallD/releases/download/v0.0.8/longcallD-v0.0.8.tar.gz
-tar -zxvf longcallD-v0.0.8.tar.gz
-cd longcallD-v0.0.8; make
+wget https://github.com/yangao07/longcallD/releases/download/v0.0.9/longcallD-v0.0.9.tar.gz
+tar -zxvf longcallD-v0.0.9.tar.gz
+cd longcallD-v0.0.9; make
 ```
 
 ## Usage

src/bam_utils.c

@@ -230,6 +230,7 @@ void update_read_var_profile_with_allele(int var_i, int allele_i, int alt_qi, re
 
 
 int get_digar_ave_qual(digar1_t *digar, const uint8_t *qual) {
+    if (digar->is_low_qual) return 0;
     if (digar->qi < 0) return 0;
     int q_start, q_end;
     if (digar->type == BAM_CDEL) {
@@ -1327,9 +1328,13 @@ int collect_digar_from_ref_seq(bam_chunk_t *chunk, int read_i, const struct call
     return skip == 1 ? -1 : 0;
 }
 
-int bam_chunk_init0(bam_chunk_t *chunk, int n_reads, int n_bam) {
+int bam_chunk_init0(bam_chunk_t *chunk, const struct call_var_opt_t *opt, int n_reads, int n_bam) {
     // input
     chunk->n_reads = 0; chunk->m_reads = n_reads; chunk->ordered_read_ids = (int*)malloc(n_reads * sizeof(int));
+    if (opt->output_sv_rnames || opt->output_somatic_sv_rnames) {
+        chunk->read_names = (char**)malloc(n_reads * sizeof(char*));
+        for (int i = 0; i < n_reads; i++) chunk->read_names[i] = NULL;
+    }
     chunk->n_bam = n_bam;
     chunk->n_up_ovlp_reads = (int*)calloc(n_bam, sizeof(int)); chunk->n_up_ovlp_skip_reads = (int*)calloc(n_bam, sizeof(int));
     chunk->up_ovlp_read_i = (int**)malloc(n_bam * sizeof(int*));
@@ -1368,9 +1373,10 @@ int bam_chunk_init0(bam_chunk_t *chunk, int n_reads, int n_bam) {
     return 0;
 }
 
-int bam_chunk_realloc(bam_chunk_t *chunk) {
+int bam_chunk_realloc(bam_chunk_t *chunk, const struct call_var_opt_t *opt) {
     int m_reads = chunk->m_reads * 2;
     chunk->reads = (bam1_t**)realloc(chunk->reads, m_reads * sizeof(bam1_t*));
+    if (opt->output_sv_rnames || opt->output_somatic_sv_rnames) chunk->read_names = (char**)realloc(chunk->read_names, m_reads * sizeof(char*));
     chunk->ordered_read_ids = (int*)realloc(chunk->ordered_read_ids, m_reads * sizeof(int));
     for (int i = 0; i < chunk->n_bam; ++i) {
         chunk->up_ovlp_read_i[i] = (int*)realloc(chunk->up_ovlp_read_i[i], m_reads * sizeof(int));
@@ -1480,6 +1486,12 @@ void bam_chunk_post_free(bam_chunk_t *chunk, const struct call_var_opt_t *opt) {
         free(chunk->reads);
     }
     free(chunk->ordered_read_ids);
+    // save read_var_profile for output_sv_rnames or output_somatic_sv_rnames, which will be used in make_variants
+    if (chunk->read_var_profile != NULL) free_read_var_profile(chunk->read_var_profile, chunk->n_reads);
+    if (opt->output_sv_rnames || opt->output_somatic_sv_rnames) {
+        for (int i = 0; i < chunk->n_reads; i++) free(chunk->read_names[i]); 
+        free(chunk->read_names);
+    }
 }
 
 // free variables that will not be used in stitch & make variants
@@ -1493,7 +1505,6 @@ void bam_chunk_mid_free(bam_chunk_t *chunk, const struct call_var_opt_t *opt) {
         } free(chunk->noisy_reg_to_reads);
     }
     if (chunk->noisy_reg_to_n_reads != NULL) free(chunk->noisy_reg_to_n_reads);
-    if (chunk->read_var_profile != NULL) free_read_var_profile(chunk->read_var_profile, chunk->n_reads);
     if (chunk->read_var_cr != NULL) cr_destroy(chunk->read_var_cr);
 }
 
@@ -1622,7 +1633,7 @@ int collect_ref_seq_bam_main(const struct call_var_pl_t *pl, struct call_var_io_
     call_var_opt_t *opt = pl->opt;
 
     int min_mapq = opt->min_mq;
-    bam_chunk_init0(chunk, 4096, io_aux->n_bam);
+    bam_chunk_init0(chunk, opt, 4096, io_aux->n_bam);
     chunk->reg_chunk_i = reg_chunk_i; chunk->reg_i = reg_i;
     chunk->tid = tid; chunk->tname = io_aux->headers[0]->target_name[tid]; chunk->reg_beg = reg_beg; chunk->reg_end = reg_end;
     hts_pos_t min_read_beg = reg_beg, max_read_end = reg_end;
@@ -1655,8 +1666,11 @@ int collect_ref_seq_bam_main(const struct call_var_pl_t *pl, struct call_var_io_
             if (chunk->reads[chunk->n_reads]->core.pos+1 < min_read_beg) min_read_beg = chunk->reads[chunk->n_reads]->core.pos+1;
             if (bam_endpos(chunk->reads[chunk->n_reads]) > max_read_end) max_read_end = bam_endpos(chunk->reads[chunk->n_reads]);
             // check if read is overlapping with next region
+            if (opt->output_sv_rnames || opt->output_somatic_sv_rnames) {
+                chunk->read_names[chunk->n_reads] = strdup(bam_get_qname(chunk->reads[chunk->n_reads]));
+            }
             chunk->n_reads++;
-            if (chunk->n_reads == chunk->m_reads) bam_chunk_realloc(chunk);
+            if (chunk->n_reads == chunk->m_reads) bam_chunk_realloc(chunk, opt);
         }
         bam_itr_destroy(iter);
     }

src/bam_utils.h

@@ -62,6 +62,7 @@ typedef struct bam_chunk_t {
     // should always be 1 region XXX
     cgranges_t *low_comp_cr; // tandem_rep_cr;
     int n_reads, m_reads; int *ordered_read_ids; // size: m_reads, for multiple input bams, merge sort reads by pos, end, NM, name
+    char **read_names; // size: m_reads, for output
     int n_bam; // for each input bam, record the number of region-overlapping reads
     int *n_up_ovlp_reads, *n_down_ovlp_reads, *n_up_ovlp_skip_reads, *n_down_ovlp_skip_reads; // number of reads overlapping with up/downstream bam chunk
     int **up_ovlp_read_i, **down_ovlp_read_i;

src/call_var_main.c

@@ -37,6 +37,8 @@ const struct option call_var_opt [] = {
     { "autosome", 0, NULL, 0},
     { "mosaic", 0, NULL, 0},
     { "refine-aln", 0, NULL, 0},
+    { "out-sv-rnames", 0, NULL, 0},
+    { "out-som-sv-rnames", 0, NULL, 0},
     // { "alpha", 1, NULL, 0},
     // { "beta", 1, NULL, 0},
     { "max-somvar", 1, NULL, 0},
@@ -206,7 +208,7 @@ call_var_opt_t *call_var_init_para(void) {
     opt->min_sv_len = LONGCALLD_MIN_SV_LEN;
     opt->min_tsd_len = LONGCALLD_MIN_TSD_LEN; opt->max_tsd_len = LONGCALLD_MAX_TSD_LEN;
     opt->min_polya_len = LONGCALLD_MIN_POLYA_LEN; opt->min_polya_ratio = LONGCALLD_MIN_POLYA_RATIO;
-    opt->output_sv_reads = 0;
+    opt->output_sv_rnames = 0; opt->output_somatic_sv_rnames = 0;
 
     initialize_lgamma_cache(opt);
     opt->te_seq_fn = NULL; opt->n_te_seqs = 0; opt->te_kmer_len = 15; opt->te_for_h = NULL; opt->te_rev_h = NULL;
@@ -297,6 +299,7 @@ void var_free(var_t *v) {
             if (v->vars[i].tsd_len > 0) {
                 free(v->vars[i].tsd_seq);
             }
+            if (v->vars[i].is_sv && v->vars[i].AD[1] > 0) free(v->vars[i].alt_read_i);
         }
         free(v->vars);
     }
@@ -793,7 +796,7 @@ static void *call_var_worker_pipeline(void *shared, int step, void *in) { // kt_
         int n_out_vars = 0, n_out_reads = 0;
         for (int i = 0; i < s->n_chunks; ++i) {
             bam_chunk_t *c = s->chunks + i;
-            n_out_vars += write_var_to_vcf(s->vars+i, pl->opt, c->tname);
+            n_out_vars += write_var_to_vcf(s->vars+i, pl->opt, c);
             if (pl->opt->out_aln_fp != NULL) n_out_reads += write_read_to_bam(c, pl->opt, pl->io_aux);
             var_free(s->vars + i);  // free output
         }
@@ -858,6 +861,8 @@ static void call_var_usage(void) {//main usage
     fprintf(stderr, "                          note: multiple input BAM/CRAM files will be merged in SAM/BAM/CRAM output\n");
     fprintf(stderr, "    --refine-aln          refine alignment in SAM/BAM/CRAM output\n");
     fprintf(stderr, "                          note: output SAM/BAM/CRAM may be unsorted when --refine-aln is set\n");
+    fprintf(stderr, "    --out-sv-rnames          output names of all reads supporting alternative allele in INFO field of VCF output for all SVs [False]\n");
+    fprintf(stderr, "    --out-som-sv-rnames      output names of all reads supporting alternative allele in INFO field of VCF output for mosaic/somatic SVs [False]\n");
     // fprintf(stderr, "\n");
     fprintf(stderr, "  Variant calling:\n");
     fprintf(stderr, "    -c --min-cov     INT  min. total read coverage for candidate variant [%d]\n", LONGCALLD_MIN_CAND_DP);
@@ -941,6 +946,8 @@ int call_var_main(int argc, char *argv[]) {
                         opt->somatic_win_max_vars = strtol(optarg, &s, 10); 
                         if (*s == ',') opt->somatic_win = strtol(s+1, &s, 10);
                     } else if (strcmp(call_var_opt[op_idx].name, "refine-aln") == 0) opt->refine_bam = 1;
+                    else if (strcmp(call_var_opt[op_idx].name, "out-sv-rnames") == 0) opt->output_sv_rnames = 1;
+                    else if (strcmp(call_var_opt[op_idx].name, "out-som-sv-rnames") == 0) opt->output_somatic_sv_rnames = 1;
                     break;
             case 's': opt->out_somatic = 1; break;
             case 'm': opt->out_methylation = 1; break;

src/call_var_main.h

@@ -109,14 +109,15 @@ 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;
+    uint8_t is_somatic, is_sv;
     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
 
     int n_alt_allele; // currently, only one alt allele
     uint8_t **alt_bases; int *alt_len;
     int DP, AD[2]; uint8_t GT[2]; // DP/AD/GT
     int QUAL, GQ, PL[6]; // phred-scaled, QUAL/FILTER/GQ/PL
+    int *alt_read_i;
 } var1_t;
 
 typedef struct var_t {
@@ -163,7 +164,7 @@ typedef struct call_var_opt_t {
     // Alu/L1/SVA sequences
     char *te_seq_fn; char **te_seq_names; int n_te_seqs;
     int te_kmer_len; kmer32_hash_t **te_for_h, **te_rev_h;
-    int output_sv_reads; // output supporting read IDs for each SV
+    int output_sv_rnames, output_somatic_sv_rnames; // output supporting read IDs for each SV
     // general
     // int max_ploidy;
     int pl_threads, n_threads;

src/collect_var.c

@@ -1350,7 +1350,7 @@ int make_variants(const call_var_opt_t *opt, bam_chunk_t *chunk, var_t **_var) {
     if (n_cand_vars <= 0) return 0;
     char *ref_seq = chunk->ref_seq; hts_pos_t ref_beg = chunk->ref_beg;
     hts_pos_t active_reg_beg = chunk->reg_beg, active_reg_end = chunk->reg_end;
-    cand_var_t *cand_vars = chunk->cand_vars;
+    cand_var_t *cand_vars = chunk->cand_vars; read_var_profile_t *p = chunk->read_var_profile;
     (*_var) = (var_t*)malloc(sizeof(var_t));
     var_t *var = *_var;
     var->n = 0; var->m = n_cand_vars;
@@ -1411,6 +1411,7 @@ int make_variants(const call_var_opt_t *opt, bam_chunk_t *chunk, var_t **_var) {
             var->vars[i].alt_bases = (uint8_t**)malloc(2 * sizeof(uint8_t*));
         }
         var->vars[i].n_alt_allele = 0;
+        var->vars[i].is_sv = 0;
         for (int hap=1; hap <= 2; ++hap) {
             int hap_alle = hap == 1 ? hap1_alle : hap2_alle;
             if (hap_alle != 0) { // alt allele
@@ -1433,12 +1434,36 @@ int make_variants(const call_var_opt_t *opt, bam_chunk_t *chunk, var_t **_var) {
                     }
                 }
                 var->vars[i].alt_len[var->vars[i].n_alt_allele] = alt_len;
+                if (abs(alt_len - var->vars[i].ref_len) >= opt->min_sv_len) var->vars[i].is_sv = 1;
                 var->vars[i].GT[hap-1] = ++var->vars[i].n_alt_allele;
                 if (is_hom) hom_alt_is_set = 1;
             } else var->vars[i].GT[hap-1] = 0;
         }
         var->vars[i].DP = cand_vars[cand_i].total_cov;
         for (int j = 0; j < cand_vars[cand_i].n_uniq_alles; ++j) var->vars[i].AD[j] = cand_vars[cand_i].alle_covs[j];
+        if (var->vars[i].is_sv && var->vars[i].AD[1] > 0) { // collect alt_read_i from read_var_profile
+            var->vars[i].alt_read_i = (int*)malloc(var->vars[i].AD[1] * sizeof(int));
+            int alt_read_i_idx = 0;
+            for (int k = 0; k < chunk->n_reads; ++k) {
+                int read_i = chunk->ordered_read_ids[k];
+                if (chunk->is_skipped[read_i]) continue;
+                read_var_profile_t *p1 = p + read_i;
+                if (p1->start_var_idx < 0 || p1->end_var_idx < 0) continue;
+                if (cand_i < p1->start_var_idx || cand_i > p1->end_var_idx) continue;
+                int allele = p1->alleles[cand_i - p1->start_var_idx];
+                if (allele == 1) {
+                    // fprintf(stderr, "Collect alt_read_i for var: %s:%" PRIi64 " read: %s allele: %d\n", chunk->tname, var->vars[i].pos, chunk->read_names[read_i], allele);
+                    if (alt_read_i_idx >= var->vars[i].AD[1]) {
+                        _err_error_exit("Error: alt_read_i_idx: %d exceeds AD[1]: %d for var: %s:%" PRIi64 " %d-%c-%d\n", alt_read_i_idx, var->vars[i].AD[1], chunk->tname, var->vars[i].pos, var->vars[i].ref_len, BAM_CIGAR_STR[var->vars[i].type], var->vars[i].alt_len[0]);
+                    }
+                    var->vars[i].alt_read_i[alt_read_i_idx++] = read_i;
+                }
+            }
+            if (alt_read_i_idx != var->vars[i].AD[1]) {
+                _err_error("Error: alt_read_i_idx: %d does not match AD[1]: %d for var: %s:%" PRIi64 " %d-%c-%d\n", alt_read_i_idx, var->vars[i].AD[1], chunk->tname, var->vars[i].pos, var->vars[i].ref_len, BAM_CIGAR_STR[var->vars[i].type], var->vars[i].alt_len[0]);
+                var->vars[i].AD[1] = alt_read_i_idx;
+            }
+        } else var->vars[i].alt_read_i = NULL;
         var->vars[i].QUAL = cal_var_QUAL1(var->vars[i].AD[0], var->vars[i].AD[1], opt->log_p, opt->log_1p, opt->max_qual);
         var->vars[i].GQ = cal_sample_GQ(var->vars[i].AD[0], var->vars[i].AD[1], opt->log_p, opt->log_1p, opt->log_2, opt->max_gq);
         i++;

src/main.c

@@ -9,7 +9,7 @@
 const char PROG[20] = "longcallD";
 const char DESCRIP[100] = "local-haplotagging-based small and structural variant calling";
 #ifndef LONGCALLD_VERSION
-const char LONGCALLD_VERSION[20] = "0.0.8";
+const char LONGCALLD_VERSION[20] = "0.0.9";
 #endif
 const char CONTACT[50] = "yangao@ds.dfci.harvard.edu";
 int LONGCALLD_VERBOSE = 0;