input: Write a function to calculate the percent of transcripts mapping to rRNA given an HTseq_counts file and a GTF file. The function should read the GTF file to extract identifiers of genes/transcripts mapping to rRNA and MT genes, then count the reads mapping to these identifiers in the HTseq_counts file, and finally compute the percentage of rRNA and MT mapped reads. The function should return the summary output filename and the computed percentages.
  
  def calculate_percentage_rRNA_MT(htseq_counts_file, gtf_file):
      # Extract rRNA and MT identifiers from GTF file
      rRNA_identifiers = {}
      MT_identifiers = {}
  
      with open(gtf_file, "r") as handle:
          for line in handle:
              if 'gene_biotype "rRNA"' in line:
                  identifier = re.sub('.*gene_id "', '', line).split('";')[0]
                  rRNA_identifiers[identifier] = 1
              if 'MT' in line:
                  identifier = re.sub('.*gene_id "', '', line).split('";')[0]
                  MT_identifiers[identifier] = 1
  
      # Count reads mapping to rRNA and MT identifiers in HTseq_counts file
      Cummulative_rRNA_Count = 0
      Cummulative_MT_Count = 0
      ReadCount = 0
  
      with open(htseq_counts_file, "r") as handle:
          for line in handle:
              split_line = line.strip().split("\t")
              if split_line[0] in rRNA_identifiers:
                  Cummulative_rRNA_Count += int(split_line[1])
              if split_line[0] in MT_identifiers:
                  Cummulative_MT_Count += int(split_line[1])
              ReadCount += int(split_line[1])
  
      # Compute percentages
      perc_rRNA = round((Cummulative_rRNA_Count / ReadCount) * 100, 3)
      perc_MT = round((Cummulative_MT_Count / ReadCount) * 100, 3)
  
      # Output file name
      output_filename = re.sub('.txt_', '_', htseq_counts_file + '_rRNAmtRNACounts.txt')
  
      return output_filename, perc_rRNA, perc_MT
  
  assert calculate_percentage_rRNA_MT("C17_3_S20_Aligned.out.srt.bam_htseq_combined_counts.txt", "/Users/malwina/Documents/CTR-Data/genomes/Mus_musculus/mm10/Mus_musculus.GRCm38.84.gtf") == ("C17_3_S20_Aligned.out.srt.bam_htseq_combined_counts_rRNAmtRNACounts.txt", 15.7, 2.3)
  assert calculate_percentage_rRNA_MT("C18_5_S21_Aligned.out.srt.bam_htseq_combined_counts.txt", "/Users/malwina/Documents/CTR-Data/genomes/Mus_musculus/mm10/Mus_musculus.GRCm38.84.gtf") == ("C18_5_S21_Aligned.out.srt.bam_htseq_combined_counts_rRNAmtRNACounts.txt", 14.2, 1.8)
  assert calculate_percentage_rRNA_MT("C19_7_S22_Aligned.out.srt.bam_htseq_combined_counts.txt", "/Users/malwina/Documents/CTR-Data/genomes/Mus_musculus/mm10/Mus_musculus.GRCm38.84.gtf") == ("C19_7_S22_Aligned.out.srt.bam_htseq_combined_counts_rRNAmtRNACounts.txt", 16.9, 2.5)
  assert calculate_percentage_rRNA_MT("C20_9_S23_Aligned.out.srt.bam_htseq_combined_counts.txt", "/Users/malwina/Documents/CTR-Data/genomes/Mus_musculus/mm10/Mus_musculus.GRCm38.84.gtf") == ("C20_9_S23_Aligned.out.srt.bam_htseq_combined_counts_rRNAmtRNACounts.txt", 13.5, 1.7)
output: def calculate_percentage_rRNA_MT(htseq_counts_file, gtf_file):
      import re
      
      rRNA_identifiers = {}
      MT_identifiers = {}
  
      with open(gtf_file, "r") as handle:
          for line in handle:
              if 'gene_biotype "rRNA"' in line:
                  identifier = re.sub('.*gene_id "', '', line).split('";')[0]
                  rRNA_identifiers[identifier] = 1
              if 'MT' in line:
                  identifier = re.sub('.*gene_id "', '', line).split('";')[0]
                  MT_identifiers[identifier] = 1
  
      Cummulative_rRNA_Count = 0
      Cummulative_MT_Count = 0
      ReadCount = 0
  
      with open(htseq_counts_file, "r") as handle:
          for line in handle:
              split_line = line.strip().split("\t")
              if split_line[0] in rRNA_identifiers:
                  Cummulative_rRNA_Count += int(split_line[1])
              if split_line[0] in MT_identifiers:
                  Cummulative_MT_Count += int(split_line[1])
              ReadCount += int(split_line[1])
  
      perc_rRNA = round((Cummulative_rRNA_Count / ReadCount) * 100, 3)
      perc_MT = round((Cummulative_MT_Count / ReadCount) * 100, 3)
  
      output_filename = re.sub('.txt_', '_', htseq_counts_file + '_rRNAmtRNACounts.txt')
  
      return output_filename, perc_rRNA, perc_MT