Dsl2/filter sv blacklist

.travis.yml

@@ -14,7 +14,7 @@ before_install:
   - chmod 777 nextflow
   # to change the test-data for travis, please download using the following command, extract, make changes, tarball again with gzip, and upload to google drive. 
   # you will have to change the link below as well. Click to share the link, making it so anyone with the link can access, then extract the id in the link and put it here after "id="
-  - wget -O test-data.tar.gz --no-check-certificate 'https://docs.google.com/uc?export=download&confirm=no_antivirus&id=1XuitIb02GZ5v928Do2Vh4hnQNkbSAT0V'
+  - wget -O test-data.tar.gz --no-check-certificate 'https://docs.google.com/uc?export=download&confirm=no_antivirus&id=13zUVw4BZ0_5QAW7zmdCsZ_CZDHbNxyMV'
   - tar -xzvf test-data.tar.gz 
 
 script:

conf/references.config

@@ -60,6 +60,10 @@ params {
       hlaDat = "${params.reference_base}/hla/hla.dat"
       neoantigenCDNA = "${params.reference_base}/neoantigen/Homo_sapiens.GRCh37.75.cdna.all.fa.gz"
       neoantigenCDS = "${params.reference_base}/neoantigen/Homo_sapiens.GRCh37.75.cds.all.fa.gz"
+      svBlacklistBed   = "${params.genome_base}/sv_calling/pcawg6_blacklist.slop.bed.gz"
+      svBlacklistBedpe = "${params.genome_base}/sv_calling/pcawg6_blacklist.slop.trunc.bedpe.gz"
+      svBlacklistFoldbackBedpe = "${params.genome_base}/sv_calling/pcawg6_blacklist_foldback_artefacts.slop.trunc.bedpe.gz"
+      svBlacklistTEBedpe       = "${params.genome_base}/sv_calling/pcawg6_blacklist_TE_pseudogene.bedpe.gz"
     }
     'GRCh37' {
       acLoci      = "${params.genome_base}/Annotation/ASCAT/1000G_phase3_20130502_SNP_maf0.3.loci"
@@ -101,6 +105,10 @@ params {
       spliceSites = "${params.reference_base}/mskcc-igenomes/grch37/splice_sites/splice_sites.bed"
       brassRefDir      = "${params.reference_base}/mskcc-igenomes/grch37/brass"
       vagrentRefDir    = "${params.reference_base}/mskcc-igenomes/grch37/vagrent"
+      svBlacklistBed   = "${params.reference_base}/mskcc-igenomes/grch37/sv_calling/pcawg6_blacklist.slop.bed"
+      svBlacklistBedpe = "${params.reference_base}/mskcc-igenomes/grch37/sv_calling/pcawg6_blacklist.slop.bedpe"
+      svBlacklistFoldbackBedpe = "${params.reference_base}/mskcc-igenomes/grch37/sv_calling/pcawg6_blacklist_foldback_artefacts.bedpe"
+      svBlacklistTEBedpe       = "${params.reference_base}/mskcc-igenomes/grch37/sv_calling/pcawg6_blacklist_TE_pseudogene.bedpe"
     } 
     'GRCh38' {
       acLoci           = "${params.genome_base}/Annotation/ASCAT/1000G_phase3_GRCh38_maf0.3.loci"

conf/resources_juno.config

@@ -113,6 +113,10 @@
     cpus = { task.attempt < 4 ? 1 : 2 } 
     memory = { task.attempt < 3 ? 1.GB * task.attempt : 5.Gb * task.attempt } 
   }
+  withName:SomaticAnnotateSVBedpe {
+    cpus = { task.attempt < 4 ? 1 : 2 }
+    memory = { 8.GB * task.attempt }
+  }
 
 //------------- Germline pipeline
 
@@ -156,6 +160,10 @@
     cpus = { task.attempt < 4 ? 1 : 2 } 
     memory = { task.attempt < 3 ? 1.GB * task.attempt : 5.Gb * task.attempt } 
   }
+  withName:GermlineAnnotateSVBedpe {
+    cpus = { task.attempt < 4 ? 1 : 2 }
+    memory = { 8.GB * task.attempt }
+  }
 
 //------------- Quality Control
 

conf/resources_juno_genome.config

@@ -152,6 +152,10 @@
     cpus = { task.attempt < 4 ? 1 : 2 }
     memory = { 4.GB * task.attempt } 
   }
