geneAnchoring.smk

rule mapHomologousRegions:
	message: " mapping homologous regions of both references using ISBPs markers as anchors"
	input: 
		#closestMarkers=config['results']+'/2.mergedClosestMarkers.txt',
		marker5prime=config['results']+'/2.closestbed_split/{chrom}.upstream.txt',
		marker3prime=config['results']+'/2.closestbed_split/{chrom}.downstream.txt',
		refQuery=config['queryFasta'],
		faiQuery=config['queryFasta']+'.fai',
		refTarget=config['targetFasta'],
		faiTarget=config['targetFasta']+'.fai',
		markersOnTarget=config['results']+"/1.filteredISBPs.sorted.{chrom}.bed"
	output: dir=directory(config['results']+'/3.mapping/{chrom}')
	log: config['results']+'/3.mapping/{chrom}/microMappingPipeline.log'
	shell: 
		"""
		bin/microMappingPipeline.py {input.marker5prime} {input.marker3prime} {output.dir} {input.refQuery} {input.refTarget} {input.markersOnTarget} &> {log}
		"""

rule getFailedMarkerPairs:
	message: "extract gene info not anchored with a valid marker pair"
	input: config['results']+'/3.mapping/{chrom}/microMappingPipeline.log'
	output: config['results']+'/3.mapping/{chrom}/noMarkerPairs.txt'
	log: config['results']+'/3.mapping/{chrom}/noMarkerPairs.log'
	shell:
		"""
		fgrep -h 'Could not identify a pair of marker to use to anchor' {input} |cut -d ' ' -f 15|cut -d '.' -f 1|sort 1> {output} 2> {log}'
		"""

rule gmapRescue:
	message: " Rescue anchoring of failed genes with gmap"
	input: anchoringDir=config['results']+'/3.mapping/{chrom}'
	output: config['results']+'/4.gmapRescue/{chrom}.txt'
	params: blastdb=config['blastdb']
	log: config['results']+'/4.gmapRescue/{chrom}.log'
	shell:
		"""
		bin/gmapRescue.sh {input.anchoringDir} {params.blastdb} 1> {output} 2> {log}
		"""