number of snps are very low in vcf file

[Waqasuddin Khan]

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.

[jvolanen]

Are you mixing old and new reference genomes:

REFERENCE=new_newhg38.fa
--known 1000G_phase1.indels.hg19.sites.vcf

[jvolanen]

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.

[Waqasuddin Khan]

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???

[jvolanen]

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