+  withName:SomaticAnnotateSVBedpe {
+    cpus = { task.attempt < 3 ? 2 : 4 }
+    memory = { 8.GB * task.attempt }
+  }
   withName:HRDetect {
     cpus = { 2 }
     memory = { 4.GB * task.attempt }
@@ -203,6 +207,10 @@
     cpus = { task.attempt < 4 ? 1 : 2 } 
     memory = { task.attempt < 3 ? 1.GB * task.attempt : 5.Gb * task.attempt } 
   }
+  withName:GermlineAnnotateSVBedpe {
+    cpus = { task.attempt < 3 ? 2 : 4 }
+    memory = { 8.GB * task.attempt }
+  }
 
 //------------- Quality Control
 

containers/iannotatesv/detect_cdna.py

@@ -51,6 +51,8 @@ def prep_intersection(intersect=pd.DataFrame(),bedpe_header_list=[]):
 def main():
 	args = usage()
 
+	print("Detecting possible cdna contamination in {}".format(os.path.basename(args.bedpe)))
+
 	[meta_header, bedpe_header_list, bedpe_df] = parse_svtools_bedpe_file(args.bedpe)
 
 	bedpe_bt = pybedtools.BedTool.from_dataframe(bedpe_df)

containers/iannotatesv/filter_regions_bedpe.py

@@ -0,0 +1,109 @@
+import argparse, sys, os
+import numpy as np
+import pandas as pd
+import pybedtools
+from utils import *
+
+def usage():
+	parser = argparse.ArgumentParser()
+	parser.add_argument('--blacklist-regions', dest="regions", help = 'regions bed/bedpe file', required = True)
+	parser.add_argument('--bedpe', help = 'bedpe file', required = True)
+	parser.add_argument('--tag', help = 'string to update FILTER column', required = True)
+	parser.add_argument('--output',dest="outfile",help = 'output file' , required = True)
+	parser.add_argument('--match-type',dest="type",help = 'both/either/notboth/neither' , required = True)
+	parser.add_argument('--ignore-strand',dest="ignore_strand", default=False, action="store_true", help = 'Use flag to ignore strand when running bedtools. Default False' )
+	return parser.parse_args()
+
+def determine_regions_filetype(filepath):
+	"""
+	Determine if the regions file is bed or bedpe, and raise error if neither
+	"""
+	if not any([filepath.endswith(i) for i in "bed|bedpe|bed.gz|bedpe.gz".split("|")]):
+		raise ValueError("--bedpe requires bed, bedpe, bed.gz or bedpe.gz file.")
+	elif any([filepath.endswith(i) for i in "bed|bed.gz".split("|")]):
+		return True
+	else: return False
+
+def validate_bedpe_input(filepath):
+	"""
+	Raise error if filetype is not bedpe
+	"""
+	if not any([filepath.endswith(i) for i in "bedpe|bedpe.gz".split("|")]):
+		raise ValueError("--bedpe requires bedpe or bedpe.gz file.")
+
+def validate_match_type(type):
+	"""
+	Raise error if type is not acceptable for use in pybedtools pair_to_pair or pair_to_bed
+	"""
+	if type not in overlap_type.keys():
+		raise ValueError( "--match-type must be {}.".format(" or ".join(overlap_type.keys())) )
+
+def main():
+	"""
+	1. validate inputs
+	2. read input beds/bedpes as pybedtools.BedTool objects
+	3. extract variant IDs from overlapping regions
+	4. update the FILTER tag for identified variants
+	"""
+	args = usage()
+
+	print("Filtering {} with {} regions".format(os.path.basename(args.bedpe), args.tag))
+
+	# parse inputs
+	try:
+		validate_bedpe_input(args.bedpe)
+		[meta_header, bedpe_header_list, bedpe_df] = parse_svtools_bedpe_file(args.bedpe)
+	except Exception as e:
+		print(e)
+		sys.exit("Unable to parse --bedpe")
+
+	try:
+		is_bed = determine_regions_filetype(args.regions)
+		regions_bt = pybedtools.BedTool(args.regions)
+	except Exception as e:
+		print(e)
+		sys.exit("Unable to parse --blacklist-regions")
+
+	# annotate bedpe
+	try:
+		bedpe_bt = pybedtools.BedTool.from_dataframe(bedpe_df)
+		ids_df = find_overlapped_ids(bedpe_bt, regions_bt, args.type, args.ignore_strand, is_bed)
+		bedpe_df = add_filter_by_id(bedpe_df,ids_df,args.tag)
+	except Exception as e:
+		print(e)
+		sys.exit("Filtering failed with inputs {} and {}".format(args.bedpe, args.regions))
+
+	# write result
+	with open(args.outfile, "w") as fw:
+		fw.write("".join(meta_header))
+	bedpe_df.to_csv(args.outfile, header=True, index=False, sep="\t", mode="a")
+
+def find_overlapped_ids(bedpe_bt,regions_bt,match_type,ignore_strand,is_bed):
+	# determine validity of match_type
+	validate_match_type(match_type)
+
+	# run pair_to_bed if regions file is bed, pair_to_pair if regions file is bedpe
+	if is_bed:
+		intersect = run_pair_to_bed(bedpe_bt,regions_bt,match_type)
+	else:
+		intersect = run_pair_to_pair(bedpe_bt,regions_bt,match_type, ignore_strand=ignore_strand)
+
+	# extract ids
+	try:
+		ids_df = pd.DataFrame({"ID":list( intersect.to_dataframe(header=None)[6].drop_duplicates() )})
+	except:
+		ids_df = pd.DataFrame(columns=["ID"])
+	return ids_df
+
+def add_filter_by_id(bedpe_df,ids_df,tag):
+	# annotate bedpe based on matching IDs
+	ids_df.columns = ["ID"]
+	ids_df["FILTER_NEW"] = tag
+	filtered_bedpe_df = bedpe_df.merge(ids_df, on="ID", how="left")
+	filtered_bedpe_df = filtered_bedpe_df.apply(lambda x: update_filter(x,x["FILTER_NEW"]) if not pd.isna(x["FILTER_NEW"]) else x, axis=1)
+	filtered_bedpe_df = filtered_bedpe_df.drop(['FILTER_NEW'], axis=1)
+
+	return filtered_bedpe_df
+
+if __name__ == "__main__":
+	main()

