<reponame>roysomak4/cosmic-preprocess<filename>process_cosmic_data.py
from collections import Counter
import sys
import os
import gzip
import sh
import progressbar


def process_cosmic_db(args):
    if len(args) < 2:
        print ('Insufficient arguments.')
        print('Expected: python process_cosmic_data.py <cosmic_coding.vcf.gz> <cosmic_mutant_export.tsv.gz><truncate large variant [true or false]>')
        print ('Truncate large variant is optional. Default = False. If true, removes variants with ref or alt fields >1000bp in size')
        sys.exit(1)

    # COSMIC coding mut VCF file
    cosmic_coding_vcf = args[0]

    # COSMIC mut export file
    cosmic_export = args[1]

    # Truncate large variants beyond a defined size (default = 1000bp)
    truncate_large_var = False
    max_allele_len = 1000
    if len(args) > 2:
        truncate_var = args[2].lower()
        if truncate_var == 'true':
            truncate_large_var = True
        elif truncate_var == 'false':
            truncate_large_var = False
        else:
            print ('Invalid truncate variant argument.')
            print('Expected: python process_cosmic_data.py <cosmic_coding.vcf.gz> <cosmic_mutant_export.tsv.gz><truncate large variant [true or false]>')
            print ('Truncate large variant is optional. Default = False. If true, removes variants with ref or alt fields >1000bp in size')
            sys.exit(1)


    # cosmic processed file
    cosmic_output_file = "cosmic_processed.txt.gz"

    # process vcf file first
    unique_variants = {}
    process_vcf(cosmic_coding_vcf, unique_variants, truncate_large_var, max_allele_len)

    # process mut export file
    sites = {}
    process_mut_export(cosmic_export, sites)

    # merge unique variants and sites count data and write to output file
    get_sites_per_variant(unique_variants, sites, cosmic_output_file)


def process_vcf(vcf_file, unique_variants, truncate_large_var, max_allele_len):
    gz_decomp = get_gzip_app()
    print('Calculating file size...')
    row_count = int(sh.bash("-c", f"{gz_decomp} {vcf_file} | wc -l"))
    print('Parsing coding cosmic VCF file...')
    counter = 0
    with progressbar.ProgressBar(max_value=row_count) as pbar:
        for vcf_entry in read_vcf(vcf_file):
            temp_arr = vcf_entry.split('\t')
            variant = temp_arr[0] + ':' + temp_arr[1] + ':' + temp_arr[3] + ':' + temp_arr[4]
            cosmic_id = temp_arr[2]
            # if truncation is activated filter variants with 
            # ref or alt fields > max allele length
            if truncate_large_var:
                if len(temp_arr[3]) <= max_allele_len and len(temp_arr[4]) <= max_allele_len:
                    unique_variants.setdefault(variant, []).append(cosmic_id)
            else:
                unique_variants.setdefault(variant, []).append(cosmic_id)
            counter += 1
            pbar.update(counter)
    print('Parsing complete.')
    

def process_mut_export(filename, sites):
    # run a shell command to create intermediate file
    sites_file = 'sites.tmp'
    print('Preprocessing cosmic mutant export file...')
    # check operating system
    gz_decomp = get_gzip_app()
    sh.bash("-c", f"{gz_decomp} {filename} | cut -f 7,8,17 >{sites_file}")
    print('Calculating file size...')
    row_count = int(sh.bash("-c", f"cat {sites_file} | wc -l"))
    
    # parse the intermediate file
    print ('Parsing intermediate file...')
    counter = 0
    with progressbar.ProgressBar(max_value=row_count) as pbar:
        for sites_item in read_sites_file(sites_file):
            temp_arr = sites_item.strip('\n').split('\t')
            cosmic_id = temp_arr[2]
            site = temp_arr[1]
            tumor_id = temp_arr[0]
            if cosmic_id in sites:
                if site in sites[cosmic_id]:
                    sites[cosmic_id][site].add(tumor_id)
                else:
                    sites[cosmic_id][site] = {tumor_id}
            else:
                temp_dict = {}
                temp_dict[site] = {tumor_id}
                sites[cosmic_id] = temp_dict
            counter += 1
            pbar.update(counter)
    print('Parsing complete.')
    print('Removing temporary file...')
    # sh.rm(sites_file)
    

