EricChu:LabNotesMDA/2015-10-28

From ZhangLabWiki
Jump to navigation Jump to search

next:

compare VCF to reference filter low quality score base remove homozygous SNV to hopefully remove all shared SNV in P1 and P2.

  • Use intersect on the identified P1 and P2 bam.
  • Try chamber 15 and chamber 18 because it has a P2 overlap at the HLA region. (see slides 2015-10-05)
  • However Chr6 was not found in any of the fragments!
ericchu@genome-miner2:/media/LTS_15T/Eric_LTS/PGP1_21_BWA_MEM/HaplotypeBAC$ /home/kunzhang/softwares/bedtools-2.20/bin/bedtools intersect -a PGP1_21_ch15.bac.hap2.txt -b PGP1_21_ch18.bac.hap2.txt > intersectCh15Ch18hap2.bed
ericchu@genome-miner2:/media/LTS_15T/Eric_LTS/PGP1_21_BWA_MEM/HaplotypeBAC$ ^C
ericchu@genome-miner2:/media/LTS_15T/Eric_LTS/PGP1_21_BWA_MEM/HaplotypeBAC$ more intersectCh15Ch18hap2.bed
chr1    237404608       237460113       161629  chamber15
chr2    26188503        26477616        699901  chamber15
chr3    107905302       108184164       536794  chamber15
chr3    179860303       179913717       108510  chamber15
chr4    101657850       101762603       321532  chamber15
chr4    148586980       148648700       451863  chamber15
chr5    130833747       130886458       109538  chamber15
chr8    104223130       104328780       105650  chamber15
chr10   43876151        44003179        627043  chamber15
chr10   63610223        63770883        743520  chamber15
chr12   43466622        43734216        813903  chamber15
chr15   48096184        48202050        105866  chamber15
chr16   33843974        34220262        1331685 chamber15
chr17   17635371        17849296        383109  chamber15
ericchu@genome-miner2:/media/LTS_15T/Eric_LTS/PGP1_21_BWA_MEM/HaplotypeBAC$ /home/kunzhang/softwares/bedtools-2.20/bin/bedtools intersect -a PGP1_21_ch15.bac.hap1.txt -b PGP1_21_ch18.bac.hap1.txt > intersectCh15Ch18hap1.bed
ericchu@genome-miner2:/media/LTS_15T/Eric_LTS/PGP1_21_BWA_MEM/HaplotypeBAC$ more intersectCh15Ch18hap1.bed
chr1    197761398       197823581       115021  chamber15
chr2    165163762       165373813       210051  chamber15
chr5    31336252        31388732        104143  chamber15
chr9    36063511        36174513        436066  chamber15
chr11   58154650        58324188        272973  chamber15
chr11   100516817       100739773       222956  chamber15
chr12   13382275        13602142        271208  chamber15
ericchu@genome-miner2:/media/LTS_15T/Eric_LTS/PGP1_21_BWA_MEM/HaplotypeBAC$ /home/kunzhang/softwares/bedtools-2.20/bin/bedtools intersect -a PGP1_21_ch15.bac.hap0.txt -b PGP1_21_ch18.bac.hap0.txt > intersectCh15Ch18hap0.bed
ericchu@genome-miner2:/media/LTS_15T/Eric_LTS/PGP1_21_BWA_MEM/HaplotypeBAC$ more intersectCh15Ch18hap0.bed
chr7    122340726       122563605       2336551 chamber15
chr11   37172028        37335873        2992837 chamber15
chr18   10407458        10713105        568646  chamber15
  • Unfortunately one fragment was found in the haplotype and the other one in the undefined.
