include clip reads, def:p1a2

Author: yangao07

abPOA

@@ -13,8 +13,8 @@
 #define LONGCALLD_GAP_LEFT_ALN 1 // default: put gap at the left-most position
 #define LONGCALLD_GAP_RIGHT_ALN 2
 
-#define LONGCALLD_EXT_ALN_RIGHT 1
-#define LONGCALLD_EXT_ALN_LEFT  2
+#define LONGCALLD_EXT_ALN_LEFT_TO_RIGHT 1
+#define LONGCALLD_EXT_ALN_RIGHT_TO_LEFT  2
 
 #ifdef __cplusplus
 extern "C" {
@@ -23,13 +23,6 @@ extern "C" {
 struct bam_chunk_t;
 struct aln_str_t;
 
-int test_wfa(char *pattern, char *text);
-int test_abpoa(uint8_t **bseqs, int n_seqs, int *seq_lens);
-int test_edlib(char *pattern, char *text);
-int test_ksw2(char *pattern, char *text, const call_var_opt_t *opt);
-
-int wfa_heuristic_aln(uint8_t *pattern, int plen, uint8_t *text, int tlen, int *n_eq, int *n_xid);
-
 // int wfa_aln(int gap_pos, char *pattern, int plen, char *text, int tlen, uint32_t **cigar_buf);
 int end2end_aln(const call_var_opt_t *opt, char *pattern, int plen, uint8_t *text, int tlen, uint32_t **cigar_buf);
 int collect_reg_ref_cseq(bam_chunk_t *chunk, hts_pos_t reg_beg, hts_pos_t reg_end, char **ref_cseq);

src/align.c

@@ -105,6 +105,7 @@ call_var_opt_t *call_var_init_para(void) {
     opt->min_noisy_reg_ratio = LONGCALLD_NOISY_REG_RATIO;
     opt->max_noisy_frac_per_read = LONGCALLD_MAX_NOISY_FRAC_PER_READ;
     opt->min_hap_full_reads = LONGCALLD_MIN_HAP_FULL_READS;
+    opt->min_hap_reads = LONGCALLD_MIN_HAP_READS;
     opt->min_no_hap_full_reads = LONGCALLD_MIN_NO_HAP_FULL_READS;
 
     opt->min_somatic_af = LONGCALLD_MIN_SOMATIC_AF;
@@ -123,6 +124,8 @@ call_var_opt_t *call_var_init_para(void) {
     opt->pl_threads = MIN_OF_TWO(CALL_VAR_PL_THREAD_N, get_num_processors());
     opt->n_threads = MIN_OF_TWO(CALL_VAR_THREAD_N, get_num_processors());
 
+    opt->p_error = 0.001; opt->log_p = -3.0; opt->log_1p = log10(1-opt->p_error); opt->log_2 = 0.301023;
+    opt->max_gq = 60; opt->max_qual = 60;
     opt->out_vcf = NULL; opt->no_vcf_header = 0; opt->out_amb_base = 0;
     opt->out_bam = NULL; opt->no_bam_header = 0;
     // opt->verbose = 0;

src/align.h

@@ -36,9 +36,10 @@
 #define LONGCALLD_NOISY_REG_READS 5 // >= 5 reads supporting noisy region
 #define LONGCALLD_NOISY_REG_RATIO 0.20 // >= 25% reads supporting noisy region
 
-#define LONGCALLD_MIN_HAP_FULL_READS 2 // full read supporting each haplotype
-#define LONGCALLD_MIN_HAP_READS 3 // XXX >= 3 reads supporting each haplotype, including partial/clipped reads
-#define LONGCALLD_MIN_NO_HAP_FULL_READS 10 // total full reads in noisy region
+#define LONGCALLD_MIN_HAP_FULL_READS 1 // >= 1full read supporting each haplotype
+#define LONGCALLD_MIN_HAP_READS 2 // >= 3 reads supporting each haplotype, including partial/clipped reads, call consensus from >= 3 reads
+#define LONGCALLD_MIN_NO_HAP_FULL_READS 10 // >10 total full reads in noisy region
+#define LONGCALLD_MIN_READ_TO_HAP_CONS_SIM 0.9 // for reads with >= 90% equal bases, assign haplotype
 // for sdust
 #define LONGCALLD_SDUST_T 5
 #define LONGCALLD_SDUST_W 20
@@ -98,18 +99,20 @@ typedef struct call_var_opt_t {
     // filters for noisy region, i.e., coverage/ratio
     int max_noisy_reg_reads, max_noisy_reg_len, min_noisy_reg_reads; 
     float max_noisy_frac_per_read, min_noisy_reg_ratio;
-    int min_hap_full_reads, min_no_hap_full_reads;
+    int min_hap_full_reads, min_hap_reads, min_no_hap_full_reads;
     // alignment
     int match, mismatch, gap_open1, gap_ext1, gap_open2, gap_ext2;
     int gap_aln; // default: 1: left (minimap2, abpoa), 2: right (wfa2)
+    float min_read_to_hap_cons_sim;
     int disable_read_realign; // disable re-alignment/MSA, only use variant calling from consensus sequence
                               // phasing is based on read cluster of the current haplotype, so not error-robust
     // general
     // int max_ploidy;
     int pl_threads, n_threads;
     // output
     htsFile *out_bam; // phased bam
-    FILE *out_vcf;
+    FILE *out_vcf; 
+    double p_error, log_p, log_1p, log_2; int max_gq; int max_qual;
     int8_t no_vcf_header, out_amb_base, no_bam_header; // phased vcf
 } call_var_opt_t;
 

src/call_var_main.c

@@ -23,8 +23,10 @@
 #define LONGCALLD_VAR_CATE_TYPE(var_cate) LONGCALLD_VAR_CATE_STR[(int)(log2(var_cate))]
 
 #define LONGCALLD_REF_CONS_ALN_STR(clu_aln_strs) clu_aln_strs
-#define LONGCALLD_CONS_READ_ALN_STR(clu_aln_strs, read_i) clu_aln_strs+(read_i+1)*2-1
-#define LONGCALLD_REF_READ_ALN_STR(clu_aln_strs, read_i) clu_aln_strs+(read_i+1)*2
+#define LONGCALLD_CONS_READ_ALN_STR(clu_aln_strs, read_i) clu_aln_strs+read_i+1
+#define LONGCALLD_REF_READ_ALN_STR(clu_aln_strs, read_i) clu_aln_strs+read_i+1
+// #define LONGCALLD_CONS_READ_ALN_STR(clu_aln_strs, read_i) clu_aln_strs+(read_i+1)*2-1
+// #define LONGCALLD_REF_READ_ALN_STR(clu_aln_strs, read_i) clu_aln_strs+(read_i+1)*2
 
 #ifdef __cplusplus
 extern "C" {
@@ -73,6 +75,8 @@ typedef struct aln_str_t {
     uint8_t *target_aln;
     uint8_t *query_aln;
     int aln_len;
+    // for partially aligned reads
+    int target_beg, target_end, query_beg, query_end; // 0..aln_len-1, including '-'/5
 } aln_str_t;
 
 struct bam_chunk_t;

src/call_var_main.h

@@ -59,6 +59,7 @@ int write_var_to_vcf(var_t *vars, const struct call_var_opt_t *opt, char *chrom)
     for (int i = 0; i < n_vars; i++) {
         var1_t var = vars->vars[i];
         if (var.n_alt_allele == 0) continue;
+        if (var.DP < opt->min_dp || var.AD[1] < opt->min_alt_dp) continue;
         if (opt->out_amb_base == 0) {
             uint8_t skip = 0;
             for (int j = 0; j < var.ref_len; j++) {
@@ -111,7 +112,12 @@ int write_var_to_vcf(var_t *vars, const struct call_var_opt_t *opt, char *chrom)
         fprintf(out_vcf, "\t%d\tPASS\tEND=%" PRId64 "", var.QUAL, var.pos + var.ref_len - 1);
         if (is_sv) fprintf(out_vcf, ";%s;%s\t", SVTYPE, SVLEN);
         else fprintf(out_vcf, "\t");
-        fprintf(out_vcf, "GT:PS\t%d|%d:%" PRId64 "\n", var.GT[0], var.GT[1], var.PS);
+        fprintf(out_vcf, "GT:DP:AD:GQ:PS\t%d|%d:%d:", var.GT[0], var.GT[1], var.DP);
+        for (int j = 0; j < 1+var.n_alt_allele; j++) {
+            if (j > 0) fprintf(out_vcf, ",");
+            fprintf(out_vcf, "%d", var.AD[j]);
+        }
+        fprintf(out_vcf, ":%d:%" PRId64 "\n", var.GQ, var.PS);
     }
     return ret;
 }