library(tidyverse)
library(optparse)

#Jia Jinbu. 2020.09.15

parser <- OptionParser()
parser <- add_option(parser, c("-i", "--inrelpos"), type="character", 
						help="Intron splicing and Pol II position file generated by prepare_data_for_splice_kinetics.py")
parser <- add_option(parser, c("-a", "--inreadinfo"), type="character",
						help="Input Read info file generated by merge_read_info.R. This file was used to extract the read_core_id which type is elongating")
parser <- add_option(parser, c("-r", "--outrna"), type="character", 
						help="Output File mRNA incompletely spliced ratio")
parser <- add_option(parser, c("-p", "--outintron"), type="character", 
						help="Output File intron unspliced ratio")

opt <- parse_args(parser)
file_rel_pos = opt$inrelpos
file_read_info = opt$inreadinfo
fileout_rna_stat = opt$outrna
fileout_intron_stat = opt$outintron

rel_pos <- read_tsv(file_rel_pos)
read_info <- read_tsv(file_read_info, na="NA")

select_read_info <- filter(read_info, type == "polya")
select_rel_pos <- filter(rel_pos, read_core_id %in% select_read_info$read_core_id)

rna_stat <- select_read_info %>%
				group_by(mRNA) %>%
				summarise(total_read_num=n(), 
				          ir_read_num = sum(retention_introns != ""),
				          ir_ratio = ir_read_num/total_read_num)

ir_stat <- select_rel_pos %>% 
				group_by(intron_id) %>% 
				summarise(total_read_num=n(), 
						  ir_read_num=sum(retention), 
						  ir_ratio=ir_read_num/total_read_num)

write_tsv(rna_stat, fileout_rna_stat)
write_tsv(ir_stat, fileout_intron_stat)