ericchu@genome-miner2:/media/LTS_15T/Eric_LTS/PGP1_21_BWA_MEM/HaplotypeBAC$ /home/kunzhang/softwares/bedtools-2.20/bin/bedtools intersect -a PGP1_21_ch15.bac.hap0.txt -b PGP1_21_ch18.bac.hap1.txt > intersectCh15hap0Ch18hap1.bed
ericchu@genome-miner2:/media/LTS_15T/Eric_LTS/PGP1_21_BWA_MEM/HaplotypeBAC$ more intersectCh15hap0Ch18hap1.bed
chr2    120685215       120917453       776418  chamber15
chr2    177847392       177899850       814223  chamber15
chr4    55095268        55212283        117015  chamber15
chr5    124392176       124713972       479615  chamber15
chr6    30521122        31340586        1062492 chamber15
chr7    121796946       122075571       2336551 chamber15
chr11   38786387        38952685        2992837 chamber15
chr19   19730163        19788627        440962  chamber15
chr20   3456456 3684368 227912  chamber15
  • Specially generate the pileup for Chamber15Haplotype0(undefined) for the next step.
/media/LTS_15T/Eric_LTS/PGP1_21_BWA_MEM$ /home/kunzhang/softwares/samtools-0.1.19/samtools view -b -L PGP1_21_ch15.hap0.txt PGP1_21_ch15.sorted.fragment.bam > PGP1_21_ch15.sorted.haplotype0.bam &
/media/LTS_15T/Eric_LTS/PGP1_21_BWA_MEM$ /home/kunzhang/softwares/samtools-0.1.8/samtools pileup -Ac -f /GenomeDB/Homo_sapiens/UCSC/hg19/Sequence/WholeGenomeFasta/genome.fa PGP1_21_ch15.sorted.haplotype0.bam > PGP1_21_ch15.haplotype0.pileup &
  • Just look into the pileup (or VCF) file /media/LTS_15T/Eric_LTS/PGP1_21_BWA_MEM
/home/kunzhang/softwares/samtools-1.2/misc/sam2vcf.pl -s < PGP1_21_ch18.haplotype1.pileup > PGP1_21_ch18.haplotype1.vcf
/home/kunzhang/softwares/samtools-1.2/misc/sam2vcf.pl -s < PGP1_21_ch15.haplotype0.pileup > PGP1_21_ch15.haplotype0.vcf
  • I don't know why suddenly the fragment is only found in haplotype 2 in ch15 and ch18.
  • Extract the region chr6: 30521122-31340586
bgzip PGP1_21_ch18.haplotype2.vcf
tabix -p vcf PGP1_21_ch18.haplotype2.vcf.gz
/media/LTS_15T/Eric_LTS/PGP1_21_BWA_MEM$ tabix PGP1_21_ch18.haplotype2.vcf.gz chr6:30521122-31340586 > PGP1_21_ch18.haplotype2.HLA.vcf
or tabix -fB my.vcf.gz reg.bed
bgzip PGP1_21_ch15.haplotype2.vcf
tabix -p vcf PGP1_21_ch15.haplotype2.vcf.gz
/media/LTS_15T/Eric_LTS/PGP1_21_BWA_MEM$ tabix PGP1_21_ch15.haplotype2.vcf.gz chr6:30521122-31340586 > PGP1_21_ch15.haplotype2.HLA.vcf
  • Need to add header back to vcf, then bgzip, then tabix
  • Then we can compare the vcf files
/home/kunzhang/softwares/vcftools_0.1.12b/bin/vcf-compare -d PGP1_21_ch15.haplotype2.HLA.vcf.gz PGP1_21_ch18.haplotype2.HLA.vcf.gz > compareCh15Ch18HLAvcfd.txt
ericchu@genome-miner2:/media/LTS_15T/Eric_LTS/PGP1_21_BWA_MEM$ more compareCh15Ch18HLAvcfd.txt
# This file was generated by vcf-compare.
# The command line was: vcf-compare(r953) PGP1_21_ch15.haplotype2.HLA.vcf.gz PGP1_21_ch18.haplotype2.HLA.vcf.gz
#
#VN 'Venn-Diagram Numbers'. Use `grep ^VN | cut -f 2-` to extract this part.
#VN The columns are:
#VN        1  .. number of sites unique to this particular combination of files
#VN        2- .. combination of files and space-separated number, a fraction of sites in the file
VN      1328    PGP1_21_ch15.haplotype2.HLA.vcf.gz (38.9%)      PGP1_21_ch18.haplotype2.HLA.vcf.gz (42.8%)
VN      1777    PGP1_21_ch18.haplotype2.HLA.vcf.gz (57.2%)
VN      2085    PGP1_21_ch15.haplotype2.HLA.vcf.gz (61.1%)
#SN Summary Numbers. Use `grep ^SN | cut -f 2-` to extract this part.
SN      Number of REF matches:  1328
SN      Number of ALT matches:  1315
SN      Number of REF mismatches:       0
SN      Number of ALT mismatches:       13
SN      Number of samples in GT comparison:     0
  • 13 ALT mismatches are listed below. Note I haven't removed the calls with 0 quality.
