EricChu:LabNotesMDA/2015-11-10

From ZhangLabWiki
Revision as of 23:57, 10 November 2015 by >Ericchu
Jump to navigation Jump to search

Pileup PGP1#21 Unique Reads

/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

$ /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

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

$ 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)

$ 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

  • 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

Export Common SNV in Chamber 15/18 HLA fragment

$ 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

$ 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

  • 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.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.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)

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

Then compare with pileup in the alternative chamber