SEQanswers

Go Back   SEQanswers > Bioinformatics > Bioinformatics



Similar Threads
Thread Thread Starter Forum Replies Last Post
combining fastq files mmmm Bioinformatics 9 02-21-2014 02:54 AM
combining htseq-count files into one mbk0asis General 6 06-17-2013 07:07 AM
Combining pileup and gff3 files saemi Bioinformatics 0 05-16-2012 09:22 AM
Comparing and combining vcf files gavin.oliver Bioinformatics 3 11-30-2011 02:16 AM
Combining *.gtf files? jhb1980 Bioinformatics 2 11-01-2011 12:43 AM

Reply
 
Thread Tools
Old 05-21-2014, 11:45 AM   #1
prs321
Member
 
Location: US

Join Date: Jun 2013
Posts: 96
Default Combining VCF files using GATK?

Here is my task. I have 4 treatment groups containing 5 replicates in each group.

For each treatment group, I am looking for at least 2 replicates that have the same base substitution.

So far I have generated a vcf file for each replicate in each treatment group.

1. What is the best way to combined these files so that I have 4 vcf files, 1 for each treatment group that contains 5 replicates?

2. What is the best way to keep track of the position of sites where at least 2 replicates showed the same base substitution so that I can produce results once the gene annotation is complete?
prs321 is offline   Reply With Quote
Old 05-21-2014, 11:56 AM   #2
lethalfang
Member
 
Location: San Francisco, CA

Join Date: Aug 2011
Posts: 91
Default

http://www.broadinstitute.org/gatk/g...eVariants.html

If you want to only show variants in 2+ duplicates, use "--minimumN 2" option.
lethalfang is offline   Reply With Quote
Old 05-21-2014, 11:58 AM   #3
prs321
Member
 
Location: US

Join Date: Jun 2013
Posts: 96
Default

The issue seems to be that when I combine all 5 vcf files for the replicates, it shows a replicate at 1 site, but shows a '.' for the other 4.

I'm not really sure how to go about this.
prs321 is offline   Reply With Quote
Old 05-21-2014, 12:01 PM   #4
lethalfang
Member
 
Location: San Francisco, CA

Join Date: Aug 2011
Posts: 91
Default

Quote:
Originally Posted by prs321 View Post
The issue seems to be that when I combine all 5 vcf files for the replicates, it shows a replicate at 1 site, but shows a '.' for the other 4.

I'm not really sure how to go about this.
Does that mean only 1 replicate has that called variant?
lethalfang is offline   Reply With Quote
Old 05-21-2014, 12:04 PM   #5
prs321
Member
 
Location: US

Join Date: Jun 2013
Posts: 96
Default

Yes so for example I have 5 alignment files: rep1.bam, rep2.bam, rep3.bam, rep4.bam, rep5.bam

All of them have been processed (markduplicates, addreadgroups etc).

Then I called snps on each alignment to produce the following: rep1.vcf, rep2.vcf, rep3.vcf, rep4.vcf, rep5.vcf
prs321 is offline   Reply With Quote
Old 05-21-2014, 12:19 PM   #6
lethalfang
Member
 
Location: San Francisco, CA

Join Date: Aug 2011
Posts: 91
Default

Quote:
Originally Posted by prs321 View Post
Yes so for example I have 5 alignment files: rep1.bam, rep2.bam, rep3.bam, rep4.bam, rep5.bam

All of them have been processed (markduplicates, addreadgroups etc).

Then I called snps on each alignment to produce the following: rep1.vcf, rep2.vcf, rep3.vcf, rep4.vcf, rep5.vcf

What will this give you?

Code:
java -Xmx2g -jar GenomeAnalysisTK.jar \
   -R ref.fasta \
   -T CombineVariants \
   -minN 2 \
   --variant rep1.vcf \
   --variant rep2.vcf \
   --variant rep3.vcf \
   --variant rep4.vcf \
   --variant rep5.vcf \
   -o combined.output.vcf
lethalfang is offline   Reply With Quote
Old 05-22-2014, 07:08 AM   #7
prs321
Member
 
Location: US

Join Date: Jun 2013
Posts: 96
Default

0/0:6:21:24,0,100:44:299:15 . . . .
prs321 is offline   Reply With Quote
Old 05-22-2014, 07:29 AM   #8
lethalfang
Member
 
Location: San Francisco, CA

Join Date: Aug 2011
Posts: 91
Default

There is only one such column? If you combine 5 samples (replicates), there should be 5 such columns.
If not, make sure the header for the sample info are different in each of the 5 vcf files. Otherwise GATK thinks it's the same sample (i.e., same replicate in this case).
0/0 means reference calls in both copies.

Last edited by lethalfang; 05-22-2014 at 07:33 AM.
lethalfang is offline   Reply With Quote
Old 05-22-2014, 09:52 AM   #9
prs321
Member
 
Location: US

Join Date: Jun 2013
Posts: 96
Default

