VcfCalculatePRS: improved documentation and error message.

Author: marc-sturm

README.md

@@ -16,24 +16,21 @@ Binaries of *ngs-bits* are available via Bioconda. Alternatively, *ngs-bits* can
 
 Changes already implemented in GIT master for next release:
 
-* GSvar: improved cfDNA sample handling (Batchimport of cfDNA panels, queue anlysis, ...)
+* none
+
+
+Changes in release 2021_09:  
+
+* GenePrioritization: Performs gene prioritization based on list of known disease genes and a PPI graph.
+* GraphStringDb: Creates simple representation of String-DB interaction graph.
+* VcfCheck: improved checking of empty INFO column
 * VcfAnnotateFromBed: added multithread support
-* MappingQC: added support for cfDNA samples
-
-Changes in release 2021_06:
-
-* General: Improved GRCh38 support in several tools.
-* General: Using BGZIP for compressed VCFs now to allow indexing them with tabix. 
-* VcfAnnotateFromBed: Made separator configurable; Added check for separator in source BED file; Fixed broken output VCF if input has no FORMAT column.
-* VcfAnnotateFromVcf: Fixed crash in VCF header parser.
-* NGSDExportSamples: Added ancestry column.
-* SampleAncestry: Improved runtime and memory use.
-* SampleGender: Improved runtime for algorithm 'hetx'.
-* SomaticQC: Added support for mutect2.
+* VcfAnnotateFromVcf: improved memory useage
+* MappingQC: added support for cfDNA samples, improved support for RNA
 * NGSD:
-	* Added disease status 'Unclear' to table 'sample'.
-	* Added table 'processed_sample_ancestry'.
-	* Added percent occupied to 'runqc_lane' (for Illumina NovaSeq).
+	* removed 'gene' and 'variant\_type' columns from 'variant' table
+	* added 'germline\_het' and 'germline\_hom' columns to 'variant' table
+	* added method 'shallow WGS' to 'variant\_validation' table
 
 For older releases see the [releases page](https://github.com/imgag/ngs-bits/releases).
 
@@ -134,6 +131,7 @@ The default output format of the quality control tools is [qcML](https://pubmed.
 * [VcfAnnotateFromBed](doc/tools/VcfAnnotateFromBed.md) - Annotates the INFO column of a VCF with data from a BED file.
 * [VcfAnnotateFromVcf](doc/tools/VcfAnnotateFromVcf.md) - Annotates the INFO column of a VCF with data from another VCF file (or multiple VCF files if config file is provided)
 * [VcfBreakMulti](doc/tools/VcfBreakMulti.md) - Breaks multi-allelic variants into several lines, making sure that allele-specific INFO/SAMPLE fields are still valid.
+* [VcfCalculatePRS](doc/tools/VcfCalculatePRS.md) - Calculates the Polgenic Risk Score(s) for a sample.
 * [VcfCheck](doc/tools/VcfCheck.md) - Checks a VCF file for errors.
 * [VcfExtractSamples](doc/tools/VcfExtractSamples.md) - Extract one or several samples from a VCF file.
 * [VcfFilter](doc/tools/VcfFilter.md) - Filters a VCF based on the given criteria.
@@ -154,6 +152,8 @@ The default output format of the quality control tools is [qcML](https://pubmed.
 
 ### Gene handling tools
 
+* [GenePrioritization](doc/tools/GenePrioritization.md): Performs gene prioritization based on list of known disease genes and a PPI graph (see also GraphStringDb).
+* [GraphStringDb](doc/tools/GraphStringDb.md): Creates simple representation of String-DB interaction graph.
 * [GenesToApproved](doc/tools/GenesToApproved.md) - Replaces gene symbols by approved symbols using the HGNC database (needs [NGSD](doc/install_ngsd.md)).
 * [GenesToBed](doc/tools/GenesToBed.md) - Converts a text file with gene names to a BED file (needs [NGSD](doc/install_ngsd.md)).
 * [NGSDExportGenes](doc/tools/NGSDExportGenes.md) - Lists genes from NGSD (needs [NGSD](doc/install_ngsd.md)).

doc/tools/VcfCalculatePRS.md

@@ -1,10 +1,14 @@
 ### VcfCalculatePRS tool help
-	VcfCalculatePRS (2020_06-40-g7bb4fdf2)
+	VcfCalculatePRS (2021_06-89-gbbd16264)
 
-	Calculates the Polgenic Risk Score for a given set of PRS VCFs
+	Calculates the Polgenic Risk Score(s) for a sample.
+	
+	The PRS VCF files have to contain a WEIGHT entry in the INFO column.
+	Additionally some information about the PRS score is required in the VCF header.
+	An example VCF file can be found at https://github.com/imgag/ngs-bits/blob/master/src/tools-TEST/data_in/VcfCalculatePRS_prs2.vcf
 
 	Mandatory parameters:
-	  -in <file>      Tabix indexed VCF.GZ file of the sample.
+	  -in <file>      Tabix indexed VCF.GZ file of a sample.
 	  -prs <filelist> List of PRS VCFs.
 	  -out <file>     Output TSV file containing Scores and PRS details
 
@@ -15,7 +19,7 @@
 	  --tdx           Writes a Tool Definition Xml file. The file name is the application name with the suffix '.tdx'.
 
 ### VcfCalculatePRS changelog
-	VcfCalculatePRS 2020_06-40-g7bb4fdf2
+	VcfCalculatePRS 2021_06-89-gbbd16264
 
 	2020-07-22 Initial version of this tool.
 [back to ngs-bits](https://github.com/imgag/ngs-bits)

src/VcfCalculatePRS/main.cpp

@@ -20,8 +20,10 @@ class ConcreteTool
 
 	virtual void setup()
 	{
-		setDescription("Calculates the Polgenic Risk Score for a given set of PRS VCFs");
-		addInfile("in", "Tabix indexed VCF.GZ file of the sample.", false);
+		setDescription("Calculates the Polgenic Risk Score(s) for a sample.");
+		setExtendedDescription(QStringList() << "The PRS VCF files have to contain a WEIGHT entry in the INFO column." << "Additionally some information about the PRS score is required in the VCF header." << "An example VCF file can be found at https://github.com/imgag/ngs-bits/blob/master/src/tools-TEST/data_in/VcfCalculatePRS_prs2.vcf");
+
+		addInfile("in", "Tabix indexed VCF.GZ file of a sample.", false);
 		addInfileList("prs", "List of PRS VCFs.", false);
 		addOutfile("out", "Output TSV file containing Scores and PRS details", false);
 
@@ -52,15 +54,14 @@ class ConcreteTool
 			//does not support multi sample
 			if(prs_variant_list.sampleIDs().count() > 1)
 			{
-				THROW(FileParseException, "VCF file can only contain one sample for Polgenic Risk Score calculation.");
+				THROW(FileParseException, "PRS VCF file must not contain more than one sample: " + prs_file_path);
 			}
 
 			//define PRS
 			double prs = 0;
 			QVector<double> percentiles;
 			QHash<QByteArray,QByteArray> column_entries;
 
-
 			//parse comment lines
 			foreach(const VcfHeaderLine& comment_line, prs_variant_list.vcfHeader().comments())
 			{
@@ -71,7 +72,7 @@ class ConcreteTool
 					{
 						if(column_entries.contains(column_name))
 						{
-							THROW(FileParseException, "Comment section of PRS VCF '" + prs_file_path + "' contains more than one entry for '" +column_name + "'!");
+							THROW(FileParseException, "Comment section of PRS VCF  file contains more than one entry for '" +column_name + "': " + prs_file_path);
 						}
 						column_entries[column_name] = comment_line.value.trimmed();
 						break;
@@ -80,10 +81,9 @@ class ConcreteTool
 
 				if(comment_line.key.startsWith("percentiles"))
 				{
-					if (percentiles.size() != 0) THROW(FileParseException, "Percentiles for PRS VCF '" + prs_file_path + "' are given twice!");
+					if (percentiles.size() != 0) THROW(FileParseException, "Percentiles in PRS VCF file given twice: " + prs_file_path);
 					QByteArrayList percentile_string = comment_line.value.trimmed().split(',');
-					if (percentile_string.size() != 100) THROW(FileParseException, "Invalid number of percentiles given (required: 100, given: "
-															   + QByteArray::number(percentile_string.size()) + "!");
+					if (percentile_string.size() != 100) THROW(FileParseException, "Invalid number of percentiles given (required: 100, given: "  + QByteArray::number(percentile_string.size()) + ": " + prs_file_path);
 					foreach (const QByteArray& value_string, percentile_string)
 					{
 						percentiles.append(Helper::toDouble(value_string, "Percentile"));
@@ -98,7 +98,7 @@ class ConcreteTool
 				if(col_entries_not_in_header.contains(key)) continue;
 				if(!column_entries.contains(key))
 				{
-					THROW(FileParseException, "Comment section of PRS VCF '" + prs_file_path + "' misses the entry for '" + key + "'!");
+					THROW(FileParseException, "Comment section of PRS VCFs does not contsin an entry for '" + key + "': " + prs_file_path);
 				}
 			}
 
@@ -110,7 +110,7 @@ class ConcreteTool
 				//does not support multi-allelic variants
 				if(prs_variant.isMultiAllelic())
 				{
-					THROW(FileParseException, "Does not support multi-allelic variants for Polgenic Risk Score calculation.");
+					THROW(FileParseException, "Multi-allelic variants in PRS VCF files are not supported: " + prs_variant.variantToString());
 				}
 
 				int allele_count = 0;
@@ -119,13 +119,17 @@ class ConcreteTool
 				QByteArrayList matching_variants;
 				foreach(const QByteArray& line, matching_lines)
 				{
-					// check if variant has same ref/alternative base(s)
-					if((Sequence(line.split('\t')[3]) == prs_variant.ref()) && (Sequence(line.split('\t')[4]) == prs_variant.alt(0))) matching_variants.append(line);
+					// check if overlapping variant is actually the one we are looking for
+					QByteArrayList parts = line.split('\t');
+					if(parts[1].toInt()==prs_variant.start() && parts[3]==prs_variant.ref() && parts[4]==prs_variant.alt(0))
+					{
+						matching_variants.append(line);
+					}
 				}
 
 				if(matching_variants.size() > 1)
 				{
-					THROW(FileParseException, "Variant '" + prs_variant.variantToString() + "' occures multiple times in sample VCF!");
+					THROW(FileParseException, "Variant occures multiple times in sample VCF: " +  prs_variant.variantToString());
 				}
 
 				if(matching_variants.size() == 1)
@@ -135,12 +139,12 @@ class ConcreteTool
 					QByteArrayList format_header_items = split_line[8].split(':');
 					QByteArrayList format_value_items = split_line[9].split(':');
 					int genotype_idx = format_header_items.indexOf("GT");
-					if(genotype_idx < 0) THROW(FileParseException, "Genotype information is missing for variant '" + prs_variant.variantToString() + "'!");
-					QByteArray genotype = format_value_items[genotype_idx].trimmed();
+					if(genotype_idx < 0) THROW(FileParseException, "Genotype information is missing for sample variant: " + matching_variants[0]);
+					QByteArray genotype = format_value_items[genotype_idx].trimmed().replace("|", "/").replace(".", "0");
 
 					if(genotype == "0/1") allele_count = 1;
 					else if(genotype == "1/1") allele_count = 2;
-					else THROW(FileParseException, "Invalid genotype '" + genotype + "' in variant " + prs_variant.variantToString() + "'!");
+					else THROW(FileParseException, "Invalid genotype '" + genotype + "' in sample variant: " + matching_variants[0]);
 
 					//calculate PRS part
 					double weight = Helper::toDouble(prs_variant.info("WEIGHT"), "PRS weight");