# how to count the number of each allele of CNVs?

I have a bunch of files (vcf format) each containing site and sample annotations of a copy number variant (following example is part of one of my files, I am showing 11 individuals but they are nearly 900 individuals)

16  32328771    esv3638416;esv3638417;esv3638418;esv3638419;esv3638420  A   <CN0>,<CN2>,<CN3>,<CN4>,<CN5>   100 PASS    AC=27,84,3220,20,26;AF=0.00539137,0.0167732,0.642971,0.00399361,0.00519169;AN=5008;CS=DUP_gs;END=32370968;NS=2504;SVTYPE=CNV;DP=26381;EAS_AF=0.0099,0.0188,0.5427,0.004,0.0089;AMR_AF=0.0101,0.0058,0.4063,0,0.0014;AFR_AF=0.0015,0.0371,0.8374,0.003,0.003;EUR_AF=0.003,0.0089,0.5189,0.005,0.005;SAS_AF=0.0051,0.0031,0.7791,0.0072,0.0072;VT=SV  GT  3|3 3|1 0|3 3|3 3|3 3|0 0|3 3|3 0|0 3|3 1|3


I want to count each allele of CNV in each individual(sample) if there are 0, 1 or 2 copies! this is my desired output (rows are individuals):

CN0 CN2 CN3 CN4 CN5
0   0   2   0   0
1   0   1   0   0
0   0   1   0   0
0   0   2   0   0
0   0   2   0   0
0   0   1   0   0
0   0   1   0   0
0   0   2   0   0
0   0   0   0   0
0   0   2   0   0
1   0   1   0   0


Any idea how to do that?

• Cold you explain your file format a bit? How do we get from the VCF line you show to the table of your desired output? Jan 2 '18 at 10:04
• just edited my post to clarify that it is a vcf form file. I used vcftools --poition option to extract the position that I needed! and how to create the table is exactly the question that I have .....! Jan 3 '18 at 0:20

Here's a pure python solution to process your example. It should be straightforward to loop over all the rows of your vcf files

my_line='16 32328771    esv3638416;esv3638417;esv3638418;esv3638419;esv3638420  A   <CN0>,<CN2>,<CN3>,<CN4>,<CN5>   100 PASS    AC=27,84,3220,20,26;AF=0.00539137,0.0167732,0.642971,0.00399361,0.00519169;AN=5008;CS=DUP_gs;END=32370968;NS=2504;SVTYPE=CNV;DP=26381;EAS_AF=0.0099,0.0188,0.5427,0.004,0.0089;AMR_AF=0.0101,0.0058,0.4063,0,0.0014;AFR_AF=0.0015,0.0371,0.8374,0.003,0.003;EUR_AF=0.003,0.0089,0.5189,0.005,0.005;SAS_AF=0.0051,0.0031,0.7791,0.0072,0.0072;VT=SV  GT  3|3 3|1 0|3 3|3 3|3 3|0 0|3 3|3 0|0 3|3 1|3 '
my_line=my_line.strip() # remove trailing whitespace
my_fields=my_line.split()
names=my_fields[4].split(',') # extract CNV names
names=[x[1:-1] for x in names] # remove '<>' symbols
print('\t'.join(names))
for sample in my_fields[9:]:
genotypes=[0]*(len(names)+1) # initialize array of zero counts, including reference
for allele in sample.split('|'):
genotypes[int(allele)]+=1
print('\t'.join([str(x) for x in genotypes[1:]])) # print all columns except reference

• Thanks @heathobrien, this may work but in the example above I just showed genotypes for a few individuals, I have nearly 900 individuals and It is not convenient to copy all genotypes in your code in "my_line" part! Do you have any alternative to be able to fix this? Jan 4 '18 at 20:45
• You'll probably want to include the code I wrote as part of a larger script that reads your genotypes file and writes the output to another file rather than printing it to the screen. This isn't a script writing service, so take a look at an python I/O tutorial or the documentation for the csv module (or pysam) Jan 4 '18 at 21:23