SEQanswers

Go Back   SEQanswers > Bioinformatics > Bioinformatics



Similar Threads
Thread Thread Starter Forum Replies Last Post
how to extract snps from vcf file cli Bioinformatics 0 08-23-2014 08:18 AM
merging neighbouring SNPs in vcf file VeBeKay Bioinformatics 1 05-13-2014 05:58 AM
How and what to use to filter out any known SNPs (with rs#) from the vcf file? BhariD Bioinformatics 0 05-16-2013 02:55 PM
How to sort out cancer-specific snps from vcf file? LI Jia Bioinformatics 0 01-21-2013 05:26 AM
Predicting true SNPs from .vcf file swbarnes2 Bioinformatics 1 04-06-2011 04:29 PM

Reply
 
Thread Tools
Old 08-27-2015, 02:54 AM   #1
Waqasuddin Khan
Junior Member
 
Location: Karachi

Join Date: Nov 2013
Posts: 7
Default number of snps are very low in vcf file

Hi,

Dear All,

I have WGS data of one healthy individual. I want to make a VCF file. I am running the following pipeline but the number of SNPs are low, only in thousands. I have a big bam file (256GB) with depth of 38X.

java -Xmx32g -Djava.io.tmpdir=java_tmp -XX:MaxPermSize=512m -XX:-UseGCOverheadLimit -jar /home/lifescope/Desktop/picard-1.133/dist/picard.jar MarkDuplicates MAX_FILE_HANDLES_FOR_READ_ENDS_MAP=16000 REMOVE_DUPLICATES=true INPUT=gcxall_sort_mapped.bam OUTPUT=gcxall_sort_mapped_dedup.bam METRICS_FILE=gcx.mat ASSUME_SORTED=true

nohup java -Xmx32g -Djava.io.tmpdir=java_tmp -XX:MaxPermSize=512m -XX:-UseGCOverheadLimit -jar /home/lifescope/Desktop/picard-1.133/dist/picard.jar AddOrReplaceReadGroups INPUT=gcxall_sort_mapped_dedup.bam OUTPUT=gcxall_sort_mapped_dedup_crctd.bam SORT_ORDER=coordinate RGLB=MP RGPL=SOLID RGPU=NOBARCODE RGSM=GCXI

java -Xmx32g -Djava.io.tmpdir=java_tmp -XX:MaxPermSize=512m -XX:-UseGCOverheadLimit -jar /home/lifescope/Desktop/picard-1.133/dist/picard.jar CreateSequenceDictionary REFERENCE=new_newhg38.fa OUTPUT=new_newhg38.dict

java -Xmx32g -Djava.io.tmpdir=java_tmp -XX:MaxPermSize=512m -XX:-UseGCOverheadLimit -jar /home/lifescope/Desktop/picard-1.133/dist/picard.jar ReorderSam I=gcxall_sort_mapped_dedup_crctd.bam O=gcxall_sort_mapped_dedup_crctd_rordr.bam REFERENCE=new_newhg38.fa

/home/lifescope/Desktop/usr/java/jdk1.7.0_79/bin/java -Xmx32g -Djava.io.tmpdir=java_tmp -XX:MaxPermSize=512m -XX:-UseGCOverheadLimit -jar /home/lifescope/Desktop/GenomeAnalysisTK.jar -T RealignerTargetCreator -R /scratch/gcxi_all_bam/new_newhg38.fa -I gcxall_sort_mapped_dedup_crctd_rordr.bam --known Mills_and_1000G_gold_standard.indels.hg19.sites.vcf --known 1000G_phase1.indels.hg19.sites.vcf -o gcxall_sort_mapped_dedup_crctd_rordr.intervals --defaultBaseQualities 35

/home/lifescope/Desktop/usr/java/jdk1.7.0_79/bin/java -Xmx32g -Djava.io.tmpdir=java_tmp -XX:MaxPermSize=512m -XX:-UseGCOverheadLimit -jar /home/lifescope/Desktop/GenomeAnalysisTK.jar -T IndelRealigner -R /scratch/gcxi_all_bam/new_newhg38.fa -I gcxall_sort_mapped_dedup_crctd_rordr.bam -targetIntervals gcxall_sort_mapped_dedup_crctd_rordr.intervals -o gcxall_sort_mapped_dedup_crctd_rordr_interval.bam -known Mills_and_1000G_gold_standard.indels.hg19.sites.vcf -known 1000G_phase1.indels.hg19.sites.vcf -l INFO --defaultBaseQualities 35

/home/lifescope/Desktop/usr/java/jdk1.7.0_79/bin/java -Xmx32g -Djava.io.tmpdir=java_tmp -XX:MaxPermSize=512m -XX:-UseGCOverheadLimit -jar /home/lifescope/Desktop/GenomeAnalysisTK.jar -T BaseRecalibrator -R /scratch/gcxi_all_bam/new_newhg38.fa -I gcxall_sort_mapped_dedup_crctd_rordr_interval.bam -knownSites 1000G_phase1.indels.hg19.sites.vcf -knownSites Mills_and_1000G_gold_standard.indels.hg19.sites.vcf -knownSites dbsnp_138.hg19_out.vcf -o gcxi_recal_data.table --covariate QualityScoreCovariate --covariate ReadGroupCovariate --covariate ContextCovariate --covariate CycleCovariate --solid_nocall_strategy PURGE_READ --solid_recal_mode SET_Q_ZERO_BASE_N

/home/lifescope/Desktop/usr/java/jdk1.7.0_79/bin/java -Xmx32g -Djava.io.tmpdir=java_tmp -XX:MaxPermSize=512m -XX:-UseGCOverheadLimit -jar /home/lifescope/Desktop/GenomeAnalysisTK.jar -T PrintReads -R /scratch/gcxi_all_bam/new_newhg38.fa -I gcxall_sort_mapped_dedup_crctd_rordr_interval.bam -BQSR gcxi_recal_data.table -o gcxall_sort_mapped_dedup_crctd_rordr_interval_recal.bam -allowPotentiallyMisencodedQuals

/home/lifescope/Desktop/usr/java/jdk1.7.0_79/bin/java -Xmx32g -Djava.io.tmpdir=java_tmp -XX:MaxPermSize=512m -XX:-UseGCOverheadLimit -jar /home/lifescope/Desktop/GenomeAnalysisTK.jar -T UnifiedGenotyper -glm BOTH -R /scratch/gcxi_all_bam/new_newhg38.fa -I gcxall_sort_mapped_dedup_crctd_rordr_interval_recal.bam -o raw_variants.vcf -D dbsnp_138.hg19_out.vcf -allowPotentiallyMisencodedQuals

Kindly help me in this regard, what I am missing..,,!!!??

Thanks.
Waqasuddin Khan is offline   Reply With Quote
Old 08-27-2015, 08:24 AM   #2
jvolanen
Junior Member
 
Location: Espoo, Finland

Join Date: Apr 2015
Posts: 4
Default

Are you mixing old and new reference genomes:
Quote:
REFERENCE=new_newhg38.fa
--known 1000G_phase1.indels.hg19.sites.vcf
jvolanen is offline   Reply With Quote
Old 08-27-2015, 08:39 AM   #3
jvolanen
Junior Member
 
Location: Espoo, Finland

Join Date: Apr 2015
Posts: 4
Default

Quote:
Originally Posted by Waqasuddin Khan View Post
Hi,
I have a big bam file (256GB) with depth of 38X.
The file size is quite high to be 38X. You can compare your quality metrics to other samples with omnomicsQ (requires registration, free trial) to see if your bam has any other issues.
jvolanen is offline   Reply With Quote
Old 08-27-2015, 11:30 AM   #4
Waqasuddin Khan
Junior Member
 
Location: Karachi

Join Date: Nov 2013
Posts: 7
Default

Hi jvolanen,

Throughout I used only one reference genome file (new_newhg38.fa, thats just my reference genome filename).

Is there GATK pipeline issue, UnifiedGenotyper???
Waqasuddin Khan is offline   Reply With Quote
Old 08-27-2015, 10:46 PM   #5
jvolanen
Junior Member
 
Location: Espoo, Finland

Join Date: Apr 2015
Posts: 4
Default

Quote:
Originally Posted by Waqasuddin Khan View Post
Throughout I used only one reference genome file (new_newhg38.fa, thats just my reference genome filename).

Is there GATK pipeline issue, UnifiedGenotyper???
Based on the file names you are using hg38 reference genome to align but you are using known indels and dbsnp files using hg19 coordinates. That does not work. Either align your bam again against hg19 or change known indel and snp files to hg38 version.

I do not know how compatible GATK pipeline is with hg38. I would use SpeedSeq instead of GATK because it is faster, much easier and the accuracy is about the same.

Jussi
jvolanen 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 07:02 PM.


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