other/pipeline/trunk/mut_stats.py

Code
Comments
Other
Rev Date Author Line
5767 02 Dec 19 nicklas 1 #!/usr/bin/python
5767 02 Dec 19 nicklas 2 """
5767 02 Dec 19 nicklas 3 $Id $
5767 02 Dec 19 nicklas 4
5767 02 Dec 19 nicklas 5 Script that calculates statistics for variants in a list of VCF files.
5771 03 Dec 19 nicklas 6 The script need three parameters:
5767 02 Dec 19 nicklas 7
5769 03 Dec 19 nicklas 8 1: Path to text file with list of VCF files to process (se below)
5769 03 Dec 19 nicklas 9 2: Path to progress reporting file
5771 03 Dec 19 nicklas 10 3: Path to a log file
5769 03 Dec 19 nicklas 11
5769 03 Dec 19 nicklas 12 The text file given as the first parameter should be a tab-separated 
5769 03 Dec 19 nicklas 13 text file. Each line should have three columns:
5769 03 Dec 19 nicklas 14
5767 02 Dec 19 nicklas 15 1: Patient identifier
5767 02 Dec 19 nicklas 16 2: Path to VCF file
5769 03 Dec 19 nicklas 17 3: Name of alignment (used for progress reporting)
5767 02 Dec 19 nicklas 18
5767 02 Dec 19 nicklas 19 Counts and frequencies for all variants will calculated for the total and
5767 02 Dec 19 nicklas 20 per-patient. 
5767 02 Dec 19 nicklas 21
5767 02 Dec 19 nicklas 22 The output is a VCF file with all variants and the following annotations:
5767 02 Dec 19 nicklas 23 CNT_ALL: The number of VCF files the variants was seen in
5767 02 Dec 19 nicklas 24 MF_ALL: The mutation frequency among all VCF files
5767 02 Dec 19 nicklas 25 CNT_NR: The number of patients the variants was seen at least one time
5767 02 Dec 19 nicklas 26 MF_NR: The mutation frequency among all patients
5767 02 Dec 19 nicklas 27
5767 02 Dec 19 nicklas 28 """
5767 02 Dec 19 nicklas 29
5767 02 Dec 19 nicklas 30 import sys
5767 02 Dec 19 nicklas 31 import gzip
5767 02 Dec 19 nicklas 32 import datetime
5769 03 Dec 19 nicklas 33 import time
5767 02 Dec 19 nicklas 34
5767 02 Dec 19 nicklas 35 # Store all variants that we have seen in the VCF files we load
5767 02 Dec 19 nicklas 36 # The 'key' is a combination of chromosome, position, id, ref and alt columns
5767 02 Dec 19 nicklas 37 # The 'value' is an array that is counting the number of times we have seen the variant
5767 02 Dec 19 nicklas 38 # [0] = Total number of times including replicates
5767 02 Dec 19 nicklas 39 # [1] = Number of times in unique patients
5767 02 Dec 19 nicklas 40 allVariants = {}
5767 02 Dec 19 nicklas 41
5767 02 Dec 19 nicklas 42 # The total number of samples, 0=including replicates, 1=only non-replicates
5767 02 Dec 19 nicklas 43 tots = [0, 0]
5767 02 Dec 19 nicklas 44
5767 02 Dec 19 nicklas 45 # Store all variants that we have seen for the current patient
5767 02 Dec 19 nicklas 46 # The 'key' is the same as for the allVarints but we only store
5767 02 Dec 19 nicklas 47 # single count value
5767 02 Dec 19 nicklas 48 currentVariants = {}
5767 02 Dec 19 nicklas 49 currentVariants['name'] = ''
5767 02 Dec 19 nicklas 50
5767 02 Dec 19 nicklas 51 def loadVcf(path, patient):
5767 02 Dec 19 nicklas 52     tots[0] += 1
5767 02 Dec 19 nicklas 53     
5767 02 Dec 19 nicklas 54     # If this is the first time we see the patient
5767 02 Dec 19 nicklas 55     # the current information is reset
5767 02 Dec 19 nicklas 56     if patient != currentVariants['name']:
5767 02 Dec 19 nicklas 57         currentVariants.clear()
5767 02 Dec 19 nicklas 58         currentVariants['name'] = patient
5767 02 Dec 19 nicklas 59         tots[1] += 1
5767 02 Dec 19 nicklas 60     
5767 02 Dec 19 nicklas 61     try:
5767 02 Dec 19 nicklas 62         # Open the VCF file and read it line by line
5767 02 Dec 19 nicklas 63         vcf = gzip.open(path, 'rt')
5767 02 Dec 19 nicklas 64         line = vcf.readline()
5767 02 Dec 19 nicklas 65         while line:
5767 02 Dec 19 nicklas 66             if not line.startswith("#"):
5767 02 Dec 19 nicklas 67                 cols = line.rstrip().split('\t')
5767 02 Dec 19 nicklas 68                 # Key is <chr>|<pos>|<id>|<ref>|<alt>
5767 02 Dec 19 nicklas 69                 key = '|'.join(cols[i] for i in [0,1,2,3,4])
5767 02 Dec 19 nicklas 70                 
5767 02 Dec 19 nicklas 71                 # Update the total count
5767 02 Dec 19 nicklas 72                 if key in allVariants:
5767 02 Dec 19 nicklas 73                     allVariants[key][0] += 1
5767 02 Dec 19 nicklas 74                 else:
5767 02 Dec 19 nicklas 75                     allVariants[key] = [1, 0]
5767 02 Dec 19 nicklas 76                     
5767 02 Dec 19 nicklas 77                 # Update count for this patient
5767 02 Dec 19 nicklas 78                 if key in currentVariants:
5767 02 Dec 19 nicklas 79                     currentVariants[key] += 1
5767 02 Dec 19 nicklas 80                 else:
5767 02 Dec 19 nicklas 81                     currentVariants[key] = 1
5767 02 Dec 19 nicklas 82                     allVariants[key][1] += 1
5767 02 Dec 19 nicklas 83                 
5767 02 Dec 19 nicklas 84             line = vcf.readline()
5767 02 Dec 19 nicklas 85     finally:
5767 02 Dec 19 nicklas 86         vcf.close()
5767 02 Dec 19 nicklas 87
5769 03 Dec 19 nicklas 88 # Report progress to the progressFile (1-90%)
5769 03 Dec 19 nicklas 89 def progressReport(progressFile, alignment, current, total):
5769 03 Dec 19 nicklas 90     percent = 1+current * 89 /  total
5769 03 Dec 19 nicklas 91     with open(progressFile, 'w') as p:
5769 03 Dec 19 nicklas 92         p.write("{0} Reading VCF file for '{1}' ({2} of {3})".format(percent, alignment, current, total))
5769 03 Dec 19 nicklas 93
5767 02 Dec 19 nicklas 94 # Read the patient/VCF list
5767 02 Dec 19 nicklas 95 # Each line should be tab-separated <PAT>\t<PATH-TO-VCF>
5767 02 Dec 19 nicklas 96 vcfList = sys.argv[1]
5767 02 Dec 19 nicklas 97 with open(vcfList) as f:
5767 02 Dec 19 nicklas 98     lines = f.read().splitlines()
5767 02 Dec 19 nicklas 99
5767 02 Dec 19 nicklas 100 # Sort lines so that entries for the same patient is after each other
5767 02 Dec 19 nicklas 101 # The counting algorithm depends on it
5767 02 Dec 19 nicklas 102 lines.sort()
5767 02 Dec 19 nicklas 103
5769 03 Dec 19 nicklas 104 progressFile = sys.argv[2]
5771 03 Dec 19 nicklas 105 logFile = sys.argv[3]
5769 03 Dec 19 nicklas 106
5767 02 Dec 19 nicklas 107 # Load the VCF files and count the variants in them
5769 03 Dec 19 nicklas 108 vcfCount = 0
5769 03 Dec 19 nicklas 109 totalVcf = len(lines)
5769 03 Dec 19 nicklas 110 nextProgressReport = time.time()+15
5769 03 Dec 19 nicklas 111
5767 02 Dec 19 nicklas 112 for line in lines:
5769 03 Dec 19 nicklas 113     cols = line.split('\t')    
5769 03 Dec 19 nicklas 114     vcfCount += 1
5769 03 Dec 19 nicklas 115     # Report progress every 15 seconds
5769 03 Dec 19 nicklas 116     if time.time() > nextProgressReport:
5769 03 Dec 19 nicklas 117         progressReport(progressFile, cols[2], vcfCount, totalVcf)
5769 03 Dec 19 nicklas 118         nextProgressReport = time.time()+15
5767 02 Dec 19 nicklas 119     loadVcf(cols[1], cols[0])
5767 02 Dec 19 nicklas 120
5771 03 Dec 19 nicklas 121 # Write the log file
5771 03 Dec 19 nicklas 122 with open(logFile, 'w') as log:
5771 03 Dec 19 nicklas 123     log.write("NumberOfVCF: {0}\n".format(tots[0]))
5771 03 Dec 19 nicklas 124     log.write("NumberOfPatients: {0}\n".format(tots[1]))
5771 03 Dec 19 nicklas 125     log.write("NumberOfVariants: {0}\n".format(len(allVariants)))
5771 03 Dec 19 nicklas 126
5767 02 Dec 19 nicklas 127 # Generate header for the output VCF
5767 02 Dec 19 nicklas 128 header = '''##fileformat=VCFv4.2
5767 02 Dec 19 nicklas 129 ##build={0}
5767 02 Dec 19 nicklas 130 '''.format(datetime.datetime.today().strftime('%Y-%m-%d'))
5767 02 Dec 19 nicklas 131 header += '''##INFO=<ID=CNT_ALL,Number=1,Type=Integer,Description="Number of samples with this mutation among ALL ({all}) samples">
5767 02 Dec 19 nicklas 132 ##INFO=<ID=CNT_NR,Number=1,Type=Integer,Description="Number of NON-replicate ({nonrep}) samples with this mutation">
5767 02 Dec 19 nicklas 133 ##INFO=<ID=MF_ALL,Number=1,Type=Float,Description="Mutant frequency among ALL ({all}) samples">
5767 02 Dec 19 nicklas 134 ##INFO=<ID=MF_NR,Number=1,Type=Float,Description="Mutant frequency among NON-replicate ({nonrep}) samples">
5767 02 Dec 19 nicklas 135 #CHROM  POS     ID      REF     ALT     QUAL    FILTER  INFO'''.format(all=tots[0], nonrep=tots[1])
5767 02 Dec 19 nicklas 136
5767 02 Dec 19 nicklas 137 print(header)
5767 02 Dec 19 nicklas 138
5767 02 Dec 19 nicklas 139 # Print out all variants that was found
5767 02 Dec 19 nicklas 140 for key in allVariants.keys():
5767 02 Dec 19 nicklas 141     cols = key.split('|')
5767 02 Dec 19 nicklas 142     vals = (cols + ['.', '.'] + allVariants[key] + [ float(num) / tot for num, tot in zip(allVariants[key], tots) ] )
5767 02 Dec 19 nicklas 143     print("{0}\t{1}\t{2}\t{3}\t{4}\t{5}\t{6}\tCNT_ALL={7};CNT_NR={8};MF_ALL={9:.4};MF_NR={10:.4}".format(*vals))
5767 02 Dec 19 nicklas 144

Generated on Tue Nov 21 03:40:16 2023 (UTC) by svndigest 0.11.2pre (dev build 1672)