How to aggregate values over a bigger than RAM gzip'ed csv file? <p>For starters I am new to bioinformatics and especially to programming, but I have built a script that will go through a so-called VCF file (only the individuals are included, one clumn = one individual), and uses a search string to find out for every variant (line) whether the individual is homozygous or heterozygous.</p> <p>This script works, at least on small subsets, but I know it stores everything in memory. I would like to do this on very large zipped files (of even whole genomes), but I do not know how to transform this script into a script that does everything line by line (because I want to count through whole columns I just do not see how to solve that).</p> <p>So the output is 5 things per individual (total variants, number homozygote, number heterozygote, and proportions of homo- and heterozygotes). See the code below:</p> <pre><code>#!usr/bin/env python import re import gzip subset_cols = 'subset_cols_chr18.vcf.gz' #nuc_div = 'nuc_div_chr18.txt' gz_infile = gzip.GzipFile(subset_cols, "r") #gz_outfile = gzip.GzipFile(nuc_div, "w") # make a dictionary of the header line for easy retrieval of elements later on headers = gz_infile.readline().rstrip().split('\t') print headers column_dict = {} for header in headers: column_dict[header] = [] for line in gz_infile: columns = line.rstrip().split('\t') for i in range(len(columns)): c_header=headers[i] column_dict[c_header].append(columns[i]) #print column_dict for key in column_dict: number_homozygotes = 0 number_heterozygotes = 0 for values in column_dict[key]: SearchStr = '(\d)/(\d):\d+,\d+:\d+:\d+:\d+,\d+,\d+' #this search string contains the regexp (this regexp was tested) Result = re.search(SearchStr,values) if Result: #here, it will skip the missing genoytypes ./. variant_one = int(Result.group(1)) variant_two = int(Result.group(2)) if variant_one == 0 and variant_two == 0: continue elif variant_one == variant_two: #count +1 in case variant one and two are equal (so 0/0, 1/1, etc.) number_homozygotes += 1 elif variant_one != variant_two: #count +1 in case variant one is not equal to variant two (so 1/0, 0/1, etc.) number_heterozygotes += 1 print "%s homozygotes %s" % (number_homozygotes, key) print "%s heterozygotes %s" % (number_heterozygotes,key) variants = number_homozygotes + number_heterozygotes print "%s variants" % variants prop_homozygotes = (1.0*number_homozygotes/variants)*100 prop_heterozygotes = (1.0*number_heterozygotes/variants)*100 print "%s %% homozygous %s" % (prop_homozygotes, key) print "%s %% heterozygous %s" % (prop_heterozygotes, key) </code></pre> <p>Any help will be much appreciated so I can go on investigating large datasets, thank you :)</p> <p>The VCF file by the way looks something like this: INDIVIDUAL_1 INDIVIDUAL_2 INDIVIDUAL_3 0/0:9,0:9:24:0,24,221 1/0:5,4:9:25:25,0,26 1/1:0,13:13:33:347,33,0</p> <p>This is then the header line with the individual ID names (I have in total 33 individuals with more complicated ID tags, I simplified here) and then I have a lot of these information lines with the same specific pattern. I am only interested in the first part with the slash so hence the regular epxression.</p>