There are actually 5 columns, but only the first column is listed with the genotype. The other 4 have a '.' under them.
prs321 is offline   Reply With Quote
Old 05-22-2014, 09:53 AM   #10
lethalfang
Member
 
Location: San Francisco, CA

Join Date: Aug 2011
Posts: 91
Default

Quote:
Originally Posted by prs321 View Post
There are actually 5 columns, but only the first column is listed with the genotype. The other 4 have a '.' under them.
Can you show me a few rows of the output?
lethalfang is offline   Reply With Quote
Old 05-22-2014, 10:06 AM   #11
prs321
Member
 
Location: US

Join Date: Jun 2013
Posts: 96
Default

#CHROM POS ID REF ALT QUAL FILTER INFO FORMAT cs1 cs2 cs3 cs4 cs5
NODE_4_length_282_cov_38.510639 66 . G A 169.04 . AC=1;AF=0.500;AN=2;CIGAR=1X;DP=66;DPRA=0;GTI=0;LEN=1;MQM=70;MQMR=70;NS=1;NUMALT=1;PAIREDR=1;PAO=0;PQA=0;PQR=0;PRO=0;RUN=1;TYPE=snp;set=Intersection GT:AOP:PL:QA:QR:RO 0/1:9:19:100,0
,100:271:306:10 . . . .
NODE_4_length_282_cov_38.510639 99 . C T 214.21 . AC=1;AF=0.500;AN=2;CIGAR=1X;DP=92;DPRA=0;GTI=0;LEN=1;MQM=70;MQMR=70;NS=1;NUMALT=1;PAIREDR=1;PAO=0;PQA=0;PQR=0;PRO=0;RUN=1;TYPE=snp;set=Intersection GT:AOP:PL:QA:QR:RO 0/1:12:24:100,
0,100:319:364:11 . . . .
NODE_4_length_282_cov_38.510639 114 . A G 346.87 . AC=1;AF=0.500;AN=2;CIGAR=1X;DP=102;DPRA=0;GTI=0;LEN=1;MEANALT=1;MQM=70;MQMR=70;NS=1;NUMALT=1;PAIREDR=1;PAO=0;PQA=0;PQR=0;PRO=0;RUN=1;TYPE=snp;set=Intersection GT:AOP:PL:QA:QR:RO 0/1:14
:26:100,0,100:475:413:12 . . . .
NODE_4_length_282_cov_38.510639 154 . G A 452.29 . AC=1;AF=0.500;AN=2;CIGAR=1X;DP=123;DPRA=0;GTI=0;LEN=1;MEANALT=1;MQM=70;MQMR=70;NS=1;NUMALT=1;PAIREDR=1;PAO=0;PQA=0;PQR=0;PRO=0;RUN=1;TYPE=snp;set=Intersection GT:AOP:PL:QA:QR:RO 0/1:18
:31:100,0,100:578:467:13 . . . .
NODE_4_length_282_cov_38.510639 247 . T G 355.03 . AC=1;AF=0.500;AN=2;CIGAR=1X;DP=131;DPRA=0;GTI=0;LEN=1;MEANALT=1;MQM=70;MQMR=70;NS=1;NUMALT=1;PAIRED=1;PAIREDR=1;PAO=0;PQA=0;PQR=0;PRO=0;RUN=1;TYPE=snp;set=Intersection GT:AOP:PL:QA:QR:RO
0/1:14:32:100,0,100:484:576:18 . . . .
NODE_9_length_16257_cov_15.910501 5514 . C A 248.89 . AB=0;ABP=0;AC=2;AF=1.00;AN=2;CIGAR=1X;DP=68;DPRA=0;EPPR=0;GTI=0;LEN=1;MEANALT=1;MQM=70;MQMR=0;NS=1;NUMALT=1;PAIRED=1;PAIREDR=0;PAO=0;PQA=0;PQR=0;PRO=0;QR=0;RO=0;RPPR=0;RUN=1;SRF=0;SR
P=0;SRR=0;TYPE=snp;set=Intersection GT:AOP:PL:QA:QR:RO 1/1:11:11:100,33,0:341:0:0 . . . .
NODE_9_length_16257_cov_15.910501 9753 . T G 0 . AB=0;ABP=0;AC=0;AF=0.00;AN=2;CIGAR=1X;DP=45;DPRA=0;GTI=0;LEN=1;MEANALT=1;MQM=70;MQMR=70;NS=1;NUMALT=1;PAIRED=1;PAIREDR=1;PAO=0;PQA=0;PQR=0;PRO=0;RUN=1;SAR=0;TYPE=snp;set=cs1-cs4
GT:AOP:PL:QA:QR:RO 0/0:6:28:7,0,100:39:587:22 . . . .
NODE_9_length_16257_cov_15.910501 12869 . T G 0 . AB=0;ABP=0;AC=0;AF=0.00;AN=2;CIGAR=1X;DP=85;DPRA=0;GTI=0;LEN=1;MEANALT=1;MQM=70;MQMR=70;NS=1;NUMALT=1;PAIRED=1;PAIREDR=1;PAO=0;PQA=0;PQR=0;PRO=0;RUN=1;TYPE=snp;set=cs1-cs2-cs3 GT:AO:
DP:PL:QA:QR:RO 0/0:8:28:34,0,100:59:485:20 . . . .
NODE_9_length_16257_cov_15.910501 14192 . T G 0 . AB=0;ABP=0;AC=0;AF=0.00;AN=2;CIGAR=1X;DP=37;DPRA=0;GTI=0;LEN=1;MQM=70;MQMR=70;NS=1;NUMALT=1;PAIRED=1;PAIREDR=1;PAO=0;PQA=0;PQR=0;PRO=0;RUN=1;SAR=0;TYPE=snp;set=cs2-cs4 GT:AOP:PL:QA
:QR:RO . 0/0:4:19:1,0,100:24:385:15 . . .
NODE_9_length_16257_cov_15.910501 14201 . A G 0 . AB=0;ABP=0;AC=0;AF=0.00;AN=2;CIGAR=1X;DP=33;DPRA=0;GTI=0;LEN=1;MEANALT=1;MQM=70;MQMR=70;NS=1;NUMALT=1;PAIRED=1;PAIREDR=1;PAO=0;PQA=0;PQR=0;PRO=0;RUN=1;SAR=0;TYPE=snp;set=cs2-cs3
GT:AOP:PL:QA:QR:RO . 0/0:4:17:5,0,100:24:385:13 . . .
NODE_9_length_16257_cov_15.910501 14240 . T G 0 . AB=0;ABP=0;AC=0;AF=0.00;AN=2;AO=4;CIGAR=1X;DP=36;DPB=18;DPRA=0;EPP=11.6962;EPPR=3.63072;GTI=0;LEN=1;MEANALT=1;MQM=70;MQMR=70;NS=1;NUMALT=1;PAIRED=1;PAIREDR=1;PAO=0;PQA=0;PQR=0;PRO=0;
RO=14;RPP=11.6962;RUN=1;SAF=4;SAP=11.6962;SAR=0;TYPE=snp;set=cs2-cs4 GT:AOP:PL:QA:QR:RO . 0/0:4:18:3,0,100:24:383:14 . . .
NODE_9_length_16257_cov_15.910501 14261 . T G 0 . AB=0;ABP=0;AC=0;AF=0.00;AN=2;CIGAR=1X;DP=36;DPRA=0;GTI=0;LEN=1;MEANALT=1;MQM=70;MQMR=70;NS=1;NUMALT=1;PAIRED=1;PAIREDR=1;PAO=0;PQA=0;PQR=0;PRO=0;QA=44;RUN=1;SAR=1;TYPE=snp;set=cs1-cs
3 GT:AOP:PL:QA:QR:RO 0/0:6:21:24,0,100:44:299:15 . . . .
NODE_9_length_16257_cov_15.910501 14263 . A G 0.01 . AB=0;ABP=0;AC=0;AF=0.00;AN=2;CIGAR=1X;DP=27;DPRA=0;GTI=0;LEN=1;MEANALT=1;MQM=70;MQMR=70;NS=1;NUMALT=1;PAIRED=1;PAIREDR=1;PAO=0;PQA=0;PQR=0;PRO=0;RUN=1;SAF=1;TYPE=snp;set=cs3-cs4
GT:AOP:PL:QA:QR:RO . . 0/0:5:15:22,0,61:35:78:10 . .
prs321 is offline   Reply With Quote
Old 05-22-2014, 10:09 AM   #12
lethalfang
Member
 
Location: San Francisco, CA

Join Date: Aug 2011
Posts: 91
Default

You may have mostly positions that are only called in one single replicates.

Try this as a test, just to extract only those positions that exist in all 5 replicates, to see if things are working as intended.

Code:
java -Xmx2g -jar GenomeAnalysisTK.jar \
   -R ref.fasta \
   -T CombineVariants \
   -minN 5 \
   --variant rep1.vcf \
   --variant rep2.vcf \
   --variant rep3.vcf \
   --variant rep4.vcf \
   --variant rep5.vcf \
   -o combined.output.vcf
lethalfang is offline   Reply With Quote
Old 05-22-2014, 10:20 AM   #13
prs321
Member
 
Location: US

Join Date: Jun 2013
Posts: 96
Default

Thank you, I will try it out. Appreciate it.
prs321 is offline   Reply With Quote
Reply

Thread Tools

Posting Rules
You may not post new threads
You may not post replies
You may not post attachments
You may not edit your posts

BB code is On
Smilies are On
[IMG] code is On
HTML code is Off




All times are GMT -8. The time now is 05:48 AM.


Powered by vBulletin® Version 3.8.9
Copyright ©2000 - 2020, vBulletin Solutions, Inc.
Single Sign On provided by vBSSO