kkonganti@105: #!/usr/bin/env python
kkonganti@105: 
kkonganti@105: import os
kkonganti@105: import argparse
kkonganti@105: import logging as log
kkonganti@105: import pandas as pd
kkonganti@105: import numpy as np
kkonganti@105: from Bio import SeqIO
kkonganti@105: 
kkonganti@105: 
kkonganti@105: def main():
kkonganti@105:     # READ IN ARGUMENTS
kkonganti@105:     desc = """
kkonganti@105:     This script is part of the centriflaken pipeline: It processes centrifuge
kkonganti@105:     output and produces either a filtered FASTQ or a text file of FASTQ IDs based
kkonganti@105:     on the supplied taxa/bug
kkonganti@105:     """
kkonganti@105:     parser = argparse.ArgumentParser(prog='process_centrifuge_output.py', description=desc)
kkonganti@105:     parser.add_argument("-v", dest='verbose', action="store_true", help="For more verbose output")
kkonganti@105:     parser.add_argument("-i", dest='input_fastq', required=False,
kkonganti@105:         help="Path to input FASTQ file (same as input to centrifuge). If not mentioned, \
kkonganti@105:             a text file of sequence IDs are produced instead of a FASTQ file")
kkonganti@105:     parser.add_argument("-t", dest='taxa_filtered_fastq_file', required=True,
kkonganti@105:         help="Path to output FASTQ or output text file filtered by the taxa specified")
kkonganti@105:     parser.add_argument("-r", dest='cent_report', required=True, help="Path to centrifuge report")
kkonganti@105:     parser.add_argument("-o", dest='cent_output', required=True, help="Path to centrifuge output")
kkonganti@105:     parser.add_argument("-b", dest='bug', required=True,
kkonganti@105:         help="Name or fragment of name of the bug by which reads are extracted")
kkonganti@105:     args = parser.parse_args()
kkonganti@105: 
kkonganti@105:     # MORE INFO IF VERBOSE
kkonganti@105:     if args.verbose:
kkonganti@105:         log.basicConfig(format="%(levelname)s: %(message)s", level=log.DEBUG)
kkonganti@105:     else:
kkonganti@105:         log.basicConfig(format="%(levelname)s: %(message)s")
kkonganti@105: 
kkonganti@105:     # ASSIGN VARIABLES
kkonganti@105:     input_fastq = args.input_fastq
kkonganti@105:     taxa_filtered_fastq_file = args.taxa_filtered_fastq_file
kkonganti@105:     cent_report = args.cent_report
kkonganti@105:     cent_output = args.cent_output
kkonganti@105:     bug = args.bug
kkonganti@105:     report_col_list = ["name", "taxID"]
kkonganti@105:     output_col_list = ["taxID", "readID"]
kkonganti@105: 
kkonganti@105:     # Match and filter taxa names and ids from centrifuge report file
kkonganti@105:     report_df = pd.read_csv(cent_report, delimiter="\t", usecols=report_col_list)
kkonganti@105:     report_df['name'] = report_df['name'].str.lower()
kkonganti@105:     filt_report_df = report_df[report_df['name'].str.contains(bug.lower())]
kkonganti@105:     #print("\nMatching taxa names and ids:\n",filt_report_df)
kkonganti@105:     taxID_list = filt_report_df['taxID']
kkonganti@105: 
kkonganti@105:     # Match the above tax ids to read ids from centrifuge output file and deduplicate
kkonganti@105:     output_df = pd.read_csv(cent_output, delimiter="\t", usecols=output_col_list)
kkonganti@105:     filt_output_df = output_df.loc[output_df['taxID'].isin(taxID_list)]
kkonganti@105:     readID_list = filt_output_df['readID']
kkonganti@105:     readID_dedup_list = np.unique(readID_list)
kkonganti@105:     TF=open(taxa_filtered_fastq_file, "w")
kkonganti@105: 
kkonganti@105:     if (not input_fastq):
kkonganti@105:         # print("\nFILTERED READ ID LIST:\n", readID_dedup_list)
kkonganti@105:         for ID in readID_dedup_list:
kkonganti@105:             TF.write(f"{ID}\n")
kkonganti@105:     else:
kkonganti@105:         # Extract filtered reads from input fastq and write to output fastq
kkonganti@105:         print ("Indexing reads..")
kkonganti@105:         rec = SeqIO.index(input_fastq,"fastq")
kkonganti@105:         for i in readID_dedup_list:
kkonganti@105:             if i in rec:
kkonganti@105:                 SeqIO.write(rec[i], TF, "fastq")
kkonganti@105: 
kkonganti@105:     TF.close()
kkonganti@105: 
kkonganti@105: if __name__ == "__main__":
kkonganti@105:     main()