# TODO: Add comment
# 
# Author: E.Korsching  2024
###############################################################################



# 1) a tool to extract complete entries in a fasta file - in two flavours - small RNA focus
# 2) a tool to extract all ">" description lines in a fasta file


# 1)
# a)  for smaller files
fasta_keep <- function(fin="", fout="", tp1=">", tp2="hsa"){
	# read fasta in and keep tp1,tp2 and following seq. line
	atext <- readLines(con=fin)
	nr <- length(atext)
	out <- NULL
	keep <- F
	for(i in 1:nr){
		tmp1 <- grepl(pattern=tp1, x=atext[i], ignore.case=F, perl=F, fixed=F)
		tmp2 <- grepl(pattern=paste(tp1,tp2,sep=""), x=atext[i], ignore.case=F, perl=F, fixed=F)
		if(tmp1 & tmp2){
			keep <- T
			out <- c(out, atext[i])
		} else if(!tmp1 & !tmp2 & keep){
			out <- c(out, atext[i])
		} else if(tmp1 & !tmp2){
			keep <- F
		}
		#cat("\n tmp1 ",tmp1," tmp2 ",tmp2)
	}
	if(is.null(out)){ stop("nothing found") }
	writeLines(text=out, con=fout, sep="\n")
	nr2 <- length(out)
	cat("\n rows in ",nr," rows out ",nr2,"\n")
	return()
}


#fasta_keep(fin="/home/korschi/new_ref_mir/mirdb/hairpin.fa", fout="/home/korschi/new_ref_mir/mirdb/hairpin_hsa.fa", tp1=">", tp2="hsa")
# rows in  120116 rows out  5644
#fasta_keep(fin="/home/korschi/new_ref_mir/mirdb/mature.fa", fout="/home/korschi/new_ref_mir/mirdb/mature_hsa.fa", tp1=">", tp2="hsa")
# rows in  97770  rows out  5312



# b)  for large files
fasta_keep2 <- function(fin="", fout="", tp1=">", tp2="hsa"){
	# read fasta in and keep tp1,tp2 and following seq. lines
	cat("", file=fout, append=F)
	nr <- 0; nr2 <- 0
	keep <- F
	con1 <- file(fin, "r")
	while(T){
		line1 <- readLines(con=con1, n=1)
		if(length(line1) == 0){ break }
		#print(line1)
		nr <- nr+1
		tmp1 <- grepl(pattern=tp1, x=line1, ignore.case=F, perl=T, fixed=F)
		tmp2 <- grepl(pattern=paste(tp1,tp2,sep=""), x=line1, ignore.case=F, perl=T, fixed=F)
		if(tmp1 & tmp2){
			keep <- T
			nr2 <- nr2+1
			cat(paste(line1,"\n",sep=""), file=fout, append=T)
		} else if(!tmp1 & !tmp2 & keep){
			nr2 <- nr2+1
			cat(paste(line1,"\n",sep=""), file=fout, append=T)
		} else if(tmp1 & !tmp2){
			keep <- F
		}
	}
	close(con1)
	cat("\n rows in ",nr," rows out ",nr2,"\n")
	return()
}


#fasta_keep2(fin="/home/korschi/new_ref_mir/RNAcentral/test.fasta", fout="/home/korschi/new_ref_mir/RNAcentral/test_out.fasta", tp1=">", tp2="URS.*_9606 ")

#fasta_keep2(fin="/home/korschi/on1/rnacentral_species_specific_ids.fasta", fout="/home/korschi/on1/rnacentral_spec_speci_ids_HS.fasta", tp1=">", tp2="URS.*_9606 ")
# rows in   503,672,999   rows out  14,416,176
# 29.8GB : 2h 45min : 864MB
# mem allocation a bit faulty over time (~1GB)




# 2)  for smaller files
fasta_keep3 <- function(fin="", fout="", tp1=">"){
	# read fasta in and keep tp1 lines
	atext <- readLines(con=fin)
	nr <- length(atext)
	out <- NULL
	keep <- F
	for(i in 1:nr){
		tmp1 <- grepl(pattern=tp1, x=atext[i], ignore.case=F, perl=F, fixed=F)
		if(tmp1){
			out <- c(out, atext[i])
		}
	}
	if(is.null(out)){ stop("nothing found") }
	writeLines(text=out, con=fout, sep="\n")
	nr2 <- length(out)
	cat("\n rows in ",nr," rows out ",nr2,"\n")
	return()
}


#fasta_keep3(fin="/home/korschi/1/t_hg38_tRNAs.fa", fout="/home/korschi/1/t_hg38_tRNAs_info.txt", tp1=">")
#fasta_keep3(fin="/home/korschi/1/m_tRNA.fasta", fout="/home/korschi/1/m_tRNA_info.txt", tp1=">")