AM      chr6    31026712        A       T
AM      chr6    31076199        T       C,T
AM      chr6    31119716        T       C,T
AM      chr6    31190931        G       C,G
AM      chr6    31191665        T       G
AM      chr6    31235197        G       A
AM      chr6    31240312        T       A
AM      chr6    31248930        A       C
AM      chr6    31253771        C,T     T
AM      chr6    31314044        C       A
AM      chr6    31317555        G       A
AM      chr6    31317556        C       A
AM      chr6    31325341        C       A
  • position 31026712 quality for A is 0 while quality for T is 99.
$ grep 31026712 PGP1_21_ch15.haplotype2.HLA.vcf 
chr6    31026712        .       C       A       0       0       DP=1    GT:GQ:DP        1/1:0:1
$ grep 31026712 PGP1_21_ch18.haplotype2.HLA.vcf
chr6    31026712        .       C       T       99      0       DP=24   GT:GQ:DP        1/1:99:24
  • position 31076199 both have good quality score. But one of them has two possible ALTs. So it's deduced to a T.
$ grep 31076199 PGP1_21_ch15.haplotype2.HLA.vcf
chr6    31076199        .       A       T       48      0       DP=7    GT:GQ:DP        1/1:48:7
$ grep 31076199 PGP1_21_ch18.haplotype2.HLA.vcf
chr6    31076199        .       A       C,T     55      0       DP=9    GT:GQ:DP        1/2:7:9
  • position 31191665 both have none 0 quality score. But one of them has only 1 time coverage. I'll look into the pileup file.
$ grep 31191665 PGP1_21_ch15.haplotype2.HLA.vcf
chr6    31191665        .       A       T       57      0       DP=11   GT:GQ:DP        1/1:42:11
$ grep 31191665 PGP1_21_ch18.haplotype2.HLA.vcf
chr6    31191665        .       A       G       6       0       DP=1    GT:GQ:DP        1/1:6:1
  • I looked into the pileup file for this position 31191665. Very likely the ALT is a T.
$ grep 31191665 PGP1_21_ch15.haplotype2.pileup
chr6    31191665        a       T       42      57      53      11      GttTtTTTtTt     'FI'I<0'FFF
$ grep 31191665 PGP1_21_ch18.haplotype2.pileup
chr6    31191665        a       G       6       6       60      1       G       '
  • I wrote a program to remove the 0 quality entries in VCF
$ perl remove0qualityVCF.pl PGP1_21_ch15.haplotype2.HLA.vcf PGP1_21_ch15.haplotype2.HLA0.vcf
$ perl remove0qualityVCF.pl PGP1_21_ch18.haplotype2.HLA.vcf PGP1_21_ch18.haplotype2.HLA0.vcf
  • Add header
bgzip PGP1_21_ch18.haplotype2.HLA0.vcf
tabix -p vcf PGP1_21_ch18.haplotype2.HLA0.vcf.gz
bgzip PGP1_21_ch15.haplotype2.HLA0.vcf
tabix -p vcf PGP1_21_ch15.haplotype2.HLA0.vcf.gz
/home/kunzhang/softwares/vcftools_0.1.12b/bin/vcf-compare -d PGP1_21_ch15.haplotype2.HLA0.vcf.gz PGP1_21_ch18.haplotype2.HLA0.vcf.gz > compareCh15Ch18HLA0vcfd.txt
#VN The columns are:
#VN        1  .. number of sites unique to this particular combination of files
#VN        2- .. combination of files and space-separated number, a fraction of sites in the file
VN      1304    PGP1_21_ch15.haplotype2.HLA0.vcf.gz (39.9%)     PGP1_21_ch18.haplotype2.HLA0.vcf.gz (44.5%)
VN      1624    PGP1_21_ch18.haplotype2.HLA0.vcf.gz (55.5%)
VN      1962    PGP1_21_ch15.haplotype2.HLA0.vcf.gz (60.1%)
SN      Number of REF matches:  1304
SN      Number of ALT matches:  1298
SN      Number of REF mismatches:       0
SN      Number of ALT mismatches:       6
SN      Number of samples in GT comparison:     0

