EricChu:LabNotesMDA/2015-10-28
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