EricChu:LabNotesMDA/2015-11-10

From ZhangLabWiki
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

File:VennDiagramHLAP2&CompleteGenomicsbyPosition.jpg