gatk-3.8/python/expandedSummaryToGeli.py

#!/usr/bin/env python
import farm_commands
import os.path
import sys
import re


def grepFile(string, file):
    return grep(string, file.readlines())

def grep(string, list):
    expr = re.compile(string)
    return [elem for elem in list if expr.match(elem)]

callFile = open(sys.argv[1])
pileupFileLines = open(sys.argv[2])
callFileWithLodScores = open(sys.argv[3])
pileupFileLines = pileupFileLines.readlines()
lodFileLines = callFileWithLodScores.readlines()

for line in callFile:
# note: file is a .csv; so comma delimited with headers
# chr:position,alleles,gene,type,rs#,base_chg,annot,dist,pop_nr_freq,f_nr_case,f_nr_con,chisqstat_assoc,lrt_stat,f+,f-,SLOD,p_val,min_snp_dist,min_amp_dist,num_sig_pools,mean_worst_dir_nonref, ...
# ..., max_worst_dir_nonref,mean_diff_dir_nonref,max_diff_dir_nonref,mean_other/minor,max_other/minor,mean_frac_non_concord,max_frac_non_concord,mean_combined_lod,max_combined_lod,mean_cov_tot, ...
# ..., max_cov_tot,mean_cov_diff,max_cov_diff,mean_lod,max_lod,mean_lod_diff,max_lod_diff,mean_lod_diff_norm,max_lod_diff_norm,target_name,sig_pools

# call file assumed already to have been split by pool

# pileupFile is a .bam.coverage file from the pooled pipeline

# call file isn't sorted by chr:position 
    if line.startswith("chr:pos"):
        continue
    elif not line.startswith("chr"):
        continue
    else:
        #print(line)
        entries = line.split(",")
        chrompos = entries.pop(0) # (now 'alleles' is 0-element)
        variant = entries.pop(0) # (now 'gene' is 0-element)
        #print(chrompos)
        #print(variant)
        #change format of "chr:pos" to "chr pos" for lookup
        #somehow this next line don't work
        #chrompos.replace(":"," ")
        g = chrompos.split(":");
        chromposspc = g.pop(0)+" "+g.pop(0)
        #print(chrompos)
        pileupLine = grep(chromposspc,pileupFileLines)
        lodLine = grep(chrompos,lodFileLines)
        #print(pileupLine)
        # line is
        # chr pos ref num_A num_C num_G num_T
        if not pileupLine:
            continue
        else:
            pileupLine = pileupLine.pop(0)
            #print(pileupLine)
            pileupList = pileupLine.split(" ")
            ref = pileupList.pop(2) # now num_A is 2
            num_A = int(pileupList.pop(2))
            num_C = int(pileupList.pop(2))
            num_G = int(pileupList.pop(2))
            num_T = int(pileupList.pop(2))
            depth = num_A+num_C+num_G+num_T
            lodLine = lodLine.pop(0)
            lodLine = lodLine.split(" ")
            # print(lodLine)
            lod = max(float(lodLine.pop(21)),float(lodLine.pop(19)))
            # output is
            # chr pos ref depth mapping call btr btnb AA AC AG AT CC CG CT GG GT TT
            outStr = chrompos+" "+ref+" "+str(depth)+" 5 "+variant+" "+str(lod)+" -2 -3 -4 -5 -6 -7 -8 -9 -10 -11 -12"
            print(outStr)