EricChu:LabNotesMDA/2015-11-10: Difference between revisions
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