Addition of scripts that couts the number of reads in the CDS in the process data pipeline.

snakemake/process_data/Snakefile

 configfile: "config.yaml"
 #from snakemake.utils import listfiles
 
-localrules: create_output_and_log_directories, remove_multimappers, read_length_histogram, finish
+localrules: create_output_and_log_directories, remove_multimappers, read_length_histogram, count_reads, sfinish
 
 #################################################################################
 ### Finish rule
@@ -9,7 +9,8 @@ localrules: create_output_and_log_directories, remove_multimappers, read_length_
 
 rule finish:
 	input:
-		sam = expand(os.path.join(config["output_dir"], "{sample}/read_length/read_length_histogram.pdf"), sample=config["sample"])
+		pdf = expand(os.path.join(config["output_dir"], "{sample}/read_length/read_length_histogram.pdf"), sample=config["sample"]),
+		counts = expand(os.path.join(config["output_dir"], "{sample}/counts.tsv"), sample=config["sample"])
 
 #################################################################################
 ### Create output and log directories
@@ -246,3 +247,22 @@ rule read_length_histogram:
 		"(python scripts/plot_read_lengths.py \
 		--sam {input.sam} \
 		--outdir {params.dir}) & > {log}"
+
+################################################################################
+### Count reads
+################################################################################
+
+rule count_reads:
+	input:
+		sam = os.path.join(config["output_dir"], "{sample}/transcripts.mapped.unique.sam"),
+		transcript_id_gene_id_CDS = config["transcript_id_gene_id_CDS"]
+	output:
+		counts = os.path.join(config["output_dir"], "{sample}/counts.tsv")
+	threads:	1
+	singularity:
+		"docker://perl:5.24-slim"
+	shell:
+		"(perl scripts/xp-count-reads-ribseq.pl \
+		{input.sam} \
+		{input.transcript_id_gene_id_CDS} \
+		> {output.counts})"

snakemake/process_data/config.yaml

@@ -6,6 +6,7 @@
   other_RNAs_index: "../prepare_annotation/results/other_RNAs_sequence.idx"
   transcripts_sequence: "../prepare_annotation/results/longest_pc_transcript_per_gene.fa"
   transcripts_index: "../prepare_annotation/results/longest_pc_transcript_per_gene.idx"
+  transcript_id_gene_id_CDS: "../prepare_annotation/results/transcript_id_gene_id_CDS.tsv"
   ##############################################################################
   ### Output and log directory
   ##############################################################################