containers/iannotatesv/pair_to_pair_annot.py

@@ -0,0 +1,55 @@
+import argparse, sys, os
+import numpy as np
+import pandas as pd
+import pybedtools
+from utils import *
+
+def usage():
+	parser = argparse.ArgumentParser()
+	parser.add_argument('--blacklist-regions', dest="regions", help = 'regions bed/bedpe file', required = True)
+	parser.add_argument('--bedpe', help = 'bedpe file', required = True)
+	parser.add_argument('--tag', help = 'string to update FILTER column', required = True)
+	parser.add_argument('--output',dest="outfile",help = 'output file' , required = True)
+	parser.add_argument('--match-type',dest="type",help = 'both/either' , required = True)
+	parser.add_argument('--ignore-strand',dest="ignore_strand", default=False, action="store_true", help = 'Use flag to ignore strand when running bedtools. Default False' )
+	return parser.parse_args()
+
+
+def main():
+	args = usage()
+
+	if args.type not in ["both","either"]:
+		sys.exit("--match-type must be both or either")
+
+	match_type = args.type	
+	# determine if regions file is bed or bedpe
+	filename=os.path.basename(args.regions)
+	is_bed = False
+	if filename.endswith("bed") or filename.endswith("bed.gz"):
+		is_bed = True
+	if filename.endswith("bedpe") or filename.endswith("bedpe.gz"):
+		pass
+	else:
+		sys.exit("--blacklist-regions must have .bed/.bed.gz/.bedpe/.bedpe.gz extension")
+
+	[meta_header, bedpe_header_list, bedpe_df] = parse_svtools_bedpe_file(args.bedpe)
+	bedpe_bt = pybedtools.BedTool.from_dataframe(bedpe_df)
+	regions_bt = pybedtools.BedTool(args.regions)
+
+	if is_bed:
+		intersect = run_pair_to_bed(bedpe_bt,regions_bt,match_type)
+	else:
+		intersect = run_pair_to_pair(bedpe_bt,regions_bt,match_type, ignore_strand=args.ignore_strand)
+
+	intersect_df = bedtool_to_df(intersect,bedpe_header_list)[["ID"]]
+	intersect_df["FILTER_NEW"] = args.tag
+	filtered_bedpe_df = bedpe_df.merge(intersect_df, on="ID", how="left")
+	filtered_bedpe_df = filtered_bedpe_df.apply(lambda x: update_filter(x,x["FILTER_NEW"]) if not pd.isna(x["FILTER_NEW"]) else x, axis=1)
+	filtered_bed_df = filtered_bedpe_df.drop(['FILTER_NEW'], axis=1)
+
+	with open(args.outfile, "w") as fw:
+		fw.write("".join(meta_header))
+	filtered_bed_df.to_csv(args.outfile, header=True, index=False, sep="\t", mode="a")
+
+if __name__ == "__main__":
+	main()