def get_sites_per_variant(unique_variants, sites, output_file):
    print ('Processing unique variants and associated sites...')
    with gzip.open(output_file, 'wb') as outf:
        # output buffer before writing to file
        outbuffer = []
        buffer_cap = 10000
        row_count = 0
        # parse the unique variant list
        counter = 0
        var_count = len(unique_variants)
        with progressbar.ProgressBar(max_value=var_count) as pbar:
            for variant, ids in unique_variants.items():
                # hold column element of each row in a list
                var_line_ele = variant.split(':')
                var_line_ele.append(';'.join(ids))
                # parse each of the cosmic ids to extract the site counts
                temp_site_info = {}
                print(variant)
                print(ids)
                for cosmic_id in ids:
                    print(cosmic_id)
                    for site, tumor_ids in sites[cosmic_id].items():
                        if site in temp_site_info:
                            temp_site_info[site] = temp_site_info[site] | tumor_ids
                        else:
                            temp_site_info[site] = tumor_ids
                # parse the temp_site_info
                site_list = []
                for site,tumor_ids in temp_site_info.items():
                    site_list.append(site + '=' + str(len(tumor_ids)))
                var_line_ele.append(';'.join(site_list))
                # check if outbuffer is full
                if row_count > buffer_cap:
                    # write to file
                    outf.write('\n'.join(outbuffer) + '\n')
                    # clear the buffer
                    del outbuffer[:]
                    # reset row_count
                    row_count = 0
                # add current line to buffer
                outbuffer.append('\t'.join(var_line_ele))
                row_count += 1
                counter += 1
                pbar.update(counter)
                # print 'Progress: {0:.2f}%      \r'.format((counter/var_count) * 100),
        # write the last set of lines to the file
        outf.write('\n'.join(outbuffer) + '\n')
        # clear the buffer
        del outbuffer[:]
        print(f'Cosmic files processed. Cosmic mutation and site counts written to {output_file}.')
        

def get_gzip_app():
    gz_app = ''
    if sys.platform.startswith("linux"):
        gz_app = "zcat"
    elif sys.platform.startswith("darwin"):
        gz_app = "gunzip -c"
    else:
        pass
    return gz_app


def read_vcf(vcf_file):
    with gzip.open(vcf_file, 'rb') as vcf:
        for line in vcf:
            linestr = line.decode('UTF-8')
            if not linestr.startswith('#'):
                yield linestr


def read_sites_file(filename):
    with open(filename, 'r') as input_file:
        for line in input_file:
            if not 'Primary' in line:
                yield line


if __name__ == '__main__':
    process_cosmic_db(sys.argv[1:])

from collections import Counter
import sys
import os
import gzip
import sh
import progressbar


def process_cosmic_db(args):
    if len(args) < 2:
        print ('Insufficient arguments.')
        print('Expected: python process_cosmic_data.py <cosmic_coding.vcf.gz> <cosmic_mutant_export.tsv.gz><truncate large variant [true or false]>')
        print ('Truncate large variant is optional. Default = False. If true, removes variants with ref or alt fields >1000bp in size')
        sys.exit(1)

    # COSMIC coding mut VCF file
    cosmic_coding_vcf = args[0]

    # COSMIC mut export file
    cosmic_export = args[1]

    # Truncate large variants beyond a defined size (default = 1000bp)
    truncate_large_var = False
    max_allele_len = 1000
    if len(args) > 2:
        truncate_var = args[2].lower()
        if truncate_var == 'true':
            truncate_large_var = True
        elif truncate_var == 'false':
            truncate_large_var = False
        else:
            print ('Invalid truncate variant argument.')
            print('Expected: python process_cosmic_data.py <cosmic_coding.vcf.gz> <cosmic_mutant_export.tsv.gz><truncate large variant [true or false]>')
            print ('Truncate large variant is optional. Default = False. If true, removes variants with ref or alt fields >1000bp in size')
            sys.exit(1)


    # cosmic processed file
    cosmic_output_file = "cosmic_processed.txt.gz"

    # process vcf file first
    unique_variants = {}
    process_vcf(cosmic_coding_vcf, unique_variants, truncate_large_var, max_allele_len)

    # process mut export file
    sites = {}
    process_mut_export(cosmic_export, sites)

    # merge unique variants and sites count data and write to output file
    get_sites_per_variant(unique_variants, sites, cosmic_output_file)


def process_vcf(vcf_file, unique_variants, truncate_large_var, max_allele_len):
    gz_decomp = get_gzip_app()
    print('Calculating file size...')
    row_count = int(sh.bash("-c", f"{gz_decomp} {vcf_file} | wc -l"))
    print('Parsing coding cosmic VCF file...')
    counter = 0
    with progressbar.ProgressBar(max_value=row_count) as pbar:
        for vcf_entry in read_vcf(vcf_file):
            temp_arr = vcf_entry.split('\t')
            variant = temp_arr[0] + ':' + temp_arr[1] + ':' + temp_arr[3] + ':' + temp_arr[4]
            cosmic_id = temp_arr[2]
            # if truncation is activated filter variants with 
            # ref or alt fields > max allele length
            if truncate_large_var:
                if len(temp_arr[3]) <= max_allele_len and len(temp_arr[4]) <= max_allele_len:
                    unique_variants.setdefault(variant, []).append(cosmic_id)
            else:
                unique_variants.setdefault(variant, []).append(cosmic_id)
            counter += 1
            pbar.update(counter)
    print('Parsing complete.')
    

