Noncpg

src/genread.c

@@ -204,7 +204,7 @@ static inline int get_rlen(core_t *core, int tid){
     return len;
 }
 
-static void methylate_dna(core_t *core, int32_t ref_len, int32_t ref_pos, int32_t rlen, char c, char *seq, int seq_i, int tid){
+void methylate_dna(core_t *core, int32_t ref_len, int32_t ref_pos, int32_t rlen, char c, char *seq, int seq_i, int tid){
     if(core->ref->ref_meth[seq_i]){
         for(int i=0; i<rlen; i++){
             int rpos = ref_pos+i;
@@ -240,6 +240,40 @@ static void methylate_dna(core_t *core, int32_t ref_len, int32_t ref_pos, int32_
     return;
 }
 
+
+
+void methylate_all_c_dna(core_t *core, int32_t ref_len, int32_t ref_pos, int32_t rlen, char c, char *seq, int seq_i, int tid){
+    if(core->ref->ref_meth[seq_i]){
+        for(int i=0; i<rlen; i++){
+            int rpos = ref_pos+i;
+            if(rpos < ref_len && i<rlen){
+
+                char base = core->ref->ref_seq[seq_i][ref_pos+i];
+                int p = -1;
+                if(c == '+' && base == 'C'){
+                    p = i;
+                } else if (c == '-' && base == 'G') {
+                    p = rlen - i - 2;
+                }
+
+                if(p != -1){
+                    if(!(p>=0 && p<rlen)){
+                        fprintf(stderr,"p: %d, i: %d, rlen: %d, c: %c, ref_pos: %d, ref_len: %d\n",p,i,rlen,c,ref_pos,ref_len);
+                        assert(0);
+                    }
+                    assert(seq[p] == 'C' || seq[p] == 'c');
+                    int methr = rng(&core->rand_meth[tid]) * 254;
+                    if( methr <= core->ref->ref_meth[seq_i][ref_pos+i]){
+                        seq[p] = 'M';
+                    }
+                }
+            }
+        }
+    }
+    return;
+}
+
+
 static char *gen_read_dna(core_t *core, char **ref_id, int32_t *ref_len, int32_t *ref_pos, int32_t *rlen, char *c, int tid){
 
     char *seq = NULL;
@@ -274,7 +308,11 @@ static char *gen_read_dna(core_t *core, char **ref_id, int32_t *ref_len, int32_t
     }
 
     if(core->opt.meth_freq){
-        methylate_dna(core, *ref_len, *ref_pos, *rlen, *c, seq, seq_i, tid);
+        if(core->opt.flag & SQ_ALL_CTX){
+            methylate_all_c_dna(core, *ref_len, *ref_pos, *rlen, *c, seq, seq_i, tid);
+        } else {
+            methylate_dna(core, *ref_len, *ref_pos, *rlen, *c, seq, seq_i, tid);
+        }
     }
 
     return seq;

src/sim.c

@@ -512,7 +512,8 @@ void fake_uuid(char *read_id, int64_t num){
 
 char *gen_read(core_t *core, char **ref_id, int32_t *ref_len, int32_t *ref_pos, int32_t *rlen, char *c, int8_t rna, int tid);
 int16_t *gen_sig(core_t *core, const char *read, int32_t len, double *offset, double *median_before, int64_t *len_raw_signal, int8_t rna, int tid, aln_t *aln);
-
+void methylate_dna(core_t *core, int32_t ref_len, int32_t ref_pos, int32_t rlen, char c, char *seq, int seq_i, int tid);
+void methylate_all_c_dna(core_t *core, int32_t ref_len, int32_t ref_pos, int32_t rlen, char c, char *seq, int seq_i, int tid);
 
 /* process the ith read in the batch db */
 void work_per_single_read(core_t* core,db_t* db, int32_t i, int tid) {
@@ -551,6 +552,13 @@ void work_per_single_read(core_t* core,db_t* db, int32_t i, int tid) {
         strand = '+';
         ref_pos_st = 0;
         ref_pos_end = rlen;
+        if(core->opt.meth_freq){
+            if (opt.flag & SQ_ALL_CTX) {
+                methylate_all_c_dna(core, rlen, ref_pos_st, rlen, strand, seq, core->total_reads+i, tid);
+            } else{
+                methylate_dna(core, rlen, ref_pos_st, rlen, strand, seq, core->total_reads+i, tid);
+            }
+        }
     } else {
         seq=gen_read(core, &rid, &ref_len, &ref_pos_st, &rlen, &strand, rna, tid);
         ref_pos_end = ref_pos_st+rlen;
@@ -703,7 +711,8 @@ static struct option long_options[] = {
     {"cdna", no_argument, 0, 0 },                   //34 cdna
     {"ont-friendly", required_argument, 0, 0},             //35 ont-friendly
     {"meth-freq", required_argument, 0, 0 },                   //36 meth-freq
-    {"meth-model", required_argument, 0, 0 },                   //37 meth-model
+    {"meth-model", required_argument, 0, 0 },                  //37 meth-model
+    {"meth-all-ctx", required_argument, 0, 0 },                 //38 meth-all-ctx
     {0, 0, 0, 0}};
 
 
@@ -768,6 +777,7 @@ static void print_help(FILE *fp_help, opt_t opt, profile_t p, int64_t nreads) {
     fprintf(fp_help,"   --median-before-std FLOAT  Median before standard deviation [%.1f]\n",p.median_before_std);
     fprintf(fp_help,"   --kmer-model FILE          custom nucleotide k-mer model file (format similar to f5c models)\n");
     fprintf(fp_help,"   --meth-model FILE          custom methylation k-mer model file (format similar to f5c models)\n");
+    fprintf(fp_help,"   --meth-all-ctx=yes/no      relax CpG condition (for custom meth-model)\n");
     fprintf(fp_help,"\n");
     fprintf(fp_help,"See the manual page on GitHub for more details, options and the format of input/output files.\n");
 
@@ -1001,6 +1011,9 @@ int sim_main(int argc, char* argv[], double realtime0) {
         } else if (c == 0 && longindex == 37){ //meth model
             opt.meth_model_file = optarg;
             //WARNING("%s","Option --meth-model is experimental. Please report any issues.")
+        } else if (c == 0 && longindex == 38){ //all context
+            WARNING("%s","Option --meth-all-ctx is experimental. Please report any issues.")
+            yes_or_no(&opt, SQ_ALL_CTX, longindex, optarg, 1);
         } else if (c == '?'){
             exit(EXIT_FAILURE);
         } else {

src/sq.h

@@ -41,6 +41,7 @@
 #define SQ_CDNA 0x200 //CDNA
 #define SQ_ONT 0x400 //ont friendly
 #define SQ_PROM 0x800 //is a promethion
+#define SQ_ALL_CTX 0x1000 //all context instead of CpG
 
 #define WORK_STEAL 1 //simple work stealing enabled or not (no work stealing mean no load balancing)
 #define STEAL_THRESH 1 //stealing threshold