count.smk

rule index :
	input:
		unpack(contigs_input)
	output:
		expand("work/index/{{sample}}.{ext}.bt2l", ext=["1", "2", "3", "4", "rev.1", "rev.2"])
	params:
		index = "work/index/{sample}"
	shell:
		"bowtie2-build "
		"--large-index "
		"{input} "
		"{params.index}"

rule mapping:
	input:
		R1 = "DATA/trim/{reads}_R1.fastq.gz",
		R2 = "DATA/trim/{reads}_R2.fastq.gz",
		i = expand("work/index/{{contigs}}.{ext}.bt2l", ext=["1", "2", "3", "4", "rev.1", "rev.2"])
	output:
		sam = temp("work/bowtie/align_{reads}_{contigs}.sam"),
		bam = "work/bowtie/align_{reads}_{contigs}.bam",
		bai = "work/bowtie/align_{reads}_{contigs}.bam.bai"
	params:
		index = "work/index/{contigs}"
	threads:
		config["THREADS"]
	shell:
		"bowtie2 "
		"--threads {threads} "
		"-x {params.index} "
		"-1 {input.R1} "
		"-2 {input.R2} "
		"-S {output.sam} "
		" ; "
		"samtools view "
		"--threads {threads} "
		"-b "
		"-o {output.bam} "
		"{output.sam} "
		" ; "
		"samtools sort "
		"--threads {threads} "
		"{output.bam} "
		"-o {output.bam} "
		" ; "
		"samtools index "
		"{output.bam}"

rule count_contig:
	input:
		bam = "work/bowtie/align_{reads}_{contigs}.bam",
		bai = "work/bowtie/align_{reads}_{contigs}.bam.bai"
	output:
		"report/count-contigs-{reads}-{contigs}.tsv"
	shell:
		"samtools idxstats "
		"{input.bam} "
		" > "
		"{output}"

rule count_gene:
	input:
		bam = "work/bowtie/align_{reads}_{contigs}.bam",
		bai = "work/bowtie/align_{reads}_{contigs}.bam.bai",
		gff = "work/{gff_tool}/{contigs}_{gff_tool}.gff"
	output:
		"report/count-genes-{gff_tool}-{reads}-{contigs}.tsv"
	shell:
		"htseq-count "
		"--format bam "