EricChu:LabNotesMDA/2015-11-10
Jump to navigation
Jump to search
Pileup PGP1#21 Unique Reads[edit]
/home/kunzhang/softwares/samtools-0.1.19/samtools view -b -L PGP1_21_FragmentBoundaryChXX.bed PGP1_21_chXX.unique.chr.bam > PGP1_21_chXX.unique.fragment.bam /home/kunzhang/softwares/samtools-0.1.8/samtools pileup -Ac -f /GenomeDB/Homo_sapiens/UCSC/hg19/Sequence/WholeGenomeFasta/genome.fa PGP1_21_chXX.unique.fragment.bam > PGP1_21_chXX.unique.fragment.txt awk -F"\t" 'FILENAME=="PGP1_21_chXX.unique.upper.fragment.txt"{A[$1$2$3$4]=$1$2$3$4} FILENAME=="haplotype.combine.txt"{if(A[$4$5$6$7]){print}}' PGP1_21_chXX.unique.upper.fragment.txt haplotype.combine.txt > PGP1_21_chXX.unique.Phasedupper.fragment.txt perl removeSNPflags.pl PGP1_21_chXX.unique.Phasedupper.fragment.txt PGP1_21_chXX.unique.heteroSNP.fragment.txt perl splithaplotype.pl PGP1_21_FragmentBoundaryChXX.bed PGP1_21_chXX.unique.heteroSNP.fragment.txt PGP1_21_chXX.hap1.txt PGP1_21_chXX.hap2.txt PGP1_21_chXX.hap0.txt /home/kunzhang/softwares/samtools-0.1.19/samtools view -b -L PGP1_21_chXX.hap1.txt PGP1_21_chXX.unique.fragment.bam > PGP1_21_chXX.unique.haplotype1.bam /home/kunzhang/softwares/bedtools-2.20/bin/bedtools bamtobed -i PGP1_21_chXX.unique.haplotype1.bam > PGP1_21_chXX.unique.haplotype1.bed
Convert Pileup to VCF[edit]
$ /home/kunzhang/softwares/samtools-0.1.8/samtools pileup -Ac -f /GenomeDB/Homo_sapiens/UCSC/hg19/Sequence/WholeGenomeFasta/genome.fa PGP1_21_ch15.unique.haplotype2.bam > PGP1_21_ch15.unique.haplotype2.pileup & $ /home/kunzhang/softwares/samtools-0.1.8/samtools pileup -Ac -f /GenomeDB/Homo_sapiens/UCSC/hg19/Sequence/WholeGenomeFasta/genome.fa PGP1_21_ch18.unique.haplotype2.bam > PGP1_21_ch18.unique.haplotype2.pileup &
/home/kunzhang/softwares/samtools-1.2/misc/sam2vcf.pl -s < PGP1_21_ch15.unique.haplotype2.pileup > PGP1_21_ch15.unique.haplotype2.vcf /home/kunzhang/softwares/samtools-1.2/misc/sam2vcf.pl -s < PGP1_21_ch18.unique.haplotype2.pileup > PGP1_21_ch18.unique.haplotype2.vcf
Extract Fragment in HLA[edit]
bgzip PGP1_21_ch18.unique.haplotype2.vcf tabix -p vcf PGP1_21_ch18.unique.haplotype2.vcf.gz tabix PGP1_21_ch18.unique.haplotype2.vcf.gz chr6:30521122-31340586 > PGP1_21_ch18.unique.haplotype2.HLA.vcf bgzip PGP1_21_ch15.unique.haplotype2.vcf tabix -p vcf PGP1_21_ch15.unique.haplotype2.vcf.gz tabix PGP1_21_ch15.unique.haplotype2.vcf.gz chr6:30521122-31340586 > PGP1_21_ch15.unique.haplotype2.HLA.vcf
Remove Zero Quality[edit]
$ perl remove0qualityVCF.pl PGP1_21_ch15.unique.haplotype2.HLA.vcf PGP1_21_ch15.unique.haplotype2.HLA0.vcf $ perl remove0qualityVCF.pl PGP1_21_ch18.unique.haplotype2.HLA.vcf PGP1_21_ch18.unique.haplotype2.HLA0.vcf $ wc -l PGP1_21_ch15.unique.haplotype2.HLA0.vcf 3255 PGP1_21_ch15.unique.haplotype2.HLA0.vcf $ wc -l PGP1_21_ch18.unique.haplotype2.HLA0.vcf 2921 PGP1_21_ch18.unique.haplotype2.HLA0.vcf $ wc -l PGP1_21_ch15.unique.haplotype2.HLA.vcf 3405 PGP1_21_ch15.unique.haplotype2.HLA.vcf $ wc -l PGP1_21_ch18.unique.haplotype2.HLA.vcf 3100 PGP1_21_ch18.unique.haplotype2.HLA.vcf
Remove Low Quality (BQ<30)[edit]
$ perl /home/ericchu/softwares/removelowquality30.pl PGP1_21_ch18.unique.haplotype2.HLA.vcf PGP1_21_ch18.unique.haplotype2.quality.HLA.vcf $ perl /home/ericchu/softwares/removelowquality30.pl PGP1_21_ch15.unique.haplotype2.HLA.vcf PGP1_21_ch15.unique.haplotype2.quality.HLA.vcf $ wc -l PGP1_21_ch15.unique.haplotype2.quality.HLA.vcf 2237 PGP1_21_ch15.unique.haplotype2.quality.HLA.vcf $ wc -l PGP1_21_ch18.unique.haplotype2.quality.HLA.vcf 1870 PGP1_21_ch18.unique.haplotype2.quality.HLA.vcf
- In comparison to using sorted BAM, where 2239 SNV and 1885 SNV were found in chamber 15 and 18 respectively.
VCF Compare[edit]
- Add header.
bgzip PGP1_21_ch15.unique.haplotype2.quality.HLA.vcf bgzip PGP1_21_ch18.unique.haplotype2.quality.HLA.vcf tabix -p vcf PGP1_21_ch18.unique.haplotype2.quality.HLA.vcf.gz tabix -p vcf PGP1_21_ch15.unique.haplotype2.quality.HLA.vcf.gz $ /home/kunzhang/softwares/vcftools_0.1.12b/bin/vcf-compare -d PGP1_21_ch15.unique.haplotype2.quality.HLA.vcf.gz PGP1_21_ch18.unique.haplotype2.quality.HLA.vcf.gz > comparech15ch18uniqueHLAvcfdquality.txt VN 692 PGP1_21_ch18.unique.haplotype2.quality.HLA.vcf.gz (37.0%) VN 1059 PGP1_21_ch15.unique.haplotype2.quality.HLA.vcf.gz (47.3%) VN 1178 PGP1_21_ch15.unique.haplotype2.quality.HLA.vcf.gz (52.7%) PGP1_21_ch18.unique.haplotype2.quality.HLA.vcf.gz (63.0%) AM chr6 31076199 T C,T AM chr6 31119716 T C,T AM chr6 31190931 G C,G AM chr6 31253771 C,T T SN Number of REF matches: 1178 SN Number of ALT matches: 1174 SN Number of REF mismatches: 0 SN Number of ALT mismatches: 4
$ grep 31076199 pgp1f_hg19_vcf_fixedCompleteGenomics_withHeaders_final.HLA.vcf chr6 31076199 rs10947135 A T -1 PASS COL8 GT:GQ:HQ 1/1:-1:-1,-1 $ grep 31119716 pgp1f_hg19_vcf_fixedCompleteGenomics_withHeaders_final.HLA.vcf chr6 31119716 rs9263769,rs114597112 G T -1 PASS COL8 GT:GQ:HQ 0/1:-1:-1,-1 $ grep 31190931 pgp1f_hg19_vcf_fixedCompleteGenomics_withHeaders_final.HLA.vcf chr6 31190931 rs28361007,rs113067399,rs116072838 A G -1 PASS COL8 GT:GQ:HQ 0/1:-1:-1,-1 $ grep 31253771 pgp1f_hg19_vcf_fixedCompleteGenomics_withHeaders_final.HLA.vcf chr6 31253771 rs2844600,rs115502452 G T -1 PASS COL8 GT:GQ:HQ 0/1:-1:-1,-1
Export Common SNV in Chamber 15/18 HLA fragment[edit]
$ awk -F"\t" 'FILENAME=="PGP1_21_ch15.unique.haplotype2.quality.HLA.vcf"{A[$1$2]=$1$2} FILENAME=="PGP1_21_ch18.unique.haplotype2.quality.HLA.vcf"{if(A[$1$2]){print}}' PGP1_21_ch15.unique.haplotype2.quality.HLA.vcf PGP1_21_ch18.unique.haplotype2.quality.HLA.vcf > PGP1_21_ch15ch18.unique.HLA.common.vcf $ wc -l PGP1_21_ch15ch18.unique.HLA.common.vcf 1184 PGP1_21_ch15ch18.unique.HLA.common.vcf (6 line header)
$ perl /home/ericchu/softwares/removeduplicateline.pl PGP1_21_ch15.unique.haplotype2.quality.HLA.vcf PGP1_21_ch15ch18.unique.HLA.common.vcf PGP1_21_ch15.unique.HLA.mismatch18.vcf $ wc -l PGP1_21_ch15.unique.HLA.mismatch18.vcf 1059 PGP1_21_ch15.unique.HLA.mismatch18.vcf $ perl /home/ericchu/softwares/removeduplicateline.pl PGP1_21_ch18.unique.haplotype2.quality.HLA.vcf PGP1_21_ch15ch18.unique.HLA.common.vcf PGP1_21_ch18.unique.HLA.mismatch15.vcf $ wc -l PGP1_21_ch18.unique.HLA.mismatch15.vcf 692 PGP1_21_ch18.unique.HLA.mismatch15.vcf
Compare Common SNV in Chamber 15/18 and Complete Genomics Data[edit]
$ awk -F"\t" 'FILENAME=="PGP1_21_ch15.unique.haplotype2.quality.HLA.vcf"{A[$1$2$4$5]=$1$2$4$5} FILENAME=="pgp1f_hg19_vcf_fixedCompleteGenomics_withHeaders_final.vcf"{if(A[$1$2$4$5]){print}}' PGP1_21_ch15.unique.haplotype2.quality.HLA.vcf pgp1f_hg19_vcf_fixedCompleteGenomics_withHeaders_final.vcf > PGP1_21_ch15CompleteGenomics.unique.HLA.common.vcf &
$ awk -F"\t" 'FILENAME=="PGP1_21_ch18.unique.haplotype2.quality.HLA.vcf"{A[$1$2$4$5]=$1$2$4$5} FILENAME=="pgp1f_hg19_vcf_fixedCompleteGenomics_withHeaders_final.vcf"{if(A[$1$2$4$5]){print}}' PGP1_21_ch18.unique.haplotype2.quality.HLA.vcf pgp1f_hg19_vcf_fixedCompleteGenomics_withHeaders_final.vcf > PGP1_21_ch18CompleteGenomics.unique.HLA.common.vcf &
$ wc -l PGP1_21_ch15CompleteGenomics.unique.HLA.common.vcf 1588 PGP1_21_ch15CompleteGenomics.unique.HLA.common.vcf (1 line header) $ wc -l PGP1_21_ch18CompleteGenomics.unique.HLA.common.vcf 1368 PGP1_21_ch18CompleteGenomics.unique.HLA.common.vcf (1 line header)
$ awk -F"\t" 'FILENAME=="PGP1_21_ch15CompleteGenomics.unique.HLA.common.vcf"{A[$1$2$4$5]=$1$2$4$5} FILENAME=="PGP1_21_ch18CompleteGenomics.unique.HLA.common.vcf"{if(A[$1$2$4$5]){print}}' PGP1_21_ch15CompleteGenomics.unique.HLA.common.vcf PGP1_21_ch18CompleteGenomics.unique.HLA.common.vcf > PGP1_21_ch15ch18CompleteGenomics.unique.HLA.common.vcf $ wc -l PGP1_21_ch15ch18CompleteGenomics.unique.HLA.common.vcf 1026 PGP1_21_ch15ch18CompleteGenomics.unique.HLA.common.vcf (1 line header)
Export Mismatch from Chamber 15 to chamber 18 and CGI data[edit]
- remove SNV that are found in CGI and report mismatch to both alternative chamber and CGI
$ perl /home/ericchu/softwares/removeduplicateline.pl PGP1_21_ch15.unique.HLA.mismatch18.vcf pgp1f_hg19_vcf_fixedCompleteGenomics_withHeaders_final.HLA.vcf PGP1_21_ch15.unique.HLA.mismatch18nCGI.vcf & $ wc -l PGP1_21_ch15.unique.HLA.mismatch18nCGI.vcf 499 PGP1_21_ch15.unique.HLA.mismatch18nCGI.vcf
$ perl /home/ericchu/softwares/removeduplicateline.pl PGP1_21_ch18.unique.HLA.mismatch15.vcf pgp1f_hg19_vcf_fixedCompleteGenomics_withHeaders_final.HLA.vcf PGP1_21_ch18.unique.HLA.mismatch15nCGI.vcf & $ wc -l PGP1_21_ch18.unique.HLA.mismatch15nCGI.vcf 350 PGP1_21_ch18.unique.HLA.mismatch15nCGI.vcf
- Extracting those mismatch SNV that were covered by the other chamber
$ perl /home/ericchu/softwares/countVCFinbed2.pl PGP1_21_ch18.sorted.haplotype2.HLA.bed PGP1_21_ch15.unique.HLA.mismatch18nCGI.vcf PGP1_21_ch15vcfuniqueinch18bed.txt & $ awk '!seen[$0]++' PGP1_21_ch15vcfuniqueinch18bed.txt > PGP1_21_ch15vcfuniqueinch18bedall.txt $ wc -l PGP1_21_ch15vcfuniqueinch18bedall.txt 251 PGP1_21_ch15vcfuniqueinch18bedall.txt $ perl /home/ericchu/softwares/countVCFinbed2.pl PGP1_21_ch15.sorted.haplotype2.HLA.bed PGP1_21_ch18.unique.HLA.mismatch15nCGI.vcf PGP1_21_ch18vcfuniqueinch15bed.txt & $ awk '!seen[$0]++' PGP1_21_ch18vcfuniqueinch15bed.txt > PGP1_21_ch18vcfuniqueinch15bedall.txt $ wc -l PGP1_21_ch18vcfuniqueinch15bedall.txt 216 PGP1_21_ch18vcfuniqueinch15bedall.txt
/home/kunzhang/softwares/samtools-0.1.19/samtools mpileup -BQ0 -d1000000 -f /GenomeDB/Homo_sapiens/UCSC/hg19/Sequence/WholeGenomeFasta/genome.fa PGP1_21_ch15.unique.fragment.bam -l PGP1_21_ch15vcfuniqueinch18bed.txt > PileupCh15HLAmismatchCh18.txt /home/kunzhang/softwares/samtools-0.1.19/samtools mpileup -BQ0 -d1000000 -f /GenomeDB/Homo_sapiens/UCSC/hg19/Sequence/WholeGenomeFasta/genome.fa PGP1_21_ch18.unique.fragment.bam -l PGP1_21_ch15vcfuniqueinch18bed.txt > PileupCh15HLAmismatchCh18complement.txt
perl /home/ericchu/softwares/countbase.pl PileupCh15HLAmismatchCh18.txt AlleleCountCh15HLAmismatchCh18.txt perl /home/ericchu/softwares/countbase.pl PileupCh15HLAmismatchCh18complement.txt AlleleCountCh15HLAmismatchCh18complement.txt
- Example first 20 lines in PileupCh15HLAmismatchCh18.txt, where these chamber 15 SNV were found called reference in chamber 18.
chr6 30543198 A 1 ^]G B chr6 30543443 A 1 g I chr6 30552647 g 1 c$ B chr6 30555378 c 1 ^]G B chr6 30555477 t 1 G B chr6 30571649 A 4 .Gg, FFFI chr6 30590696 C 53 t$,.,t.T.,.,t.tT.,.ttTTT.T,,t,,t,..,...,,,.,.,,,.,,.., BBBFFBF0FFF<FFFFBBFF7FFFFFIIFFFFBBF<IBFFFIhFBBBFFBFFB chr6 30591764 c 45 .$,,,.,,.,T,,,..,.,t,,..,T..,,,,.T..tt,tt,t.T. FBIFFFIBIFIFF<IFFIIIFIIFBFIFBFFIBBFBB<BFFBFFB chr6 30611816 A 4 G.g. BFFI chr6 30611991 A 1 g B chr6 30612004 C 1 t B chr6 30612008 C 1 t F chr6 30612929 C 6 T.t.,. FFIFBF chr6 30614744 T 2 Cc FF chr6 30614938 G 1 A W chr6 30615658 A 4 T.,t FBBB chr6 30616279 a 1 G B chr6 30616288 t 1 G B chr6 30616292 c 1 T F chr6 30617172 G 1 a B
- Example first 20 lines in PileupCh15HLAmismatchCh18complement.txt, where the reference was called in chamber 18. We can see that the confidence of these reference calls were much higher. Therefore, we only look at those AlleleCount in chamber 18 (AlleleCountCh15HLAmismatchCh18complement.txt) that were not high in confidence (where alternate ALT was called.)
chr6 30543198 A 9 ..,,,,... 0BIIFFIFI chr6 30543443 A 14 ,.....,,...,., BFF<FIIBFBIFFF chr6 30552647 g 6 ..,,.^]. B<IFFB chr6 30555378 c 13 ,$..,,,,..,... BFFIIIFIIBFIF chr6 30555477 t 14 .,.,.,,.,,,,,^]. <FBFFIFFFFFBBB chr6 30571649 A 1 . F chr6 30590696 C 20 ..,,..,.,.,...,.,,^].^]. FFBBBFFIFIFFfIF0FBB< chr6 30591764 c 33 ..,..,,.,,,,,.,..,,,.,.,.,.,.,,,. FFFFFIIFIFFIIFFFFFIIIBFBFBI0IBFfB chr6 30611816 A 56 .$......,,,.....,,,,.,.C.,,,,...,,..,,,.,,,..,,..,..,,,,, FBBBFFBIII<BIFFFBFIIII'IIBIFIIFFFIIFBFIFFFIIFBIIBIFBFBFF chr6 30611991 A 35 ,,,...,,.,,,...,.....,,,,,....,..., BBFBFFIIFIFIFBIIFI7IIFFFF'IIIIB<FFf chr6 30612004 C 30 .,,.,,,...,.....,,,,,....,..., FFBFIFIBBFIFFBIIBFBF'IIIF<0BIf chr6 30612008 C 30 ,$,.,,,...,.....,,,,,....,...,, BFFFBIFFFIFF<FFIFFI<IIIFBFFIfB chr6 30612929 C 34 g,.....,..,....,.,..,,..........., BFBFFFIFIFIFIFIII0III<IFIIIIIIFIF7 chr6 30614744 T 10 ,..,,,,,,, <BFB<BFBF< chr6 30614938 G 2 ,. FI chr6 30615658 A 17 ,,...,.,,.....,,^], <FBF<FBFB<IBFF<B< chr6 30616279 a 26 .....,,.,,....,,,,.,,.,.,, 0FF<BIIFIIFIIIBBF0F<BF<BB' chr6 30616288 t 31 ....,,.,,....,,,,.,,.,.,,.,,^\.^],^], FFBFIFBIFFFIFB<FBFBBIBFB<F'0B0< chr6 30616292 c 31 .$..,,.,,....,,,,.,,.,.,,.,,.,,^]. F<<FFFIFFFIFFBFFIBFIBIB'FFB0B chr6 30617172 G 36 .$.$.....,,,,.,..,,,,,,,,..,.,..,..,^].^], B0FBBBFIIIFFFFIIFIFIIFFII0FBIIB<FBB< chr6 30619551 a 181 ,,,,..,...,,,.,,,,,...,,,,.......,.,,,,..$.,,,,,,,.,,...,,,..,,..,,,,,,.....,...,,,..,.....,,,.,,,,,....,,...,,.,...,,,...,,,..,,,,..,,..,,...,,,,.,,,.,,,...g..,,,,,.,,g.,,,,,,,,..^], B0BF<FFBFBFFFBFFFFFFBBFF<FBBF<F0FIFIIII<F<IIIFIII<FF<FFIIIIIIIIIIBFIIBIFFIFIFIFIFF0FBIIIIIBFFIIIIFfIII_IFIIFIIIf<IIIIFBIFIFIIIIBIIIFIFFFIIIIIFFeFFFBFFFFFIIF<FFBBBBFFFF0FFFFFBBFBFBf chr6 30621413 g 66 aa...........,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,...,,.,,,,,...,,^].^], FFIIBFFIIIIFhIIIIBFIFIFIIBF<IFFIFIFIIBIIFFFIIIIFIIIFFFBFfBBFFFFfBB
- These are the ones that are low confidence in Chamber 18 (AlleleCountCh15HLAmismatchCh18complement.txt). The criteria of confidence here was ALT base frequency > REF base.
chr6 31106499 G A Depth = 1 A1, C0, G0, T0 chr6 31239752 C T Depth = 1 A0, C0, G0, T1 chr6 31240473 T C Depth = 3 A0, C2, G0, T1 chr6 31240489 A G Depth = 3 A1, C0, G2, T0 chr6 31240490 G C Depth = 3 A0, C2, G1, T0 chr6 31240496 G C Depth = 3 A0, C2, G1, T0 chr6 31240497 G A Depth = 3 A2, C0, G1, T0 chr6 31240498 T C Depth = 3 A0, C2, G0, T1 chr6 31240592 T C Depth = 4 A0, C3, G0, T1 chr6 31240593 G A Depth = 4 A3, C0, G1, T0 chr6 31240759 A G Depth = 4 A1, C0, G3, T0 chr6 31240789 G A Depth = 7 A6, C0, G1, T0 chr6 31240854 C A Depth = 3 A2, C1, G0, T0 chr6 31240864 C T Depth = 3 A0, C1, G0, T2 chr6 31241260 T C Depth = 2 A0, C2, G0, T0 chr6 31242174 G A Depth = 1 A1, C0, G0, T0 chr6 31323766 T C Depth = 1 A0, C1, G0, T0 chr6 31325526 G A Depth = 2 A2, C0, G0, T0 chr6 31325880 G A Depth = 2 A2, C0, G0, T0 chr6 31325881 A G Depth = 2 A0, C0, G2, T0
- The first SNV was dropped in chamber 18 and called a reference in concensus because of low (1) base quality. But we can actually see that both chambers were calling the same SNV.
In PileupCh15HLAmismatchCh18.txt, chr6 31106499 G 11 AAAAaAAAaAa B'BFBIFIBB7 In PileupCh15HLAmismatchCh18complement.txt, chr6 31106499 G 1 a 7
- The second SNV was also dropped in chamber 18 for the same reason
chr6 31239752 C 72 TttTTTtTTTttttTttTTTTTtttTTTTtttttTTTttttTtttttttttttttttTtT.TTTttTTtttT 7FFF7<B777FBFF'FI<7FB7BII7BFBIFIIFB7BFFFFFFFBBBBBBBBFBF<FF<IBIBIBBFFB'BB chr6 31239752 C 1 T 7
- The next 6 mismatch have a pattern. This is definitely a mapping issue! Two very similar reads were mapped to same position.
chamber 15 chr6 31240473 T 54 ....,...,...,,,,..,,.C,,,.cc,cCccc,c.CCCCCCCCC,C,,CC.C <B<<BB0<IF7<FIBFBBFI<FIIIFFFIFIFFFFFIBIIFFFFFFdFFFF<<B chr6 31240489 A 56 .,...,,,,..,,.G,,,.gg,gGggg,g.GGGGGGGGG,G,,GG.Ggg,,...^].^], FBFFFBFFIIIFIIFIIFIIIFIFBIIIIIBFIFIFIFBfFFFIFIIFFFFFFFIF chr6 31240490 G 56 .$,...,,,,..,,.C,,,.tc,cCccc,c.CCCCCCCCC,C,,CC.Ccc,,...., FBFFF<F<FIIFFIFIIFIIIFFF<IIFII7IIFIFIF<dIFFFBIIFFFFFFBIF chr6 31240496 G 58 ...,,,,..,,.C,,,.cc,cCccc,c.CCCCCCCCC,C,,CC.Ccc,,....,,,,^\, FBFFIIIFIIIIIIIIIIFFFIIIFIII'IBFIIIIIfIFIIIIIFFBBII0IFFFFB chr6 31240497 G 58 ...,,,,..,,.A,,,.aa,aAaaa,a.AAAAAAAAA,A,,AA.Aaa,,....,,,,, FFFBIFIFIIIIFIIFFIFIFIIIFFII0IFFIIIFFdIIIIFIIFF<<IIBFFFFBB chr6 31240498 T 58 ...,,,,..,,.C,,,.cc,cCccc,c.CCCCCCCCCcC,,CC.Ccc,,....,,,,, B<BBFFI<<IIFIIIFFIIBIIIIIIIF'IFIIIIIIhIFIIFFIIFB'BF<FFFFBB
chamber 18 chr6 31240473 T 3 .cc FBF chr6 31240489 A 3 .gg FFF chr6 31240490 G 3 .cc BBF chr6 31240496 G 3 .cc B<F chr6 31240497 G 3 .aa BBF chr6 31240498 T 3 .cc BFI
- Next two. Again mapping issue.
chamber 15 chr6 31240592 T 68 ,,,,.c$c$c$.C.CCc,,c,,CC.cccc,C,cccc,cc.c,.CC.C.,...,.,,..,..,,.,,...., BBBFFBBB<FF<FBIIIIIFFFFIF7III<FIIFIIFIIIBFIIFIIFBIIIIFIBFFBFIFBFBBBB chr6 31240593 G 66 ,$,,,..A.AAa,,a,,AA.aaaa,A,aaaa,aa.a,.$AA.A.,...,.,,..,..,,.,,....,^]. BBBBFBFFBFIIFFFIFFFFIIFIIIIFFIBIIFIIIFFIIFIIIFIIIIIIBFFFFIFFFFBBFB chamber 18 chr6 31240592 T 4 ccc, IFIF chr6 31240593 G 4 aaa, IFIF
- Then now, the others with higher confidence reference call in Chamber 18 may also have the same mapping issue.
/home/kunzhang/softwares/samtools-0.1.19/samtools mpileup -BQ0 -d1000000 -f /GenomeDB/Homo_sapiens/UCSC/hg19/Sequence/WholeGenomeFasta/genome.fa PGP1_21_ch18.unique.fragment.bam -l PGP1_21_ch18vcfuniqueinch15bed.txt > PileupCh18HLAmismatchCh15.txt & /home/kunzhang/softwares/samtools-0.1.19/samtools mpileup -BQ0 -d1000000 -f /GenomeDB/Homo_sapiens/UCSC/hg19/Sequence/WholeGenomeFasta/genome.fa PGP1_21_ch15.unique.fragment.bam -l PGP1_21_ch18vcfuniqueinch15bed.txt > PileupCh18HLAmismatchCh15complement.txt
perl /home/ericchu/softwares/countbase.pl PileupCh18HLAmismatchCh15.txt AlleleCountCh18HLAmismatchCh15.txt perl /home/ericchu/softwares/countbase.pl PileupCh18HLAmismatchCh15complement.txt AlleleCountCh18HLAmismatchCh15complement.txt
Redo From Removing Zero Quality and Kept Low Quality[edit]
$ bgzip PGP1_21_ch15.unique.haplotype2.HLA0.vcf $ bgzip PGP1_21_ch18.unique.haplotype2.HLA0.vcf $ tabix -p vcf PGP1_21_ch15.unique.haplotype2.HLA0.vcf.gz $ tabix -p vcf PGP1_21_ch18.unique.haplotype2.HLA0.vcf.gz $ /home/kunzhang/softwares/vcftools_0.1.12b/bin/vcf-compare -d PGP1_21_ch15.unique.haplotype2.HLA0.vcf.gz PGP1_21_ch18.unique.haplotype2.HLA0.vcf.gz > comparech15ch18uniqueHLAvcfd0.txt VN 1304 PGP1_21_ch15.unique.haplotype2.HLA0.vcf.gz (40.1%) PGP1_21_ch18.unique.haplotype2.HLA0.vcf.gz (44.6%) VN 1617 PGP1_21_ch18.unique.haplotype2.HLA0.vcf.gz (55.4%) VN 1951 PGP1_21_ch15.unique.haplotype2.HLA0.vcf.gz (59.9%) SN Number of REF matches: 1304 SN Number of ALT matches: 1297 SN Number of REF mismatches: 0 SN Number of ALT mismatches: 7
AM chr6 31076199 T C,T AM chr6 31119716 T C,T AM chr6 31149314 T G AM chr6 31190931 G C,G AM chr6 31191665 T G AM chr6 31253771 C,T T AM chr6 31314044 C A
- Problem with low quality base HERE. Three poor quality base SNV caused this mismatch problem.
Export Common SNV in Chamber 15/18 HLA fragment[edit]
$ awk -F"\t" 'FILENAME=="PGP1_21_ch15.unique.haplotype2.HLA0.vcf"{A[$1$2$4$5]=$1$2$4$5} FILENAME=="PGP1_21_ch18.unique.haplotype2.HLA0.vcf"{if(A[$1$2$4$5]){print}}' PGP1_21_ch15.unique.haplotype2.HLA0.vcf PGP1_21_ch18.unique.haplotype2.HLA0.vcf > PGP1_21_ch15ch18.unique.HLA0.common.vcf $ wc -l PGP1_21_ch15ch18.unique.HLA0.common.vcf 1303 PGP1_21_ch15ch18.unique.HLA0.common.vcf (6 line header)
$ perl /home/ericchu/softwares/removeduplicateline.pl PGP1_21_ch15.unique.haplotype2.HLA0.vcf PGP1_21_ch15ch18.unique.HLA0.common.vcf PGP1_21_ch15.unique.HLA0.mismatch18.vcf $ wc -l PGP1_21_ch15.unique.HLA0.mismatch18.vcf 1958 PGP1_21_ch15.unique.HLA0.mismatch18.vcf
$ perl /home/ericchu/softwares/removeduplicateline.pl PGP1_21_ch18.unique.haplotype2.HLA0.vcf PGP1_21_ch15ch18.unique.HLA0.common.vcf PGP1_21_ch18.unique.HLA0.mismatch15.vcf $ wc -l PGP1_21_ch18.unique.HLA0.mismatch15.vcf 1624 PGP1_21_ch18.unique.HLA0.mismatch15.vcf
Compare Common SNV in Chamber 15/18 and Complete Genomics Data[edit]
$ awk -F"\t" 'FILENAME=="PGP1_21_ch15.unique.haplotype2.HLA0.vcf"{A[$1$2$4$5]=$1$2$4$5} FILENAME=="pgp1f_hg19_vcf_fixedCompleteGenomics_withHeaders_final.vcf"{if(A[$1$2$4$5]){print}}' PGP1_21_ch15.unique.haplotype2.HLA0.vcf pgp1f_hg19_vcf_fixedCompleteGenomics_withHeaders_final.vcf > PGP1_21_ch15CompleteGenomics.unique.HLA0.common.vcf & $ wc -l PGP1_21_ch15CompleteGenomics.unique.HLA0.common.vcf 1648 PGP1_21_ch15CompleteGenomics.unique.HLA0.common.vcf
$ awk -F"\t" 'FILENAME=="PGP1_21_ch18.unique.haplotype2.HLA0.vcf"{A[$1$2$4$5]=$1$2$4$5} FILENAME=="pgp1f_hg19_vcf_fixedCompleteGenomics_withHeaders_final.vcf"{if(A[$1$2$4$5]){print}}' PGP1_21_ch18.unique.haplotype2.HLA0.vcf pgp1f_hg19_vcf_fixedCompleteGenomics_withHeaders_final.vcf > PGP1_21_ch18CompleteGenomics.unique.HLA0.common.vcf & $ wc -l PGP1_21_ch18CompleteGenomics.unique.HLA0.common.vcf 1428 PGP1_21_ch18CompleteGenomics.unique.HLA0.common.vcf
$ awk -F"\t" 'FILENAME=="PGP1_21_ch15CompleteGenomics.unique.HLA0.common.vcf"{A[$1$2$4$5]=$1$2$4$5} FILENAME=="PGP1_21_ch18CompleteGenomics.unique.HLA0.common.vcf"{if(A[$1$2$4$5]){print}}' PGP1_21_ch15CompleteGenomics.unique.HLA0.common.vcf PGP1_21_ch18CompleteGenomics.unique.HLA0.common.vcf > PGP1_21_ch15ch18CompleteGenomics.unique.HLA0.common.vcf $ wc -l PGP1_21_ch15ch18CompleteGenomics.unique.HLA0.common.vcf 1091 PGP1_21_ch15ch18CompleteGenomics.unique.HLA0.common.vcf
Export Mismatch from Chamber 15 to chamber 18 and CGI data[edit]
- remove SNV that are found in CGI and report mismatch to both alternative chamber and CGI
$ perl /home/ericchu/softwares/removeduplicateline.pl PGP1_21_ch15.unique.HLA0.mismatch18.vcf pgp1f_hg19_vcf_fixedCompleteGenomics_withHeaders_final.HLA.vcf PGP1_21_ch15.unique.HLA0.mismatch18nCGI.vcf & $ wc -l PGP1_21_ch15.unique.HLA0.mismatch18nCGI.vcf 1399 PGP1_21_ch15.unique.HLA0.mismatch18nCGI.vcf $ perl /home/ericchu/softwares/removeduplicateline.pl PGP1_21_ch18.unique.HLA0.mismatch15.vcf pgp1f_hg19_vcf_fixedCompleteGenomics_withHeaders_final.HLA.vcf PGP1_21_ch18.unique.HLA0.mismatch15nCGI.vcf & $ wc -l PGP1_21_ch18.unique.HLA0.mismatch15nCGI.vcf 1281 PGP1_21_ch18.unique.HLA0.mismatch15nCGI.vcf
- Extracting those mismatch SNV that were covered by the other chamber
$ perl /home/ericchu/softwares/countVCFinbed2.pl PGP1_21_ch18.sorted.haplotype2.HLA.bed PGP1_21_ch15.unique.HLA0.mismatch18nCGI.vcf PGP1_21_ch15vcfunique0inch18bed.txt & $ awk '!seen[$0]++' PGP1_21_ch15vcfunique0inch18bed.txt > PGP1_21_ch15vcfunique0inch18bedall.txt $ wc -l PGP1_21_ch15vcfunique0inch18bedall.txt 877 PGP1_21_ch15vcfunique0inch18bedall.txt
$ perl /home/ericchu/softwares/countVCFinbed2.pl PGP1_21_ch15.sorted.haplotype2.HLA.bed PGP1_21_ch18.unique.HLA0.mismatch15nCGI.vcf PGP1_21_ch18vcfunique0inch15bed.txt & $ awk '!seen[$0]++' PGP1_21_ch18vcfunique0inch15bed.txt > PGP1_21_ch18vcfunique0inch15bedall.txt $ wc -l PGP1_21_ch18vcfunique0inch15bedall.txt 881 PGP1_21_ch18vcfunique0inch15bedall.txt
/home/kunzhang/softwares/samtools-0.1.19/samtools mpileup -BQ0 -d1000000 -f /GenomeDB/Homo_sapiens/UCSC/hg19/Sequence/WholeGenomeFasta/genome.fa PGP1_21_ch15.unique.fragment.bam -l PGP1_21_ch15vcfunique0inch18bed.txt > PileupCh15HLA0mismatchCh18.txt $ wc -l PileupCh15HLA0mismatchCh18.txt 877 PileupCh15HLA0mismatchCh18.txt
/home/kunzhang/softwares/samtools-0.1.19/samtools mpileup -BQ0 -d1000000 -f /GenomeDB/Homo_sapiens/UCSC/hg19/Sequence/WholeGenomeFasta/genome.fa PGP1_21_ch18.unique.fragment.bam -l PGP1_21_ch15vcfunique0inch18bed.txt > PileupCh15HLA0mismatchCh18complement.txt $ wc -l PileupCh15HLA0mismatchCh18complement.txt 868 PileupCh15HLA0mismatchCh18complement.txt
$ /home/kunzhang/softwares/samtools-0.1.19/samtools mpileup -BQ0 -d1000000 -f /GenomeDB/Homo_sapiens/UCSC/hg19/Sequence/WholeGenomeFasta/genome.fa PGP1_21_ch18.unique.fragment.bam -l PGP1_21_ch18vcfunique0inch15bed.txt > PileupCh18HLA0mismatchCh15.txt $ wc -l PileupCh18HLA0mismatchCh15.txt 881 PileupCh18HLA0mismatchCh15.txt
$ /home/kunzhang/softwares/samtools-0.1.19/samtools mpileup -BQ0 -d1000000 -f /GenomeDB/Homo_sapiens/UCSC/hg19/Sequence/WholeGenomeFasta/genome.fa PGP1_21_ch15.unique.fragment.bam -l PGP1_21_ch18vcfunique0inch15bed.txt > PileupCh18HLA0mismatchCh15complement.txt $ wc -l PileupCh18HLA0mismatchCh15complement.txt 874 PileupCh18HLA0mismatchCh15complement.txt
- Example first 20 lines in PileupCh15HLA0mismatchCh18.txt, where these chamber 15 SNV were found called reference in chamber 18.
chr6 30528469 a 1 T F chr6 30528475 a 1 T F chr6 30536395 t 1 ^!c 7 chr6 30536413 t 1 a B chr6 30536414 c 1 t < chr6 30536415 a 1 g < chr6 30536416 c 1 t B chr6 30536421 c 2 g^!G BB chr6 30536444 c 2 tT FF chr6 30536453 c 2 aA <B chr6 30537007 c 1 a ^ chr6 30543198 A 1 ^]G B chr6 30543443 A 1 g I chr6 30552647 g 1 c$ B chr6 30555378 c 1 ^]G B chr6 30555438 c 1 A ' chr6 30555477 t 1 G B chr6 30555740 t 4 .,.C <IBI chr6 30560986 g 1 t F chr6 30561010 t 1 a F
- Example first 20 lines in PileupCh15HLAmismatchCh18complement.txt, where the reference was called in chamber 18. We can see that the confidence of these reference calls were much higher.
chr6 30528469 a 15 ,$..,,,,,.,*,... BB7FIFFFFIIFIFB chr6 30528475 a 16 ..,,,,,.,,,...,^(g F0BFFFIFBB<IFFB< chr6 30536395 t 67 ,$.$,,..,,.,...,,,,,,.,.,.........,$,.....,..,,,,,,,.,.,,..,,,g,.,,... BFFFFFFFFF<BFIIFIIB<FFIIIIIIIIIIhIFIIIIWFIFFBFBFBF0IFcFhFFB'F0BBFFB chr6 30536413 t 66 .$.,,,,,,.,.,.........,.......,,,,,,,.,.,,..,,, ,,.,,.....,,,,,.,,., 7BFF<IBI<FBFFBFFFFFFFIIIIFIFIIFIF0F<F'IFYFcBFB<BBBBIFIIFFB<B'FB<FF chr6 30536414 c 68 .$,,,,,,.,.,.........,.......,,,,,,,.,.,,..,,,,,.,,.....,,,,,.,,.,^;.^;,^9, FFFBFFI'BBFFBBFBFFFFIIIFBIIIIFIF7FFI'IFWBdBF7BB<FFIFIIIFBBB7FBFFBBB0 chr6 30536415 a 67 ,,,,,,.,.,.........,.......,,,,,,,.,.,g..,,,,,.,,.....,,,,,.,,.,.,, FF0FFI'<FFFBBBFFFFFIBIFBIIII7FI'BII0FFKFd<B7BF'BFFFIIFFB<<'F<BFBB7' chr6 30536416 c 69 ,,,,,,.t.,.........,.......,,,,,,,.,.,,..,,,,,.,,.....,,,,,.,,.,.,,^:,^:, FF<FFF'<BIF0BBBFFFFIBFBFIIII<II0FII0FFWFcFFB<B<BFBBIIFFFBF77BF<FBB7'B chr6 30536421 c 71 ,,,,.g.,.........,.......,,,,,,,.,.,,..,,,,,.,,.....,,,,,.,,.,.,,,,,... <FFF0FBFFFBFFFFBBIFFFFIIIIIII0IFF7IFYFdBFF<FBFFIFFIBFFBB<IBFFFFB<0BFFFF chr6 30536444 c 62 ,.......,,,,,,,.,.,.,,,,,.,,.....,,,,,.,,.,.,,,,,...,.,,,.,,., FFFFFFFFIBIIFIIFBBIFIBF'F7FIIFIIFBIIF7IFFFFIF<BFFIIFFF<BFF<BFF chr6 30536453 c 59 ....,,,,,,,.,.,.,,,,,.,,.....,,,,,.,,.,.,,,,,...,.,,,.,,.,, 77B<FFFFBFI0F<F7BFFFF'FIB'<B<FIIFBBFIIIIB<<FFFFBF'BFFIFFFBB chr6 30537007 c 20 ....,.,.,...,,,,,,,, BFFFBIFFFFIIFFFBF<FF chr6 30543198 A 9 ..,,,,... 0BIIFFIFI chr6 30543443 A 14 ,.....,,...,., BFF<FIIBFBIFFF chr6 30552647 g 6 ..,,.^]. B<IFFB chr6 30555378 c 13 ,$..,,,,..,... BFFIIIFIIBFIF chr6 30555438 c 13 .,....,.,.,,. <I0'<'F<FBFFF chr6 30555477 t 14 .,.,.,,.,,,,,^]. <FBFFIFFFFFBBB chr6 30555740 t 16 ,C.,,...,,,,..,. F'0FFBBFBBBFBB7< chr6 30560986 g 2 ,, BF chr6 30561010 t 2 ,,
Resolved the Majority of SNV Calling Error[edit]
- These are the ones that are low confidence in Chamber 18 (AlleleCountCh15HLA0mismatchCh18complement.txt). The high quality ones were removed from list. The criteria of confidence here was ALT base frequency > REF base. It is much better than previously removing low quality base, which both chamber were called the same ALT. Some issue like the pattern misalignment up there were only masked but not removed.
$ perl /home/ericchu/softwares/countbase.pl PileupCh15HLA0mismatchCh18complement.txt AlleleCountCh15HLA0mismatchCh18complement.txt $ more AlleleCountCh15HLA0mismatchCh18complement.txt chr6 31240473 T C Depth = 3 A0, C2, G0, T1 chr6 31240489 A G Depth = 3 A1, C0, G2, T0 chr6 31240759 A G Depth = 4 A1, C0, G3, T0 chr6 31240854 C A Depth = 3 A2, C1, G0, T0 chr6 31240864 C T Depth = 3 A0, C1, G0, T2
$ perl /home/ericchu/softwares/countbase.pl PileupCh18HLA0mismatchCh15complement.txt AlleleCountCh18HLA0mismatchCh15complement.txt $ more AlleleCountCh18HLA0mismatchCh15complement.txt chr6 31134118 C A Depth = 91 A66, C-107, G132, T0 $ wc -l AlleleCountCh18HLA0mismatchCh15complement.txt 1 AlleleCountCh18HLA0mismatchCh15complement.txt $ grep 31134118 PileupCh18HLA0mismatchCh15complement.txt chr6 31134118 C 91 ..$,,,,$,.....$.,,,.+3GAG.+3GAG.+3GAG.+3GAG.+3GAG,+3gag,.+3GAG.+3GAG,+3gag,+3gag,$.$,+3gag.$.+3GAG.+3GAG..+3GAG,+3gag.+3GAG,+3gag,+3gag,+3gag,+3gag.+3GAG.+3GAG.+3GAG,+3gag.+3GAG.+3GAG.+3GAG,+3gag.+3AAG.+3GAG,+3gag.+3GAG.+3GAG.+3GAG.,+3gag,+3gag,+3gag,+3gag,+3gag,+3gag,+3gag.+3GAG.+3GAG.+3GAG.+3GAG,+3gag.+3GAG.+3GAG,+3gag,+3gag,+3gag,+3gag,+3gag.+3GAG.+3GAG,+3gag,+3gag.+3GAG.+3GAG.+3GAG,+3gag,+3gag,+3gag.+3GAG..+3GAG^]G^]G^]g BB<BB0FFFFFBFFFBFFBBFFFFFFIFFIFBBBBBFIIIIIIFIBFFIFIFFFIfFFFFBFeFIFIFFIFBB'FFI7BIIIF0BBBFBB<
Extract All Unmatched SNV after Ref Correction[edit]
$ perl /home/ericchu/softwares/removeduplicateline.pl PGP1_21_ch15.unique.HLA0.mismatch18nCGI.vcf PileupCh15HLA0mismatchCh18complement.txt PGP1_21_ch15vcfunmatched.txt $ wc -l PGP1_21_ch15vcfunmatched.txt 531 PGP1_21_ch15vcfunmatched.txt (Still 5 uncertain from alignment)
- remove low quality
$ perl /home/ericchu/softwares/removelowquality30.pl PGP1_21_ch15vcfunmatched.txt PGP1_21_ch15vcfunmatched.quality.txt $ wc -l PGP1_21_ch15vcfunmatched.quality.txt 252 PGP1_21_ch15vcfunmatched.quality.txt
$ perl /home/ericchu/softwares/removeduplicateline.pl PGP1_21_ch18.unique.HLA0.mismatch15nCGI.vcf PileupCh18HLA0mismatchCh15complement.txt PGP1_21_ch18vcfunmatched.txt $ wc -l PGP1_21_ch18vcfunmatched.txt 407 PGP1_21_ch18vcfunmatched.txt
- remove low quality
$ perl /home/ericchu/softwares/removelowquality30.pl PGP1_21_ch18vcfunmatched.txt PGP1_21_ch18vcfunmatched.quality.txt $ wc -l PGP1_21_ch18vcfunmatched.quality.txt 135 PGP1_21_ch18vcfunmatched.quality.txt
- However a lot of them have good base quality but very low depth (DP)
chr6 30523430 . C G 30 0 DP=1 GT:GQ:DP 1/1:30:1 chr6 30538637 . T A 75 0 DP=11 GT:GQ:DP 0/1:75:11 chr6 30544307 . C T 114 0 DP=16 GT:GQ:DP 0/1:64:16 chr6 30546029 . A G 30 0 DP=1 GT:GQ:DP 1/1:30:1 chr6 30546679 . A G 30 0 DP=1 GT:GQ:DP 1/1:30:1 chr6 30547425 . G A 30 0 DP=1 GT:GQ:DP 1/1:30:1 chr6 30547443 . C T 30 0 DP=1 GT:GQ:DP 1/1:30:1 chr6 30547523 . A T 33 0 DP=2 GT:GQ:DP 1/1:33:2 chr6 30548506 . A T 30 0 DP=1 GT:GQ:DP 1/1:30:1 chr6 30548776 . A G 30 0 DP=1 GT:GQ:DP 1/1:30:1 chr6 30548786 . T C 30 0 DP=1 GT:GQ:DP 1/1:30:1 chr6 30548795 . C T 30 0 DP=1 GT:GQ:DP 1/1:30:1 chr6 30551806 . A G 34 0 DP=3 GT:GQ:DP 0/1:1:3 chr6 30551807 . A G 34 0 DP=3 GT:GQ:DP 0/1:1:3 chr6 30555369 . G A 63 0 DP=13 GT:GQ:DP 1/1:26:13 chr6 30556424 . C T 35 0 DP=5 GT:GQ:DP 0/1:35:5 chr6 30556814 . C T 33 0 DP=2 GT:GQ:DP 1/1:33:2 ...
- If DP equals or large than 5, then unmatched in chamber 15 and 18 would drop to 149 SNVs and 36 SNVs respectively.
File:VennDiagramHLAP2&CompleteGenomics3.jpg
Compare Common SNV in Chamber 15/18 and Complete Genomics Data (by Position Only)[edit]
$ awk -F"\t" 'FILENAME=="PGP1_21_ch15.unique.haplotype2.HLA0.vcf"{A[$1$2$4]=$1$2$4} FILENAME=="PGP1_21_ch18.unique.haplotype2.HLA0.vcf"{if(A[$1$2$4]){print}}' PGP1_21_ch15.unique.haplotype2.HLA0.vcf PGP1_21_ch18.unique.haplotype2.HLA0.vcf > PGP1_21_ch15ch18.unique.HLA0.commonpos.vcf $ wc -l PGP1_21_ch15ch18.unique.HLA0.commonpos.vcf 1310 PGP1_21_ch15ch18.unique.HLA0.commonpos.vcf
$ awk -F"\t" 'FILENAME=="PGP1_21_ch15.unique.haplotype2.HLA0.vcf"{A[$1$2$4]=$1$2$4} FILENAME=="pgp1f_hg19_vcf_fixedCompleteGenomics_withHeaders_final.vcf"{if(A[$1$2$4]){print}}' PGP1_21_ch15.unique.haplotype2.HLA0.vcf pgp1f_hg19_vcf_fixedCompleteGenomics_withHeaders_final.vcf > PGP1_21_ch15CompleteGenomics.unique.HLA0.commonpos.vcf & $ wc -l PGP1_21_ch15CompleteGenomics.unique.HLA0.commonpos.vcf 1651 PGP1_21_ch15CompleteGenomics.unique.HLA0.commonpos.vcf
$ awk -F"\t" 'FILENAME=="PGP1_21_ch18.unique.haplotype2.HLA0.vcf"{A[$1$2$4]=$1$2$4} FILENAME=="pgp1f_hg19_vcf_fixedCompleteGenomics_withHeaders_final.vcf"{if(A[$1$2$4]){print}}' PGP1_21_ch18.unique.haplotype2.HLA0.vcf pgp1f_hg19_vcf_fixedCompleteGenomics_withHeaders_final.vcf > PGP1_21_ch18CompleteGenomics.unique.HLA0.commonpos.vcf & $ wc -l PGP1_21_ch18CompleteGenomics.unique.HLA0.commonpos.vcf 1435 PGP1_21_ch18CompleteGenomics.unique.HLA0.commonpos.vcf
$ awk -F"\t" 'FILENAME=="PGP1_21_ch15CompleteGenomics.unique.HLA0.commonpos.vcf"{A[$1$2$4]=$1$2$4} FILENAME=="PGP1_21_ch18CompleteGenomics.unique.HLA0.commonpos.vcf"{if(A[$1$2$4]){print}}' PGP1_21_ch15CompleteGenomics.unique.HLA0.commonpos.vcf PGP1_21_ch18CompleteGenomics.unique.HLA0.commonpos.vcf > PGP1_21_ch15ch18CompleteGenomics.unique.HLA0.commonpos.vcf & $ wc -l PGP1_21_ch15ch18CompleteGenomics.unique.HLA0.commonpos.vcf 1097 PGP1_21_ch15ch18CompleteGenomics.unique.HLA0.commonpos.vcf
- Number of False Positives = Difference of SNV Calling by just position and actual SNV Calling