AM      chr6    31076199        T       C,T
AM      chr6    31119716        T       C,T
AM      chr6    31190931        G       C,G
AM      chr6    31191665        T       G
AM      chr6    31253771        C,T     T
AM      chr6    31314044        C       A

  • After looking at pileup, all of the mismatching ALT are resolved.
AM      chr6    31076199        T       C,T
Ch15 chr6    31076199        a       T       48      48      59      7       TttTTtT IIIIIFI
Ch18 chr6    31076199        a       Y       7       55      55      9       c$c$ttTtTTT     BBIFFIIIB
--> T
--> no found SNV from HiC+BAC vcf
AM      chr6    31119716        T       C,T
ch15 chr6    31119716        g       T       33      33      60      2       tt      II
ch18 chr6    31119716        g       Y       1       37      60      4       tTCt    FBIF
--> T
--> found in Hap2 of HiC+BAC vcf
AM      chr6    31190931        G       C,G
ch15 chr6    31190931        A       G       217     217     60      63      gGggggggggGGgggGGGGggggggGgGGGggGggggGggggGGGGGggGGGgGgGGGGGggg 
ch18 chr6    31190931        A       S       51      114     60      17      GccCCCcgCCGCgGcc^]c     FFFFFIF0BIII0FFFB
--> G
--> found in Hap2 of HiC+BAC vcf
AM      chr6    31191665        T       G
ch15 chr6    31191665        a       T       42      57      53      11      GttTtTTTtTt     'FI'I<0'FFF
ch18 chr6    31191665        a       G       6       6       60      1       G       '
--> T
--> found in Hap2 of HiC+BAC vcf
AM      chr6    31253771        C,T     T
ch15 chr6    31253771        G       Y       1       34      60      3       CcT     FIF
ch18 chr6    31253771        G       T       120     120     60      31      TTTTTttTTTTtTTttttTtTTTtTttTtt^]t 
--> T
--> found in Hap2 of HiC+BAC vcf
AM      chr6    31314044        C       A
ch15 chr6    31314044        g       C       4       11      7       2       ac      FF
ch18 chr6    31314044        g       A       43      43      5       17      *AAaaaaAAAAaaAAaa
--> A
--> no found SNV from HiC+BAC vcf

In /media/LTS_15T/Eric_LTS/PGP1_21_BWA_MEM

$ /home/kunzhang/softwares/samtools-0.1.19/samtools view -b PGP1_21_ch18.sorted.chr.bam chr6:30521122-31340586 > PGP1_21_ch18.HLAfrag.bam
$ /home/kunzhang/softwares/samtools-0.1.19/samtools view -b PGP1_21_ch15.sorted.chr.bam chr6:30521122-31340586 > PGP1_21_ch15.HLAfrag.bam
$ /home/kunzhang/softwares/bedtools-2.20/bin/bedtools bamtobed -i PGP1_21_ch18.HLAfrag.bam > PGP1_21_ch18.HLAfrag.bed
$ /home/kunzhang/softwares/bedtools-2.20/bin/bedtools bamtobed -i PGP1_21_ch15.HLAfrag.bam > PGP1_21_ch15.HLAfrag.bed
$ cat PGP1_21_ch15.HLAfrag.bed PGP1_21_ch18.HLAfrag.bed > PGP1_21_ch15ch18.HLAfrag.bed
$ sort -k1,1 -k2,2n PGP1_21_ch15ch18.HLAfrag.bed > PGP1_21_ch15ch18.HLAfrag.sorted.bed
$ /home/kunzhang/softwares/bedtools-2.20/bin/bedtools intersect -a PGP1_21_ch15ch18.HLAfrag.merge.bed -b HLAfrag.bed > intersectCh15Ch18HLAfrag.bed
$ awk '{s+=($3-$2)}END{print s}' intersectCh15Ch18HLAfrag.bed
750308
$ wc -l intersectCh15Ch18HLAfrag.bed
423 intersectCh15Ch18HLAfrag.bed
  • then coverage = 750308+423 = 750731bp
  • in the region of 819464bp


  • Let's find out if I have the complimentary (same as REF) call for those SNV calls
  • Quality score filters?!
  • Here's the low quality score for most of the reads in BAC pileup.
