calculateMRD.py: added alternative MRD calculation

Test/*: added test for calculateMRD.py

.gitignore

@@ -16,3 +16,5 @@ umiVar_personalized
 /Test/chr10.fa
 /Test/chr10.fa.fai
 /Test/postFilter01
+/Test/test_output/
+/Test/filter_variants.log

README.md

@@ -132,3 +132,31 @@ $umiVar_R_binary="[PATH_TO_R]/RScript"
 $umiVar_samtools_binary="[PATH_TO_SAMTOOLS]/samtools"
 ```
 
+### Alternative MRD calculation
+
+The Script `calculateMRD.py` provides an alternative approach to calculate the minimum residual disease based on the background noise. For that it uses pileup files created by `umiVar` and calculates the MRD with a Fisher's Exact Test:
+```
+usage: calculateMRD.py [-h] [--max_af MAX_AF] [--keep_gonosomes] [--blacklist BLACKLIST] 
+                        [--remove_off_target] [--keep_indels] folder panel output variants
+
+Script to calculate MRD based on background noise.
+
+positional arguments:
+  folder                umiVar folder containing dedup TSV files.
+  panel                 VCF containing the cfDNA panel.(commas seperated if multiple)
+  output                Output TSV containing MRD value.
+  variants              Output TSV containing read counts and p-values for each monitoring variants.
+
+optional arguments:
+  -h, --help            show this help message and exit
+  --max_af MAX_AF       Maximum allele frequency for variants to include in the computation. (>1 to disable)
+  --keep_gonosomes      Do not remove gonosomes
+  --blacklist BLACKLIST
+                        VCF of variants which should be excluded for the background error rate
+  --remove_off_target   Remove tumor off-target variants
+  --keep_indels         Do not remove InDels.
+
+```
+
+
+

Test/Makefile

@@ -1,7 +1,7 @@
 all: clean CFDNA_01 CFDNA_02 select_monitoring post_filtering
 
 ##Test case for CI
-ci: clean download_test_reference CFDNA_02 select_monitoring post_filtering
+ci: clean download_test_reference CFDNA_02 select_monitoring post_filtering mrd_calculation_01
 
 clean:
 	rm -rf umiVar_CFDNA_01 umiVar_CFDNA_02 select_monitoring_01 select_monitoring_02 select_monitoring_03 postFilter01
@@ -93,4 +93,21 @@ check_post_filter_01:
 	diff -u CFDNA05_01.GSvar.out postFilter01/CFDNA05_01.GSvar
 	diff -u CFDNA06_01.GSvar.out postFilter01/CFDNA06_01.GSvar
 
+mrd_calculation_01: run_mrd_calculation_01 check_mrd_calculation_01
+
+run_mrd_calculation_01:
+	mkdir -p test_output
+	mkdir -p test_output/mrdCalculation01
+	python3 ../calculateMRD.py data_input/mrdCalculation01 data_input/mrdCalculation01/cfdna_panel.vcf test_output/mrdCalculation01/mrdCalculation01_bg.mrd \
+		test_output/mrdCalculation01/mrdCalculation01_monitoring_counts.tsv
+	python3 ../calculateMRD.py data_input/mrdCalculation01 data_input/mrdCalculation01/cfdna_panel.vcf test_output/mrdCalculation01/mrdCalculation01_bg_unfiltered.mrd \
+		test_output/mrdCalculation01/mrdCalculation01_monitoring_counts_unfiltered.tsv --max_af 1.0 --keep_gonosomes --keep_indels
+
+check_mrd_calculation_01:
+	diff -u data_output/mrdCalculation01/mrdCalculation01_bg.mrd test_output/mrdCalculation01/mrdCalculation01_bg.mrd
+	diff -u data_output/mrdCalculation01/mrdCalculation01_monitoring_counts.tsv test_output/mrdCalculation01/mrdCalculation01_monitoring_counts.tsv
+	diff -u data_output/mrdCalculation01/mrdCalculation01_bg_unfiltered.mrd test_output/mrdCalculation01/mrdCalculation01_bg_unfiltered.mrd
+	diff -u data_output/mrdCalculation01/mrdCalculation01_monitoring_counts_unfiltered.tsv test_output/mrdCalculation01/mrdCalculation01_monitoring_counts_unfiltered.tsv
+
+
 

Test/data_input/mrdCalculation01/cfdna_panel.vcf

@@ -0,0 +1,36 @@
+##fileformat=VCFv4.2
+##GENOME_BUILD=GRCh38
+##CREATION_DATE=2023-08-02
+##INFO=<ID=gene_change,Number=.,Type=String,Description="Affected gene and location.">
+##FORMAT=<ID=GT,Number=1,Type=String,Description="Genotype">
+#CHROM	POS	ID	REF	ALT	QUAL	FILTER	INFO	FORMAT	SAMPLE
+chr1	26300987	ID	N	N	.	.	.	GT	./.
+chr1	114713909	M	G	T	.	.	gene_change=NRAS_Q61K	GT	./.
+chr1	114713915	M	C	T	.	.	gene_change=NRAS_A59T	GT	./.
+chr2	24918669	ID	N	N	.	.	.	GT	./.
+chr2	68655336	ID	N	N	.	.	.	GT	./.
+chr3	179218303	M	G	A	.	.	gene_change=PIK3CA_E545K	GT	./.
+chr7	55174771	M	AGGAATTAAGAGAAGC	A	.	.	gene_change=EGFR_delE746-A750	GT	./.
+chr7	55181316	M	GG	GGCCAGCGTGG	.	.	gene_change=EGFR_V769-D770insASV	GT	./.
+chr7	55181378	M	C	T	.	.	gene_change=EGFR_T790M	GT	./.
+chr7	55191822	M	T	G	.	.	gene_change=EGFR_L858R	GT	./.
+chr7	100577850	ID	N	N	.	.	.	GT	./.
+chr7	150742412	ID	N	N	.	.	.	GT	./.
+chr8	106770107	ID	N	N	.	.	.	GT	./.
+chr8	143984818	ID	N	N	.	.	.	GT	./.
+chr9	128261226	ID	N	N	.	.	.	GT	./.
+chr9	128424504	ID	N	N	.	.	.	GT	./.
+chr10	17644392	ID	N	N	.	.	.	GT	./.
+chr10	71362188	ID	N	N	.	.	.	GT	./.
+chr11	47324594	ID	N	N	.	.	.	GT	./.
+chr12	25245350	M	C	T	.	.	gene_change=KRAS_G12D	GT	./.
+chr14	74286575	ID	N	N	.	.	.	GT	./.
+chr16	66913480	ID	N	N	.	.	.	GT	./.
+chr17	7513364	ID	N	N	.	.	.	GT	./.
+chr17	7847039	ID	N	N	.	.	.	GT	./.
+chr17	19093611	ID	N	N	.	.	.	GT	./.
+chr17	50356597	ID	N	N	.	.	.	GT	./.
+chr19	46674656	ID	N	N	.	.	.	GT	./.
+chr20	7982394	ID	N	N	.	.	.	GT	./.
+chr20	63540928	ID	N	N	.	.	.	GT	./.
+chr22	19132425	ID	N	N	.	.	.	GT	./.