def process_mut_export(filename, sites):
    # run a shell command to create intermediate file
    sites_file = 'sites.tmp'
    print('Preprocessing cosmic mutant export file...')
    # check operating system
    gz_decomp = get_gzip_app()
    sh.bash("-c", f"{gz_decomp} {filename} | cut -f 7,8,17 >{sites_file}")
    print('Calculating file size...')
    row_count = int(sh.bash("-c", f"cat {sites_file} | wc -l"))
    
    # parse the intermediate file
    print ('Parsing intermediate file...')
    counter = 0
    with progressbar.ProgressBar(max_value=row_count) as pbar:
        for sites_item in read_sites_file(sites_file):
            temp_arr = sites_item.strip('\n').split('\t')
            cosmic_id = temp_arr[2]
            site = temp_arr[1]
            tumor_id = temp_arr[0]
            if cosmic_id in sites:
                if site in sites[cosmic_id]:
                    sites[cosmic_id][site].add(tumor_id)
                else:
                    sites[cosmic_id][site] = {tumor_id}
            else:
                temp_dict = {}
                temp_dict[site] = {tumor_id}
                sites[cosmic_id] = temp_dict
            counter += 1
            pbar.update(counter)
    print('Parsing complete.')
    print('Removing temporary file...')
    # sh.rm(sites_file)
    

def get_sites_per_variant(unique_variants, sites, output_file):
    print ('Processing unique variants and associated sites...')
    with gzip.open(output_file, 'wb') as outf:
        # output buffer before writing to file
        outbuffer = []
        buffer_cap = 10000
        row_count = 0
        # parse the unique variant list
        counter = 0
        var_count = len(unique_variants)
        with progressbar.ProgressBar(max_value=var_count) as pbar:
            for variant, ids in unique_variants.items():
                # hold column element of each row in a list
                var_line_ele = variant.split(':')
                var_line_ele.append(';'.join(ids))
                # parse each of the cosmic ids to extract the site counts
                temp_site_info = {}
                print(variant)
                print(ids)
                for cosmic_id in ids:
                    print(cosmic_id)
                    for site, tumor_ids in sites[cosmic_id].items():
                        if site in temp_site_info:
                            temp_site_info[site] = temp_site_info[site] | tumor_ids
                        else:
                            temp_site_info[site] = tumor_ids
                # parse the temp_site_info
                site_list = []
                for site,tumor_ids in temp_site_info.items():
                    site_list.append(site + '=' + str(len(tumor_ids)))
                var_line_ele.append(';'.join(site_list))
                # check if outbuffer is full
                if row_count > buffer_cap:
                    # write to file
                    outf.write('\n'.join(outbuffer) + '\n')
                    # clear the buffer
                    del outbuffer[:]
                    # reset row_count
                    row_count = 0
                # add current line to buffer
                outbuffer.append('\t'.join(var_line_ele))
                row_count += 1
                counter += 1
                pbar.update(counter)
                # print 'Progress: {0:.2f}%      \r'.format((counter/var_count) * 100),
        # write the last set of lines to the file
        outf.write('\n'.join(outbuffer) + '\n')
        # clear the buffer
        del outbuffer[:]
        print(f'Cosmic files processed. Cosmic mutation and site counts written to {output_file}.')
        

def get_gzip_app():
    gz_app = ''
    if sys.platform.startswith("linux"):
        gz_app = "zcat"
    elif sys.platform.startswith("darwin"):
        gz_app = "gunzip -c"
    else:
        pass
    return gz_app


def read_vcf(vcf_file):
    with gzip.open(vcf_file, 'rb') as vcf:
        for line in vcf:
            linestr = line.decode('UTF-8')
            if not linestr.startswith('#'):
                yield linestr


def read_sites_file(filename):
    with open(filename, 'r') as input_file:
        for line in input_file:
            if not 'Primary' in line:
                yield line


if __name__ == '__main__':
    process_cosmic_db(sys.argv[1:])
# COSMIC coding mut VCF file # COSMIC mut export file # Truncate large variants beyond a defined size (default = 1000bp) # cosmic processed file # process vcf file first # process mut export file # merge unique variants and sites count data and write to output file # if truncation is activated filter variants with # ref or alt fields > max allele length # run a shell command to create intermediate file # check operating system # parse the intermediate file # sh.rm(sites_file) # output buffer before writing to file # parse the unique variant list # hold column element of each row in a list # parse each of the cosmic ids to extract the site counts # parse the temp_site_info # check if outbuffer is full # write to file # clear the buffer # reset row_count # add current line to buffer # print 'Progress: {0:.2f}%      \r'.format((counter/var_count) * 100), # write the last set of lines to the file # clear the buffer