/media/LTS_15T/Eric_LTS/BACPoolPileup$ more Indx86.bac.pileup.vcf
#CHROM POS ID REF       ALT     QUAL    FILTER  INFO    FORMAT  data
chr1 53126 . A  T       4       0       DP=1    GT:GQ:DP        1/1:4:1
chr1 53196 . G  C       4       0       DP=1    GT:GQ:DP        1/1:4:1
chr1 54857 . T  C       4       0       DP=1    GT:GQ:DP        1/1:4:1
chr1 60013 . A  G       4       0       DP=1    GT:GQ:DP        1/1:4:1
chr1 60901 . A  G       8       0       DP=1    GT:GQ:DP        1/1:8:1
chr1 61579 . G  A       4       0       DP=1    GT:GQ:DP        1/1:4:1
chr1 63200 . T  C       8       0       DP=1    GT:GQ:DP        1/1:8:1
chr1 63950 . A  G       4       0       DP=1    GT:GQ:DP        1/1:4:1
chr1 67834 . T  C       4       0       DP=1    GT:GQ:DP        1/1:4:1
chr1 68659 . G  A       4       0       DP=1    GT:GQ:DP        1/1:4:1
chr1 68675 . T  C       4       0       DP=1    GT:GQ:DP        1/1:4:1
chr1 69428 . T  G       7       0       DP=1    GT:GQ:DP        1/1:7:1
chr1 72220 . G  A       4       0       DP=1    GT:GQ:DP        1/1:4:1
chr1 72627 . A  G       4       0       DP=1    GT:GQ:DP        1/1:4:1
chr1 74262 . T  C       4       0       DP=1    GT:GQ:DP        1/1:4:1
chr1 74822 . T  C       4       0       DP=1    GT:GQ:DP        1/1:4:1
chr1 74891 . T  C       4       0       DP=1    GT:GQ:DP        1/1:4:1
chr1 74893 . T  G       4       0       DP=1    GT:GQ:DP        1/1:4:1
  • compare to haplotype 1 to see if there are matches
$ perl remove0qualityVCF.pl PGP1_21_ch24.haplotype1.vcf PGP1_21_ch24.haplotype1.0.vcf
$ perl remove0qualityVCF.pl PGP1_21_ch23.haplotype1.vcf PGP1_21_ch23.haplotype1.0.vcf
$ bgzip PGP1_21_ch24.haplotype1.0.vcf
$ bgzip PGP1_21_ch23.haplotype1.0.vcf
$ tabix -p vcf PGP1_21_ch23.haplotype1.0.vcf.gz
$ tabix -p vcf PGP1_21_ch24.haplotype1.0.vcf.gz
$ /home/kunzhang/softwares/vcftools_0.1.12b/bin/vcf-compare -d PGP1_21_ch15.haplotype2.HLA0.vcf.gz PGP1_21_ch18.haplotype2.HLA0.vcf.gz PGP1_21_ch23.haplotype1.HLA.vcf.gz PGP1_21_ch24.haplotype1.HLA.vcf.gz > compareCh15Ch18Ch23Ch24HLA0vcfd.txt
#VN The columns are:
#VN        1  .. number of sites unique to this particular combination of files
#VN        2- .. combination of files and space-separated number, a fraction of sites in the file
VN      18      ch15.haplotype2.HLA (0.6%) ch23.haplotype1.HLA (3.0%)
VN      23      ch15.haplotype2.HLA (0.7%) ch18.haplotype2.HLA (0.8%) ch23.haplotype1.HLA (3.9%)
VN      28      ch18.haplotype2.HLA (1.0%) ch23.haplotype1.HLA. (4.7%)
VN      44      ch18.haplotype2.HLA (1.5%) ch24.haplotype1.HLA. (3.9%)
VN      116     ch15.haplotype2.HLA (3.6%) ch24.haplotype1.HLA (10.2%)
VN      258     ch15.haplotype2.HLA (7.9%) ch18.haplotype2.HLA (8.8%) ch24.haplotype1.HLA (22.6%)
VN      523     ch23.haplotype1.HLA (88.3%)
VN      723     ch24.haplotype1.HLA (63.4%)
VN      1023    ch15.haplotype2.HLA (31.3%) ch18.haplotype2.HLA (34.9%)
VN      1552    ch18.haplotype2.HLA (53.0%)
VN      1828    ch15.haplotype2.HLA (56.0%)
  • compare one P1 and one P2 fragments
