EricChu:LabNotesMDA/2015-11-10: Difference between revisions

From ZhangLabWiki
Jump to navigation Jump to search
>Ericchu
(Created page with "==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.uniqu...")
 
>Ericchu
No edit summary
Line 22: Line 22:
  tabix -p vcf PGP1_21_ch15.unique.haplotype2.vcf.gz
  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
  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)
==Export Mismatch from Chamber 15 to chamber 18 and CGI data==
$ 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 &
$ 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 &
Then compare with pileup in the alternative chamber

Revision as of 20:04, 10 November 2015

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)

Export Mismatch from Chamber 15 to chamber 18 and CGI data

$ 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 &
$ 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 &
Then compare with pileup in the alternative chamber