|
| 1 | +#!/usr/bin/env python |
| 2 | +import logging |
| 3 | +import argparse |
| 4 | +import glob |
| 5 | +import os |
| 6 | +from collections import Counter, defaultdict, OrderedDict |
| 7 | +from collections.abc import Set |
| 8 | +from typing import Dict |
| 9 | + |
| 10 | +# Configure logging |
| 11 | +logging.basicConfig(format="%(name)s - %(asctime)s %(levelname)s: %(message)s") |
| 12 | +logger = logging.getLogger(__name__) |
| 13 | +logger.setLevel(logging.INFO) |
| 14 | + |
| 15 | + |
| 16 | +def read_top_transcripts(quant_dir: str, file_pattern: str) -> Set[str]: |
| 17 | + """ |
| 18 | + Read the top 100 transcripts from the quantification file. |
| 19 | +
|
| 20 | + Parameters: |
| 21 | + quant_dir (str): Directory where quantification files are located. |
| 22 | + file_pattern (str): Pattern to match quantification files. |
| 23 | +
|
| 24 | + Returns: |
| 25 | + set: A set containing the top 100 transcripts. |
| 26 | + """ |
| 27 | + try: |
| 28 | + # Find the quantification file within the directory |
| 29 | + quant_file_path = glob.glob(os.path.join(quant_dir, "*", file_pattern))[0] |
| 30 | + with open(quant_file_path, "r") as file_handle: |
| 31 | + # Read the file and extract the top 100 transcripts |
| 32 | + return {line.split()[0] for i, line in enumerate(file_handle) if i > 0 and i <= 100} |
| 33 | + except IndexError: |
| 34 | + # Log an error and raise a FileNotFoundError if the quant file does not exist |
| 35 | + logger.error("No quantification files found.") |
| 36 | + raise FileNotFoundError("Quantification file not found.") |
| 37 | + |
| 38 | + |
| 39 | +def discover_transcript_attribute(gtf_file: str, transcripts: Set[str]) -> str: |
| 40 | + """ |
| 41 | + Discover the attribute in the GTF that corresponds to transcripts, prioritizing 'transcript_id'. |
| 42 | +
|
| 43 | + Parameters: |
| 44 | + gtf_file (str): Path to the GTF file. |
| 45 | + transcripts (Set[str]): A set of transcripts to match in the GTF file. |
| 46 | +
|
| 47 | + Returns: |
| 48 | + str: The attribute name that corresponds to transcripts in the GTF file. |
| 49 | + """ |
| 50 | + votes = Counter() |
| 51 | + with open(gtf_file) as inh: |
| 52 | + # Read GTF file, skipping header lines |
| 53 | + for line in filter(lambda x: not x.startswith("#"), inh): |
| 54 | + cols = line.split("\t") |
| 55 | + # Parse attribute column and update votes for each attribute found |
| 56 | + attributes = dict(item.strip().split(" ", 1) for item in cols[8].split(";") if item.strip()) |
| 57 | + votes.update(key for key, value in attributes.items() if value.strip('"') in transcripts) |
| 58 | + |
| 59 | + if not votes: |
| 60 | + # Log a warning if no matching attribute is found |
| 61 | + logger.warning("No attribute in GTF matching transcripts") |
| 62 | + return "" |
| 63 | + |
| 64 | + # Check if 'transcript_id' is among the attributes with the highest votes |
| 65 | + if "transcript_id" in votes and votes["transcript_id"] == max(votes.values()): |
| 66 | + logger.info("Attribute 'transcript_id' corresponds to transcripts.") |
| 67 | + return "transcript_id" |
| 68 | + |
| 69 | + # If 'transcript_id' isn't the highest, determine the most common attribute that matches the transcripts |
| 70 | + attribute, _ = votes.most_common(1)[0] |
| 71 | + logger.info(f"Attribute '{attribute}' corresponds to transcripts.") |
| 72 | + return attribute |
| 73 | + |
| 74 | + |
| 75 | +def parse_attributes(attributes_text: str) -> Dict[str, str]: |
| 76 | + """ |
| 77 | + Parse the attributes column of a GTF file. |
| 78 | +
|
| 79 | + :param attributes_text: The attributes column as a string. |
| 80 | + :return: A dictionary of the attributes. |
| 81 | + """ |
| 82 | + # Split the attributes string by semicolon and strip whitespace |
| 83 | + attributes = attributes_text.strip().split(";") |
| 84 | + attr_dict = OrderedDict() |
| 85 | + |
| 86 | + # Iterate over each attribute pair |
| 87 | + for attribute in attributes: |
| 88 | + # Split the attribute into key and value, ensuring there are two parts |
| 89 | + parts = attribute.strip().split(" ", 1) |
| 90 | + if len(parts) == 2: |
| 91 | + key, value = parts |
| 92 | + # Remove any double quotes from the value |
| 93 | + value = value.replace('"', "") |
| 94 | + attr_dict[key] = value |
| 95 | + |
| 96 | + return attr_dict |
| 97 | + |
| 98 | + |
| 99 | +def map_transcripts_to_gene( |
| 100 | + quant_type: str, gtf_file: str, quant_dir: str, gene_id: str, extra_id_field: str, output_file: str |
| 101 | +) -> bool: |
| 102 | + """ |
| 103 | + Map transcripts to gene names and write the output to a file. |
| 104 | +
|
| 105 | + Parameters: |
| 106 | + quant_type (str): The quantification method used (e.g., 'salmon'). |
| 107 | + gtf_file (str): Path to the GTF file. |
| 108 | + quant_dir (str): Directory where quantification files are located. |
| 109 | + gene_id (str): The gene ID attribute in the GTF file. |
| 110 | + extra_id_field (str): Additional ID field in the GTF file. |
| 111 | + output_file (str): The output file path. |
| 112 | +
|
| 113 | + Returns: |
| 114 | + bool: True if the operation was successful, False otherwise. |
| 115 | + """ |
| 116 | + # Read the top transcripts based on quantification type |
| 117 | + transcripts = read_top_transcripts(quant_dir, "quant.sf" if quant_type == "salmon" else "abundance.tsv") |
| 118 | + # Discover the attribute that corresponds to transcripts in the GTF |
| 119 | + transcript_attribute = discover_transcript_attribute(gtf_file, transcripts) |
| 120 | + |
| 121 | + if not transcript_attribute: |
| 122 | + # If no attribute is found, return False |
| 123 | + return False |
| 124 | + |
| 125 | + # Open GTF and output file to write the mappings |
| 126 | + # Initialize the set to track seen combinations |
| 127 | + seen = set() |
| 128 | + |
| 129 | + with open(gtf_file) as inh, open(output_file, "w") as output_handle: |
| 130 | + # Parse each line of the GTF, mapping transcripts to genes |
| 131 | + for line in filter(lambda x: not x.startswith("#"), inh): |
| 132 | + cols = line.split("\t") |
| 133 | + attr_dict = parse_attributes(cols[8]) |
| 134 | + if gene_id in attr_dict and transcript_attribute in attr_dict: |
| 135 | + # Create a unique identifier for the transcript-gene combination |
| 136 | + transcript_gene_pair = (attr_dict[transcript_attribute], attr_dict[gene_id]) |
| 137 | + |
| 138 | + # Check if the combination has already been seen |
| 139 | + if transcript_gene_pair not in seen: |
| 140 | + # If it's a new combination, write it to the output and add to the seen set |
| 141 | + extra_id = attr_dict.get(extra_id_field, attr_dict[gene_id]) |
| 142 | + output_handle.write(f"{attr_dict[transcript_attribute]}\t{attr_dict[gene_id]}\t{extra_id}\n") |
| 143 | + seen.add(transcript_gene_pair) |
| 144 | + |
| 145 | + return True |
| 146 | + |
| 147 | + |
| 148 | +# Main function to parse arguments and call the mapping function |
| 149 | +if __name__ == "__main__": |
| 150 | + parser = argparse.ArgumentParser(description="Map transcripts to gene names for tximport.") |
| 151 | + parser.add_argument("--quant_type", type=str, help="Quantification type", default="salmon") |
| 152 | + parser.add_argument("--gtf", type=str, help="GTF file", required=True) |
| 153 | + parser.add_argument("--quants", type=str, help="Output of quantification", required=True) |
| 154 | + parser.add_argument("--id", type=str, help="Gene ID in the GTF file", required=True) |
| 155 | + parser.add_argument("--extra", type=str, help="Extra ID in the GTF file") |
| 156 | + parser.add_argument("-o", "--output", dest="output", default="tx2gene.tsv", type=str, help="File with output") |
| 157 | + |
| 158 | + args = parser.parse_args() |
| 159 | + if not map_transcripts_to_gene(args.quant_type, args.gtf, args.quants, args.id, args.extra, args.output): |
| 160 | + logger.error("Failed to map transcripts to genes.") |
0 commit comments