#VN The columns are:
#VN        1  .. number of sites unique to this particular combination of files
#VN        2- .. combination of files and space-separated number, a fraction of sites in the file
VN      374     PGP1_21_ch15.haplotype2.HLA0.vcf.gz (11.5%)     PGP1_21_ch24.haplotype1.HLA.vcf.gz (32.8%)
VN      767     PGP1_21_ch24.haplotype1.HLA.vcf.gz (67.2%)
VN      2892    PGP1_21_ch15.haplotype2.HLA0.vcf.gz (88.5%)
#SN Summary Numbers. Use `grep ^SN | cut -f 2-` to extract this part.
SN      Number of REF matches:  374
SN      Number of ALT matches:  372
SN      Number of REF mismatches:       0
SN      Number of ALT mismatches:       2
SN      Number of samples in GT comparison:     0
#RM REF mismatches. The columns are:
#RM        1 .. chromosome
#RM        2 .. position
#RM        3 .. REF in the first file
#RM        4 .. differing REF
AM      chr6    31239101        T       A
AM      chr6    31266952        A       G
chr6    31239101        G       T       135     135     44      36      t$ttttttttTtttttttTtttTtTtttttTtTTttt   BBFFFFFFFFIIFFIBIBIFB0FFFFBF<FBBB<0f
chr6    31239101        G       A       33      33      49      2       aa      FF
--> no SNV from HiC+BAC vcf
chr6    31266952        t       A       171     202     54      59      AAAaaAaaaAAAAAAAAAAAaA,$AAAAaAAaAAAAAaaAAAaaaaAAAAaAAAAaa^[a^Ya^[a    BFFFF<IIIF'IIIFIIBFIIIBIIFFIFFIIIIIIBBIIhFFFFIIFIFIFIF'F0FB
chr6    31266952        t       G       117     117     25      30      GGGGGGGGGgGgGGGgGGGGGGGGggGg^<g^:g      BFFBFFIIIIIIIIFIIIIIIFFFBFFFFF
--> no SNV from HiC+BAC vcf
#VN The columns are:
#VN        1  .. number of sites unique to this particular combination of files
#VN        2- .. combination of files and space-separated number, a fraction of sites in the file
VN      41      PGP1_21_ch15.haplotype2.HLA0.vcf.gz (1.3%)      PGP1_21_ch23.haplotype1.HLA.vcf.gz (6.9%)
VN      551     PGP1_21_ch23.haplotype1.HLA.vcf.gz (93.1%)
VN      3225    PGP1_21_ch15.haplotype2.HLA0.vcf.gz (98.7%)
#SN Summary Numbers. Use `grep ^SN | cut -f 2-` to extract this part.
SN      Number of REF matches:  41
SN      Number of ALT matches:  40
SN      Number of REF mismatches:       0
SN      Number of ALT mismatches:       1
SN      Number of samples in GT comparison:     0
AM      chr6    30928757        T       A,T
--> no SNV from HiC+BAC vcf