2
$\begingroup$

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?

Improve this question
$\endgroup$
2
  • $\begingroup$ 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? $\endgroup$
    – terdon ♦
    Jan 2, 2018 at 10:04
  • $\begingroup$ 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 .....! $\endgroup$
    – Anna1364
    Jan 3, 2018 at 0:20

1 Answer 1

Reset to default
1
$\begingroup$

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
Improve this answer
$\endgroup$
2
  • $\begingroup$ 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? $\endgroup$
    – Anna1364
    Jan 4, 2018 at 20:45
  • $\begingroup$ 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) $\endgroup$
    – heathobrien
    Jan 4, 2018 at 21:23

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge that you have read and understand our privacy policy and code of conduct.

Not the answer you're looking for? Browse other questions tagged or ask your own question.