Benchmarking and v0.2 regions progress

Made a new version of the tr regions that's a single, consolidated file
Working on the benchmarking (this will need to be redone)
updating a lot of .gitignores because I have data littered throughout as I work on this

Author: ACEnglish

realignment/.gitignore

@@ -0,0 +1,2 @@
+simulation/
+v0.1variants_v0.2regions_realignment/

regions/.gitignore

@@ -1,3 +1,4 @@
+adotto_TRannotations_v0.2.bed.gz
 *.tbi
 data/merged.slop25.bed.gz
 data/tr_regions.bed.gz

regions/README.md

@@ -119,7 +119,31 @@ bcftools view -r chr:start-end some.vcf.gz # Query variants
 Note `samtools faidx` is not the same coordinate system as `tabix`. But the bed files are 0-based and vcfs are 1-based,
 both half-open(?).
 
-Next steps
-==========
+Variants
+========
 After creating variants, go to `intersection/` to use the variants to perform more filtering and analysis on the
 tr_regions/annotations.
+
+Creating the final annotation
+=============================
+As of right now, we have tr regions (simple chr:start-end) and annotations (coordinates + other info). Since there can
+be multiple annotations per-region, it makes sense to combine these two files into one to facilitate analysis. We'll
+design this file for use with `truvari anno trf`, but it can be adopted to other uses. The final file format will be a
+tab-delimited file with columns:
+
+* chrom - tr_region chromosome
+* start - tr_region 0-based start position
+* end - tr_region 0-based end position
+* annos - json of annotations
+
+The json of annotations will be a list of dicts with key:values listed above in "Defining Repeats".
+
+To create this file, simply run:
+
+```bash
+python scripts/tr_reganno_maker.py intersection/data/tr_regions.bed intersection/data/tr_annotations.bed.gz \
+    | bedtools sort -i - | bgzip > adotto_TRannotations_v0.2.bed.gz
+```
+
+The tr_regions can be `.bed` or `.bed.gz` but the tr_annotations need to be `.bed.gz` with a `.tbi` index.
+

regions/intersection/.gitignore

@@ -0,0 +1,3 @@
+counts_tr_annotations_vinter.bed.gz
+data/
+filtered_tr_annotations_vinter.bed

regions/scripts/tr_reganno_maker.py

@@ -0,0 +1,29 @@
+import sys
+import json
+import tabix
+import truvari
+
+def main(reg_fn, anno_fn):
+    """
+    Combine regions and annotations into a file
+    """
+    header = [('chrom', str), ('start', int), ('end', int), ('period', float),
+              ('copies', float), ('score', int), ('entropy', float), ('repeat', str)]
+
+    tb = tabix.open(anno_fn)
+    for line in truvari.opt_gz_open(reg_fn):
+        chrom, start, end = line.strip().split('\t')[:3]
+        start = int(start)
+        end = int(end)
+        m_data = []
+        for i in tb.query(chrom, start, end):
+            try:
+                m_data.append({fmt[0]: fmt[1](x) for x, fmt in zip(i, header)})
+            except tabix.TabixError as e:
+                pass # allow empty regions
+        data_str = json.dumps(m_data)
+        print(f"{chrom}\t{start}\t{end}\t{data_str}")
+
+if __name__ == '__main__':
+    reg, anno = sys.argv[1:]
+    main(reg, anno)

variants/benchmarking/.gitignore

@@ -1,4 +1,5 @@
 inputs/merged.vcf.gz
+inputs/merged.vcf.gz.tbi
 inputs/GRCh38.variants.sqoff.vcf.gz
 inputs/HG002_comp.vcf.gz
 inputs/HG002_comp.vcf.gz.tbi

variants/benchmarking/results/.gitignore

@@ -0,0 +1,12 @@
+HG002_comp.jl
+NA24385_comp.jl
+base_cnt.jl
+base_sv.jl
+by_size_performance_maker.py
+deliver.tgz
+li:NA24385_comp.jl
+multi_vcfeval/
+old_regions/
+test.jl
+test2.jl
+x.sh

variants/benchmarking/scripts/run_bench.sh

@@ -10,6 +10,8 @@ comp_vcfs[eichler]=/users/u233287/scratch/code/adotto/variants/benchmarking/inpu
 bdir=/users/u233287/scratch/code/adotto/variants/benchmarking/
 thfa_vcf=$bdir/truthsets/HPRC-cur.20211005-align2-GRCh38.dip.singlealleles.vcf.gz
 thfa_bed=$bdir/truthsets/HPRC-cur.20211005-align2-GRCh38.dip.bed
+thfa_sm_bed=$bdir/truthsets/temp_hc_regions/GRCh38_HG2-HPRC-20211005_dipcall-z2k.smallvar.benchmark.bed
+thfa_lg_bed=$bdir/truthsets/temp_hc_regions/GRCh38_HG2-HPRC-20211005_dipcall-z2k.SV.benchmark.bed
 cmrg_vcf=$bdir/truthsets/HG002_GRCh38_CMRG_SV_v1.00.vcf.gz
 cmrg_bed=$bdir/truthsets/HG002_GRCh38_CMRG_SV_v1.00.bed
 
@@ -21,13 +23,13 @@ mkdir -p results/
 for samp_name in "${!comp_vcfs[@]}"
 do
     m_vcf="${comp_vcfs[$samp_name]}"
-    #$rtg vcfeval -t $rtg_ref --squash-ploidy -e $cmrg_sm_bed \
-                 #-b $cmrg_sm_vcf -c $m_vcf -o results/rtg_cmrg_${samp_name}
-    #$rtg vcfeval -t $rtg_ref --squash-ploidy -e $thfa_bed \
-                 #-b $thfa_ma_vcf -c $m_vcf -o results/rtg_thfa_${samp_name}
+    $rtg vcfeval -t $rtg_ref --squash-ploidy -e $cmrg_sm_bed \
+                 -b $cmrg_sm_vcf -c $m_vcf -o results/rtg_cmrg_${samp_name}
+    $rtg vcfeval -t $rtg_ref --squash-ploidy -e $thfa_sm_bed \
+                 -b $thfa_vcf -c $m_vcf -o results/rtg_thfa_${samp_name}
 
     truvari bench -b $cmrg_vcf -c $m_vcf -f $ref --includebed $cmrg_bed \
                   --passonly -o results/truvari_cmrg_${samp_name}
-    truvari bench -b $thfa_vcf -c $m_vcf -f $ref --includebed $thfa_bed \
+    truvari bench -b $thfa_vcf -c $m_vcf -f $ref --includebed $thfa_lg_bed \
                   --passonly -o results/truvari_thfa_${samp_name}
 done

variants/data/.gitignore

@@ -0,0 +1,4 @@
+GRCh38.variants.sqoff.vcf.gz
+GRCh38.variants.sqoff.vcf.gz.tbi
+adotto_variants.grch38.sqoff.bcf.gz
+adotto_variants.grch38.sqoff.bcf.gz.csi

variants/diversity/.gitignore

@@ -0,0 +1 @@
+results/